OGR_CT: handle point motion operations (PROJ >= 9.4)

Author: rouault

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

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``

autotest/proj_grids/ca_nrc_NAD83v70VG.tif

@@ -1470,24 +1470,35 @@ int OGRProjCT::Initialize(const OGRSpatialReference *poSourceIn,
         bSourceLatLong = CPL_TO_BOOL(poSRSSource->IsGeographic());
         bSourceIsDynamicCRS = poSRSSource->IsDynamic();
         dfSourceCoordinateEpoch = poSRSSource->GetCoordinateEpoch();
+        if (!bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0)
+        {
+            bSourceIsDynamicCRS = poSRSSource->HasPointMotionOperation();
+        }
         poSRSSource->GetAxis(nullptr, 0, &m_eSourceFirstAxisOrient);
     }
     if (poSRSTarget)
     {
         bTargetLatLong = CPL_TO_BOOL(poSRSTarget->IsGeographic());
         bTargetIsDynamicCRS = poSRSTarget->IsDynamic();
         dfTargetCoordinateEpoch = poSRSTarget->GetCoordinateEpoch();
+        if (!bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0)
+        {
+            bTargetIsDynamicCRS = poSRSTarget->HasPointMotionOperation();
+        }
         poSRSTarget->GetAxis(nullptr, 0, &m_eTargetFirstAxisOrient);
     }
 
+#if PROJ_VERSION_MAJOR < 9 ||                                                  \
+    (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR < 4)
     if (bSourceIsDynamicCRS && bTargetIsDynamicCRS &&
         dfSourceCoordinateEpoch > 0 && dfTargetCoordinateEpoch > 0 &&
         dfSourceCoordinateEpoch != dfTargetCoordinateEpoch)
     {
         CPLError(CE_Warning, CPLE_AppDefined,
-                 "Coordinate transformation between different epochs are "
-                 "not currently supported");
+                 "Coordinate transformation between different epochs only"
+                 "supported since PROJ 9.4");
     }
+#endif
 
     /* -------------------------------------------------------------------- */
     /*      Setup source and target translations to radians for lat/long    */
@@ -1659,6 +1670,24 @@ int OGRProjCT::Initialize(const OGRSpatialReference *poSourceIn,
                                         options.d->bOnlyBest ? "YES" : "NO");
             }
 #endif
+
+#if PROJ_VERSION_MAJOR > 9 ||                                                  \
+    (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR >= 4)
+            if (bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0 &&
+                bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0)
+            {
+                auto srcCM = proj_coordinate_metadata_create(
+                    ctx, srcCRS, dfSourceCoordinateEpoch);
+                proj_destroy(srcCRS);
+                srcCRS = srcCM;
+
+                auto targetCM = proj_coordinate_metadata_create(
+                    ctx, targetCRS, dfTargetCoordinateEpoch);
+                proj_destroy(targetCRS);
+                targetCRS = targetCM;
+            }
+#endif
+
             m_pj = proj_create_crs_to_crs_from_pj(ctx, srcCRS, targetCRS, area,
                                                   aosOptions.List());
             proj_destroy(srcCRS);
@@ -1694,7 +1723,9 @@ int OGRProjCT::Initialize(const OGRSpatialReference *poSourceIn,
 #endif
     }
 
-    if (options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget)
+    if (options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
+        (dfSourceCoordinateEpoch == 0 || dfTargetCoordinateEpoch == 0 ||
+         dfSourceCoordinateEpoch == dfTargetCoordinateEpoch))
     {
         // Determine if we can skip the transformation completely.
         const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",