Skip to content

Commit

Permalink
Implement 'Geographic3D Offset by velocity grid (NRCan byn)' transfor…
Browse files Browse the repository at this point in the history 
…mation method (e.g for NAD83(CSRS)v7 to NAD83(CSRS)v3)
  • Loading branch information
@rouault
rouault committed Aug 29, 2023
1 parent 788b131 commit c9dff55
Show file tree
Hide file tree
Showing 5 changed files with 245 additions and 2 deletions.
22 changes: 22 additions & 0 deletions data/sql/grid_transformation.sql
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1440,6 +1440,22 @@ 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','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);
INSERT INTO "usage" VALUES('EPSG','18210','grid_transformation','EPSG','10120','EPSG','1061','EPSG','1274');
INSERT INTO "grid_transformation" VALUES('EPSG','10121','NAD83(CSRS)v6 to NAD83(CSRS)v2 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8251','EPSG','8235',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8252','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','18203','grid_transformation','EPSG','10121','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10122','NAD83(CSRS)v6 to NAD83(CSRS)v3 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8251','EPSG','8239',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8252','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','18351','grid_transformation','EPSG','10122','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10123','NAD83(CSRS)v6 to NAD83(CSRS)v4 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8251','EPSG','8244',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8252','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','18259','grid_transformation','EPSG','10123','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10124','NAD83(CSRS)v7 to NAD83(CSRS)v2 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8254','EPSG','8235',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8255','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','18206','grid_transformation','EPSG','10124','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10125','NAD83(CSRS)v7 to NAD83(CSRS)v3 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8254','EPSG','8239',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8255','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','18207','grid_transformation','EPSG','10125','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10126','NAD83(CSRS)v7 to NAD83(CSRS)v4 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','8254','EPSG','8244',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8255','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','18350','grid_transformation','EPSG','10126','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10128','NAD83(CSRS)v4 to NAD83(CSRS)v4 + CGVD2013a(2002) height (1)','Reversible alternative to NAD83(CSRS)v4 to CGVD2013a(2002) height (1) (code 10110).','EPSG','1090','Geog3D to Geog2D+GravityRelatedHeight (NRCan byn)','EPSG','8244','EPSG','20037',0.05,'EPSG','8666','Geoid (height correction) model file','CGG2013an83.byn',NULL,NULL,NULL,NULL,'EPSG','8246','NR-Can CGG2013a 2002',0);
INSERT INTO "usage" VALUES('EPSG','18287','grid_transformation','EPSG','10128','EPSG','1061','EPSG','1270');
INSERT INTO "grid_transformation" VALUES('EPSG','10129','NAD83(CSRS)v3 to NAD83(CSRS)v3 + CGVD2013a(1997) height (1)','Reversible alternative to NAD83(CSRS)v3 to CGVD2013a(1997) height (1) (code 10111).','EPSG','1090','Geog3D to Geog2D+GravityRelatedHeight (NRCan byn)','EPSG','8239','EPSG','20038',0.05,'EPSG','8666','Geoid (height correction) model file','CGG2013an83.byn',NULL,NULL,NULL,NULL,'EPSG','8240','NR-Can CGG2013a 1997',0);
Expand Down Expand Up @@ -1546,6 +1562,12 @@ INSERT INTO "grid_transformation" VALUES('EPSG','10417','NAD83(CSRS)v8 to CGVD20
INSERT INTO "usage" VALUES('EPSG','20252','grid_transformation','EPSG','10417','EPSG','1061','EPSG','1133');
INSERT INTO "grid_transformation" VALUES('EPSG','10418','NAD83(CSRS)v8 to CGVD28 height (1)','Valid at epoch 2010.0. Hybrid geoid model, grid derived at epoch 1997.0 through NAD83(CSRS)v3 (CT code 9983) and then modified to include correction for propagation of height between 1997.0 and 2010.0 derived from the Canada velocity grid v7.0.','EPSG','1060','Geographic3D to GravityRelatedHeight (NRCan byn)','EPSG','10413','EPSG','5713',0.05,'EPSG','8666','Geoid (height correction) model file','HT2_2010v70.byn',NULL,NULL,NULL,NULL,NULL,NULL,'NRC-Can CGG2000 2010',0);
INSERT INTO "usage" VALUES('EPSG','20137','grid_transformation','EPSG','10418','EPSG','1289','EPSG','1133');
INSERT INTO "grid_transformation" VALUES('EPSG','10421','NAD83(CSRS)v8 to NAD83(CSRS)v2 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','10413','EPSG','8235',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','10414','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','20188','grid_transformation','EPSG','10421','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10422','NAD83(CSRS)v8 to NAD83(CSRS)v3 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','10413','EPSG','8239',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','10414','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','20189','grid_transformation','EPSG','10422','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10423','NAD83(CSRS)v8 to NAD83(CSRS)v4 (1)','','EPSG','1114','Geographic3D Offset by velocity grid (NRCan byn)','EPSG','10413','EPSG','8244',0.05,'EPSG','1050','Point motion velocity grid file','NAD83v70VG.gvb',NULL,NULL,NULL,NULL,'EPSG','8255','EPSG-Can cvg70',0);
INSERT INTO "usage" VALUES('EPSG','20190','grid_transformation','EPSG','10423','EPSG','1061','EPSG','1026');
INSERT INTO "grid_transformation" VALUES('EPSG','10469','ETRS89 to COV23-IRF (1)','In conjunction with the COV23-TM map projection (code 10470) applied to COV23-IRF (code 10468), emulates the COV23 Snake projection. Applied to ETRS89 (as realized through the OSNet v2009) defines COV23-IRF hence is errorless. ','EPSG','9615','NTv2','EPSG','4258','EPSG','10468',0.0,'EPSG','8656','Latitude and longitude difference file','TN15-ETRS89-to-COV23-IRF.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'CCC-Gbr COV23 OSNet2009',0);
INSERT INTO "usage" VALUES('EPSG','20323','grid_transformation','EPSG','10469','EPSG','4743','EPSG','1141');
INSERT INTO "grid_transformation" VALUES('EPSG','15486','CH1903 to CH1903+ (1)','For improved accuracy (0.01m) use CHENyx06 interpolation programme FINELTRA. File CHENyx06 replaced by CHENyx06a; there is a small area at the border of the data where some more real data has been introduced. swisstopo consider the change insignificant.','EPSG','9615','NTv2','EPSG','4149','EPSG','4150',0.2,'EPSG','8656','Latitude and longitude difference file','CHENyx06a.gsb',NULL,NULL,NULL,NULL,NULL,NULL,'BfL-Che',0);
Expand Down
5 changes: 3 additions & 2 deletions scripts/build_db.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)'))")
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)'))")
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
Expand Down Expand Up @@ -775,12 +775,13 @@ 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)
# 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, 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, 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
Expand Down
183 changes: 183 additions & 0 deletions src/iso19111/operation/singleoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1963,6 +1963,34 @@ _getGeocentricTranslationFilename(const SingleOperation *op,

// ---------------------------------------------------------------------------

//! @cond Doxygen_Suppress
static const std::string &
_getGeographic3DOffsetByVelocityGridFilename(const SingleOperation *op,
bool allowInverse) {

const auto &l_method = op->method();
const auto &methodName = l_method->nameStr();
if (l_method->getEPSGCode() ==
EPSG_CODE_METHOD_GEOGRAPHIC3D_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) {
Expand Down Expand Up @@ -2410,6 +2438,45 @@ TransformationNNPtr SingleOperation::substitutePROJAlternativeGridNames(
}
}

const auto &geographic3DOffsetByVelocityGridFilename =
_getGeographic3DOffsetByVelocityGridFilename(this, false);
if (!geographic3DOffsetByVelocityGridFilename.empty()) {
if (databaseContext->lookForGridAlternative(
geographic3DOffsetByVelocityGridFilename, projFilename,
projGridFormat, inverseDirection)) {

if (inverseDirection) {
throw util::UnsupportedOperationException(
"Inverse direction for "
"GeocentricTranslation not supported");
}

if (geographic3DOffsetByVelocityGridFilename == 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());
}
}

if (methodEPSGCode == EPSG_CODE_METHOD_VERTCON ||
isRegularVerticalGridMethod(methodEPSGCode)) {
auto fileParameter =
Expand Down Expand Up @@ -3837,6 +3904,122 @@ bool SingleOperation::exportToPROJStringGeneric(
return true;
}

const auto &geographic3DOffsetByVelocityGridFilename =
_getGeographic3DOffsetByVelocityGridFilename(this, true);
if (!geographic3DOffsetByVelocityGridFilename.empty()) {
auto sourceCRSGeog =
dynamic_cast<const crs::GeographicCRS *>(sourceCRS().get());
if (!sourceCRSGeog) {
throw io::FormattingException(
concat("Can apply ", methodName, " only to GeographicCRS"));
}

auto targetCRSGeog =
dynamic_cast<const crs::GeographicCRS *>(targetCRS().get());
if (!targetCRSGeog) {
throw io::FormattingException(
concat("Can apply ", methodName, " only to GeographicCRS"));
}

const auto &interpCRS = interpolationCRS();
if (!interpCRS) {
throw io::FormattingException(
"InterpolationCRS required "
"for"
" " EPSG_NAME_METHOD_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN);
}
const bool interpIsSrc = interpCRS->_isEquivalentTo(
sourceCRS()->demoteTo2D(std::string(), nullptr).get(),
util::IComparable::Criterion::EQUIVALENT);
const bool interpIsTarget = interpCRS->_isEquivalentTo(
targetCRS()->demoteTo2D(std::string(), nullptr).get(),
util::IComparable::Criterion::EQUIVALENT);
if (!interpIsSrc && !interpIsTarget) {
throw io::FormattingException(
"For"
" " EPSG_NAME_METHOD_GEOGRAPHIC3D_OFFSET_BY_VELOCITY_GRID_NRCAN
", interpolation CRS should be the source or target CRS");
}

formatter->startInversion();
sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
formatter->stopInversion();

if (isMethodInverseOf) {
formatter->startInversion();
}

const bool addPushPopV3 =
((sourceCRSGeog &&
sourceCRSGeog->coordinateSystem()->axisList().size() == 2) ||
(targetCRSGeog &&
targetCRSGeog->coordinateSystem()->axisList().size() == 2));

if (addPushPopV3) {
formatter->addStep("push");
formatter->addParam("v_3");
}

formatter->addStep("cart");
sourceCRSGeog->ellipsoid()->_exportToPROJString(formatter);

formatter->addStep("deformation");
auto srcName = sourceCRS()->nameStr();
auto dstName = targetCRS()->nameStr();
const struct {
const char *name;
double epoch;
} realizationEpochs[] = {
{"NAD83(CSRS)v2", 1997.0}, {"NAD83(CSRS)v3", 1997.0},
{"NAD83(CSRS)v4", 2002.0}, {"NAD83(CSRS)v5", 2006.0},
{"NAD83(CSRS)v6", 2010.0}, {"NAD83(CSRS)v7", 2010.0},
{"NAD83(CSRS)v8", 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_GEOGRAPHIC3D_OFFSET_BY_VELOCITY_GRID_NRCAN
", missing epoch for source CRS");
}
if (targetYear == 0.0) {
throw io::FormattingException(
"For"
" " EPSG_NAME_METHOD_GEOGRAPHIC3D_OFFSET_BY_VELOCITY_GRID_NRCAN
", missing epoch for target CRS");
}
formatter->addParam("dt", targetYear - sourceYear);
formatter->addParam("grids", geographic3DOffsetByVelocityGridFilename);
(interpIsTarget ? targetCRSGeog : sourceCRSGeog)
->ellipsoid()
->_exportToPROJString(formatter);

formatter->startInversion();
formatter->addStep("cart");
targetCRSGeog->ellipsoid()->_exportToPROJString(formatter);
formatter->stopInversion();

if (addPushPopV3) {
formatter->addStep("pop");
formatter->addParam("v_3");
}

if (isMethodInverseOf) {
formatter->stopInversion();
}

targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);

return true;
}

const auto &heightFilename = _getHeightToGeographic3DFilename(this, true);
if (!heightFilename.empty()) {
auto l_targetCRS = targetCRS();
Expand Down
Loading

0 comments on commit c9dff55

Please sign in to comment.