Merge pull request #4038 from rouault/gridshift_projected

+proj=gridshift: enhance to support grids referenced in projected CRS, and with easting_offset/northing_offset corrections

docs/source/operations/transformations/gridshift.rst

@@ -13,19 +13,21 @@ Translation of geodetic coordinates using a grid shift.
 +-----------------+-------------------------------------------------------------------+
 | **Domain**      | 2D and 3D                                                         |
 +-----------------+-------------------------------------------------------------------+
-| **Input type**  | Geodetic coordinates (horizontal), meters (vertical)              |
+| **Input type**  | Geodetic or projected coordinates (horizontal), meters (vertical) |
 +-----------------+-------------------------------------------------------------------+
-| **Output type** | Geodetic coordinates (horizontal), meters (vertical)              |
+| **Output type** | Geodetic or projected coordinates (horizontal), meters (vertical) |
 +-----------------+-------------------------------------------------------------------+
 
-The transformation may apply horizontal geodetic offsetting and/or vertical
+The transformation may apply horizontal geodetic or projected 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.
 
+.. note:: Support for grids referenced in a projected CRS has been added in PROJ 9.4.0
+
 
 Example
 -------------------------------------------------------------------------------

docs/source/specifications/geodetictiffgrids.rst

@@ -212,10 +212,16 @@ is an easy way to inspect such grid files:
     Values recognized by PROJ currently are:
 
     - ``HORIZONTAL_OFFSET``: implies the presence of at least two samples.
-      The first sample must contain the latitude offset and the second
-      sample must contain the longitude offset. The offset may also be expressed
-      as a speed per year for temporal gridshifting.
-      Corresponds to PROJ :ref:`hgridshift` method.
+      For grids referenced in geographic coordinates, the first sample must
+      contain the latitude offset and the second sample must contain the
+      longitude offset.
+      For grids referenced in projected coordinates (supported since PROJ 9.4),
+      the first sample must contain the easting offset and the second sample
+      must contain the northing offset.
+      The offset may also be expressed as a speed per year for temporal gridshifting.
+      Corresponds to PROJ :ref:`hgridshift` (only for grids referenced in
+      geographic coordinates) and :ref:`gridshift` methods (for grids referenced
+      both in geographic and projected coordinates)
 
     - ``GEOGRAPHIC_3D_OFFSET``: implies the presence of at least 3 samples.
       The first sample must contain the latitude offset, the second
@@ -278,6 +284,14 @@ is an easy way to inspect such grid files:
       the value to add a longitude expressed in the CRS encoded in the GeoKeys
       to obtain a longitude value expressed in the target CRS.
 
+    + ``easting_offset``: valid for TYPE=HORIZONTAL_OFFSET. Sample values should be
+      the value to add a easting expressed in the CRS encoded in the GeoKeys
+      to obtain a easting value expressed in the target CRS.
+
+    + ``northing_offset``: valid for TYPE=HORIZONTAL_OFFSET. Sample values should be
+      the value to add a northing expressed in the CRS encoded in the GeoKeys
+      to obtain a northing 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
@@ -383,6 +397,19 @@ is an easy way to inspect such grid files:
         <Item name="UNITTYPE" sample="0" role="unittype">arc-second</Item>
         <Item name="UNITTYPE" sample="1" role="unittype">arc-second</Item>
 
+  * For grids with TYPE=HORIZONTAL_OFFSET and with a ``easting_offset`` and
+    ``northing_offset`` channel, an extra offset to apply after the grid correction
+    in the forward direction of the transformation can be specified with a
+    ``constant_offset`` metadata item, expressed in the same units as the grid
+    values (only metre supported at the moment).
+
+    Example:
+
+    .. code-block:: xml
+
+        <Item name="constant_offset" sample="0">-5000000</Item>
+        <Item name="constant_offset" sample="1">-5000000</Item>
+
   * For TYPE=DEFORMATION_MODEL, the type of the displacement must be specified
     with a `Item` whose ``name`` is set to ``DISPLACEMENT_TYPE``.
 

scripts/reference_exported_symbols.txt

@@ -322,7 +322,7 @@ osgeo::proj::GenericShiftGridSet::gridAt(std::string const&, double, double) con
 osgeo::proj::GenericShiftGridSet::open(pj_ctx*, std::string const&)
 osgeo::proj::GenericShiftGridSet::reassign_context(pj_ctx*)
 osgeo::proj::GenericShiftGridSet::reopen(pj_ctx*)
