16.41. mapproject

The tool mapproject is used to orthorectify (mapproject) a camera image onto a DEM or datum. ASP is able to use mapprojected images to run stereo, see Section 6.1.7.

The mapproject program can be run using multiple processes and can be distributed over multiple machines (options --nodes-list and --processes).

This is particularly useful for ISIS cameras, as in that case any single process must use only one thread due to the limitations of ISIS. The tool splits the image up into tiles, distributes the tiles to sub-processes, and then merges the tiles into the requested output image. If the input image is small but takes a while to process, smaller tiles can be used to start more simultaneous processes (use the parameters --tile-size and --processes).

It is important to note that processing more tiles at a time may actually slow things down, if all processes write to the same disk and if processing each tile is dominated by the speed of writing to disk. Hence some benchmarking may be necessary for your camera type and storage setup.

The grid size, that is the dimension of pixels on the ground, set via the --tr option, should be in units as expected by the projection string obtained either from the DEM to project onto, or, if specified, from the --t_srs option. If the grid size is not set, it will be estimated as the mean ground sampling distance (GSD). See the --tr option for how this affects the extent of the output image.

If the resulting mapprojected images are used for stereo, all mapprojected images should have the same grid size and projection (Section 6.1.7.4).

16.41.1. Determination of projection

Use the option --t_srs to set a desired output projection. The projection should be local to the area of interest, in units of meter.

If this is not set, the projection from the DEM will be used, unless it is the longlat projection. In that case a good output projection is auto-determined.

For Earth, with the WGS84 datum, the auto-determined projection is UTM with an auto-computed zone, except for latitudes above 84 degrees North and below 80 degrees South, where the NSDIC polar stereographic projections are used.

For other Earth datums and other planetary bodies, the automatic determination produces a local stereographic projection. The projection center is found by a median calculation based on a sample of image pixels. Or consider using the cylindrical equal area projection.

To ensure the automatic projection determination is always invoked, overriding all other cases from above, use --t_srs auto.

All mapprojected images passed to stereo should use the same projection and grid size (Section 6.1.7).

16.41.2. Examples

16.41.2.1. Earth image with auto projection

Mapproject an image with a pinhole camera (Section 20.1), with a grid size of 2 meters, for Earth (WGS84):

mapproject --tr 2.0 DEM.tif image.tif camera.tsai output.tif

If the DEM has a longlat projection, a projection in meters is found first (Section 16.41.1).

16.41.2.2. Moon image with a custom projection

Mapproject a .cub file (Section 4.2). Such a file has both image and camera information. The planetary body is the Moon. Use a custom stereographic projection:

proj="+proj=stere +lat_0=-85.364 +lon_0=31.238 +R=1737400 +units=m +no_defs"

The grid size is set to 1 meter/pixel:

mapproject --tr 1.0 --t_srs "$proj" DEM.tif image.cub output.tif

16.41.2.3. RPC camera with bundle adjustment

Mapproject an image file with an RPC camera model (Section 8.22) in XML format. Use bundle-adjusted cameras (Section 16.5):

mapproject -t rpc --bundle-adjust-prefix ba/run \
  DEM.tif image.tif camera.xml output.tif

Here, the grid size is auto-determined.

See Section 8.22 for other ways of specifying the RPC camera model.

16.41.2.4. CSM camera

Mapproject with the CSM camera model (Section 8.12):

mapproject -t csm DEM.tif image.cub camera.json output.tif

16.41.2.5. Preexisting projection and grid size

The projection and grid size of a given mapprojected image can be borrowed when mapprojecting another image:

mapproject -t rpc                \
  --ref-map image1_map.tif       \
  DEM.tif image2.tif camera2.xml \
  image2_map.tif

This becomes important for stereo, when the two input mapprojected images must share these attributes (Section 6.1.7).

16.41.2.6. Multiple camera models

A DigitalGlobe / Maxar camera file has both an exact linescan model and an approximate RPC model. The RPC model is somewhat faster to use.

To choose between these with mapproject, invoke it either with -t dg or -t rpc. See Section 5 for more information.

16.41.2.7. Mapproject with no DEM

Mapproject onto the surface of zero height above a datum:

mapproject -t rpc WGS84 image.tif image.xml output.tif

Valid datum names include WGS84, NAD83, NAD27, D_MOON, D_MARS, and MOLA.

16.41.3. Saved metadata

The output image will have the following metadata saved to its geoheader:

  • INPUT_IMAGE_FILE, the input image name.

  • BUNDLE_ADJUST_PREFIX, the bundle adjustment prefix. Set to NONE if not present.

  • CAMERA_MODEL_TYPE, this is the session name, such as set with -t rpc.

  • CAMERA_FILE, the camera file used on input. Can be empty if the camera is contained within the input image.

  • DEM_FILE, the DEM used in mapprojection.

These metadata values are used to undo the mapprojection in stereo triangulation (Section 6.1.7.8). The geoheader can be inspected with gdalinfo (Section 16.25).

In addition, if the cameras have been bundle-adjusted, the translation and quaternion rotation from the .adjust file will be saved to the fields ADJUSTMENT_TRANSLATION and ADJUSTMENT_QUATERNION. This is useful for having mapprojection be reproducible if the separately stored .adjust files are not available.

