Merge pull request OSGeo#8340 from rouault/pointmotionoperation

Coordinate transformation: enhacements related to PointMotionOperation and CoordinateMetadata
ojwb · Sep 15, 2023 · 566204a · 566204a
2 parents c73205d + c4f5b92
commit 566204a
Show file tree

Hide file tree

Showing 17 changed files with 321 additions and 40 deletions.
diff --git a/alg/gdaltransformer.cpp b/alg/gdaltransformer.cpp
@@ -2259,8 +2259,12 @@ void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
     const bool bMayInsertCenterLong =
         (bCanUseSrcGeoTransform && !oSrcSRS.IsEmpty() && hSrcDS &&
          CPLFetchBool(papszOptions, "INSERT_CENTER_LONG", true));
+    const char *pszSrcCoordEpoch =
+        CSLFetchNameValue(papszOptions, "SRC_COORDINATE_EPOCH");
+    const char *pszDstCoordEpoch =
+        CSLFetchNameValue(papszOptions, "DST_COORDINATE_EPOCH");
     if ((!oSrcSRS.IsEmpty() && !oDstSRS.IsEmpty() &&
-         (!oSrcSRS.IsSame(&oDstSRS) ||
+         (pszSrcCoordEpoch || pszDstCoordEpoch || !oSrcSRS.IsSame(&oDstSRS) ||
           (oSrcSRS.IsGeographic() && bMayInsertCenterLong))) ||
         pszCO)
     {
@@ -2298,16 +2302,12 @@ void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
             aosOptions.SetNameValue("COORDINATE_EPOCH", pszCoordEpoch);
         }
 
-        const char *pszSrcCoordEpoch =
-            CSLFetchNameValue(papszOptions, "SRC_COORDINATE_EPOCH");
         if (pszSrcCoordEpoch)
         {
             aosOptions.SetNameValue("SRC_COORDINATE_EPOCH", pszSrcCoordEpoch);
             oSrcSRS.SetCoordinateEpoch(CPLAtof(pszSrcCoordEpoch));
         }
 
-        const char *pszDstCoordEpoch =
-            CSLFetchNameValue(papszOptions, "DST_COORDINATE_EPOCH");
         if (pszDstCoordEpoch)
         {
             aosOptions.SetNameValue("DST_COORDINATE_EPOCH", pszDstCoordEpoch);

diff --git a/apps/gdaltransform.cpp b/apps/gdaltransform.cpp
@@ -61,6 +61,7 @@ static void Usage(const char *pszErrorMsg = nullptr)
 {
     printf("Usage: gdaltransform [--help] [--help-general]\n"
            "    [-i] [-s_srs srs_def] [-t_srs srs_def] [-to \"NAME=VALUE\"]\n"
+           "    [-s_coord_epoch epoch] [-t_coord_epoch epoch]\n"
            "    [-ct proj_string] [-order n] [-tps] [-rpc] [-geoloc] \n"
            "    [-gcp pixel line easting northing [elevation]]* [-output_xy]\n"
            "    [srcfile [dstfile]]\n"
@@ -178,6 +179,18 @@ MAIN_START(argc, argv)
                 exit(1);
             aosTO.SetNameValue("SRC_SRS", pszSRS);
         }
+        else if (EQUAL(argv[i], "-s_coord_epoch"))
+        {
+            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+            const char *pszCoordinateEpoch = argv[++i];
+            aosTO.SetNameValue("SRC_COORDINATE_EPOCH", pszCoordinateEpoch);
+        }
+        else if (EQUAL(argv[i], "-t_coord_epoch"))
+        {
+            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
+            const char *pszCoordinateEpoch = argv[++i];
+            aosTO.SetNameValue("DST_COORDINATE_EPOCH", pszCoordinateEpoch);
+        }
         else if (EQUAL(argv[i], "-ct"))
         {
             CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);

diff --git a/apps/gdalwarp_bin.cpp b/apps/gdalwarp_bin.cpp
@@ -70,7 +70,7 @@ static void Usage(const char *pszErrorMsg = nullptr)
         "    [-b|-srcband n]* [-dstband n]*\n"
         "    [-s_srs srs_def] [-t_srs srs_def] [-ct string]\n"
         "    [-to \"NAME=VALUE\"]* [-vshift | -novshift]\n"
-        "    [[-s_coord_epoch epoch] | [-t_coord_epoch epoch]]\n"
+        "    [-s_coord_epoch epoch] [-t_coord_epoch epoch]\n"
         "    [-order n | -tps | -rpc | -geoloc] [-et err_threshold]\n"
         "    [-refine_gcps tolerance [minimum_gcps]]\n"
         "    [-te xmin ymin xmax ymax] [-te_srs srs_def]\n"

diff --git a/apps/ogr2ogr_bin.cpp b/apps/ogr2ogr_bin.cpp
@@ -111,8 +111,8 @@ static void Usage(const char *pszAdditionalMsg = nullptr, bool bShort = true)
         "               [-explodecollections] [-zfield field_name]\n"
         "               [-gcp ungeoref_x ungeoref_y georef_x georef_y "
         "[elevation]]* [-order n | -tps]\n"
-        "               [[-s_coord_epoch epoch] | [-t_coord_epoch epoch] | "
-        "[-a_coord_epoch epoch]]\n"
+        "               [-s_coord_epoch epoch] [-t_coord_epoch epoch] "
+        "[-a_coord_epoch epoch]\n"
         "               [-nomd] [-mo \"META-TAG=VALUE\"]* [-noNativeData]\n");
 
     if (bShort)

diff --git a/autotest/osr/osr_basic.py b/autotest/osr/osr_basic.py
@@ -2472,3 +2472,43 @@ def test_osr_basic_default_axis_mapping_strategy():
 
     assert crs1.GetAxisMappingStrategy() == osr.OAMS_TRADITIONAL_GIS_ORDER
     assert crs2.GetAxisMappingStrategy() == osr.OAMS_AUTHORITY_COMPLIANT
+
+
+###############################################################################
+
+
+@pytest.mark.require_proj(9, 4)
+def test_osr_basic_urn_coordinateMetadata():
+
+    srs = osr.SpatialReference()
+    assert (
+        srs.SetFromUserInput(
+            "urn:ogc:def:coordinateMetadata:NRCAN::NAD83_CSRS_1997_MTM7_HT2_1997"
+        )
+        == ogr.OGRERR_NONE
+    )
+    assert srs.GetName() == "NAD83(CSRS)v7 / MTM zone 7 + CGVD28 height"
+    assert srs.GetCoordinateEpoch() == 1997.0
+
+
+###############################################################################
+
+
+@pytest.mark.require_proj(9, 4)
+def test_osr_basic_set_from_user_input_EPSG_8254_at_something():
+
+    srs = osr.SpatialReference()
+    srs.SetFromUserInput("EPSG:8255 @ 2002.5")  # NAD83(CSRS)v7
+    assert srs.GetName() == "NAD83(CSRS)v7"
+    assert srs.GetCoordinateEpoch() == 2002.5
+
+
+###############################################################################
+
+
+@pytest.mark.require_proj(9, 4)
+def test_osr_basic_has_point_motion_operation():
+
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(8255)  # NAD83(CSRS)v7
+    assert srs.HasPointMotionOperation()
diff --git a/autotest/osr/osr_ct.py b/autotest/osr/osr_ct.py
@@ -768,3 +768,35 @@ def test_osr_ct_OGR_CT_PREFER_OFFICIAL_SRS_DEF():
         x, y, _ = ct.TransformPoint(826158.063, 2405844.125, 0)
         assert abs(x - 9.867) < 0.001, x
         assert abs(y - 71.125) < 0.001, y
+
+
+###############################################################################
+# Test NAD83(CSRS)v7 change of epoch
+
+
+@pytest.mark.require_proj(9, 4)
+def test_osr_ct_point_motion_operation():
+
+    s = osr.SpatialReference()
+    s.ImportFromEPSG(8254)  # NAD83(CSRS)v7 3D
+    s.SetCoordinateEpoch(2002)
+
+    t = osr.SpatialReference()
+    t.ImportFromEPSG(8254)  # NAD83(CSRS)v7 3D
+    t.SetCoordinateEpoch(2010)
+
+    ct = osr.CoordinateTransformation(s, t)
+    x, y, z = ct.TransformPoint(60.5, -79.5)
+    assert abs(x - 60.49999994) < 1e-8, x
+    assert abs(y - -79.49999963) < 1e-8, y
+    assert abs(z - 0.060) < 1e-3, z
+
+    t = osr.SpatialReference()
+    t.ImportFromEPSG(22717)  # "NAD83(CSRS)v7 / UTM zone 17N"
+    t.SetCoordinateEpoch(2010)
+
+    ct = osr.CoordinateTransformation(s, t)
+    x, y, z = ct.TransformPoint(60.5, -79.5)
+    assert abs(x - 582395.993) < 1e-3, x
+    assert abs(y - 6708035.973) < 1e-3, y
+    assert abs(z - 0.060) < 1e-3, z
diff --git a/autotest/proj_grids/README.txt b/autotest/proj_grids/README.txt
@@ -1,2 +1,3 @@
 * ``conus`` is the "official" NADCON4 conus NAD27->NAD84 file in CTable2 format
 * ``egm96_15_extract.gtx`` is an extract generated with ``gdal_translate /usr/share/proj/egm96_15.gtx ../autotest/proj_grids/egm96_15_extract.gtx -srcwin 234 233 3 3``
+* ``ca_nrc_NAD83v70VG.tif`` is an extract generated with ``gdal_translate ~/proj/PROJ-data/ca_nrc/ca_nrc_NAD83v70VG.tif ca_nrc_NAD83v70VG.tif -b 1 -b 2 -b 3 -co compress=deflate  -projwin -80 61 -79 60``
diff --git a/autotest/proj_grids/ca_nrc_NAD83v70VG.tif b/autotest/proj_grids/ca_nrc_NAD83v70VG.tif
diff --git a/autotest/utilities/test_gdaltransform.py b/autotest/utilities/test_gdaltransform.py
@@ -28,11 +28,14 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
+import sys
 
 import gdaltest
 import pytest
 import test_cli_utilities
 
+from osgeo import osr
+
 pytestmark = pytest.mark.skipif(
     test_cli_utilities.get_gdaltransform_path() is None,
     reason="gdaltransform not available",
@@ -225,3 +228,60 @@ def test_gdaltransform_ct_4D(gdaltransform_path):
     assert values[0] == pytest.approx(2.0000005420366, abs=1e-10), ret
     assert values[1] == pytest.approx(49.0000003766711, abs=1e-10), ret
     assert values[2] == pytest.approx(-0.0222802283242345, abs=1e-8), ret
+
+
+###############################################################################
+# Test s_coord_epoch
+
+
+@pytest.mark.require_proj(7, 2)
+def test_gdaltransform_s_coord_epoch(gdaltransform_path):
+
+    ret = gdaltest.runexternal(
+        gdaltransform_path
+        + " -s_srs EPSG:9000 -s_coord_epoch 2030 -t_srs EPSG:7844 -coord 150 -30"
+    )
+
+    values = [float(x) for x in ret.split(" ")]
+    assert len(values) == 3, ret
+    assert abs(values[0] - 150) > 1e-8, ret
+    assert abs(values[1] - -30) > 1e-8, ret
+
+
+###############################################################################
+# Test t_coord_epoch
+
+
+@pytest.mark.require_proj(7, 2)
+def test_gdaltransform_t_coord_epoch(gdaltransform_path):
+
+    ret = gdaltest.runexternal(
+        gdaltransform_path
+        + " -s_srs EPSG:7844 -t_srs EPSG:9000 -t_coord_epoch 2030 -coord 150 -30"
+    )
+
+    values = [float(x) for x in ret.split(" ")]
+    assert len(values) == 3, ret
+    assert abs(values[0] - 150) > 1e-8, ret
+    assert abs(values[1] - -30) > 1e-8, ret
+
+
+###############################################################################
+# Test s_coord_epoch and t_coord_epoch
+
+
+@pytest.mark.require_proj(9, 4)
+def test_gdaltransform_s_t_coord_epoch(gdaltransform_path):
+
+    sep = ";" if sys.platform == "win32" else ":"
+    PROJ_DATA = sep.join(osr.GetPROJSearchPaths())
+
+    ret = gdaltest.runexternal(
+        gdaltransform_path
+        + f' -s_srs EPSG:8254 -s_coord_epoch 2002 -t_srs EPSG:8254 -t_coord_epoch 2010 -coord -79.5 60.5 --config PROJ_DATA "{PROJ_DATA}"'
+    )
+
+    values = [float(x) for x in ret.split(" ")]
+    assert len(values) == 3, ret
+    assert abs(values[0] - -79.499999630188) < 1e-8, ret
+    assert abs(values[1] - 60.4999999378478) < 1e-8, ret
diff --git a/doc/source/programs/gdaltransform.rst b/doc/source/programs/gdaltransform.rst
@@ -17,6 +17,7 @@ Synopsis
 
     gdaltransform [--help] [--help-general]
         [-i] [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]
+        [-s_coord_epoch epoch] [-t_coord_epoch epoch]
         [-ct proj_string] [-order n] [-tps] [-rpc] [-geoloc]
         [-gcp pixel line easting northing [elevation]]* [-output_xy]
         [srcfile [dstfile]]
@@ -39,6 +40,17 @@ projection,including GCP-based transformations.
     (i.e. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prj file
     containing well known text.
 
+.. option:: -s_coord_epoch <epoch>
+
+    .. versionadded:: 3.8
+
+    Assign a coordinate epoch, linked with the source SRS. Useful when the
+    source SRS is a dynamic CRS. Only taken into account if :option:`-s_srs`
+    is used.
+
+    Before PROJ 9.4, :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
+    mutually exclusive, due to lack of support for transformations between two dynamic CRS.
+
 .. option:: -t_srs <srs_def>
 
     set target spatial reference.
@@ -47,6 +59,17 @@ projection,including GCP-based transformations.
     (i.e. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prj file
     containing well known text.
 
+.. option:: -t_coord_epoch <epoch>
+
+    .. versionadded:: 3.8
+
+    Assign a coordinate epoch, linked with the output SRS. Useful when the
+    output SRS is a dynamic CRS. Only taken into account if :option:`-t_srs`
+    is used.
+
+    Before PROJ 9.4, :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
+    mutually exclusive, due to lack of support for transformations between two dynamic CRS.
+
 .. option:: -ct <string>
 
     A PROJ string (single step operation or multiple step string

diff --git a/doc/source/programs/gdalwarp.rst b/doc/source/programs/gdalwarp.rst
@@ -19,7 +19,7 @@ Synopsis
         [-b|-srcband n]* [-dstband n]*
         [-s_srs srs_def] [-t_srs srs_def] [-ct string]
         [-to "NAME=VALUE"]* [-vshift | -novshift]
-        [[-s_coord_epoch epoch] | [-t_coord_epoch epoch]]
+        [-s_coord_epoch epoch] [-t_coord_epoch epoch]
         [-order n | -tps | -rpc | -geoloc] [-et err_threshold]
         [-refine_gcps tolerance [minimum_gcps]]
         [-te xmin ymin xmax ymax] [-te_srs srs_def]
@@ -118,7 +118,7 @@ with control information.
     source SRS is a dynamic CRS. Only taken into account if :option:`-s_srs`
     is used.
 
-    Currently :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
+    Before PROJ 9.4, :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
     mutually exclusive, due to lack of support for transformations between two dynamic CRS.
 
 .. option:: -t_srs <srs_def>
@@ -139,7 +139,7 @@ with control information.
     target SRS is a dynamic CRS. Only taken into account if :option:`-t_srs`
     is used.
 
-    Currently :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
+    Before PROJ 9.4, :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
     mutually exclusive, due to lack of support for transformations between two dynamic CRS.
 
 .. option:: -ct <string>

diff --git a/doc/source/programs/ogr2ogr.rst b/doc/source/programs/ogr2ogr.rst
@@ -49,7 +49,7 @@ Synopsis
             [-resolveDomains]
             [-explodecollections] [-zfield field_name]
             [-gcp ungeoref_x ungeoref_y georef_x georef_y [elevation]]* [-order n | -tps]
-            [[-s_coord_epoch epoch] | [-t_coord_epoch epoch] | [-a_coord_epoch epoch]]
+            [-s_coord_epoch epoch] [-t_coord_epoch epoch] [-a_coord_epoch epoch]
             [-nomd] [-mo "META-TAG=VALUE"]* [-noNativeData]
 
 Description
@@ -237,7 +237,7 @@ output coordinate system or even reprojecting the features during translation.
     output SRS is a dynamic CRS. Only taken into account if :option:`-t_srs`
     is used. It is also mutually exclusive with  :option:`-a_coord_epoch`.
 
-    Currently :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
+    Before PROJ 9.4, :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
     mutually exclusive, due to lack of support for transformations between two dynamic CRS.
 
 .. option:: -s_srs <srs_def>
@@ -256,7 +256,7 @@ output coordinate system or even reprojecting the features during translation.
     source SRS is a dynamic CRS. Only taken into account if :option:`-s_srs`
     is used.
 
-    Currently :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
+    Before PROJ 9.4, :option:`-s_coord_epoch` and :option:`-t_coord_epoch` are
     mutually exclusive, due to lack of support for transformations between two dynamic CRS.
 
 .. option:: -ct <string>

diff --git a/ogr/ogr_spatialref.h b/ogr/ogr_spatialref.h
@@ -375,6 +375,10 @@ class CPL_DLL OGRSpatialReference
     int IsProjected() const;
     int IsGeocentric() const;
     bool IsDynamic() const;
+
+    // cppcheck-suppress functionStatic
+    bool HasPointMotionOperation() const;
+
     int IsLocal() const;
     int IsVertical() const;
     int IsCompound() const;

diff --git a/ogr/ogr_srs_api.h b/ogr/ogr_srs_api.h
@@ -551,6 +551,7 @@ int CPL_DLL OSRIsCompound(OGRSpatialReferenceH);
 int CPL_DLL OSRIsGeocentric(OGRSpatialReferenceH);
 int CPL_DLL OSRIsVertical(OGRSpatialReferenceH);
 int CPL_DLL OSRIsDynamic(OGRSpatialReferenceH);
+int CPL_DLL OSRHasPointMotionOperation(OGRSpatialReferenceH);
 int CPL_DLL OSRIsSameGeogCS(OGRSpatialReferenceH, OGRSpatialReferenceH);
 int CPL_DLL OSRIsSameVertCS(OGRSpatialReferenceH, OGRSpatialReferenceH);
 int CPL_DLL OSRIsSame(OGRSpatialReferenceH, OGRSpatialReferenceH);