-osgeo::proj::GenericShiftGrid::valuesAt(int, int, int, int, int, int const*, float*) const
+osgeo::proj::GenericShiftGrid::valuesAt(int, int, int, int, int, int const*, float*, bool&) const
 osgeo::proj::Grid::~Grid()
 osgeo::proj::Grid::Grid(std::string const&, int, int, osgeo::proj::ExtentAndRes const&)
 osgeo::proj::HorizontalShiftGrid::gridAt(double, double) const

src/grids.cpp

@@ -484,7 +484,8 @@ class GTiffGrid : public Grid {
     bool valueAt(uint16_t sample, int x, int y, float &out) const;
 
     bool valuesAt(int x_start, int y_start, int x_count, int y_count,
-                  int sample_count, const int *sample_idx, float *out) const;
+                  int sample_count, const int *sample_idx, float *out,
+                  bool &nodataFound) const;
 
     bool isNodata(float val) const;
 
@@ -769,11 +770,12 @@ bool GTiffGrid::valueAt(uint16_t sample, int x, int yFromBottom,
 // ---------------------------------------------------------------------------
 
 bool GTiffGrid::valuesAt(int x_start, int y_start, int x_count, int y_count,
-                         int sample_count, const int *sample_idx,
-                         float *out) const {
+                         int sample_count, const int *sample_idx, float *out,
+                         bool &nodataFound) const {
     const auto getTIFFRow = [this](int y) {
         return m_bottomUp ? y : m_height - 1 - y;
     };
+    nodataFound = false;
     if (m_blockIs256Pixel && m_planarConfig == PLANARCONFIG_CONTIG &&
         m_dt == TIFFDataType::Float32 &&
         (x_start / 256) == (x_start + x_count - 1) / 256 &&
@@ -915,6 +917,9 @@ bool GTiffGrid::valuesAt(int x_start, int y_start, int x_count, int y_count,
                 if (!valueAt(static_cast<uint16_t>(sample_idx[isample]), x, y,
                              *out))
                     return false;
+                if (isNodata(*out)) {
+                    nodataFound = true;
+                }
                 ++out;
             }
         }
@@ -2845,8 +2850,8 @@ class GTiffGenericGrid final : public GenericShiftGrid {
     bool valueAt(int x, int y, int sample, float &out) const override;
 
     bool valuesAt(int x_start, int y_start, int x_count, int y_count,
-                  int sample_count, const int *sample_idx,
-                  float *out) const override;
+                  int sample_count, const int *sample_idx, float *out,
+                  bool &nodataFound) const override;
 
     int samplesPerPixel() const override { return m_grid->samplesPerPixel(); }
 
@@ -2917,9 +2922,10 @@ bool GTiffGenericGrid::valueAt(int x, int y, int sample, float &out) const {
 
 bool GTiffGenericGrid::valuesAt(int x_start, int y_start, int x_count,
                                 int y_count, int sample_count,
-                                const int *sample_idx, float *out) const {
+                                const int *sample_idx, float *out,
+                                bool &nodataFound) const {
     return m_grid->valuesAt(x_start, y_start, x_count, y_count, sample_count,
-                            sample_idx, out);
+                            sample_idx, out, nodataFound);
 }
 
 // ---------------------------------------------------------------------------
@@ -3047,7 +3053,9 @@ GenericShiftGrid::~GenericShiftGrid() = default;
 
 bool GenericShiftGrid::valuesAt(int x_start, int y_start, int x_count,
                                 int y_count, int sample_count,
-                                const int *sample_idx, float *out) const {
+                                const int *sample_idx, float *out,
+                                bool &nodataFound) const {
+    nodataFound = false;
     for (int y = y_start; y < y_start + y_count; ++y) {
         for (int x = x_start; x < x_start + x_count; ++x) {
             for (int isample = 0; isample < sample_count; ++isample) {

src/grids.hpp

@@ -215,8 +215,8 @@ class PROJ_GCC_DLL GenericShiftGrid : public Grid {
 
     PROJ_FOR_TEST virtual bool valuesAt(int x_start, int y_start, int x_count,
                                         int y_count, int sample_count,
-                                        const int *sample_idx,
-                                        float *out) const;
+                                        const int *sample_idx, float *out,
+                                        bool &nodataFound) const;
 
     PROJ_FOR_TEST virtual void reassign_context(PJ_CONTEXT *ctx) = 0;
 };