These fields are editable with image_calc (Section 16.33.1.7), but this is not recommended.

16.41.4. Custom metadata

Custom metadata can be added to the geoheader with the --mo option. As with the GDAL -mo option, a value is everything after the first equal sign, so it may contain spaces and further equal signs. Pass several items in one quoted string, or repeat the option, with one VAR=VALUE pair each time:

mapproject                      \
  --mo 'VAR1=value with spaces' \
  --mo 'VAR2=second value'      \
  dem.tif image.tif camera.xml  \
  output.tif

16.41.5. ISIS compatibility

ISIS is in the process of merging in logic for the cam2map program that, when called with asp_map=true, uses the same per-pixel projection algorithm as mapproject.

When the grid size is not specified, each program auto-computes it from the camera and DEM, and the results may differ slightly (under 1e-6 relative error) due to implementation details. The produced projection bounds may disagree somewhat as well.

However, given the same inputs and the same grid size, these programs will create results that agree to float numerical precision at every grid point.

16.41.5.1. Example with camera in cube file

Consider an LRO NAC (Section 8.7) image having the camera information in the .cub file and an equicylindrical DEM on the Moon. Set up the projection string:

proj="+proj=eqc +lat_ts=0 +lon_0=15.3 +x_0=0 +y_0=0 +R=1737400 +units=m +no_defs"

Run both programs with the same inputs and grid size:

mapproject --tr 1.2 --t_srs "$proj" \
  dem.tif image.cub output_asp.tif

cam2map asp_map=true useproj=true   \
  projstring="$proj"                \
  pixres=mpp resolution=1.2         \
  dem=dem.tif from=image.cub to=output_isis.cub

The --tr option in mapproject and pixres=mpp resolution= in cam2map both set the grid size (ground sample distance) in meters per pixel. The DEM must be a GeoTIFF file with pixel values representing height above datum, in meters.

16.41.5.2. Example with CSM camera model

Use a CSM (Section 8.12) camera model from an ISD file instead of the SPICE camera in the cube. Assume the projection string is defined as above.

mapproject --tr 1.2 --t_srs "$proj" \
  dem.tif image.cub image.json output_asp.tif

cam2map asp_map=true useproj=true   \
  projstring="$proj"                \
  pixres=mpp resolution=1.2         \
  dem=dem.tif isd=image.json        \
  from=image.cub to=output_isis.cub

The CSM and SPICE camera models differ somewhat in implementation. When comparing map-projected images obtained with the two, the pixel difference may be on the order of 1e-4, and the extents may differ too.

However, mapproject and cam2map asp_map=true will still agree to float numerical precision at every grid point, when the same CSM camera model is used as input to both.

The isd parameter requires CSM plugins to be installed, which is the default in recent versions of ISIS.

16.41.5.3. Example with preexisting map

If an existing georeferenced image or .cub file is available, say named ref.cub, its projection, extent, and grid size can be reused:

cam2map asp_map=true dem=dem.tif     \
  map=ref.cub matchmap=true          \
  from=image.cub to=output_isis.cub

A PVL .map file with these projection parameters can also be passed to the map argument. Using a georeferenced image supports a wider range of projections, since GDAL reads the projection directly from the file.

16.41.5.4. Validation

The two results can be compared with geodiff:

geodiff output_isis.cub output_asp.tif -o run
gdalinfo -stats run-diff.tif

The expected pixel difference should be on the order of 1e-7 or less.

The gdalinfo command can also be employed to verify that both outputs have the same projection, extent, and grid size.

16.41.5.5. Running stereo

Images mapprojected with cam2map asp_map=true do not retain the original camera information in the cube header. The original cameras must be provided as additional arguments to parallel_stereo (Section 16.51). The first two arguments are the mapprojected images and the next two are the cameras:

parallel_stereo                          \
  left_cam2map.cub right_cam2map.cub     \
  left.cub right.cub                     \
  run/run                                \
  --dem dem.tif                          \
  --stereo-algorithm asp_mgm             \
  --subpixel-mode 9

point2dem --tr 1.2 run/run-PC.tif

This is the same workflow as stereo with ISIS mapprojected images (Section 14.6.2), where the camera arguments are the original unprojected .cub files.

The DEMs produced with these two invocations should agree with nearly float precision if invoked with the same choices of parameters.

16.41.5.6. Saved metadata

When cam2map is run with asp_map=true, the output .cub file contains an AspMapproject PVL group with the same fields as mapproject saves in the GeoTIFF geoheader (Section 16.41.3): INPUT_IMAGE_FILE, BUNDLE_ADJUST_PREFIX, CAMERA_MODEL_TYPE, CAMERA_FILE, and DEM_FILE. This group can be inspected with gdalinfo -mdd all or catlab.

16.41.5.7. Other cam2map options

The cam2map options minlat, maxlat, minlon, maxlon are supported in asp_map mode. Since cam2map can only create mapprojected cube files with projection in meters, these input bounds will be converted to projected units for the input projection.

