Merge pull request #4688 from OSGeo/backport-4685-to-9.8

rouault · web-flow · commit c4be2f905b75 · 2026-03-02T17:52:03.000+01:00
[Backport 9.8] Make sure that epoch is set in more scenarios of time-dependent transformations
diff --git a/docs/source/operations/index.rst b/docs/source/operations/index.rst
@@ -20,3 +20,4 @@ exert a change in reference frame are called transformations.
    transformations/index
    pipeline
    operations_computation
+   time_dependent_transformations
diff --git a/docs/source/operations/time_dependent_transformations.rst b/docs/source/operations/time_dependent_transformations.rst
@@ -0,0 +1,188 @@
+.. _time_dependent_transformations:
+
+================================================================================
+Time-dependent transformations
+================================================================================
+
+Introduction
+------------
+
+PROJ offers several coordinate operation methods that can deal with time-dependent
+transformations:
+
+* 15-parameter :ref:`Helmert <helmert>` operation, with the translation rates
+  (``dx``, ``dy``, ``dz``), the scale rate (``dz``),
+  the rotation rates (``drx``, ``dry``, ``drz``) and the central epoch of
+  transformation (``t_epoch``). Such transformations are typically used between
+  a static CRS (bound to a tectonic plate) and a dynamic CRS (such as one of the
+  ITRF or WGS 84 realizations).
+
+  To transform Sydney coordinates at epoch 2026.0 from ITRF2014 to GDA2020,
+  you can use for example:
+
+  ::
+
+        echo -33.8623 151.2077 0 2026 | cs2cs -d 8 ITRF2014 GDA2020
+
+        -33.86230294 151.20769882 0.00101330 2026
+
+
+  or equivalently using the PROJ pipeline returned by ``projinfo ITRF2014 GDA2020``:
+
+  ::
+
+        echo -33.8623 151.2077 0 2026 | cct -d 8 \
+          +proj=pipeline \
+          +step +proj=axisswap +order=2,1 \
+          +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+          +step +proj=cart +ellps=GRS80 \
+          +step +proj=helmert +x=0 +y=0 +z=0 +rx=0 +ry=0 +rz=0 +s=0 +dx=0 +dy=0 +dz=0 \
+                +drx=0.00150379 +dry=0.00118346 +drz=0.00120716 +ds=0 +t_epoch=2020 \
+                +convention=coordinate_frame \
+          +step +inv +proj=cart +ellps=GRS80 \
+          +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+          +step +proj=axisswap +order=2,1
+
+        -33.86230294    151.20769882    0.00101330     2026.0000
+
+
+  and the reverse transformation:
+
+  ::
+
+        echo -33.86230294 151.20769882 0.00101330 2026 | cs2cs -d 8 GDA2020 ITRF2014
+
+        -33.86230000 151.20770000 0.00000002 2026
+
+  Note that in the reverse transformation the epoch is also specified as seemingly
+  being an input information, whereas it is the target epoch for ITRF2014. This
+  is an historic PROJ behavior before the introduction of a more generic framework
+  to deal with coordinate epochs. Since PROJ 9.4, you can also qualify a dynamic
+  CRS with a coordinate epoch using the syntax "<CRS name>@<epoch>" or "<CRS code>@<epoch>"
+
+  For example, re-using above use cases 
+
+  ::
+
+        echo -33.8623 151.2077 | cs2cs -d 8 ITRF2014@2026.0 GDA2020
+
+        -33.86230294 151.20769882 0.00101330
+
+  ::
+
+        echo -33.86230294 151.20769882 0.00101330 | cs2cs -d 8 EPSG:7844 EPSG:9000@2026.0
+
+        -33.86230000 151.20770000 0.00000002
+
+
+* Point motion operation between two epochs of a dynamic CRS using a velocity grid-based
+  transformation, with the :ref:`deformation` operator.
+
+  For example, to transform coordinates in NAD83(CSRS)v7 CRS from epoch 2010.0 to epoch 2002.0
+
+  ::
+
+        echo 60 -100 0 | PROJ_NETWORK=ON cs2cs -d 8 "NAD83(CSRS)v7 @ 2010.0" "NAD83(CSRS)v7 @ 2002.0"
+
+        60.00000011 -100.00000045 -0.06764700
+
+  or equivalently using the PROJ pipeline returned by
+  ``projinfo "NAD83(CSRS)v7 @ 2010.0" "NAD83(CSRS)v7 @ 2002.0" --spatial-test intersects``
+  
+  ::
+
+        echo 60 -100 0 | PROJ_NETWORK=ON  cct -d 8 \
+          +proj=pipeline \
+          +step +proj=set +v_4=2010 \
+          +step +proj=axisswap +order=2,1 \
+          +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m \
+          +step +proj=cart +ellps=GRS80 \
+          +step +proj=set +v_4=2010 +omit_fwd \
+          +step +proj=deformation +dt=-8 +grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80 \
+          +step +proj=set +v_4=2002 +omit_inv \
+          +step +inv +proj=cart +ellps=GRS80 \
+          +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+          +step +proj=axisswap +order=2,1 \
+          +step +proj=set +v_4=2002
+
+        60.00000011 -100.00000045 -0.06764700 2002.0000
+
+
+* Using a multi-component time-based deformation model with the :ref:`defmodel` operator.
+
+  For example, to transform coordinates from ITRF96 at epoch 2026.0 to NZGD2000
+
+  ::
+
+        $ echo -41.2924 174.7787 0 2026 | PROJ_NETWORK=ON cs2cs -d 8 ITRF96 NZGD2000
+
+        -41.29240855 174.77870601 0.00000000 2026
+
+  or (since PROJ 9.9)
+
+  ::
+
+        $ echo -41.2924 174.7787 0 | PROJ_NETWORK=ON cs2cs -d 8 ITRF96@2026 NZGD2000
+
+        -41.29240855 174.77870601 0.00000000
+
+
+  or equivalently using the PROJ pipeline returned by
+  ``projinfo ITRF96@2026 NZGD2000 --spatial-test intersects``
+  
+  ::
+
+        echo -41.2924 174.7787 0 | PROJ_NETWORK=ON cct -d 8 \
+          +proj=pipeline \
+          +step +proj=axisswap +order=2,1 \
+          +step +proj=unitconvert +xy_in=deg +xy_out=rad \
+          +step +proj=set +v_4=2026 \
+          +step +inv +proj=defmodel +model=nz_linz_nzgd2000-20180701.json \
+          +step +proj=set +v_4=2026 \
+          +step +proj=unitconvert +xy_in=rad +xy_out=deg \
+          +step +proj=axisswap +order=2,1
+
+        -41.29240855    174.77870601    0.00000000     2026.0000
+
+
+Other types of time-dependent transformations are not yet supported when
+PROJ must infer the transformation pipelines:
+
+- time-specific Helmert transformations
+- transformations between two dynamic CRS at different epochs
+- plate motion models (as available in :source_file:`data/ITRF2008`,
+  :source_file:`data/ITRF2014` and :source_file:`data/ITRF2020`).
+
+
+Through C API
+-------------
+
+Since PROJ 9.9, the most convenient way, that will work with all above mentioned
+coordinate operations, is to attach an epoch to a dynamic CRS,
+forming a :cpp:class:`osgeo::proj::coordinates::CoordinateMetadata` instance,
+either by passing a string of the form "<CRS string> @ <epoch>" to :cpp:func:`proj_create`
+or :cpp:func:`proj_create_crs_to_crs`. Or when operating on ``PJ*`` instance, to
+use :cpp:func:`proj_coordinate_metadata_create` from an existing CRS instance,
+and pass it to :cpp:func:`proj_create_crs_to_crs_from_pj` or :cpp:func:`proj_create_operations`
+
+.. code-block:: c
+
+    PJ_CONTEXT *ctxt = proj_context_create();
+    PJ *transf = proj_create_crs_to_crs(C, "ITRF2014@2026.0", "GDA2020", nullptr);
+    PROJ_COORD coord;
+    coord.xyz.y = -33.8623;  // latitude in degree
+    coord.xyz.y = 151.2077;  // longitude in degree
+    coord.xyz.z = 0;         // ellipsoidal height in meter
+    PROJ_COORD coord_in_gda2020 = proj_trans(transf, PJ_FWD, coord);
+    proj_destroy(transf);
+    proj_context_destroy(ctxt);
+
+
+When using the :cpp:func:`proj_create_operations` functionality,
+:cpp:func:`proj_coordoperation_requires_per_coordinate_input_time` can be used on
+the returned operations to check if they have a time dependency.
+
+Note that with PROJ 9.3 or earlier for 15-parameter Helmert or, PROJ 9.8 or
+earlier, for multi-component time-based deformation model transformations, you
+must set the epoch in the ``PROJ_COORD::xyzt::t`` or ``PROJ_COORD::v[3]`` member
+since CoordinateMetadata is not available.
diff --git a/src/apps/projapps_lib.h b/src/apps/projapps_lib.h
@@ -55,7 +55,7 @@ typedef void (*projinfo_cb_t)(PJ_PROJINFO_LOG_LEVEL level, const char *msg,
  * @param argv list of char* with the command parameters of projinfo.
  *  It does not contain the program name as first parameter.
  * @param cb callback that to get the output of projinfo.
- *  It can be very fragmented, no necesarily by lines.
+ *  It can be very fragmented, no necessarily by lines.
  * @param user_data pointer for data passed to the callback.
  */
 int PROJ_DLL projinfo(PJ_CONTEXT *ctx, int argc, char **argv, projinfo_cb_t cb,
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
@@ -8930,6 +8930,11 @@ const std::string &PROJStringFormatter::toString() const {
             step.inverted = false;
             continue;
         }
+
+        // set does the same in forward and inverse paths
+        if (step.name == "set") {
+            step.inverted = false;
+        }
     }
 
     {
@@ -9398,6 +9403,34 @@ const std::string &PROJStringFormatter::toString() const {
                 }
             }
 
+            // "+proj=set +v_4=X" followed by "+proj=set +v_4=X +omit_fwd"
+            // can be optimized as "+proj=set +v_4=X"
+            if (curStep.name == "set" && prevStep.name == "set" &&
+                !curStep.inverted && !prevStep.inverted &&
+                curStepParamCount == 2 && prevStepParamCount == 1 &&
+                curStep.paramValues[0].keyEquals("v_4") &&
+                prevStep.paramValues[0].keyEquals("v_4") &&
+                curStep.paramValues[1].keyEquals("omit_fwd") &&
+                curStep.paramValues[0].value == prevStep.paramValues[0].value) {
+
+                iterCur = steps.erase(iterCur);
+                continue;
+            }
+
+            // "+proj=set +v_4=X +omit_inv" followed by "+proj=set +v_4=X"
+            // can be optimized as "+proj=set +v_4=X"
+            if (curStep.name == "set" && prevStep.name == "set" &&
+                !curStep.inverted && !prevStep.inverted &&
+                curStepParamCount == 1 && prevStepParamCount == 2 &&
+                curStep.paramValues[0].keyEquals("v_4") &&
+                prevStep.paramValues[0].keyEquals("v_4") &&
+                prevStep.paramValues[1].keyEquals("omit_inv") &&
+                curStep.paramValues[0].value == prevStep.paramValues[0].value) {
+
+                deletePrevIter();
+                continue;
+            }
+
             // for practical purposes WGS84 and GRS80 ellipsoids are
             // equivalents (cartesian transform between both lead to differences
             // of the order of 1e-14 deg..).
diff --git a/src/iso19111/operation/singleoperation.cpp b/src/iso19111/operation/singleoperation.cpp
@@ -301,6 +301,12 @@ void CoordinateOperation::setSourceCoordinateEpoch(
     const util::optional<common::DataEpoch> &epoch) {
     d->sourceCoordinateEpoch_ =
         std::make_shared<util::optional<common::DataEpoch>>(epoch);
+
+    if (epoch) {
+        auto invOp = dynamic_cast<InverseCoordinateOperation *>(this);
+        if (invOp)
+            invOp->forwardOperation()->setTargetCoordinateEpoch(epoch);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -309,6 +315,12 @@ void CoordinateOperation::setTargetCoordinateEpoch(
     const util::optional<common::DataEpoch> &epoch) {
     d->targetCoordinateEpoch_ =
         std::make_shared<util::optional<common::DataEpoch>>(epoch);
+
+    if (epoch) {
+        auto invOp = dynamic_cast<InverseCoordinateOperation *>(this);
+        if (invOp)
+            invOp->forwardOperation()->setSourceCoordinateEpoch(epoch);
+    }
 }
 
 // ---------------------------------------------------------------------------
@@ -3548,6 +3560,27 @@ bool SingleOperation::exportToPROJStringGeneric(
                                        addPushPopV3, "Helmert");
         }
 
+        double sourceYear =
+            sourceCoordinateEpoch().has_value()
+                ? getRoundedEpochInDecimalYear(
+                      sourceCoordinateEpoch()->coordinateEpoch().convertToUnit(
+                          common::UnitOfMeasure::YEAR))
+                : 0;
+        double targetYear =
+            targetCoordinateEpoch().has_value()
+                ? getRoundedEpochInDecimalYear(
+                      targetCoordinateEpoch()->coordinateEpoch().convertToUnit(
+                          common::UnitOfMeasure::YEAR))
+                : 0;
+        if (sourceYear > 0 && targetYear == 0)
+            targetYear = sourceYear;
+        else if (targetYear > 0 && sourceYear == 0)
+            sourceYear = targetYear;
+        if (sourceYear > 0) {
+            formatter->addStep("set");
+            formatter->addParam("v_4", sourceYear);
+        }
+
         formatter->addStep("helmert");
         if (fullMatrix)
             formatter->addParam("exact");
@@ -3612,6 +3645,11 @@ bool SingleOperation::exportToPROJStringGeneric(
             }
         }
 
+        if (targetYear > 0) {
+            formatter->addStep("set");
+            formatter->addParam("v_4", targetYear);
+        }
+
         if (l_targetCRS) {
             setupPROJGeodeticTargetCRS(formatter, NN_NO_CHECK(l_targetCRS),
                                        addPushPopV3, "Helmert");
@@ -4733,6 +4771,29 @@ bool SingleOperation::exportToPROJStringGeneric(
             sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
             formatter->stopInversion();
 
+            double sourceYear =
+                sourceCoordinateEpoch().has_value()
+                    ? getRoundedEpochInDecimalYear(
+                          sourceCoordinateEpoch()
+                              ->coordinateEpoch()
+                              .convertToUnit(common::UnitOfMeasure::YEAR))
+                    : 0;
+            double targetYear =
+                targetCoordinateEpoch().has_value()
+                    ? getRoundedEpochInDecimalYear(
+                          targetCoordinateEpoch()
+                              ->coordinateEpoch()
+                              .convertToUnit(common::UnitOfMeasure::YEAR))
+                    : 0;
+            if (sourceYear > 0 && targetYear == 0)
+                targetYear = sourceYear;
+            else if (targetYear > 0 && sourceYear == 0)
+                sourceYear = targetYear;
+            if (sourceYear > 0) {
+                formatter->addStep("set");
+                formatter->addParam("v_4", sourceYear);
+            }
+
             if (isMethodInverseOf) {
                 formatter->startInversion();
             }
@@ -4748,6 +4809,11 @@ bool SingleOperation::exportToPROJStringGeneric(
                 formatter->stopInversion();
             }
 
+            if (targetYear > 0) {
+                formatter->addStep("set");
+                formatter->addParam("v_4", targetYear);
+            }
+
             targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
 
             return true;
@@ -5200,6 +5266,12 @@ void InverseCoordinateOperation::setPropertiesFromForward() {
         forwardOperation_->hasBallparkTransformation());
     setRequiresPerCoordinateInputTime(
         forwardOperation_->requiresPerCoordinateInputTime());
+    if (auto sourceEpoch = forwardOperation_->sourceCoordinateEpoch()) {
+        setTargetCoordinateEpoch(sourceEpoch);
+    }
+    if (auto targetEpoch = forwardOperation_->targetCoordinateEpoch()) {
+        setSourceCoordinateEpoch(targetEpoch);
+    }
 }
 
 // ---------------------------------------------------------------------------
diff --git a/src/networkfilemanager.cpp b/src/networkfilemanager.cpp
@@ -2028,7 +2028,7 @@ static size_t pj_emscripten_read_range(PJ_CONTEXT *ctx,
     // Documentation says compiling with pthread flag is needed.
     // https://emscripten.org/docs/api_reference/fetch.html#synchronous-fetches
     // We encapsulate the code related to empscripten_fetch in a lamda.
-    // Some tests runing this lambda in a thread were partially successful.
+    // Some tests running this lambda in a thread were partially successful.
     size_t real_read = 0;
     const std::string url = handle->m_url;
     auto fetching = [&]() {
diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp
