Merge pull request #3510 from rouault/nadcon5

Add support for NADCON5 grids and transformation method
OSGeo · Dec 17, 2022 · 252c7b6 · 252c7b6
2 parents 1c1b3a3 + e732237
commit 252c7b6
Show file tree

Hide file tree

Showing 39 changed files with 2,993 additions and 162 deletions.
diff --git a/data/sql/concatenated_operation.sql b/data/sql/concatenated_operation.sql
@@ -421,6 +421,8 @@ INSERT INTO "concatenated_operation" VALUES('EPSG','8659','Kertau (RSO) to WGS 8
 INSERT INTO "usage" VALUES('EPSG','10751','concatenated_operation','EPSG','8659','EPSG','1309','EPSG','1164');
 INSERT INTO "concatenated_operation" VALUES('EPSG','9103','NAD27 to ITRF2014 (1)','For use with legacy data - see CT code 9104 for alternative for new areas. Note that steps 1 and 2 are documented in the geog2D domain, steps 3 and 4 in the geocentric domain. Steps 3 and 4 may be implemented in one operation using CT code 8970.','EPSG','4267','EPSG','7789',1.5,'IOGP-Usa GoM legacy',0);
 INSERT INTO "usage" VALUES('EPSG','10917','concatenated_operation','EPSG','9103','EPSG','3357','EPSG','1136');
+INSERT INTO "concatenated_operation" VALUES('EPSG','9104','NAD27 to ITRF2014 (2)','For use in new operations - see CT code 9103 for alternative for legacy data. Note that steps 1 to 5 are documented in the geog2D domain, steps 6 and 7 in the geocentric domain. Steps 6 and 7 may be implemented in one operation using CT code 8970.','EPSG','4267','EPSG','7789',0.2,'IOGP-Usa GoM non-legacy',0);
+INSERT INTO "usage" VALUES('EPSG','10918','concatenated_operation','EPSG','9104','EPSG','3357','EPSG','1136');
 INSERT INTO "concatenated_operation" VALUES('EPSG','9336','NAD27 to NAD83(CSRS)v4 (3)','Can be taken as an approximate transformation NAD27 to WGS 84 - see code 8585.','EPSG','4267','EPSG','8246',2.5,'EPSG-Can AB',0);
 INSERT INTO "usage" VALUES('EPSG','14011','concatenated_operation','EPSG','9336','EPSG','2376','EPSG','1151');
 INSERT INTO "concatenated_operation" VALUES('EPSG','9337','NTF (Paris) to RGF93 v1 (1)','See transformation code 7811 for an alternative which uses the NTv2 method as an emulation of the geocentric interpolation in the second step.','EPSG','4807','EPSG','4171',1.0,'IOGP-Fra',0);

diff --git a/data/sql/concatenated_operation_step.sql b/data/sql/concatenated_operation_step.sql
@@ -425,6 +425,13 @@ INSERT INTO "concatenated_operation_step" VALUES('EPSG','9103',1,'EPSG','1241');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','9103',2,'EPSG','8971');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','9103',3,'EPSG','7807');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','9103',4,'EPSG','7790');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',1,'EPSG','8555');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',2,'EPSG','8556');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',3,'EPSG','8861');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',4,'EPSG','8862');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',5,'EPSG','8559');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',6,'EPSG','7807');
+INSERT INTO "concatenated_operation_step" VALUES('EPSG','9104',7,'EPSG','7790');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','9336',1,'EPSG','1313');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','9336',2,'EPSG','9244');
 INSERT INTO "concatenated_operation_step" VALUES('EPSG','9337',1,'EPSG','1763');

diff --git a/data/sql/grid_alternatives_generated_noaa.sql b/data/sql/grid_alternatives_generated_noaa.sql
diff --git a/data/sql/grid_transformation.sql b/data/sql/grid_transformation.sql
diff --git a/data/sql/metadata.sql b/data/sql/metadata.sql
@@ -18,4 +18,4 @@ INSERT INTO "metadata" VALUES('PROJ.VERSION', '${PROJ_VERSION}');
 
 -- Version of the PROJ-data package with which this database is the most
 -- compatible.
