Add grid_tools/

1 files changed, 185 insertions, 0 deletions

diff --git a/grid_tools/README.md b/grid_tools/README.md
new file mode 100644
index 0000000..f8cab3c
--- /dev/null
+++ b/grid_tools/README.md
@@ -0,0 +1,185 @@
+# Grid tools
+
+This directory contains a set of Python scripts to convert grids to
+the [Geodetic TIFF grids (GTG)](https://github.com/OSGeo/PROJ/blob/master/docs/source/specifications/geodetictiffgrids.rst),
+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 "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 "osgeo/gdal" can also be used. It is larger, using a Ubuntu base, whereas
+"osgeo/gdal" uses a smaller Alpine Linux base)
+
+For Windows, the development version of [GisInternals builds](http://gisinternals.com/development.php)
+can also be used. Or the "gdal-dev-python" package of [OSGeo4W](https://trac.osgeo.org/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](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 column 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](https://github.com/OSGeo/PROJ/wiki/Units-of-NTv2-accuracy-samples-%3F). 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 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](vertoffset_grid_to_gtiff.py) script.
+
+Usage:
+```
+usage: vertoffset_grid_to_gtiff.py [-h] --type
+                                   {GEOGRAPHIC_TO_VERTICAL,VERTICAL_TO_VERTICAL}
+                                   --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.
+* --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 \
+    --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 osgeo/gdal:alpine-normal-latest python3 $PWD/vertoffset_grid_to_gtiff.py \
+    --type GEOGRAPHIC_TO_VERTICAL \
+    --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](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](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](https://gdal.org/programs/gdal_translate.html) or
+[gdal_edit.py](https://gdal.org/programs/gdal_edit.html)
+
+Note: this script will not compress data. This must be done in a prior step.