grid_tools

Grid tools

This directory contains a set of Python scripts to convert grids to the Geodetic TIFF grids (GTG), and check them

Requirements

GDAL >= 3.1 (currently master) is required to write geoid models with GeoTIFF 1.1 encoding of Geographic 3D CRS. Otherwise, GDAL 2.4/3.0 can be used to produce horizontal shift grids.

The "ghcr.io/osgeo/gdal:alpine-normal-latest" docker image documented at https://github.com/OSGeo/gdal/tree/master/gdal/docker can be used for that purpose (the default "ghcr.io/osgeo/gdal" can also be used. It is larger, using a Ubuntu base, whereas "ghcr.io/osgeo/gdal:alpine-normal-latest" uses a smaller Alpine Linux base)

For Windows, the development version of GisInternals builds can also be used. Or the "gdal-dev-python" package of OSGeo4W

Python 3 is suggested, but the scripts should still be compatible of Python 2.7

Converting from NTv2 to GTG

With the ntv2_to_gtiff.py script.

Usage:

usage: ntv2_to_gtiff.py [-h] --source-crs SOURCE_CRS --target-crs TARGET_CRS
                        --copyright COPYRIGHT [--description DESCRIPTION]
                        [--do-not-write-accuracy-samples]
                        [--positive-longitude-shift-value {east,west}]
                        [--uint16-encoding] [--datetime DATETIME]
                        [--accuracy-unit {arc-second,metre,unknown}]
                        [--area-of-use AREA_OF_USE]
                        source dest

Options:

• --source-crs SOURCE_CRS: the source CRS, as a "EPSG:XXXX" value or a CRS WKT string. Mandatory. Must be a Geographic 2D CRS.
• --target-crs SOURCE_CRS: the target CRS, as a "EPSG:XXXX" value or a CRS WKT string. Mandatory. Must be a Geographic 2D CRS.
• --copyright COPYRIGHT: Attribution and license to be stored in the TIFF Copyright tag. Mandatory
• --area-of-use AREA_OF_USE: To specify a textual area of use for the grid. For example "France", "Germany - Saxony". Strongly recommended.
• --description DESCRIPTION: Text describing the grid. If not provided, a default value will be automatically set. For example "NAD27 (EPSG:4267) to NAD83 (EPSG:4269). Converted from ntv2_0.gsb (last updated on 1995-07-05)"
• --do-not-write-accuracy-samples: To prevent accuracy samples from the NTv2 grid to be written in the output. Note: if it is detected that those samples are set to dummy values (negative values), they will automatically be discarded
• --positive-longitude-shift-value {east,west}: To force the convention for the longitude shift value. NTv2 uses a positive-is-west convention, that is confusing. By default, the script will negate the sign of the longitude shift values to output a positive-is-east convention. Setting this option to "west" will preserve the original convention. Not recommended
• --uint16-encoding: Whether values should be encoded on a 16-bit unsigned value, using a offset and scale floating point values. Default behaviour is to use Float32 encoding, which will preserve the binary values of the original NTv2 file
• --datetime DATETIME: to specify the value of the TIFF DateTime tag. Must be formatted as "YYYY:MM:DD HH:MM:SS" (note the use of colon as the separator for the Y-M-D part, as mandated by the TIFF specification). If not specified, the script will try to use the value from the corresponding NTv2 header
• --accuracy-unit {arc-second,metre,unknown}: to specify the unit of the accuracy samples. The NTv2 specification has been historically interpreted in different ways regarding that. Mandatory if accuracy samples are written (the script contains a few hardcoded rules for known datasets)

Example:

python3 ntv2_to_gtiff.py \
    --area-of-use "Canada" \
    --copyright "Derived from work by Natural Resources Canada. Open Government Licence - Canada: http://open.canada.ca/en/open-government-licence-canada" \
    --source-crs EPSG:4267 \
    --target-crs EPSG:4269 \
    ntv2_0.gsb ca_nrc_ntv2_0.tif

Using the Docker image:

docker run --rm -v /home:/home ghcr.io/osgeo/gdal:alpine-normal-latest python3 $PWD/ntv2_to_gtiff.py  \
    --area-of-use "Canada" \
    --copyright "Derived from work by Natural Resources Canada. Open Government Licence - Canada: http://open.canada.ca/en/open-government-licence-canada" \
    --source-crs EPSG:4267 \
    --target-crs EPSG:4269 \
    $PWD/ntv2_0.gsb $PWD/ca_nrc_ntv2_0.tif

Converting from GTX to GTG

With the vertoffset_grid_to_gtiff.py script.

Usage:

usage: vertoffset_grid_to_gtiff.py [-h] --type
                                   {GEOGRAPHIC_TO_VERTICAL,VERTICAL_TO_VERTICAL}
                                   [--parameter-name
                                   {geoid_undulation,hydroid_height,vertical_offset}]
                                   --source-crs SOURCE_CRS
                                   [--interpolation-crs INTERPOLATION_CRS]
                                   --target-crs TARGET_CRS --copyright
                                   COPYRIGHT [--description DESCRIPTION]
                                   [--encoding {float32,uint16,int32-scale-1-1000}]
                                   [--ignore-nodata] [--datetime DATETIME]
                                   [--area-of-use AREA_OF_USE]
                                   source dest

Options:

• --type {GEOGRAPHIC_TO_VERTICAL,VERTICAL_TO_VERTICAL}: specify the type of the grid. GEOGRAPHIC_TO_VERTICAL is for geoid-like grids transforming between a vertical CRS and ellipsoid heights. The value in the grid must be the offset to add to the height in the vertical CRS to get an ellipsoidal height. VERTICAL_TO_VERTICAL is for transformations between 2 vertical CRS: the value in the grid must be the offset to add to heights in the source CRS to get heights in the target CRS.
• --parameter-name {geoid_undulation,hydroid_height,vertical_offset}: semantics of the offset value. For type=GEOGRAPHIC_TO_VERTICAL, this option is mandatory, and the only valid choices are geoid_undulation and hydroid_height. For type=VERTICAL_TO_VERTICAL, vertical_offset is the only valid choice.`
• --source-crs SOURCE_CRS: the source CRS, as a "EPSG:XXXX" value or a CRS WKT string. Mandatory. For type=GEOGRAPHIC_TO_VERTICAL, this must be a Geographic 3D CRS. For type=VERTICAL_TO_VERTICAL, this must be a Vertical CRS.
• --interpolation-crs INTERPOLATION_CRS: the geographic CRS in which the grid is referenced to. This is ignored for type=GEOGRAPHIC_TO_VERTICAL (the interpolation CRS is the source CRS), but mandatory for type=VERTICAL_TO_VERTICAL
• --target-crs SOURCE_CRS: the target CRS, as a "EPSG:XXXX" value or a CRS WKT string. Mandatory. Must be a vertical CRS.
• --area-of-use AREA_OF_USE: To specify a textual area of use for the grid. For example "France", "Germany - Saxony". Strongly recommended.
• --description DESCRIPTION: Text describing the grid. If not provided, a default value will be automatically set. For example "NAD27 (EPSG:4267) to NAD83 (EPSG:4269). Converted from ntv2_0.gsb (last updated on 1995-07-05)"
• --encoding {float32,uint16,int32-scale-1-1000}: How values should be encoded. Defaults to float32. If specifying uint16, values will be encoded on a 16-bit unsigned value, using a offset and scale floating point values. If specifying int32-scale-1-1000, values will be conded on a 32-bit signed integer with a scaling factor of 1000 (this is the encoding naturally used by the .byn format of Canadian vertical shift grids)
• --ignore-nodata: whether the nodata value from the source should be encoded. This might be useful because there is an ambiguity for the GTX format. In some circumstances, the -88.8888 value must be interpreted as a nodata value, in others as a valid value. If not specifying this option, the script will print a warning if encountering a grid value that matches the nodata value, so as to call for manual verification.
• --datetime DATETIME: to specify the value of the TIFF DateTime tag. Must be formatted as "YYYY:MM:DD HH:MM:SS" (note the use of column as the separator for the Y-M-D part, as mandated by the TIFF specification). If not specified, the script will use the modification date of the source file.

Example:

python3 vertoffset_grid_to_gtiff.py \
    --type GEOGRAPHIC_TO_VERTICAL \
    --parameter-name geoid_undulation \
    --area-of-use "World" \
    --source-crs EPSG:4326 \
    --target-crs EPSG:5773 \
    --copyright "Public Domain. Derived from work at http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html" \
    egm96_15.gtx us_nga_egm96_15.tif 

Using the Docker image:

docker run --rm -v /home:/home ghcr.io/osgeo/gdal:alpine-normal-latest python3 $PWD/vertoffset_grid_to_gtiff.py \
    --type GEOGRAPHIC_TO_VERTICAL \
    --parameter-name geoid_undulation \
    --area-of-use "World" \
    --source-crs EPSG:4326 \
    --target-crs EPSG:5773 \
    --copyright "Public Domain. Derived from work at http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html" \
    $PWD/egm96_15.gtx $PWD/us_nga_egm96_15.tif

Checking compliance of GeoTIFF files with the GTG specification

With the check_gtiff_grid.py script.

Usage:

usage: check_gtiff_grid.py [-h] filename

If the script outputs anything and return the success error code, then everything is perfect.

Otherwise, it will messages among one of the three following categories:

• information: optional metadata missing, or extra metadata elements found. No corrective action is required
• warning: the file will probably be read correctly by PROJ, but corrective action may be required. Adressing those warnings is recommended
• error: the file will not be read correctly by PROJ. Corrective action is required. The return code of the script will be an error code.

Example of output on non-conformant file:

/home/even/gdal/git/gdal/gdal/byte.tif is NOT a valid PROJ GeoTIFF file.
The following errors were found (corrective action is required):
 - CRS found, but not a Geographic CRS
 - Datatype Byte of band 1 is not support

The following warnings were found (the file will probably be read correctly by PROJ, but corrective action may be required):
 - This file uses a RasterTypeGeoKey = PixelIsArea convention. While this will work properly with PROJ, a georeferencing using RasterTypeGeoKey = PixelIsPoint is more common.
 - Missing TYPE metadata item
 - GDAL area_of_use metadata item is missing. Typically used to capture plain text information about where the grid applies
 - TIFF tag ImageDescription is missing. Typically used to capture plain text information about what the grid does
 - TIFF tag Copyright is missing. Typically used to capture information on the source and licensing terms
 - TIFF tag DateTime is missing. Typically used to capture when the grid has been generated/converted
 - Image is uncompressed. Could potentially benefit from compression

The following informative message are issued (no corrective action required):
 - Metadata LAYER_TYPE=athematic on band 1 ignored

Optimize an existing GeoTIFF file

With the cloud_optimize_gtiff.py script.

Usage:

usage: cloud_optimize_gtiff.py [-h] source dest

This file takes an existing GeoTIFF file and optimize its layout, so that headers are put at the beginning of the file. This is used internally by the ntv2_to_gtiff.py and vertoffset_grid_to_gtiff.py script

This can be useful as the final step if preparing a file with the GDAL API / tools such as gdal_translate or gdal_edit.py

Note: this script will not compress data. This must be done in a prior step.