Implement 'Vertical Offset by velocity grid (NRCan byn)' method

Author: rouault

data/sql/customizations.sql

@@ -188,12 +188,19 @@ INSERT INTO usage SELECT
     scope_code
     FROM usage WHERE object_table_name = 'grid_transformation' AND object_auth_name = 'EPSG' AND object_code = '9123';
 
--- Duplicate EPSG:10015 with NAD83(CSRS)v7 as interpolation CRS
+-- Duplicate EPSG:10115, 10116, 10117 with NAD83(CSRS)v7 as interpolation CRS
 INSERT INTO "grid_transformation" VALUES
 ('PROJ','EPSG_10115_WITH_NAD83CSRSV7_INTERPOLATION','CGVD28 height to CGVD2013a(2010) height (1)','Equivalent to HT2_2010v70 hybrid geoid model minus CGG2013a geoid model (i.e. CT code 9986 minus CT 9247 = CT 9987 minus CT 10109). Specifies NAD83(CSRS)v7 (code 8255) as interpolation CRS.','EPSG','1112','Vertical Offset by Grid Interpolation (NRCan byn)','EPSG','5713','EPSG','9245',0.05,'EPSG','8732','Vertical offset file','HT2_2010v70_CGG2013a.byn',NULL,NULL,NULL,NULL,'EPSG','8255','NR-Can HT2 2010',0);
 
 INSERT INTO "usage" VALUES('PROJ','USAGE_EPSG_10115_WITH_NAD83CSRSV7_INTERPOLATION','grid_transformation','PROJ','EPSG_10115_WITH_NAD83CSRSV7_INTERPOLATION','EPSG','1061','EPSG','1133');
 
