From deb31c9d3c6ca15399d3286456d34bb917805fdc Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 16 Dec 2024 19:15:16 +0100 Subject: [PATCH 1/4] AuthorityFactory::createCoordinateOperation(): consider accuracy = '999.0' as unknown --- src/iso19111/factory.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 58effcc5c0..50e93aa55a 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -6281,7 +6281,7 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation( .set(common::IdentifiedObject::NAME_KEY, method_name); std::vector accuracies; - if (!accuracy.empty()) { + if (!accuracy.empty() && accuracy != "999.0") { accuracies.emplace_back( metadata::PositionalAccuracy::create(accuracy)); } @@ -6393,7 +6393,7 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation( .set(common::IdentifiedObject::NAME_KEY, method_name); std::vector accuracies; - if (!accuracy.empty()) { + if (!accuracy.empty() && accuracy != "999.0") { accuracies.emplace_back( metadata::PositionalAccuracy::create(accuracy)); } @@ -6531,7 +6531,7 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation( } std::vector accuracies; - if (!accuracy.empty()) { + if (!accuracy.empty() && accuracy != "999.0") { accuracies.emplace_back( metadata::PositionalAccuracy::create(accuracy)); } @@ -6646,8 +6646,10 @@ operation::CoordinateOperationNNPtr AuthorityFactory::createCoordinateOperation( std::vector accuracies; if (!accuracy.empty()) { - accuracies.emplace_back( - metadata::PositionalAccuracy::create(accuracy)); + if (accuracy != "999.0") { + accuracies.emplace_back( + metadata::PositionalAccuracy::create(accuracy)); + } } else { // Try to compute a reasonable accuracy from the members double totalAcc = -1; From e7ef01b1de5ef586576d3303476705a2c66da150 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 16 Dec 2024 19:15:18 +0100 Subject: [PATCH 2/4] Take into account new EPSG methods 'Cartesian Grid Offsets by TIN Interpolation (JSON)' and 'Vertical Offset by TIN Interpolation (JSON)' and remove custom transformations for Finland TIN base transformations that were previously under the PROJ authority --- data/CMakeLists.txt | 2 +- data/sql/consistency_checks_triggers.sql | 6 + data/sql/grid_transformation.sql | 6 + data/sql/other_transformation_custom.sql | 33 +-- data/sql/proj_db_table_defs.sql | 4 +- scripts/build_db.py | 8 +- .../operation/coordinateoperationfactory.cpp | 50 +++-- src/iso19111/operation/parammappings.cpp | 17 ++ src/iso19111/operation/singleoperation.cpp | 198 +++++++++++++----- src/proj_constants.h | 15 ++ 10 files changed, 236 insertions(+), 103 deletions(-) diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt index 2434f204fe..36efc29975 100644 --- a/data/CMakeLists.txt +++ b/data/CMakeLists.txt @@ -45,7 +45,7 @@ set(ALL_SQL_IN "${CMAKE_CURRENT_BINARY_DIR}/all.sql.in") set(PROJ_DB "${CMAKE_CURRENT_BINARY_DIR}/proj.db") include(sql_filelist.cmake) -set(PROJ_DB_SQL_EXPECTED_MD5 "179a5d10e801cb758db7aaa925b97a8f") +set(PROJ_DB_SQL_EXPECTED_MD5 "c2184d9fca631e88f21e4c0d8899d391") add_custom_command( OUTPUT ${PROJ_DB} diff --git a/data/sql/consistency_checks_triggers.sql b/data/sql/consistency_checks_triggers.sql index 6e4db5aebe..cb3fbd1d5b 100644 --- a/data/sql/consistency_checks_triggers.sql +++ b/data/sql/consistency_checks_triggers.sql @@ -346,6 +346,12 @@ FOR EACH ROW BEGIN SELECT RAISE(ABORT, 'insert on grid_transformation violates constraint: target_crs(auth_name, code) not found') WHERE NOT EXISTS (SELECT 1 FROM crs_view WHERE crs_view.auth_name = NEW.target_crs_auth_name AND crs_view.code = NEW.target_crs_code); + SELECT RAISE(ABORT, 'insert on grid_transformation violates constraint: interpolation_crs(auth_name, code) not found') + WHERE NEW.interpolation_crs_code IS NOT NULL AND NOT EXISTS (SELECT 1 FROM crs_view WHERE crs_view.auth_name = NEW.interpolation_crs_auth_name AND crs_view.code = NEW.interpolation_crs_code); + + SELECT RAISE(ABORT, 'insert on grid_transformation violates constraint: interpolation_crs must be a GeodeticCRS on non-TIN shift based files') + WHERE NEW.method_name NOT LIKE '%JSON%' AND NEW.interpolation_crs_code IS NOT NULL AND NOT EXISTS (SELECT 1 FROM geodetic_crs WHERE geodetic_crs.auth_name = NEW.interpolation_crs_auth_name AND geodetic_crs.code = NEW.interpolation_crs_code); + SELECT RAISE(ABORT, 'insert on grid_transformation violates constraint: source_crs must not be deprecated when grid_transformation is not deprecated') WHERE EXISTS(SELECT 1 FROM crs_view crs WHERE crs.auth_name = NEW.source_crs_auth_name AND crs.code = NEW.source_crs_code AND crs.deprecated != 0) AND NEW.deprecated = 0 AND NOT (NEW.auth_name = 'ESRI'); SELECT RAISE(ABORT, 'insert on grid_transformation violates constraint: target_crs must not be deprecated when grid_transformation is not deprecated') diff --git a/data/sql/grid_transformation.sql b/data/sql/grid_transformation.sql index 6d3086f5fe..e6642a498d 100644 --- a/data/sql/grid_transformation.sql +++ b/data/sql/grid_transformation.sql @@ -1762,6 +1762,12 @@ INSERT INTO "grid_transformation" VALUES('EPSG','10697','EUREF-FIN to N2000 heig INSERT INTO "usage" VALUES('EPSG','22014','grid_transformation','EPSG','10697','EPSG','3333','EPSG','1133'); INSERT INTO "grid_transformation" VALUES('EPSG','10698','EUREF-FIN to EUREF-FIN + N2000 height (2)','Replaces EUREF-FIN to EUREF-FIN + N2000 height (1) (code 10696). Reversible alternative to EUREF-FIN to N2000 height (2) (code 10697). File also available in 3 ASCII formats.','EPSG','1124','Geog3D to Geog2D+GravityRelatedHeight (gtg)','EPSG','10689','EPSG','10692',0.014,'EPSG','8666','Geoid (height correction) model file','fi_nls_fin2023n2000.tif',NULL,NULL,NULL,NULL,'EPSG','10690','NLS-Fin 2023',0); INSERT INTO "usage" VALUES('EPSG','22015','grid_transformation','EPSG','10698','EPSG','3333','EPSG','1270'); +INSERT INTO "grid_transformation" VALUES('EPSG','10703','KKJ / Finland Uniform Coordinate System to EUREF-FIN / TM35FIN(E,N) (1)','This YKJ to ETRS-TM35FIN TINshift transformation is the recommended method for KKJ to EUREF-FIN transformations, to be used in preference to transformation KKJ to ETRS89 (2) (code 10098).','EPSG','1138','Cartesian Grid Offsets by TIN Interpolation (JSON)','EPSG','2393','EPSG','3067',0.03,'EPSG','1064','TIN offset file','fi_nls_ykj_etrs35fin.json',NULL,NULL,NULL,NULL,NULL,NULL,'NLS-FIN TINshift',0); +INSERT INTO "usage" VALUES('EPSG','22055','grid_transformation','EPSG','10703','EPSG','3333','EPSG','1273'); +INSERT INTO "grid_transformation" VALUES('EPSG','10704','N43 height to N60 height (1)','Accuracy has not been determined.','EPSG','1137','Vertical Offset by TIN Interpolation (JSON)','EPSG','8675','EPSG','5717',999.0,'EPSG','1064','TIN offset file','fi_nls_n43_n60.json',NULL,NULL,NULL,NULL,'EPSG','2393','NLS-Fin TINshift',0); +INSERT INTO "usage" VALUES('EPSG','22018','grid_transformation','EPSG','10704','EPSG','4522','EPSG','1059'); +INSERT INTO "grid_transformation" VALUES('EPSG','10705','N60 height to N2000 height (1)','Recommended method for N60<>N2000 trasformations, as more accurate than concatenating geoid models.','EPSG','1137','Vertical Offset by TIN Interpolation (JSON)','EPSG','5717','EPSG','3900',0.005,'EPSG','1064','TIN offset file','fi_nls_n60_n2000.json',NULL,NULL,NULL,NULL,'EPSG','2393','NLS-Fin TINshift',0); +INSERT INTO "usage" VALUES('EPSG','22201','grid_transformation','EPSG','10705','EPSG','3333','EPSG','1059'); 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); INSERT INTO "usage" VALUES('EPSG','11497','grid_transformation','EPSG','15486','EPSG','1286','EPSG','1085'); INSERT INTO "grid_transformation" VALUES('EPSG','15488','RRAF 1991 to IGN 1988 MG height (1)','May be used for transformations from WGS 84 to IGN 1988 MG. Accuracy at each 0.025 deg x 0.025 degree grid node is given within the geoid model file.','EPSG','9664','Geographic3D to GravityRelatedHeight (IGN1997)','EPSG','4973','EPSG','5617',0.2,'EPSG','8666','Geoid (height correction) model file','ggg00_mg.txt',NULL,NULL,NULL,NULL,NULL,NULL,'IGN Glp MG',1); diff --git a/data/sql/other_transformation_custom.sql b/data/sql/other_transformation_custom.sql index 1630e3f3e5..b52ec92f89 100644 --- a/data/sql/other_transformation_custom.sql +++ b/data/sql/other_transformation_custom.sql @@ -18,17 +18,9 @@ INSERT INTO "usage" VALUES('PROJ','NGO48_TO_ETRS89NO_USAGE','other_transformatio -- Finland triangulated files -INSERT INTO other_transformation VALUES( - 'PROJ','YKJ_TO_ETRS35FIN','KKJ / Finland Uniform Coordinate System to ETRS35FIN', - 'Transformation based on a triangulated irregular network', - 'PROJ','PROJString','+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=tinshift +file=fi_nls_ykj_etrs35fin.json', - 'EPSG','2393','EPSG','3067',0.1, - NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0); -INSERT INTO "usage" VALUES('PROJ','YKJ_TO_ETRS35FIN_USAGE','other_transformation','PROJ','YKJ_TO_ETRS35FIN','EPSG','3333','EPSG','1024'); - -INSERT INTO "concatenated_operation" VALUES('PROJ','KKJ_TO_ETRS89','KKJ to ETRS89 (using PROJ:YKJ_TO_ETRS35FIN)','Transformation based on a triangulated irregular network','EPSG','4123','EPSG','4258',NULL,NULL,0); +INSERT INTO "concatenated_operation" VALUES('PROJ','KKJ_TO_ETRS89','KKJ to ETRS89 (using EPSG:10703)','Transformation based on a triangulated irregular network','EPSG','4123','EPSG','4258',NULL,NULL,0); INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_ETRS89',1,'EPSG','18193'); -INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_ETRS89',2,'PROJ','YKJ_TO_ETRS35FIN'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_ETRS89',2,'EPSG','10703'); INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_ETRS89',3,'EPSG','16065'); INSERT INTO "usage" VALUES( 'PROJ', @@ -40,9 +32,9 @@ INSERT INTO "usage" VALUES( 'EPSG','1024' -- unknown ); -INSERT INTO "concatenated_operation" VALUES('PROJ','KKJ_TO_EUREF_FIN','KKJ to EUREF-FIN (using PROJ:YKJ_TO_ETRS35FIN)','Transformation based on a triangulated irregular network','EPSG','4123','EPSG','10690',NULL,NULL,0); +INSERT INTO "concatenated_operation" VALUES('PROJ','KKJ_TO_EUREF_FIN','KKJ to EUREF-FIN (using EPSG:10703)','Transformation based on a triangulated irregular network','EPSG','4123','EPSG','10690',NULL,NULL,0); INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_EUREF_FIN',1,'EPSG','18193'); -INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_EUREF_FIN',2,'PROJ','YKJ_TO_ETRS35FIN'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_EUREF_FIN',2,'EPSG','10703'); INSERT INTO "concatenated_operation_step" VALUES('PROJ','KKJ_TO_EUREF_FIN',3,'EPSG','16065'); INSERT INTO "usage" VALUES( 'PROJ', @@ -53,20 +45,3 @@ INSERT INTO "usage" VALUES( 'EPSG','3333', -- extent 'EPSG','1024' -- unknown ); - -INSERT INTO other_transformation VALUES( - 'PROJ','N43_TO_N60','N43 height to N60 height', - 'Transformation based on a triangulated irregular network', - 'PROJ','PROJString','+proj=pipeline +step +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0 +ellps=intl +step +proj=tinshift +file=fi_nls_n43_n60.json +step +inv +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0 +ellps=intl', - 'EPSG','8675','EPSG','5717',0.01, - NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'EPSG','4123',NULL,0); -INSERT INTO "usage" VALUES('PROJ','N43_TO_N60_USAGE','other_transformation','PROJ','N43_TO_N60','EPSG','4522','EPSG','1024'); - - -INSERT INTO other_transformation VALUES( - 'PROJ','N60_TO_N2000','N60 height to N2000 height', - 'Transformation based on a triangulated irregular network', - 'PROJ','PROJString','+proj=pipeline +step +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0 +ellps=intl +step +proj=tinshift +file=fi_nls_n60_n2000.json +step +inv +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0 +ellps=intl', - 'EPSG','5717','EPSG','3900',0.01, - NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,'EPSG','4123',NULL,0); -INSERT INTO "usage" VALUES('PROJ','N60_TO_N2000_USAGE','other_transformation','PROJ','N60_TO_N2000','EPSG','3333','EPSG','1024'); diff --git a/data/sql/proj_db_table_defs.sql b/data/sql/proj_db_table_defs.sql index d7d61c367f..45393cacf9 100644 --- a/data/sql/proj_db_table_defs.sql +++ b/data/sql/proj_db_table_defs.sql @@ -833,11 +833,11 @@ CREATE TABLE grid_transformation( deprecated BOOLEAN NOT NULL CHECK (deprecated IN (0, 1)), - CONSTRAINT pk_grid_transformation PRIMARY KEY (auth_name, code), + CONSTRAINT pk_grid_transformation PRIMARY KEY (auth_name, code) --CONSTRAINT fk_grid_transformation_coordinate_operation FOREIGN KEY (auth_name, code) REFERENCES coordinate_operation(auth_name, code) ON DELETE CASCADE, --CONSTRAINT fk_grid_transformation_source_crs FOREIGN KEY (source_crs_auth_name, source_crs_code) REFERENCES crs(auth_name, code) ON DELETE CASCADE, --CONSTRAINT fk_grid_transformation_target_crs FOREIGN KEY (target_crs_auth_name, target_crs_code) REFERENCES crs(auth_name, code) ON DELETE CASCADE, - CONSTRAINT fk_grid_transformation_interpolation_crs FOREIGN KEY (interpolation_crs_auth_name, interpolation_crs_code) REFERENCES geodetic_crs(auth_name, code) ON DELETE CASCADE + -- CONSTRAINT fk_grid_transformation_interpolation_crs FOREIGN KEY (interpolation_crs_auth_name, interpolation_crs_code) REFERENCES crs_view(auth_name, code) ON DELETE CASCADE ) WITHOUT ROWID; -- Table that describe packages/archives that contain several grids diff --git a/scripts/build_db.py b/scripts/build_db.py index dd798fbce9..e099cbbcab 100755 --- a/scripts/build_db.py +++ b/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_code NOT IN (1131, 1136) 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 % by %grid%' 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 (NTv2_Vel)', 'New Zealand Deformation Model'))") + 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_code NOT IN (1131, 1136) 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 % by %grid%' 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 (NTv2_Vel)', 'New Zealand Deformation Model', 'Cartesian Grid Offsets by TIN Interpolation (JSON)', 'Vertical Offset by TIN Interpolation (JSON)'))") 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 @@ -731,7 +731,7 @@ def fill_grid_transformation(proj_db_cursor): expected_order += 1 n_params = expected_order - 1 - assert param_code[0] in (1048, 1050, 1063, 8656, 8657, 8666, 8732, 8727), (code, param_code[0]) + assert param_code[0] in (1048, 1050, 1063, 1064, 8656, 8657, 8666, 8732, 8727), (code, param_code[0]) grid2_param_auth_name = None grid2_param_code = None @@ -783,9 +783,11 @@ def fill_grid_transformation(proj_db_cursor): # 1124: Geog3D to Geog2D+GravityRelatedHeight (gtg) # 1126: Vertical change by geoid grid difference (NRCan) # 1135: Geog3D to Geog2D+GravityRelatedHeight (NGS bin) + # 1137: Vertical Offset by TIN Interpolation (JSON) + # 1138: Cartesian Grid Offsets by TIN Interpolation (JSON) # 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, 1113, 1114, 1115, 1118, 1122, 1124, 1126, 1128, 1129, 1135) 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, 1124, 1126, 1128, 1129, 1135, 1137, 1138) 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 diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index fe17a6ebc1..591286a036 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -2326,17 +2326,17 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final CoordinateOperationPtr opSrcCRSToGeogCRS{}; CoordinateOperationPtr verticalTransform{}; CoordinateOperationPtr opGeogCRStoDstCRS{}; - crs::GeographicCRSPtr interpolationGeogCRS{}; + crs::SingleCRSPtr interpolationCRS{}; MyPROJStringExportableHorizVerticalHorizPROJBased( const CoordinateOperationPtr &opSrcCRSToGeogCRSIn, const CoordinateOperationPtr &verticalTransformIn, const CoordinateOperationPtr &opGeogCRStoDstCRSIn, - const crs::GeographicCRSPtr &interpolationGeogCRSIn) + const crs::SingleCRSPtr &interpolationCRSIn) : opSrcCRSToGeogCRS(opSrcCRSToGeogCRSIn), verticalTransform(verticalTransformIn), opGeogCRStoDstCRS(opGeogCRStoDstCRSIn), - interpolationGeogCRS(interpolationGeogCRSIn) {} + interpolationCRS(interpolationCRSIn) {} ~MyPROJStringExportableHorizVerticalHorizPROJBased() override; @@ -2416,12 +2416,20 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final formatter->addParam("v_2"); } + auto interpolationCRSGeog = + dynamic_cast(interpolationCRS.get()); + auto interpolationCRSProj = + dynamic_cast(interpolationCRS.get()); + formatter->pushOmitZUnitConversion(); opSrcCRSToGeogCRS->_exportToPROJString(formatter); formatter->startInversion(); - interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter); + if (interpolationCRSGeog) + interpolationCRSGeog->addAngularUnitConvertAndAxisSwap(formatter); + else if (interpolationCRSProj) + interpolationCRSProj->addUnitConvertAndAxisSwap(formatter, false); formatter->stopInversion(); formatter->popOmitZUnitConversion(); @@ -2432,7 +2440,10 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final formatter->pushOmitZUnitConversion(); - interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter); + if (interpolationCRSGeog) + interpolationCRSGeog->addAngularUnitConvertAndAxisSwap(formatter); + else if (interpolationCRSProj) + interpolationCRSProj->addUnitConvertAndAxisSwap(formatter, false); opGeogCRStoDstCRS->_exportToPROJString(formatter); @@ -2682,12 +2693,12 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased( const operation::CoordinateOperationNNPtr &opSrcCRSToGeogCRS, const operation::CoordinateOperationNNPtr &verticalTransform, const operation::CoordinateOperationNNPtr &opGeogCRStoDstCRS, - const crs::GeographicCRSPtr &interpolationGeogCRS, bool checkExtent) { + const crs::SingleCRSPtr &interpolationCRS, bool checkExtent) { auto exportable = util::nn_make_shared( opSrcCRSToGeogCRS, verticalTransform, opGeogCRStoDstCRS, - interpolationGeogCRS); + interpolationCRS); std::vector ops; if (!(starts_with(opSrcCRSToGeogCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET) && @@ -6860,15 +6871,14 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( for (const auto &verticalTransform : verticalTransforms) { - auto interpolationGeogCRS = NN_NO_CHECK(srcGeog); + auto interpolationCRS = + NN_NO_CHECK(std::static_pointer_cast(srcGeog)); auto interpTransformCRS = verticalTransform->interpolationCRS(); if (interpTransformCRS) { - auto nn_interpTransformCRS = NN_NO_CHECK(interpTransformCRS); - if (dynamic_cast( - nn_interpTransformCRS.get())) { - interpolationGeogCRS = NN_NO_CHECK( - util::nn_dynamic_pointer_cast( - nn_interpTransformCRS)); + auto interpTransformSingleCRS = + std::static_pointer_cast(interpTransformCRS); + if (interpTransformSingleCRS) { + interpolationCRS = NN_NO_CHECK(interpTransformSingleCRS); } } else { auto compSrc0BoundCrs = @@ -6880,8 +6890,8 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( compSrc0BoundCrs->hubCRS().get()) && compSrc0BoundCrs->hubCRS()->_isEquivalentTo( compDst0BoundCrs->hubCRS().get())) { - interpolationGeogCRS = NN_NO_CHECK( - util::nn_dynamic_pointer_cast( + interpolationCRS = + NN_NO_CHECK(util::nn_dynamic_pointer_cast( compSrc0BoundCrs->hubCRS())); } } @@ -6919,7 +6929,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( verticalTransform->nameStr().find( "CGVD28 height to CGVD2013a(2002) height (1)") != std::string::npos))) { - interpolationGeogCRS = std::move(nad83CSRSv7); + interpolationCRS = std::move(nad83CSRSv7); } } catch (const std::exception &) { } @@ -6927,9 +6937,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( const auto opSrcCRSToGeogCRS = createOperations( componentsSrc[0], util::optional(), - interpolationGeogCRS, util::optional(), context); + interpolationCRS, util::optional(), context); const auto opGeogCRStoDstCRS = createOperations( - interpolationGeogCRS, util::optional(), + interpolationCRS, util::optional(), componentsDst[0], util::optional(), context); for (const auto &opSrc : opSrcCRSToGeogCRS) { for (const auto &opDst : opGeogCRStoDstCRS) { @@ -6937,7 +6947,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( try { auto op = createHorizVerticalHorizPROJBased( sourceCRS, targetCRS, opSrc, verticalTransform, opDst, - interpolationGeogCRS, true); + interpolationCRS, true); res.emplace_back(op); } catch (const InvalidOperationEmptyIntersection &) { } catch (const io::FormattingException &) { diff --git a/src/iso19111/operation/parammappings.cpp b/src/iso19111/operation/parammappings.cpp index 0abaa8b61a..faea73b089 100644 --- a/src/iso19111/operation/parammappings.cpp +++ b/src/iso19111/operation/parammappings.cpp @@ -1035,6 +1035,8 @@ const struct MethodNameCode methodNameCodesList[] = { METHOD_NAME_CODE(VERTICAL_OFFSET_AND_SLOPE), METHOD_NAME_CODE(NTV2), METHOD_NAME_CODE(NTV1), + METHOD_NAME_CODE(VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON), + METHOD_NAME_CODE(CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON), METHOD_NAME_CODE(NADCON), METHOD_NAME_CODE(NADCON5_2D), METHOD_NAME_CODE(NADCON5_3D), @@ -1486,6 +1488,13 @@ static const ParamMapping *const paramsPoleRotationNetCDFCFConvention[] = { ¶mGridNorthPoleLatitudeNetCDF, ¶mGridNorthPoleLongitudeNetCDF, ¶mNorthPoleGridLongitudeNetCDF, nullptr}; +static const ParamMapping paramTINOffsetFile = { + EPSG_NAME_PARAMETER_TIN_OFFSET_FILE, EPSG_CODE_PARAMETER_TIN_OFFSET_FILE, + nullptr, common::UnitOfMeasure::Type::NONE, nullptr}; + +static const ParamMapping *const paramsTINOffsetFile[] = {¶mTINOffsetFile, + nullptr}; + static const MethodMapping gOtherMethodMappings[] = { {EPSG_NAME_METHOD_CHANGE_VERTICAL_UNIT, EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT, nullptr, nullptr, nullptr, @@ -1646,6 +1655,14 @@ static const MethodMapping gOtherMethodMappings[] = { EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN, nullptr, nullptr, nullptr, paramsGeocentricTranslationGridInterpolationIGN}, + {EPSG_NAME_METHOD_VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON, + EPSG_CODE_METHOD_VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON, nullptr, + nullptr, nullptr, paramsTINOffsetFile}, + + {EPSG_NAME_METHOD_CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON, + EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON, nullptr, + nullptr, nullptr, paramsTINOffsetFile}, + {EPSG_NAME_METHOD_NADCON, EPSG_CODE_METHOD_NADCON, nullptr, nullptr, nullptr, paramsNADCON}, diff --git a/src/iso19111/operation/singleoperation.cpp b/src/iso19111/operation/singleoperation.cpp index 42ca9ffc63..60c280dcb4 100644 --- a/src/iso19111/operation/singleoperation.cpp +++ b/src/iso19111/operation/singleoperation.cpp @@ -2700,60 +2700,81 @@ TransformationNNPtr SingleOperation::substitutePROJAlternativeGridNames( } } - if (methodEPSGCode == EPSG_CODE_METHOD_NEW_ZEALAND_DEFORMATION_MODEL) { - auto fileParameter = - parameterValue(EPSG_NAME_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE, - EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE); - if (fileParameter && - fileParameter->type() == ParameterValue::Type::FILENAME) { + static const struct { + int methodEPSGCode; + int gridFilenameParamEPSGCode; + const char *gridFilenameParamName; + } gridTransformations[] = { + {EPSG_CODE_METHOD_NEW_ZEALAND_DEFORMATION_MODEL, + EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE, + EPSG_NAME_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE}, + {EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON, + EPSG_CODE_PARAMETER_TIN_OFFSET_FILE, + EPSG_NAME_PARAMETER_TIN_OFFSET_FILE}, + {EPSG_CODE_METHOD_VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON, + EPSG_CODE_PARAMETER_TIN_OFFSET_FILE, + EPSG_NAME_PARAMETER_TIN_OFFSET_FILE}, + }; - const auto &filename = fileParameter->valueFile(); - if (databaseContext->lookForGridAlternative( - filename, projFilename, projGridFormat, inverseDirection)) { + for (const auto &gridTransf : gridTransformations) { + if (methodEPSGCode == gridTransf.methodEPSGCode) { + auto fileParameter = + parameterValue(gridTransf.gridFilenameParamName, + gridTransf.gridFilenameParamEPSGCode); + if (fileParameter && + fileParameter->type() == ParameterValue::Type::FILENAME) { - if (filename == projFilename) { - if (inverseDirection) { - throw util::UnsupportedOperationException( - "Inverse direction for " + projFilename + - " not supported"); + const auto &filename = fileParameter->valueFile(); + if (databaseContext->lookForGridAlternative( + filename, projFilename, projGridFormat, + inverseDirection)) { + + if (filename == projFilename) { + if (inverseDirection) { + throw util::UnsupportedOperationException( + "Inverse direction for " + projFilename + + " not supported"); + } + return self; } - 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{ - createOpParamNameEPSGCode( - EPSG_CODE_PARAMETER_POINT_MOTION_VELOCITY_GRID_FILE)}; - if (inverseDirection) { - return Transformation::create( - createPropertiesForInverse( - self.as_nullable().get(), true, false), - l_targetCRS, l_sourceCRS, l_interpolationCRS, - createSimilarPropertiesMethod(method()), - parameters, - {ParameterValue::createFilename(projFilename)}, - coordinateOperationAccuracies()) - ->inverseAsTransformation(); - } else { - return Transformation::create( - createSimilarPropertiesOperation(self), l_sourceCRS, - l_targetCRS, l_interpolationCRS, - createSimilarPropertiesMethod(method()), parameters, - {ParameterValue::createFilename(projFilename)}, - coordinateOperationAccuracies()); + 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{ + createOpParamNameEPSGCode( + gridTransf.gridFilenameParamEPSGCode)}; + if (inverseDirection) { + return Transformation::create( + createPropertiesForInverse( + self.as_nullable().get(), true, false), + l_targetCRS, l_sourceCRS, l_interpolationCRS, + createSimilarPropertiesMethod(method()), + parameters, + {ParameterValue::createFilename( + projFilename)}, + coordinateOperationAccuracies()) + ->inverseAsTransformation(); + } else { + return Transformation::create( + createSimilarPropertiesOperation(self), l_sourceCRS, + l_targetCRS, l_interpolationCRS, + createSimilarPropertiesMethod(method()), parameters, + {ParameterValue::createFilename(projFilename)}, + coordinateOperationAccuracies()); + } } } + break; } } @@ -4765,6 +4786,87 @@ bool SingleOperation::exportToPROJStringGeneric( } } + if (methodEPSGCode == + EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON) { + auto sourceCRSProj = + dynamic_cast(sourceCRS().get()); + if (!sourceCRSProj) { + throw io::FormattingException( + concat("Can apply ", methodName, " only to ProjectedCRS")); + } + + auto targetCRSProj = + dynamic_cast(targetCRS().get()); + if (!targetCRSProj) { + throw io::FormattingException( + concat("Can apply ", methodName, " only to ProjectedCRS")); + } + + auto fileParameter = + parameterValue(EPSG_NAME_PARAMETER_TIN_OFFSET_FILE, + EPSG_CODE_PARAMETER_TIN_OFFSET_FILE); + if (fileParameter && + fileParameter->type() == ParameterValue::Type::FILENAME) { + + formatter->startInversion(); + sourceCRSProj->addUnitConvertAndAxisSwap(formatter, false); + formatter->stopInversion(); + + if (isMethodInverseOf) { + formatter->startInversion(); + } + + formatter->addStep("tinshift"); + formatter->addParam("file", fileParameter->valueFile()); + + if (isMethodInverseOf) { + formatter->stopInversion(); + } + + targetCRSProj->addUnitConvertAndAxisSwap(formatter, false); + + return true; + } + } + + if (methodEPSGCode == + EPSG_CODE_METHOD_VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON) { + auto sourceCRSVert = + dynamic_cast(sourceCRS().get()); + if (!sourceCRSVert) { + throw io::FormattingException( + concat("Can apply ", methodName, " only to VerticalCRS")); + } + + auto targetCRSVert = + dynamic_cast(targetCRS().get()); + if (!targetCRSVert) { + throw io::FormattingException( + concat("Can apply ", methodName, " only to VerticalCRS")); + } + + auto fileParameter = + parameterValue(EPSG_NAME_PARAMETER_TIN_OFFSET_FILE, + EPSG_CODE_PARAMETER_TIN_OFFSET_FILE); + + if (fileParameter && + fileParameter->type() == ParameterValue::Type::FILENAME) { + + if (isMethodInverseOf) { + formatter->startInversion(); + } + + formatter->addStep("tinshift"); + formatter->addParam("file", fileParameter->valueFile()); + + if (isMethodInverseOf) { + formatter->stopInversion(); + } + + return true; + } + } + const char *prefix = "PROJ-based operation method: "; if (starts_with(method()->nameStr(), prefix)) { auto projString = method()->nameStr().substr(strlen(prefix)); diff --git a/src/proj_constants.h b/src/proj_constants.h index e64ec76540..a57ae2e36d 100644 --- a/src/proj_constants.h +++ b/src/proj_constants.h @@ -742,6 +742,21 @@ /* ------------------------------------------------------------------------ */ +/* TIN-based transformations */ + +#define EPSG_NAME_METHOD_VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON \ + "Vertical Offset by TIN Interpolation (JSON)" +#define EPSG_CODE_METHOD_VERTICAL_OFFSET_BY_TIN_INTERPOLATION_JSON 1137 + +#define EPSG_NAME_METHOD_CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON \ + "Cartesian Grid Offsets by TIN Interpolation (JSON)" +#define EPSG_CODE_METHOD_CARTESIAN_GRID_OFFSETS_BY_TIN_INTERPOLATION_JSON 1138 + +#define EPSG_NAME_PARAMETER_TIN_OFFSET_FILE "TIN offset file" +#define EPSG_CODE_PARAMETER_TIN_OFFSET_FILE 1064 + +/* ------------------------------------------------------------------------ */ + #define EPSG_CODE_METHOD_CHANGE_VERTICAL_UNIT 1069 #define EPSG_NAME_METHOD_CHANGE_VERTICAL_UNIT "Change of Vertical Unit" From 40170a847c8d134f228c60232d321c80717159d2 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 17 Dec 2024 01:33:09 +0100 Subject: [PATCH 3/4] lib_proj.cmake: add a dependency to proj.db in the C23 code path --- data/CMakeLists.txt | 2 +- data/generate_proj_db.cmake | 2 ++ src/embedded_resources.c | 4 ++++ src/lib_proj.cmake | 2 ++ 4 files changed, 9 insertions(+), 1 deletion(-) diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt index 36efc29975..ede6cd5232 100644 --- a/data/CMakeLists.txt +++ b/data/CMakeLists.txt @@ -50,7 +50,7 @@ set(PROJ_DB_SQL_EXPECTED_MD5 "c2184d9fca631e88f21e4c0d8899d391") add_custom_command( OUTPUT ${PROJ_DB} COMMAND ${CMAKE_COMMAND} -E remove -f ${PROJ_DB} - COMMAND ${CMAKE_COMMAND} "-DALL_SQL_IN=${ALL_SQL_IN}" "-DEXE_SQLITE3=${EXE_SQLITE3}" "-DPROJ_DB=${PROJ_DB}" "-DPROJ_VERSION=${PROJ_VERSION}" "-DPROJ_DB_CACHE_DIR=${PROJ_DB_CACHE_DIR}" "-DPROJ_DB_SQL_EXPECTED_MD5=${PROJ_DB_SQL_EXPECTED_MD5}" + COMMAND ${CMAKE_COMMAND} "-DALL_SQL_IN=${ALL_SQL_IN}" "-DEXE_SQLITE3=${EXE_SQLITE3}" "-DPROJ_DB=${PROJ_DB}" "-DPROJ_VERSION=${PROJ_VERSION}" "-DPROJ_DB_CACHE_DIR=${PROJ_DB_CACHE_DIR}" "-DPROJ_DB_SQL_EXPECTED_MD5=${PROJ_DB_SQL_EXPECTED_MD5}" "-DDATA_BINARY_DIR=${CMAKE_CURRENT_BINARY_DIR}" -P "${CMAKE_CURRENT_SOURCE_DIR}/generate_proj_db.cmake" COMMAND ${CMAKE_COMMAND} -E copy ${PROJ_DB} ${CMAKE_CURRENT_BINARY_DIR}/for_tests DEPENDS ${SQL_FILES} "${CMAKE_CURRENT_SOURCE_DIR}/generate_proj_db.cmake" diff --git a/data/generate_proj_db.cmake b/data/generate_proj_db.cmake index 68df7ca687..25a78640f9 100644 --- a/data/generate_proj_db.cmake +++ b/data/generate_proj_db.cmake @@ -28,6 +28,8 @@ endfunction() generate_all_sql_in("${ALL_SQL_IN}" OFF PROJ_DB_SQL_MD5) +file(WRITE "${DATA_BINARY_DIR}/PROJ_DB_SQL_MD5.h" "const char* PROJ_DB_SQL_MD5=\"${PROJ_DB_SQL_MD5}\";\n") + if (NOT "${PROJ_DB_SQL_MD5}" STREQUAL "${PROJ_DB_SQL_EXPECTED_MD5}") message(WARNING "all.sql.in content has changed. Running extra validation checks when building proj.db...") diff --git a/src/embedded_resources.c b/src/embedded_resources.c index 04b721f11e..0dc93956e7 100644 --- a/src/embedded_resources.c +++ b/src/embedded_resources.c @@ -1,7 +1,11 @@ #include "embedded_resources.h" #if USE_SHARP_EMBED + +#include "PROJ_DB_SQL_MD5.h" + const unsigned char *pj_get_embedded_proj_db(unsigned int *pnSize) { + (void)PROJ_DB_SQL_MD5; static const unsigned char proj_db[] = { #embed PROJ_DB }; diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index a36633926c..3e578fdcff 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -425,6 +425,8 @@ if (EMBED_RESOURCE_FILES AND NOT IS_SHARP_EMBED_AVAILABLE_RES) target_sources(proj PRIVATE embedded_resources.c "${EMBEDDED_PROJ_DB}" "${EMBEDDED_PROJ_INI}") elseif(EMBED_RESOURCE_FILES AND IS_SHARP_EMBED_AVAILABLE_RES) add_library(proj_resources OBJECT embedded_resources.c) + target_include_directories(proj_resources PRIVATE "${PROJECT_BINARY_DIR}/data") + set_source_files_properties(embedded_resources.c OBJECT_DEPENDS ${PROJECT_BINARY_DIR}/data/proj.db) target_compile_definitions(proj_resources PRIVATE "PROJ_DB=\"${PROJECT_BINARY_DIR}/data/proj.db\"") target_compile_definitions(proj_resources PRIVATE "PROJ_INI=\"${PROJECT_SOURCE_DIR}/data/proj.ini\"") target_compile_definitions(proj_resources PRIVATE USE_SHARP_EMBED) From 4064893421b861868cf26d60cd934fdee509ac6e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 17 Dec 2024 00:55:55 +0100 Subject: [PATCH 4/4] ConcatenatedOperation::fixStepsDirection(): add even more hacks to deal with EPSG:10675 'BES2020 to Saba height' --- data/CMakeLists.txt | 2 +- data/sql/concatenated_operation.sql | 2 + data/sql/concatenated_operation_step.sql | 2 + data/sql/consistency_checks_triggers.sql | 19 +++ scripts/build_db.py | 7 -- .../operation/concatenatedoperation.cpp | 115 ++++++++++++------ test/unit/test_operation.cpp | 23 ++++ 7 files changed, 128 insertions(+), 42 deletions(-) diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt index ede6cd5232..95ad22fcdc 100644 --- a/data/CMakeLists.txt +++ b/data/CMakeLists.txt @@ -45,7 +45,7 @@ set(ALL_SQL_IN "${CMAKE_CURRENT_BINARY_DIR}/all.sql.in") set(PROJ_DB "${CMAKE_CURRENT_BINARY_DIR}/proj.db") include(sql_filelist.cmake) -set(PROJ_DB_SQL_EXPECTED_MD5 "c2184d9fca631e88f21e4c0d8899d391") +set(PROJ_DB_SQL_EXPECTED_MD5 "99e510ac698e4961c38cc2a09d07ccb0") add_custom_command( OUTPUT ${PROJ_DB} diff --git a/data/sql/concatenated_operation.sql b/data/sql/concatenated_operation.sql index ee621c82b2..bc8b04c986 100644 --- a/data/sql/concatenated_operation.sql +++ b/data/sql/concatenated_operation.sql @@ -495,3 +495,5 @@ INSERT INTO "concatenated_operation" VALUES('EPSG','10496','ETRS89 + DVR90(2013) INSERT INTO "usage" VALUES('EPSG','20516','concatenated_operation','EPSG','10496','EPSG','1080','EPSG','1273'); INSERT INTO "concatenated_operation" VALUES('EPSG','10616','SRGI2013 + INAGeoid2020 v1 height to SRGI2013 + INAGeoid v2 height (1)','In central Java INAGeoid2020 v2 height minus INAGeoid2020 v1 height is approximately +0.2m (v1 surface is above the v2 surface). This difference varies significantly across Indonesia.','EPSG','9529','EPSG','20043',0.2,'BIG-Idn 2022',0); INSERT INTO "usage" VALUES('EPSG','21321','concatenated_operation','EPSG','10616','EPSG','1122','EPSG','1178'); +INSERT INTO "concatenated_operation" VALUES('EPSG','10675','BES2020 to Saba height (1)','Although documented with source and target CRSs in the geog2D domain, step 1 is a 3D transformation in which implicit concatenated operations through geocentric Cartesian coordinates should be used. Refer to EPSG Guidance Notes 7-1 and 7-2.','EPSG','10638','EPSG','10642',0.1,'NSDI-Bes 2023',0); +INSERT INTO "usage" VALUES('EPSG','21877','concatenated_operation','EPSG','10675','EPSG','4757','EPSG','1133'); diff --git a/data/sql/concatenated_operation_step.sql b/data/sql/concatenated_operation_step.sql index 39b1f0a1c8..fd6eda7aec 100644 --- a/data/sql/concatenated_operation_step.sql +++ b/data/sql/concatenated_operation_step.sql @@ -504,3 +504,5 @@ INSERT INTO "concatenated_operation_step" VALUES('EPSG','10496',1,'EPSG','10492' INSERT INTO "concatenated_operation_step" VALUES('EPSG','10496',2,'EPSG','10494'); INSERT INTO "concatenated_operation_step" VALUES('EPSG','10616',1,'EPSG','9629'); INSERT INTO "concatenated_operation_step" VALUES('EPSG','10616',2,'EPSG','10145'); +INSERT INTO "concatenated_operation_step" VALUES('EPSG','10675',1,'EPSG','10646'); +INSERT INTO "concatenated_operation_step" VALUES('EPSG','10675',2,'EPSG','10657'); diff --git a/data/sql/consistency_checks_triggers.sql b/data/sql/consistency_checks_triggers.sql index cb3fbd1d5b..0fc43a5928 100644 --- a/data/sql/consistency_checks_triggers.sql +++ b/data/sql/consistency_checks_triggers.sql @@ -611,6 +611,25 @@ FOR EACH ROW BEGIN AND concat_op.source_crs_auth_name = step_op.target_crs_auth_name AND concat_op.source_crs_code = step_op.target_crs_code) + -- same as above, but check by CRS name, and only for geodetic CRS. + -- For example the concatenated operation EPSG:10675 ("BES2020 to Saba height (1)") + -- has EPSG:10638 "BES2020" (geographic 3D) as source CRS + -- but its first step is EPSG:10646 ("Saba to BES2020 Saba (1)") which has + -- EPSG:10639 "BES2020" (geographic 2D) as target CRS ! + AND NOT EXISTS ( + SELECT 1 FROM coordinate_operation_view step_op + LEFT JOIN concatenated_operation concat_op ON + concat_op.auth_name = NEW.operation_auth_name AND concat_op.code = NEW.operation_code + LEFT JOIN geodetic_crs concat_op_source_crs ON + concat_op_source_crs.auth_name = concat_op.source_crs_auth_name + AND concat_op_source_crs.code = concat_op.source_crs_code + LEFT JOIN geodetic_crs step_op_target_crs ON + step_op_target_crs.auth_name = step_op.target_crs_auth_name + AND step_op_target_crs.code = step_op.target_crs_code + WHERE concat_op.deprecated = 0 + AND step_op.auth_name = NEW.step_auth_name AND step_op.code = NEW.step_code + AND concat_op_source_crs.name = step_op_target_crs.name) + -- or if source_crs of step 1 is a conversion AND NOT EXISTS ( SELECT 1 FROM conversion_table step_op diff --git a/scripts/build_db.py b/scripts/build_db.py index e099cbbcab..6afbf5864e 100755 --- a/scripts/build_db.py +++ b/scripts/build_db.py @@ -930,13 +930,6 @@ def fill_concatenated_operation(proj_db_cursor): expected_order = 1 steps_code = [] - # FIXME: https://epsg.org/concatenated-operation_10675/BES2020-to-Saba-height-1.html ill defined in EPSG 11.023 - # due to first step referencing BES2020 Saba geographic 2D (EPSG:10639), but source CRS of concatenated - # operation referencing BES2020 Saba geographic 3D (EPSG:10638) - if code == 10675: - print("FIXME! Skipping EPSG:10675 'BES2020 to Saba height (1)' for now") - continue - iterator = proj_db_cursor.execute("SELECT op_path_step, single_operation_code FROM epsg_coordoperationpath WHERE concat_operation_code = ? ORDER BY op_path_step", (code,)) for (order, single_operation_code) in iterator: assert order == expected_order diff --git a/src/iso19111/operation/concatenatedoperation.cpp b/src/iso19111/operation/concatenatedoperation.cpp index 0c0a309e10..cd1ccf201b 100644 --- a/src/iso19111/operation/concatenatedoperation.cpp +++ b/src/iso19111/operation/concatenatedoperation.cpp @@ -569,36 +569,73 @@ void ConcatenatedOperation::fixStepsDirection( auto newOp(Conversion::createGeographicGeocentric( NN_NO_CHECK(prevOpTarget), NN_NO_CHECK(l_targetCRS))); operationsInOut.insert(operationsInOut.begin() + i, newOp); - // Particular case for https://github.com/OSGeo/PROJ/issues/3819 - // where the antepenultimate transformation goes to A - // (geographic 3D) // and the last transformation is a NADCON 3D - // but between A (geographic 2D) to B (geographic 2D), and the - // concatenated transformation target CRS is B (geographic 3D) - // This is due to an oddity of the EPSG database that registers - // the NADCON 3D transformation between the 2D geographic CRS - // and not the 3D ones. - } else if (i + 1 == operationsInOut.size() && - l_sourceCRS->nameStr() == prevOpTarget->nameStr() && - l_targetCRS->nameStr() == concatOpTargetCRS->nameStr() && - isGeographic(l_targetCRS.get()) && - isGeographic(concatOpTargetCRS.get()) && - isGeographic(l_sourceCRS.get()) && - isGeographic(prevOpTarget.get()) && - dynamic_cast( - prevOpTarget.get()) - ->coordinateSystem() - ->axisList() - .size() == 3 && - dynamic_cast( - l_sourceCRS.get()) - ->coordinateSystem() - ->axisList() - .size() == 2 && - dynamic_cast( - l_targetCRS.get()) - ->coordinateSystem() - ->axisList() - .size() == 2) { + } else if (i == 0 && + concatOpSourceCRS->nameStr() == l_targetCRS->nameStr() && + isGeographic(concatOpSourceCRS.get()) && + isGeographic(l_sourceCRS.get())) { + // Particular case for EPSG:10675 "BES2020 to Saba height" + op = op->inverse(); + + // This logic to convert between Geog2D <--> Geog3D could + // potentiall by generalized... + operationsInOut.insert( + operationsInOut.begin(), + NN_CHECK_ASSERT( + CoordinateOperationFactory::create()->createOperation( + concatOpSourceCRS, NN_NO_CHECK(l_targetCRS)))); + + // This logic to convert between Geog2D <--> Geog3D could + // potentiall by generalized... + if (operationsInOut.size() >= 3 && + operationsInOut[1]->targetCRS() && + operationsInOut[2]->sourceCRS() && + !areCRSMoreOrLessEquivalent( + operationsInOut[1]->targetCRS().get(), + operationsInOut[2]->sourceCRS().get()) && + operationsInOut[1]->targetCRS()->nameStr() == + operationsInOut[2]->sourceCRS()->nameStr() && + isGeographic(operationsInOut[1]->targetCRS().get()) && + isGeographic(operationsInOut[2]->sourceCRS().get())) { + operationsInOut.insert( + operationsInOut.begin() + 2, + NN_CHECK_ASSERT( + CoordinateOperationFactory::create() + ->createOperation( + NN_NO_CHECK( + operationsInOut[1]->targetCRS()), + NN_NO_CHECK( + operationsInOut[2]->sourceCRS())))); + } + } + // Particular case for https://github.com/OSGeo/PROJ/issues/3819 + // where the antepenultimate transformation goes to A + // (geographic 3D) + // and the last transformation is a NADCON 3D + // but between A (geographic 2D) to B (geographic 2D), and the + // concatenated transformation target CRS is B (geographic 3D) + // This is due to an oddity of the EPSG database that registers + // the NADCON 3D transformation between the 2D geographic CRS + // and not the 3D ones. + else if (i + 1 == operationsInOut.size() && + l_sourceCRS->nameStr() == prevOpTarget->nameStr() && + l_targetCRS->nameStr() == concatOpTargetCRS->nameStr() && + isGeographic(l_targetCRS.get()) && + isGeographic(concatOpTargetCRS.get()) && + isGeographic(l_sourceCRS.get()) && + isGeographic(prevOpTarget.get()) && + dynamic_cast( + prevOpTarget.get()) + ->coordinateSystem() + ->axisList() + .size() == 3 && + dynamic_cast(l_sourceCRS.get()) + ->coordinateSystem() + ->axisList() + .size() == 2 && + dynamic_cast(l_targetCRS.get()) + ->coordinateSystem() + ->axisList() + .size() == 2) { const auto transf = dynamic_cast(op.get()); if (transf && @@ -618,10 +655,20 @@ void ConcatenatedOperation::fixStepsDirection( auto l_sourceCRS = operationsInOut.front()->sourceCRS(); if (l_sourceCRS && !areCRSMoreOrLessEquivalent( l_sourceCRS.get(), concatOpSourceCRS.get())) { - throw InvalidOperation("The source CRS of the first step of " - "concatenated operation is not the same " - "as the source CRS of the concatenated " - "operation itself"); + if (l_sourceCRS->nameStr() == concatOpSourceCRS->nameStr() && + ((isGeographic(l_sourceCRS.get()) && + isGeocentric(concatOpSourceCRS.get())) || + (isGeocentric(l_sourceCRS.get()) && + isGeographic(concatOpSourceCRS.get())))) { + auto newOp(Conversion::createGeographicGeocentric( + NN_NO_CHECK(l_sourceCRS), concatOpSourceCRS)); + operationsInOut.push_back(newOp); + } else { + throw InvalidOperation("The source CRS of the first step of " + "concatenated operation is not the same " + "as the source CRS of the concatenated " + "operation itself"); + } } auto l_targetCRS = operationsInOut.back()->targetCRS(); diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index f2326760ba..c4ac5727fa 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6164,3 +6164,26 @@ TEST(operation, CoordinateFrameRotationFullMatrixGeog2D) { "+step +proj=unitconvert +xy_in=rad +xy_out=deg " "+step +proj=axisswap +order=2,1"); } + +// --------------------------------------------------------------------------- + +TEST(operation, ConcatenatedOperationEPSG_10675) { + auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto op = factory->createCoordinateOperation("10675", false); + auto concatOp = nn_dynamic_pointer_cast(op); + ASSERT_TRUE(concatOp != nullptr); + EXPECT_EQ(concatOp->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_3 " + "+step +proj=cart +ellps=GRS80 " + "+step +inv +proj=helmert +exact " + "+x=1138.7432 +y=-2064.4761 +z=110.7016 " + "+rx=-214.615206 +ry=479.360036 +rz=-164.703951 +s=-402.32073 " + "+convention=coordinate_frame " + "+step +inv +proj=cart +ellps=intl " + "+step +proj=pop +v_3 " + "+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m " + "+step +proj=axisswap +order=2,1"); +}