-INSERT INTO "metadata" VALUES('PROJ_DATA.VERSION', '1.12');
+INSERT INTO "metadata" VALUES('PROJ_DATA.VERSION', '1.13');
diff --git a/data/sql/nadcon5_concatenated_operations.sql b/data/sql/nadcon5_concatenated_operations.sql
diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql
@@ -1098,7 +1098,7 @@ CREATE TABLE grid_alternatives(
     proj_grid_name TEXT NOT NULL,                   -- PROJ >= 7 grid name (e.g us_nga_egm08_25.tif)
     old_proj_grid_name TEXT,                        -- PROJ < 7 grid name (e.g egm08_25.gtx)
     proj_grid_format TEXT NOT NULL,                 -- 'GTiff', 'GTX', 'NTv2', JSON
-    proj_method TEXT NOT NULL,                      -- hgridshift, vgridshift, geoid_like, geocentricoffset, tinshift or velocity_grid
+    proj_method TEXT NOT NULL,                      -- gridshift, hgridshift, vgridshift, geoid_like, geocentricoffset, tinshift or velocity_grid
     inverse_direction BOOLEAN NOT NULL CHECK (inverse_direction IN (0, 1)), -- whether the PROJ grid direction is reversed w.r.t to the authority one (TRUE in that case)
     package_name TEXT,                              -- no longer used. Must be NULL
     url TEXT,                                       -- optional URL where to download the PROJ grid
@@ -1108,7 +1108,7 @@ CREATE TABLE grid_alternatives(
 
     CONSTRAINT fk_grid_alternatives_grid_packages FOREIGN KEY (package_name) REFERENCES grid_packages(package_name) ON DELETE CASCADE,
     CONSTRAINT check_grid_alternatives_grid_fromat CHECK (proj_grid_format IN ('GTiff', 'GTX', 'NTv2', 'JSON')),
-    CONSTRAINT check_grid_alternatives_proj_method CHECK (proj_method IN ('hgridshift', 'vgridshift', 'geoid_like', 'geocentricoffset', 'tinshift', 'velocity_grid')),
+    CONSTRAINT check_grid_alternatives_proj_method CHECK (proj_method IN ('gridshift', 'hgridshift', 'vgridshift', 'geoid_like', 'geocentricoffset', 'tinshift', 'velocity_grid')),
     CONSTRAINT check_grid_alternatives_inverse_direction CHECK (NOT(proj_method = 'geoid_like' AND inverse_direction = 1)),
     CONSTRAINT check_grid_alternatives_package_name CHECK (package_name IS NULL),
     CONSTRAINT check_grid_alternatives_direct_download_url CHECK (NOT(direct_download IS NULL AND url IS NOT NULL)),

diff --git a/data/sql_filelist.cmake b/data/sql_filelist.cmake
@@ -37,6 +37,7 @@ set(SQL_FILES
   "${SQL_DIR}/iau.sql"
   "${SQL_DIR}/grid_alternatives.sql"
   "${SQL_DIR}/grid_alternatives_generated_noaa.sql"
+  "${SQL_DIR}/nadcon5_concatenated_operations.sql"
   "${SQL_DIR}/customizations.sql"
   "${SQL_DIR}/nkg_post_customizations.sql"
   "${SQL_DIR}/commit.sql"

diff --git a/data/tests/nkgrf03vel_realigned_extract_tiled_256x256.tif b/data/tests/nkgrf03vel_realigned_extract_tiled_256x256.tif
diff --git a/data/tests/us_noaa_nadcon5_nad83_2007_nad83_2011_alaska_extract.tif b/data/tests/us_noaa_nadcon5_nad83_2007_nad83_2011_alaska_extract.tif
diff --git a/data/tests/us_noaa_nadcon5_nad83_2007_nad83_2011_conus_extract.tif b/data/tests/us_noaa_nadcon5_nad83_2007_nad83_2011_conus_extract.tif
diff --git a/docs/source/operations/transformations/gridshift.rst b/docs/source/operations/transformations/gridshift.rst
@@ -0,0 +1,70 @@
+.. _gridshift:
+
+================================================================================
+General grid shift
+================================================================================
+
+.. versionadded:: 9.2.0
+
+Translation of geodetic coordinates using a grid shift.
+
++-----------------+-------------------------------------------------------------------+
+| **Alias**       | gridshift                                                         |
++-----------------+-------------------------------------------------------------------+
+| **Domain**      | 2D and 3D                                                         |
++-----------------+-------------------------------------------------------------------+
+| **Input type**  | Geodetic coordinates (horizontal), meters (vertical)              |
++-----------------+-------------------------------------------------------------------+
+| **Output type** | Geodetic coordinates (horizontal), meters (vertical)              |
++-----------------+-------------------------------------------------------------------+
+
+The transformation may apply horizontal geodetic offsetting and/or vertical
+(ellipsoidal or orthometric height) offsetting, depending on the type of the
+grid(s).
+
+This is a generalization of the :ref:`hgridshift` and :ref:`vgridshift` methods,
+that may be used in particular for US NADCON5 grids that contain both horizontal
+geodetic and ellipsoidal height offsets.
+
+
+Example
+-------------------------------------------------------------------------------
+
+Transformation of a point from NAD83(NSRS2007) to NAD83(2011) on conterminous USA
+using a NADCON5 grid that contains horizontal and ellipsoidal height offsets.
+
+::
+
+    +proj=gridshift +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif
+
+Parameters
+################################################################################
+
+Required
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. option:: +grids=<list>
+
+    Comma-separated list of grids to load. If a grid is prefixed by an ``@`` the
+    grid is considered optional and PROJ will the not complain if the grid is
+    not available.
+
+    Grids must be in GeoTIFF format (:ref:`geodetictiffgrids`) and have an
+    explicit TYPE metadata item whose value is ``HORIZONTAL_OFFSET``,
+    ``GEOGRAPHIC_3D_OFFSET``, ``ELLIPSOIDAL_HEIGHT_OFFSET``
+    ``VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL`` or ``VERTICAL_OFFSET_VERTICAL_TO_VERTICAL``.
+
+Optional
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. option:: +interpolation=bilinear/biquadratic
+
+    Default is bilinear, unless the grid contains a ``interpolation_method``
+    metadata item specifying the method.
+    Biquadratic is typically used for NADCON5 grids, and is defined in
+    `NOAA Technical Memorandum NOS NGS 84 - Biquadratic Interpolation
+    <https://geodesy.noaa.gov/library/pdfs/NOAA_TM_NOS_NGS_0084.pdf>`__
+
+.. option:: +no_z_transform
+
+    If specified, vertical coordinate transformation will be skipped.
diff --git a/docs/source/operations/transformations/hgridshift.rst b/docs/source/operations/transformations/hgridshift.rst
@@ -43,6 +43,8 @@ about all three formats can be found in the GDAL documentation and/or driver sou
 code. GDAL reads and writes all three formats. Using GDAL for construction of
 new grids is recommended.
 
+To apply as well ellipsoidal height differences sometimes present in some grids
+(such as US NADCON5 grids), use the :ref:`gridshift` method.
 
 Temporal gridshifting
 ################################################################################

diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst
@@ -14,6 +14,7 @@ systems are based on different datums.
    defmodel
    deformation
    geogoffset
+   gridshift
    helmert
    horner
    molodensky

diff --git a/docs/source/specifications/geodetictiffgrids.rst b/docs/source/specifications/geodetictiffgrids.rst
@@ -217,6 +217,21 @@ is an easy way to inspect such grid files:
       as a speed per year for temporal gridshifting.
       Corresponds to PROJ :ref:`hgridshift` method.
 
+    - ``GEOGRAPHIC_3D_OFFSET``: implies the presence of at least 3 samples.
+      The first sample must contain the latitude offset, the second
+      sample must contain the longitude offset and the third one the ellipsoidal
+      height difference. Typically used for NADCON5 grids.
+      Added in PROJ 9.2
+      Corresponds to PROJ :ref:`gridshift` method.
+
+    - ``ELLIPSOIDAL_HEIGHT_OFFSET``: implies the presence of one sample with
+      the ellipsoidal height difference. Generally used in combination with
+      another grid of type ``HORIZONTAL_OFFSET`` to perform Geographic 3D
+      offseting when the horizontal and vertical grids do not have the same
+      resolution, as found in some NADCON5 grids.
+      Added in PROJ 9.2
+      Corresponds to PROJ :ref:`gridshift` method.
+
     - ``VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL``: implies the presence of at least one sample.
       The first sample must contain the vertical adjustment. Must be used when
       the source/interpolation CRS is a Geographic CRS and the target CRS a Vertical CRS.
@@ -253,14 +268,21 @@ is an easy way to inspect such grid files:
 
     Values recognized by PROJ for this Item are currently:
 
-    + ``latitude_offset``: valid for TYPE=HORIZONTAL_OFFSET. Sample values should be
+    + ``latitude_offset``: valid for TYPE=HORIZONTAL_OFFSET or
+      GEOGRAPHIC_3D_OFFSET . Sample values should be
       the value to add a latitude expressed in the CRS encoded in the GeoKeys
       to obtain a latitude value expressed in the target CRS.
 
-    + ``longitude_offset``: valid for TYPE=HORIZONTAL_OFFSET. Sample values should be
+    + ``longitude_offset``: valid for TYPE=HORIZONTAL_OFFSET or
+      GEOGRAPHIC_3D_OFFSET . Sample values should be
       the value to add a longitude expressed in the CRS encoded in the GeoKeys
       to obtain a longitude value expressed in the target CRS.
 
+    + ``ellipsoidal_height_offset``: valid for TYPE=ELLIPSOIDAL_HEIGHT_OFFSET or
+      GEOGRAPHIC_3D_OFFSET . Sample values should be the value to add to the
+      ellipsoidal height of the source CRS to obtain the ellipsoidal height of
+      the target CRS.
+
     + ``geoid_undulation``: valid for TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL.
       For a source CRS being a geographic CRS and a target CRS being a vertical CRS,
       sample values should be the value to add to a geoid-related height (that
@@ -415,6 +437,10 @@ is an easy way to inspect such grid files:
     of this grid).
     Will be ignored by PROJ (this information can be inferred by the grids extent)
 
+  * The ``interpolation_method`` metadata item may be present to indicate the
+    interpolation method to apply. ``bilinear`` or ``biquadratic`` are supported.
+    If not specified, defaults to ``bilinear``.
+
 Example
 +++++++
 
@@ -648,5 +674,7 @@ georeferencing is contained within the first 40 KB of the file.
 Revisions
 +++++++++
 
+* v0.3: addition of TYPE=GEOGRAPHIC_3D_OFFSET, ELLIPSOIDAL_HEIGHT_OFFSET and
+        interpolation_method (PROJ 9.2)
 * v0.2: addition of "arc-seconds per year" as a valid unit (PROJ 9.1.1)
 * v0.1: initial version for PROJ 7.0
diff --git a/scripts/build_db.py b/scripts/build_db.py
@@ -688,7 +688,7 @@ def fill_helmert_transformation(proj_db_cursor):
             '?,?,?, ?, ?,?,?, ?,?, ?,?, ?, ?,?,?,?,?, ?,?,?,?,?, ?,?,?, ?,?,?,?,?, ?,?,?,?,?, ?,?,?, ?,?,?, ?,?,?,?,?, ?,?)', arg)
 
 def fill_grid_transformation(proj_db_cursor):
-    proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, coord_op_method_code, coord_op_method_name, source_crs_code, target_crs_code, coord_op_accuracy, coord_tfm_version, epsg_coordoperation.deprecated, epsg_coordoperation.remarks FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_type IN ('transformation', 'point motion operation') AND (coord_op_method_name LIKE 'Geographic3D to%' OR coord_op_method_name LIKE 'Geog3D to%' OR coord_op_method_name LIKE 'Point motion by grid%' OR coord_op_method_name LIKE 'Vertical Offset by Grid Interpolation%' OR coord_op_method_name IN ('NADCON', 'NADCON5 (2D)', 'NTv1', 'NTv2', 'VERTCON', 'Geocentric translation by Grid Interpolation (IGN)'))")
+    proj_db_cursor.execute("SELECT coord_op_code, coord_op_name, coord_op_method_code, coord_op_method_name, source_crs_code, target_crs_code, coord_op_accuracy, coord_tfm_version, epsg_coordoperation.deprecated, epsg_coordoperation.remarks FROM epsg.epsg_coordoperation LEFT JOIN epsg.epsg_coordoperationmethod USING (coord_op_method_code) WHERE coord_op_type IN ('transformation', 'point motion operation') AND (coord_op_method_name LIKE 'Geographic3D to%' OR coord_op_method_name LIKE 'Geog3D to%' OR coord_op_method_name LIKE 'Point motion by grid%' OR coord_op_method_name LIKE 'Vertical Offset by Grid Interpolation%' OR coord_op_method_name IN ('NADCON', 'NADCON5 (2D)', 'NADCON5 (3D)', 'NTv1', 'NTv2', 'VERTCON', 'Geocentric translation by Grid Interpolation (IGN)'))")
     for (code, name, method_code, method_name, source_crs_code, target_crs_code, coord_op_accuracy, coord_tfm_version, deprecated, remarks) in proj_db_cursor.fetchall():
         expected_order = 1
         max_n_params = 3 if method_name == 'Geocentric translation by Grid Interpolation (IGN)' else 2
@@ -708,7 +708,12 @@ def fill_grid_transformation(proj_db_cursor):
                 order_inc = 1
             order += order_inc
             first = False
-            assert order <= max_n_params
+            # NADCON5 lists 3 grids (lat_shift, lon_shift, ellipsoidal_height_shift). Our database
+            # can only list 2. Truncate. Not critical as we ultimately have one GeoTIFF
+            # grid for the 3 original grids
+            if method_name == "NADCON5 (3D)" and order > max_n_params:
+                break
+            assert order <= max_n_params, (code, name)
             assert order == expected_order, (code, name, method_code, method_name, param_code, param_name, order)
             if parameter_value is not None:
                 assert param_value_file_ref is None or len(param_value_file_ref) == 0, (order, parameter_code, parameter_name, parameter_value, param_value_file_ref, uom_code)
@@ -737,7 +742,7 @@ def fill_grid_transformation(proj_db_cursor):
             grid2_param_code = param_code[1]
             grid2_param_name = param_name[1]
             grid2_value = param_value[1]
-        elif method_code == 1074: # NADCON5 (2D)
+        elif method_code in (1074, 1075): # NADCON5 (2D) and NADCON5 (3D)
             assert param_code[1] == 8658, param_code[1]
             grid2_param_auth_name = EPSG_AUTHORITY
             grid2_param_code = param_code[1]

diff --git a/scripts/build_grid_alternatives_generated_noaa.py b/scripts/build_grid_alternatives_generated_noaa.py
@@ -90,7 +90,10 @@
 script_dir_name = os.path.dirname(os.path.realpath(__file__))
 sql_dir_name = os.path.join(os.path.dirname(script_dir_name), 'data', 'sql')
 
-f = open(os.path.join(sql_dir_name, 'grid_alternatives_generated') + '.sql', 'wb')
+out_filename = os.path.join(sql_dir_name, 'grid_alternatives_generated_noaa') + '.sql'
+#print(out_filename)
+
+f = open(out_filename, 'wb')
 f.write("--- This file has been generated by scripts/build_grid_alternatives_generated_noaa.py. DO NOT EDIT !\n\n".encode('UTF-8'))
 
 f.write("-- NADCON (NAD27 -> NAD83) entries\n\n".encode('UTF-8'))
@@ -163,4 +166,60 @@
                               '%s', 1, 1, NULL);""" % (las_filename, tiff_name, ntv2_name, 'https://cdn.proj.org/' + tiff_name)
     f.write((sql + '\n').encode('UTF-8'))
 
+
+
+f.write("-- NADCON5 entries\n\n".encode('UTF-8'))
+
+nadcon5 = ["as62.nad83_1993.as",
+           "gu63.nad83_1993.guamcnmi",
+           "nad27.nad83_1986.alaska",
+           "nad27.nad83_1986.conus",
+           "nad83_1986.nad83_1992.alaska",
+           "nad83_1986.nad83_1993.hawaii",
+           "nad83_1986.nad83_1993.prvi",
+           "nad83_1986.nad83_harn.conus",
+           "nad83_1992.nad83_2007.alaska",
+           "nad83_1993.nad83_1997.prvi",
+           "nad83_1993.nad83_2002.as",
+           "nad83_1993.nad83_2002.guamcnmi",
+           "nad83_1993.nad83_pa11.hawaii",
+           "nad83_1997.nad83_2002.prvi",
+           "nad83_2002.nad83_2007.prvi",
+           "nad83_2002.nad83_ma11.guamcnmi",
+           "nad83_2002.nad83_pa11.as",
+           "nad83_2007.nad83_2011.alaska",
+           "nad83_2007.nad83_2011.conus",
+           "nad83_2007.nad83_2011.prvi",
+           "nad83_fbn.nad83_2007.conus",
+           "nad83_harn.nad83_fbn.conus",
+           "ohd.nad83_1986.hawaii",
+           "pr40.nad83_1986.prvi",
+           #"sg1897.sg1952.stgeorge",
+           "sg1952.nad83_1986.stgeorge",
+           "sl1952.nad83_1986.stlawrence",
+           #"sp1897.sp1952.stpaul",
+           "sp1952.nad83_1986.stpaul",
+           #"ussd.nad27.conus",
+]
+for subdir in nadcon5:
+    tiff_name = 'us_noaa_nadcon5_' + subdir.replace('.', '_') + '.tif'
+    lat_filename = 'nadcon5.' + subdir + '.lat.trn.20160901.b'
+    sql = """INSERT INTO grid_alternatives(original_grid_name,
+                          proj_grid_name,
+                          old_proj_grid_name,
+                          proj_grid_format,
+                          proj_method,
+                          inverse_direction,
+                          package_name,
+                          url, direct_download, open_license, directory)
+                  VALUES ('%s',
+                          '%s',
+                          NULL,
+                          'GTiff',
+                          'gridshift',
+                          0,
+                          NULL,
+                          '%s', 1, 1, NULL);""" % (lat_filename, tiff_name, 'https://cdn.proj.org/' + tiff_name)
+    f.write((sql + '\n').encode('UTF-8'))
+
 f.close()