+INSERT INTO "grid_transformation" VALUES('PROJ','EPSG_10116_WITH_NAD83CSRSV7_INTERPOLATION','CGVD2013a(2010) height to CGVD2013a(2002) height (1)',' Specifies NAD83(CSRS)v7 (code 8255) as interpolation CRS. Equivalent to CGVD28 to CGVD2013a(2002) minus CGVD28 to CGVD2013a(2010) (i.e. CT code 10114 minus CT 10115).','EPSG','1113','Vertical Offset by velocity grid (NRCan byn)','EPSG','9245','EPSG','20034',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8255','EPSG-Can',0);
+INSERT INTO "usage" VALUES('PROJ','USAGE_EPSG_10116_WITH_NAD83CSRSV7_INTERPOLATION','grid_transformation','PROJ','EPSG_10116_WITH_NAD83CSRSV7_INTERPOLATION','EPSG','1061','EPSG','1026');
+
+INSERT INTO "grid_transformation" VALUES('PROJ','EPSG_10117_WITH_NAD83CSRSV7_INTERPOLATION','CGVD2013a(2010) height to CGVD2013a(1997) height (1)','Specifies NAD83(CSRS)v7 (code 8255) as interpolation CRS. Equivalent to CGVD28 to CGVD2013a(1997) minus CGVD28 to CGVD2013a(2010) (i.e. CT code 10113 minus CT 10115).','EPSG','1113','Vertical Offset by velocity grid (NRCan byn)','EPSG','9245','EPSG','20035',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8255','EPSG-Can',0);
+INSERT INTO "usage" VALUES('PROJ','USAGE_EPSG_10117_WITH_NAD83CSRSV7_INTERPOLATION','grid_transformation','PROJ','EPSG_10117_WITH_NAD83CSRSV7_INTERPOLATION','EPSG','1061','EPSG','1026');
+
+
 -- Define the allowed authorities, and their precedence, when researching a
 -- coordinate operation
 

data/sql/grid_transformation.sql

@@ -1440,6 +1440,12 @@ INSERT INTO "grid_transformation" VALUES('EPSG','10114','CGVD28 height to CGVD20
 INSERT INTO "usage" VALUES('EPSG','18322','grid_transformation','EPSG','10114','EPSG','1061','EPSG','1133');
 INSERT INTO "grid_transformation" VALUES('EPSG','10115','CGVD28 height to CGVD2013a(2010) height (1)','Equivalent to HT2_2010v70 hybrid geoid model minus CGG2013a geoid model (i.e. CT code 9986 minus CT 9247 = CT 9987 minus CT 10109). Specifies NAD83(CSRS)v6 as interpolation CRS, but NAD83(CSRS)v7 (code 8255) may equally be used.','EPSG','1112','Vertical Offset by Grid Interpolation (NRCan byn)','EPSG','5713','EPSG','9245',0.05,'EPSG','8732','Vertical offset file','HT2_2010v70_CGG2013a.byn',NULL,NULL,NULL,NULL,'EPSG','8252','NR-Can HT2 2010',0);
 INSERT INTO "usage" VALUES('EPSG','18323','grid_transformation','EPSG','10115','EPSG','1061','EPSG','1133');
+INSERT INTO "grid_transformation" VALUES('EPSG','10116','CGVD2013a(2010) height to CGVD2013a(2002) height (1)','Interpolation CRS = NAD83(CSRS) 2010, i.e. may be any one of NAD83(CSRS)v6, v7 or v8 (CRS codes 8252, 8255 or 10414). Equivalent to CGVD28 to CGVD2013a(2002) minus CGVD28 to CGVD2013a(2010) (i.e. CT code 10114 minus CT 10115).','EPSG','1113','Vertical Offset by velocity grid (NRCan byn)','EPSG','9245','EPSG','20034',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8252','EPSG-Can',0);
+INSERT INTO "usage" VALUES('EPSG','18325','grid_transformation','EPSG','10116','EPSG','1061','EPSG','1026');
+INSERT INTO "grid_transformation" VALUES('EPSG','10117','CGVD2013a(2010) height to CGVD2013a(1997) height (1)','Interpolation CRS = NAD83(CSRS) 2010, i.e. may be any one of NAD83(CSRS)v6, v7 or v8 (CRS codes 8252, 8255 or 10414). Equivalent to CGVD28 to CGVD2013a(1997) minus CGVD28 to CGVD2013a(2010) (i.e. CT code 10113 minus CT 10115).','EPSG','1113','Vertical Offset by velocity grid (NRCan byn)','EPSG','9245','EPSG','20035',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8252','EPSG-Can',0);
+INSERT INTO "usage" VALUES('EPSG','18352','grid_transformation','EPSG','10117','EPSG','1061','EPSG','1026');
+INSERT INTO "grid_transformation" VALUES('EPSG','10118','CGVD2013a(2002) height to CGVD2013a(1997) height (1)','Equivalent to CGVD28 to CGVD2013a(1997) minus CGVD28 to CGVD2013a(2002) (i.e. CT code 10113 minus CT 10114).','EPSG','1113','Vertical Offset by velocity grid (NRCan byn)','EPSG','20034','EPSG','20035',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8246','EPSG-Can',0);
+INSERT INTO "usage" VALUES('EPSG','18328','grid_transformation','EPSG','10118','EPSG','1061','EPSG','1026');
 INSERT INTO "grid_transformation" VALUES('EPSG','10119','NAD83(CSRS)v4 to NAD83(CSRS)v2 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8244','EPSG','8235',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8246','EPSG-Can cvg70',0);
 INSERT INTO "usage" VALUES('EPSG','18201','grid_transformation','EPSG','10119','EPSG','1061','EPSG','1026');
 INSERT INTO "grid_transformation" VALUES('EPSG','10120','NAD83(CSRS)v4 to NAD83(CSRS)v3 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8244','EPSG','8239',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8246','EPSG-Can cvg70',0);

scripts/build_db.py

@@ -692,7 +692,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)', 'NADCON5 (3D)', 'NTv1', 'NTv2', 'VERTCON', 'Geocentric translation by Grid Interpolation (IGN)', 'Geographic3D Offset by velocity grid (NRCan byn)'))")
+    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)', 'Geographic3D Offset by velocity grid (NRCan byn)', 'Vertical Offset by velocity grid (NRCan byn)'))")
     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
@@ -775,13 +775,14 @@ def fill_grid_transformation(proj_db_cursor):
         # 1105: Geog3D to Geog2D+GravityRelatedHeight (ITAL2005)
         # 1110: Geog3D to Geog2D+Depth (Gravsoft)
         # 1112: Vertical Offset by Grid Interpolation (NRCan byn)
+        # 1113: Vertical Offset by velocity grid (NRCan byn)
         # 1114: Geographic3D Offset by velocity grid (NRCan byn)
         # 1115: Geog3D to Geog2D+Depth (txt)
         # 1118: Geog3D to Geog2D+GravityRelatedHeight (ISG)
         # 1122: Geog3D to Geog2D+Depth (gtx)
         # WARNING: update Transformation::isGeographic3DToGravityRelatedHeight()
         # in src/iso19111/operation/singleoperation.cpp if adding new methods
-        elif method_code in (1071, 1080, 1081, 1083, 1084, 1085, 1088, 1089, 1090, 1091, 1092, 1093, 1094, 1095, 1096, 1097, 1098, 1100, 1101, 1103, 1105, 1110, 1112, 1114, 1115, 1118, 1122) and n_params == 2:
+        elif method_code in (1071, 1080, 1081, 1083, 1084, 1085, 1088, 1089, 1090, 1091, 1092, 1093, 1094, 1095, 1096, 1097, 1098, 1100, 1101, 1103, 1105, 1110, 1112, 1113, 1114, 1115, 1118, 1122) and n_params == 2:
             assert param_code[1] == 1048, (code, method_code, param_code[1])
             interpolation_crs_auth_name = EPSG_AUTHORITY
             interpolation_crs_code = str(int(param_value[1])) # needed to avoid codes like XXXX.0

src/iso19111/operation/singleoperation.cpp

@@ -1991,6 +1991,34 @@ _getGeographic3DOffsetByVelocityGridFilename(const SingleOperation *op,
 
 // ---------------------------------------------------------------------------
 
+//! @cond Doxygen_Suppress
+static const std::string &
+_getVerticalOffsetByVelocityGridFilename(const SingleOperation *op,
+                                         bool allowInverse) {
+
+    const auto &l_method = op->method();
+    const auto &methodName = l_method->nameStr();
+    if (l_method->getEPSGCode() ==
+            EPSG_CODE_METHOD_VERTICAL_OFFSET_BY_VELOCITY_GRID_NRCAN ||
+        (allowInverse &&
+         ci_equal(
+             methodName,
+             INVERSE_OF +
+                 EPSG_NAME_METHOD_GEOGRAPHIC3D_OFFSET_BY_VELOCITY_GRID_NRCAN))) {
+        const auto &fileParameter = op->parameterValue(
+            EPSG_NAME_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE,
+            EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE);
+        if (fileParameter &&
+            fileParameter->type() == ParameterValue::Type::FILENAME) {
+            return fileParameter->valueFile();
+        }
+    }
+    return nullString;
+}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
 //! @cond Doxygen_Suppress
 static const std::string &
 _getHeightToGeographic3DFilename(const SingleOperation *op, bool allowInverse) {
@@ -2485,6 +2513,45 @@ TransformationNNPtr SingleOperation::substitutePROJAlternativeGridNames(
         }
     }
 
+    const auto &verticalOffsetByVelocityGridFilename =
+        _getVerticalOffsetByVelocityGridFilename(this, false);
+    if (!verticalOffsetByVelocityGridFilename.empty()) {
+        if (databaseContext->lookForGridAlternative(
+                verticalOffsetByVelocityGridFilename, projFilename,
+                projGridFormat, inverseDirection)) {
+
+            if (inverseDirection) {
+                throw util::UnsupportedOperationException(
+                    "Inverse direction for "
+                    "VerticalOffsetByVelocityGrid not supported");
+            }
+
+            if (verticalOffsetByVelocityGridFilename == projFilename) {
+                return self;
+            }
+
+            const auto l_sourceCRSNull = sourceCRS();
+            const auto l_targetCRSNull = targetCRS();
+            if (l_sourceCRSNull == nullptr) {
+                throw util::UnsupportedOperationException("Missing sourceCRS");
+            }
+            if (l_targetCRSNull == nullptr) {
+                throw util::UnsupportedOperationException("Missing targetCRS");
+            }
+            auto l_sourceCRS = NN_NO_CHECK(l_sourceCRSNull);
+            auto l_targetCRS = NN_NO_CHECK(l_targetCRSNull);
+            auto parameters =
+                std::vector<OperationParameterNNPtr>{createOpParamNameEPSGCode(
+                    EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE)};
+            return Transformation::create(
+                createSimilarPropertiesOperation(self), l_sourceCRS,
+                l_targetCRS, l_interpolationCRS,
+                createSimilarPropertiesMethod(method()), parameters,
+                {ParameterValue::createFilename(projFilename)},
+                coordinateOperationAccuracies());
+        }
+    }
+
     bool reverseOffsetSign = false;
     if (methodEPSGCode == EPSG_CODE_METHOD_VERTCON ||
         isRegularVerticalGridMethod(methodEPSGCode, reverseOffsetSign)) {
@@ -4029,6 +4096,87 @@ bool SingleOperation::exportToPROJStringGeneric(
         return true;
     }
 
+    const auto &verticalOffsetByVelocityGridFilename =
+        _getVerticalOffsetByVelocityGridFilename(this, true);
+    if (!verticalOffsetByVelocityGridFilename.empty()) {
+
+        const auto &interpCRS = interpolationCRS();
+        if (!interpCRS) {
+            throw io::FormattingException(
+                "InterpolationCRS required "
+                "for"
+                " " EPSG_NAME_METHOD_VERTICAL_OFFSET_BY_VELOCITY_GRID_NRCAN);
+        }
+
+        auto interpCRSGeog =
+            dynamic_cast<const crs::GeographicCRS *>(interpCRS.get());
+        if (!interpCRSGeog) {
+            throw io::FormattingException(
+                concat("Can apply ", methodName,
+                       " only to a GeographicCRS interpolation CRS"));
+        }
+
+        if (isMethodInverseOf) {
+            formatter->startInversion();
+        }
+        formatter->addStep("push");
+        formatter->addParam("v_1");
+        formatter->addParam("v_2");
+
+        formatter->addStep("cart");
+        interpCRSGeog->ellipsoid()->_exportToPROJString(formatter);
+
+        formatter->addStep("deformation");
+        auto srcName = sourceCRS()->nameStr();
+        auto dstName = targetCRS()->nameStr();
+        const struct {
+            const char *name;
+            double epoch;
+        } realizationEpochs[] = {
+            {"CGVD2013a(1997) height", 1997.0},
+            {"CGVD2013a(2002) height", 2002.0},
+            {"CGVD2013a(2010) height", 2010.0},
+        };
+        double sourceYear = 0.0;
+        double targetYear = 0.0;
+        for (const auto &iter : realizationEpochs) {
+            if (iter.name == srcName)
+                sourceYear = iter.epoch;
+            if (iter.name == dstName)
+                targetYear = iter.epoch;
+        }
+        if (sourceYear == 0.0) {
+            throw io::FormattingException(
+                "For"
+                " " EPSG_NAME_METHOD_VERTICAL_OFFSET_BY_VELOCITY_GRID_NRCAN
+                ", missing epoch for source CRS");
+        }
+        if (targetYear == 0.0) {
+            throw io::FormattingException(
+                "For"
+                " " EPSG_NAME_METHOD_VERTICAL_OFFSET_BY_VELOCITY_GRID_NRCAN
+                ", missing epoch for target CRS");
+        }
+        formatter->addParam("dt", targetYear - sourceYear);
+        formatter->addParam("grids", verticalOffsetByVelocityGridFilename);
+        interpCRSGeog->ellipsoid()->_exportToPROJString(formatter);
+
+        formatter->startInversion();
+        formatter->addStep("cart");
+        interpCRSGeog->ellipsoid()->_exportToPROJString(formatter);
+        formatter->stopInversion();
+
+        formatter->addStep("pop");
+        formatter->addParam("v_1");
+        formatter->addParam("v_2");
+
+        if (isMethodInverseOf) {
+            formatter->stopInversion();
+        }
+
+        return true;
+    }
+
     const auto &heightFilename = _getHeightToGeographic3DFilename(this, true);
     if (!heightFilename.empty()) {
         auto l_targetCRS = targetCRS();

src/proj_constants.h

@@ -596,6 +596,12 @@
 
 /* ------------------------------------------------------------------------ */
 
+#define EPSG_NAME_METHOD_VERTICAL_OFFSET_BY_VELOCITY_GRID_NRCAN                \
+    "Vertical Offset by velocity grid (NRCan byn)"
+#define EPSG_CODE_METHOD_VERTICAL_OFFSET_BY_VELOCITY_GRID_NRCAN 1113
+
+/* ------------------------------------------------------------------------ */
+
 #define PROJ_WKT2_NAME_METHOD_HEIGHT_TO_GEOG3D                                 \
     "GravityRelatedHeight to Geographic3D"
 

test/unit/test_operationfactory.cpp

@@ -10003,3 +10003,43 @@ TEST(operation, createOperation_test_createOperationsWithDatumPivot_iter_1) {
     EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
               "+proj=noop");
 }
+
+// ---------------------------------------------------------------------------
+
+TEST(operation, createOperation_Vrtical_Offset_by_velocity_grid) {
+    auto dbContext = DatabaseContext::create();
+
+    auto objSrc = createFromUserInput("NAD83(CSRS)v7 + CGVD2013a(2002) height",
+                                      dbContext);
+    auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
+    ASSERT_TRUE(src != nullptr);
+
+    auto objDest = createFromUserInput("NAD83(CSRS)v7 + CGVD2013a(2010) height",
+                                       dbContext);
+    auto dst = nn_dynamic_pointer_cast<CRS>(objDest);
+    ASSERT_TRUE(dst != nullptr);
+
+    auto ctxt = CoordinateOperationContext::create(
+        AuthorityFactory::create(dbContext, std::string()), nullptr, 0);
+    ctxt->setSpatialCriterion(
+        CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
+    ctxt->setGridAvailabilityUse(
+        CoordinateOperationContext::GridAvailabilityUse::
+            IGNORE_GRID_AVAILABILITY);
+    auto list = CoordinateOperationFactory::create()->createOperations(
+        NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
+    ASSERT_GE(list.size(), 1U);
+    EXPECT_FALSE(list[0]->hasBallparkTransformation());
+    EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
+              "+proj=pipeline "
+              "+step +proj=axisswap +order=2,1 "
+              "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
+              "+step +proj=push +v_1 +v_2 "
+              "+step +proj=cart +ellps=GRS80 "
+              "+step +inv +proj=deformation +dt=-8 "
+              "+grids=ca_nrc_NAD83v70VG.tif +ellps=GRS80 "
+              "+step +inv +proj=cart +ellps=GRS80 "
+              "+step +proj=pop +v_1 +v_2 "
+              "+step +proj=unitconvert +xy_in=rad +xy_out=deg "
+              "+step +proj=axisswap +order=2,1");
+}