When the option pixres is not set, the default value is used (camera), which amounts to auto-determining the grid size from the camera and DEM. This is the same behavior as for mapproject when the --tr option is not set.

This implementation supports multi-band inputs with float values. ASP’s mapproject program can handle single-band float images and RGB images. For other multi-band inputs (including RGBA), only the first band is used. When comparing with cam2map on multi-band inputs, mapproject results correspond to the first band of the cam2map output.

The defaultrange option is partially supported: defaultrange=camera unlocks the minlat/maxlat/minlon/maxlon overrides, and defaultrange=map uses the map file’s extent. Other values are silently ignored (bounds are auto-computed from the camera footprint and DEM).

The cam2map options interp, warpalgorithm, patchsize, trim, occlusion, and lonseam are ignored in asp_map mode, which always uses per-pixel bicubic interpolation.

16.41.5.8. Grid snapping and output extent

When cam2map is run with asp_map=true, the output grid is snapped to integer multiples of the grid size, consistent with mapproject. The produced extent goes half a grid pixel beyond the snapped grid on each side, because the grid is at pixel centers.

16.41.6. Usage

mapproject [options] <dem> <camera-image> <camera-model> <output-image>

16.41.7. Command-line options

--t_srs <string (default: “”)>

Specify the output projection as a GDAL projection string (WKT, GeoJSON, or PROJ). See Section 16.41.1 for details.

--tr <float>

Set the output file resolution (ground sample distance) in target georeferenced units per pixel. This may be in meters or degrees, depending on the projection. The center of each output pixel will be at integer multiples of this grid size, unless --gdal-tap is set.

-t, --session-type <string>

Select the stereo session type to use for processing. See Section 16.51.8 for the list of types.

--t_projwin <xmin ymin xmax ymax>

Limit the mapprojected image to this region, with the corners given in georeferenced coordinates (xmin ymin xmax ymax). See also --gdal-tap.

--t_pixelwin <xmin ymin xmax ymax>

Limit the mapprojected image to this region, with the corners given in pixels (xmin ymin xmax ymax). Max is exclusive.

--gdal-tap

Ensure that the output image bounds (as printed by gdalinfo, Section 16.25) are integer multiples of the grid size (as set with --tr). When --t_projwin is set and its entries are integer multiples of the grid size, that precise extent will be produced on output. This functions as the GDAL -tap option.

--bundle-adjust-prefix <name>

Use the camera adjustment obtained by previously running bundle_adjust with this output prefix.

--ref-map <filename>

Read the projection and grid size from this mapprojected image (Section 16.41.2.5).

--processes <integer>

Number of processes to use on each node (the default is for the program to choose).

--num-processes <integer>

Same as –processes. Used for backwards compatibility.

--nodes-list

List of available computing nodes to use. If not set, use the local machine. See also Section 8.19.

--tile-size

Size of square tiles to break up processing into. Each tile is run by an individual process. The default is 1024 pixels for ISIS cameras, as then each process is single-threaded, and 5120 pixels for other cameras, as such a process is multi-threaded, and disk I/O becomes a bigger consideration.

--mpp <float>

Set the output file resolution in meters per pixel.

--ppd <float>

Set the output file resolution in pixels per degree.

--datum-offset <float>

When projecting to a datum instead of a DEM, add this elevation offset to the datum.

--ot <type (default: Float32)>

Output data type, when the input is single channel. Supported types: Byte, UInt16, Int16, UInt32, Int32, Float32. If the output type is a kind of integer, values are rounded and then clamped to the limits of that type. This option will be ignored for multi-channel images, when the output type is set to be the same as the input type.

--nearest-neighbor

Use nearest neighbor interpolation instead of bicubic interpolation. This is not recommended, as it can result in artifacts.

--mo <string>

Write metadata to the output file, as with the GDAL -mo option. The option can be repeated. See Section 16.41.4 for details and examples.

--query-projection

Display the computed projection information, estimated ground sample distance, and projected bounding box. Save the projection as a WKT file named <output image>.wkt. Quit afterwards.

--query-pixel <double double>

Trace a ray from this input image pixel (values start from 0) to the ground. Print the intersection point with the DEM as lon, lat, height, then as DEM column, row, height. Quit afterwards.

--parallel-options <string (default: “–sshdelay 0.2”)>

Options to pass directly to GNU Parallel.

--no-geoheader-info

Do not write information in the geoheader. Otherwise mapproject will write the camera model type, the bundle adjustment prefix used, the rotation and translation from the .adjust file, the DEM it mapprojected onto, and the value of the --mo option.

--nodata-value <float (default: -32768)>

No-data value to use unless specified in the input image.

--suppress-output

Suppress output from sub-processes.

--threads <integer (default: 0)>

Select the number of threads to use for each process. If 0, use the value in ~/.vwrc.

--cache-size-mb <integer (default = 1024)>

Set the system cache size, in MB, for each process.

--no-bigtiff

Tell GDAL to not create BigTiff files.

--tif-compress <None|LZW|Deflate|Packbits>

TIFF compression method.

--cog

Write a cloud-optimized GeoTIFF (COG). See Section 6.2.7.

-v, --version

Display the version of software.

-h, --help

Display the help message.