From 10d8e9af8d3b4b5219d37552fc253ba716436f4a Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 29 May 2020 16:48:36 +0200 Subject: Database: register NZGD2000 -> ITRF96 transformation for NZGD2000 database Related to https://github.com/OSGeo/PROJ-data/pull/22 An entry is added in the ``other_transformation`` table, using a raw PROJ string. Later, once the deformation model has been registered in EPSG, we'll have to add code to map the EPSG transformation method and parameters to +proj=defmodel In the meantime, this works pretty well: ``` $ src/projinfo -s EPSG:4959 -t EPSG:7907 Candidate operations found: 1 ------------------------------------- Operation No. 1: PROJ:NZGD2000-20180701, NZGD2000 to ITRF96 deformation model, unknown accuracy, New Zealand PROJ string: +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg WKT2:2019 string: COORDINATEOPERATION["NZGD2000 to ITRF96 deformation model", VERSION["20180701"], SOURCECRS[ GEOGCRS["NZGD2000", DATUM["New Zealand Geodetic Datum 2000", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,3], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], AXIS["ellipsoidal height (h)",up, ORDER[3], LENGTHUNIT["metre",1]], ID["EPSG",4959]]], TARGETCRS[ GEOGCRS["ITRF96", DATUM["International Terrestrial Reference Frame 1996", ELLIPSOID["GRS 1980",6378137,298.257222101, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,3], AXIS["geodetic latitude (Lat)",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433]], AXIS["ellipsoidal height (h)",up, ORDER[3], LENGTHUNIT["metre",1]], ID["EPSG",7907]]], METHOD["PROJ-based operation method: +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg"], USAGE[ SCOPE["unknown"], AREA["New Zealand"], BBOX[-55.95,160.6,-25.88,-171.2]], ID["PROJ","NZGD2000-20180701"], REMARK["New Zealand Deformation Model. Defines the secular model (National Deformation Model) and patches for significant deformation events since 2000"]] ``` ``` $ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cct -d 8 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg -40.99999402 172.99999938 0.00130333 2016.5000 ``` ``` $ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" EPSG:4959 EPSG:7907 -40.99999402 172.99999938 0.00130333 2016.5 ``` --- data/sql/customizations.sql | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql index efb0a4bc..1556bf62 100644 --- a/data/sql/customizations.sql +++ b/data/sql/customizations.sql @@ -25,6 +25,31 @@ INSERT INTO "geodetic_crs" VALUES('OGC','CRS84','WGS 84 (CRS84)',NULL,NULL,'geog INSERT INTO "other_transformation" VALUES('PROJ','CRS84_TO_EPSG_4326','OGC:CRS84 to WGS 84',NULL,NULL,'EPSG','9843','Axis Order Reversal (2D)','OGC','CRS84','EPSG','4326','EPSG','1262',0.0,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); +-- Temporary entry for NZGD2000 deformation model +INSERT INTO other_transformation VALUES( + 'PROJ','NZGD2000-20180701','NZGD2000 to ITRF96', + 'New Zealand Deformation Model. Defines the secular model (National Deformation Model) and patches for significant deformation events since 2000', NULL, + 'PROJ', 'PROJString', + '+proj=pipeline ' || + '+step +proj=unitconvert +xy_in=deg +xy_out=rad ' || + '+step +proj=axisswap +order=2,1 ' || + '+step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json ' || + '+step +proj=axisswap +order=2,1 ' || + '+step +proj=unitconvert +xy_in=rad +xy_out=deg', + 'EPSG','4959', + 'EPSG','7907', + 'EPSG','1175', + NULL, --accuracy + NULL,NULL,NULL,NULL,NULL,NULL, -- param1 + NULL,NULL,NULL,NULL,NULL,NULL, -- param2 + NULL,NULL,NULL,NULL,NULL,NULL, -- param3 + NULL,NULL,NULL,NULL,NULL,NULL, -- param4 + NULL,NULL,NULL,NULL,NULL,NULL, -- param5 + NULL,NULL,NULL,NULL,NULL,NULL, -- param6 + NULL,NULL,NULL,NULL,NULL,NULL, -- param7 + '20180701', -- operation version + 0); + -- alias of EPSG:3857 INSERT INTO "projected_crs" VALUES('EPSG','900913','Google Maps Global Mercator',NULL,NULL,'EPSG','4499','EPSG','4326','EPSG','3856','EPSG','3544',NULL,1); -- cgit v1.2.3 From ad9efb19d37be616d69676e732abe44e78a8305b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 29 May 2020 22:53:58 +0200 Subject: Database: add concatenated operations for NZGD2000 -> ITRF97/2000/2005/2008/2014 --- data/sql/customizations.sql | 20 ++++++++++++++++++ src/iso19111/coordinateoperation.cpp | 41 ++++++++++++++++++++++-------------- test/cli/testprojinfo | 28 ++++++++++++++++++++++++ test/cli/testprojinfo_out.dist | 18 ++++++++++++++++ 4 files changed, 91 insertions(+), 16 deletions(-) diff --git a/data/sql/customizations.sql b/data/sql/customizations.sql index 1556bf62..86baf691 100644 --- a/data/sql/customizations.sql +++ b/data/sql/customizations.sql @@ -50,6 +50,26 @@ INSERT INTO other_transformation VALUES( '20180701', -- operation version 0); +INSERT INTO "concatenated_operation" VALUES('PROJ','NZGD2000_TO_ITRF97','NZGD2000 to ITRF97','Concatenation of PROJ:NZGD2000-20180701 and 9079','','EPSG','4959','EPSG','7908','EPSG','1175',NULL,NULL,0); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF97',1,'PROJ','NZGD2000-20180701'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF97',2,'EPSG','9079'); + +INSERT INTO "concatenated_operation" VALUES('PROJ','NZGD2000_TO_ITRF2000','NZGD2000 to ITRF2000','Concatenation of PROJ:NZGD2000-20180701 and 9080','','EPSG','4959','EPSG','7909','EPSG','1175',NULL,NULL,0); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2000',1,'PROJ','NZGD2000-20180701'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2000',2,'EPSG','9080'); + +INSERT INTO "concatenated_operation" VALUES('PROJ','NZGD2000_TO_ITRF2005','NZGD2000 to ITRF2005','Concatenation of PROJ:NZGD2000-20180701 and 9081','','EPSG','4959','EPSG','7910','EPSG','1175',NULL,NULL,0); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2005',1,'PROJ','NZGD2000-20180701'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2005',2,'EPSG','9081'); + +INSERT INTO "concatenated_operation" VALUES('PROJ','NZGD2000_TO_ITRF2008','NZGD2000 to ITRF2008','Concatenation of PROJ:NZGD2000-20180701 and EPSG:9082','','EPSG','4959','EPSG','7911','EPSG','1175',NULL,NULL,0); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2008',1,'PROJ','NZGD2000-20180701'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2008',2,'EPSG','9082'); + +INSERT INTO "concatenated_operation" VALUES('PROJ','NZGD2000_TO_ITRF2014','NZGD2000 to ITRF2014','Concatenation of PROJ:NZGD2000-20180701 and EPSG:9083','','EPSG','4959','EPSG','7912','EPSG','1175',NULL,NULL,0); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2014',1,'PROJ','NZGD2000-20180701'); +INSERT INTO "concatenated_operation_step" VALUES('PROJ','NZGD2000_TO_ITRF2014',2,'EPSG','9083'); + -- alias of EPSG:3857 INSERT INTO "projected_crs" VALUES('EPSG','900913','Google Maps Global Mercator',NULL,NULL,'EPSG','4499','EPSG','4326','EPSG','3856','EPSG','3544',NULL,1); diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index bdb2ad2e..f2f674ce 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -10119,6 +10119,17 @@ void ConcatenatedOperation::fixStepsDirection( // Set of heuristics to assign CRS to steps, and possibly reverse them. + const auto isGeographic = [](const crs::CRS *crs) -> bool { + return dynamic_cast(crs) != nullptr; + }; + + const auto isGeocentric = [](const crs::CRS *crs) -> bool { + auto geodCRS = dynamic_cast(crs); + if (geodCRS && geodCRS->coordinateSystem()->axisList().size() == 3) + return true; + return false; + }; + for (size_t i = 0; i < operationsInOut.size(); ++i) { auto &op = operationsInOut[i]; auto l_sourceCRS = op->sourceCRS(); @@ -10209,18 +10220,6 @@ void ConcatenatedOperation::fixStepsDirection( } } else if (!conv && l_sourceCRS && l_targetCRS) { - const auto isGeographic = [](const crs::CRS *crs) -> bool { - return dynamic_cast(crs) != nullptr; - }; - - const auto isGeocentric = [](const crs::CRS *crs) -> bool { - auto geodCRS = dynamic_cast(crs); - if (geodCRS && - geodCRS->coordinateSystem()->axisList().size() == 3) - return true; - return false; - }; - // Transformations might be mentioned in their forward directions, // whereas we should instead use the reverse path. auto prevOpTarget = (i == 0) ? concatOpSourceCRS.as_nullable() @@ -10265,10 +10264,20 @@ void ConcatenatedOperation::fixStepsDirection( auto l_targetCRS = operationsInOut.back()->targetCRS(); if (l_targetCRS && !compareStepCRS(l_targetCRS.get(), concatOpTargetCRS.get())) { - throw InvalidOperation("The target CRS of the last step of " - "concatenated operation is not the same " - "as the target CRS of the concatenated " - "operation itself"); + if (l_targetCRS->nameStr() == concatOpTargetCRS->nameStr() && + ((isGeographic(l_targetCRS.get()) && + isGeocentric(concatOpTargetCRS.get())) || + (isGeocentric(l_targetCRS.get()) && + isGeographic(concatOpTargetCRS.get())))) { + auto newOp(Conversion::createGeographicGeocentric( + NN_NO_CHECK(l_targetCRS), concatOpTargetCRS)); + operationsInOut.push_back(newOp); + } else { + throw InvalidOperation("The target CRS of the last step of " + "concatenated operation is not the same " + "as the target CRS of the concatenated " + "operation itself"); + } } } } diff --git a/test/cli/testprojinfo b/test/cli/testprojinfo index df8cac8c..9896e23a 100755 --- a/test/cli/testprojinfo +++ b/test/cli/testprojinfo @@ -201,6 +201,34 @@ fi rm testprojinfo_out_remotedata.txt unset PROJ_NETWORK +###################### +# NZGD2000 -> ITRFxx # +###################### +echo 'Testing -s NZGD2000 -t ITRF96 -o PROJ -q' >> ${OUT} +$EXE -s NZGD2000 -t ITRF96 -o PROJ -q >>${OUT} 2>&1 +echo "" >>${OUT} + +echo 'Testing -s NZGD2000 -t ITRF97 -o PROJ -q' >> ${OUT} +$EXE -s NZGD2000 -t ITRF97 -o PROJ -q >>${OUT} 2>&1 +echo "" >>${OUT} + +echo 'Testing -s NZGD2000 -t ITRF2000 -o PROJ -q' >> ${OUT} +$EXE -s NZGD2000 -t ITRF2000 -o PROJ -q >>${OUT} 2>&1 +echo "" >>${OUT} + +echo 'Testing -s NZGD2000 -t ITRF2005 -o PROJ -q' >> ${OUT} +$EXE -s NZGD2000 -t ITRF2005 -o PROJ -q >>${OUT} 2>&1 +echo "" >>${OUT} + +echo 'Testing -s NZGD2000 -t ITRF2008 -o PROJ -q' >> ${OUT} +$EXE -s NZGD2000 -t ITRF2008 -o PROJ -q >>${OUT} 2>&1 +echo "" >>${OUT} + +echo 'Testing -s NZGD2000 -t ITRF2014 -o PROJ -q' >> ${OUT} +$EXE -s NZGD2000 -t ITRF2014 -o PROJ -q >>${OUT} 2>&1 +echo "" >>${OUT} + + # do 'diff' with distribution results echo "diff ${OUT} with testprojinfo_out.dist" diff -u ${OUT} ${TEST_CLI_DIR}/testprojinfo_out.dist diff --git a/test/cli/testprojinfo_out.dist b/test/cli/testprojinfo_out.dist index 9e03b003..f0434da9 100644 --- a/test/cli/testprojinfo_out.dist +++ b/test/cli/testprojinfo_out.dist @@ -1128,3 +1128,21 @@ DATUM["World Geodetic System 1984", Testing -k operation EPSG:8457 -o PROJ -q +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=bessel +step +proj=helmert +x=674.374 +y=15.056 +z=405.346 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 +Testing -s NZGD2000 -t ITRF96 -o PROJ -q ++proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg + +Testing -s NZGD2000 -t ITRF97 -o PROJ -q ++proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=0 +y=-0.00051 +z=0.01553 +rx=-0.00016508 +ry=0.00026897 +rz=5.984e-05 +s=-0.00151099 +dx=0.00069 +dy=-0.0001 +dz=0.00186 +drx=-1.347e-05 +dry=1.514e-05 +drz=-2.7e-07 +ds=-0.00019201 +t_epoch=2000 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 + +Testing -s NZGD2000 -t ITRF2000 -o PROJ -q ++proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=0.0067 +y=0.00379 +z=-0.00717 +rx=-0.00016508 +ry=0.00026897 +rz=0.00011984 +s=6.901e-05 +dx=0.00069 +dy=-0.0007 +dz=0.00046 +drx=-1.347e-05 +dry=1.514e-05 +drz=1.973e-05 +ds=-0.00018201 +t_epoch=2000 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 + +Testing -s NZGD2000 -t ITRF2005 -o PROJ -q ++proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=0.0068 +y=0.00299 +z=-0.01297 +rx=-0.00016508 +ry=0.00026897 +rz=0.00011984 +s=0.00046901 +dx=0.00049 +dy=-0.0006 +dz=-0.00134 +drx=-1.347e-05 +dry=1.514e-05 +drz=1.973e-05 +ds=-0.00010201 +t_epoch=2000 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 + +Testing -s NZGD2000 -t ITRF2008 -o PROJ -q ++proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=0.0048 +y=0.00209 +z=-0.01767 +rx=-0.00016508 +ry=0.00026897 +rz=0.00011984 +s=0.00140901 +dx=0.00079 +dy=-0.0006 +dz=-0.00134 +drx=-1.347e-05 +dry=1.514e-05 +drz=1.973e-05 +ds=-0.00010201 +t_epoch=2000 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 + +Testing -s NZGD2000 -t ITRF2014 -o PROJ -q ++proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=cart +ellps=GRS80 +step +inv +proj=helmert +x=0.0064 +y=0.00399 +z=-0.01427 +rx=-0.00016508 +ry=0.00026897 +rz=0.00011984 +s=0.00108901 +dx=0.00079 +dy=-0.0006 +dz=-0.00144 +drx=-1.347e-05 +dry=1.514e-05 +drz=1.973e-05 +ds=-7.201e-05 +t_epoch=2000 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 + -- cgit v1.2.3