Skip to content

Commit 1e3d0f2

Browse files
rouaultrouault
authored
Merge pull request #4577 from rouault/fix_4576
createOperations(): fix compound to compound when vertical CRS definition is equivalent but not strictly identical
2 parents 58799b3 + 9babeda commit 1e3d0f2

File tree

2 files changed

+136
-6
lines changed

2 files changed

+136
-6
lines changed

‎src/iso19111/operation/coordinateoperationfactory.cpp‎

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3714,7 +3714,8 @@ CoordinateOperationFactory::Private::createOperations(
37143714
}
37153715

37163716
if (dynamic_cast<const crs::EngineeringCRS *>(sourceCRS.get()) &&
3717-
sourceCRS->_isEquivalentTo(targetCRS.get())) {
3717+
sourceCRS->_isEquivalentTo(targetCRS.get(),
3718+
util::IComparable::Criterion::EQUIVALENT)) {
37183719
std::string name("Identity transformation from ");
37193720
name += sourceCRS->nameStr();
37203721
name += " to ";
@@ -3845,7 +3846,9 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
38453846
if (pmoDst.size() == pmoSrc.size()) {
38463847
bool ok = true;
38473848
for (size_t i = 0; i < pmoSrc.size(); ++i) {
3848-
if (pmoSrc[i]->_isEquivalentTo(pmoDst[i].get())) {
3849+
if (pmoSrc[i]->_isEquivalentTo(
3850+
pmoDst[i].get(),
3851+
util::IComparable::Criterion::EQUIVALENT)) {
38493852
auto pmo = pmoSrc[i]->cloneWithEpochs(*sourceEpoch,
38503853
*targetEpoch);
38513854
std::vector<operation::CoordinateOperationNNPtr> ops;
@@ -6629,7 +6632,9 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
66296632
const auto vertSrc = componentsSrc[1]->extractVerticalCRS();
66306633
const auto vertDst = componentsDst[1]->extractVerticalCRS();
66316634
if (vertSrc && vertDst &&
6632-
!componentsSrc[1]->_isEquivalentTo(componentsDst[1].get())) {
6635+
!componentsSrc[1]->_isEquivalentTo(
6636+
componentsDst[1].get(),
6637+
util::IComparable::Criterion::EQUIVALENT)) {
66336638
if ((!vertSrc->geoidModel().empty() ||
66346639
!vertDst->geoidModel().empty()) &&
66356640
// To be able to use "CGVD28 height to
@@ -6966,7 +6971,8 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
69666971
dynamic_cast<crs::GeographicCRS *>(
69676972
compSrc0BoundCrs->hubCRS().get()) &&
69686973
compSrc0BoundCrs->hubCRS()->_isEquivalentTo(
6969-
compDst0BoundCrs->hubCRS().get())) {
6974+
compDst0BoundCrs->hubCRS().get(),
6975+
util::IComparable::Criterion::EQUIVALENT)) {
69706976
interpolationCRS =
69716977
NN_NO_CHECK(util::nn_dynamic_pointer_cast<crs::SingleCRS>(
69726978
compSrc0BoundCrs->hubCRS()));
@@ -7155,7 +7161,8 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound(
71557161
const auto compDst0BoundCrsHubAsGeogCRSDatum =
71567162
compDst0BoundCrsHubAsGeogCRS->datumNonNull(dbContext);
71577163
if (boundSrcHubAsGeogCRSDatum->_isEquivalentTo(
7158-
compDst0BoundCrsHubAsGeogCRSDatum.get())) {
7164+
compDst0BoundCrsHubAsGeogCRSDatum.get(),
7165+
util::IComparable::Criterion::EQUIVALENT)) {
71597166
auto cs = cs::EllipsoidalCS::
71607167
createLatitudeLongitudeEllipsoidalHeight(
71617168
common::UnitOfMeasure::DEGREE,
@@ -7181,7 +7188,8 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound(
71817188
for (const auto &opDst : geog3DToTargetOps) {
71827189
if (opSrc->targetCRS() && opDst->sourceCRS() &&
71837190
!opSrc->targetCRS()->_isEquivalentTo(
7184-
opDst->sourceCRS().get())) {
7191+
opDst->sourceCRS().get(),
7192+
util::IComparable::Criterion::EQUIVALENT)) {
71857193
// Shouldn't happen normally, but typically
71867194
// one of them can be 2D and the other 3D
71877195
// due to above createOperations() not

‎test/unit/test_operationfactory.cpp‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6037,6 +6037,128 @@ TEST(operation, compoundCRS_to_compoundCRS_GDA94_AHD_to_GDA2020_AVWS) {
60376037

60386038
// ---------------------------------------------------------------------------
60396039

6040+
TEST(operation,
6041+
compoundCRS_with_horizontal_boundCRS_to_compoundCRS_and_same_vertical) {
6042+
6043+
auto dbContext = DatabaseContext::create();
6044+
auto objSrc = WKTParser().attachDatabaseContext(dbContext).createFromWKT(
6045+
"COMPOUNDCRS[\"AGD66 / AMG zone 55 + AHD height\",\n"
6046+
" BOUNDCRS[\n"
6047+
" SOURCECRS[\n"
6048+
" PROJCRS[\"AGD66 / AMG zone 55 + AHD height\",\n"
6049+
" BASEGEOGCRS[\"AGD66\",\n"
6050+
" DATUM[\"User datum (no grid)\",\n"
6051+
" ELLIPSOID[\"Australian National "
6052+
"Spheroid\",6378160,298.25,\n"
6053+
" LENGTHUNIT[\"metre\",1]]],\n"
6054+
" PRIMEM[\"Greenwich\",0,\n"
6055+
" ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
6056+
" CONVERSION[\"Australian Map Grid zone 55\",\n"
6057+
" METHOD[\"Transverse Mercator\"],\n"
6058+
" PARAMETER[\"Latitude of natural origin\",0,\n"
6059+
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
6060+
" PARAMETER[\"Longitude of natural origin\",147,\n"
6061+
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
6062+
" PARAMETER[\"Scale factor at natural "
6063+
"origin\",0.9996,\n"
6064+
" SCALEUNIT[\"unity\",1]],\n"
6065+
" PARAMETER[\"False easting\",500000,\n"
6066+
" LENGTHUNIT[\"metre\",1]],\n"
6067+
" PARAMETER[\"False northing\",10000000,\n"
6068+
" LENGTHUNIT[\"metre\",1]]],\n"
6069+
" CS[Cartesian,2],\n"
6070+
" AXIS[\"(E)\",east,\n"
6071+
" ORDER[1],\n"
6072+
" LENGTHUNIT[\"metre\",1]],\n"
6073+
" AXIS[\"(N)\",north,\n"
6074+
" ORDER[2],\n"
6075+
" LENGTHUNIT[\"metre\",1]],\n"
6076+
" USAGE[\n"
6077+
" SCOPE[\"Engineering survey, topographic "
6078+
"mapping.\"],\n"
6079+
" AREA[\"Australia - onshore and offshore between "
6080+
"144°E and 150°E. Papua New Guinea (PNG) - onshore between 144°E and "
6081+
"150°E.\"],\n"
6082+
" BBOX[-47.2,144,-1.3,150.01]]]],\n"
6083+
" TARGETCRS[\n"
6084+
" GEOGCRS[\"GDA2020\",\n"
6085+
" DATUM[\"Geocentric Datum of Australia 2020\",\n"
6086+
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
6087+
" LENGTHUNIT[\"metre\",1]]],\n"
6088+
" PRIMEM[\"Greenwich\",0,\n"
6089+
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
6090+
" CS[ellipsoidal,2],\n"
6091+
" AXIS[\"geodetic latitude (Lat)\",north,\n"
6092+
" ORDER[1],\n"
6093+
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
6094+
" AXIS[\"geodetic longitude (Lon)\",east,\n"
6095+
" ORDER[2],\n"
6096+
" ANGLEUNIT[\"degree\",0.0174532925199433]]]],\n"
6097+
" ABRIDGEDTRANSFORMATION[\"User-defined Bursa-Wolf "
6098+
"Transform\",\n"
6099+
" METHOD[\"Position Vector transformation (geog2D "
6100+
"domain)\"],\n"
6101+
" PARAMETER[\"X-axis translation\",-153.501992013804,\n"
6102+
" ID[\"EPSG\",8605]],\n"
6103+
" PARAMETER[\"Y-axis translation\",91.8979603422544,\n"
6104+
" ID[\"EPSG\",8606]],\n"
6105+
" PARAMETER[\"Z-axis translation\",-76.7203604254192,\n"
6106+
" ID[\"EPSG\",8607]],\n"
6107+
" PARAMETER[\"X-axis rotation\",-11.9539464931016,\n"
6108+
" ID[\"EPSG\",8608]],\n"
6109+
" PARAMETER[\"Y-axis rotation\",13.9315768664084,\n"
6110+
" ID[\"EPSG\",8609]],\n"
6111+
" PARAMETER[\"Z-axis rotation\",-3.45993996519333,\n"
6112+
" ID[\"EPSG\",8610]],\n"
6113+
" PARAMETER[\"Scale difference\",0.999965937122937,\n"
6114+
" ID[\"EPSG\",8611]]]],\n"
6115+
" VERTCRS[\"AHD height\",\n"
6116+
" VDATUM[\"Australian Height Datum\"],\n"
6117+
" CS[vertical,1],\n"
6118+
" AXIS[\"gravity-related height (H)\",up,\n"
6119+
" LENGTHUNIT[\"metre\",1]],\n"
6120+
" USAGE[\n"
6121+
" SCOPE[\"Cadastre, engineering surveying applications over "
6122+
"distances up to 10km.\"],\n"
6123+
" AREA[\"Australia - Australian Capital Territory, New "
6124+
"South Wales, Northern Territory, Queensland, South Australia, "
6125+
"Tasmania, Western Australia and Victoria - onshore. Christmas Island "
6126+
"- onshore. Cocos and Keeling Islands - onshore.\"],\n"
6127+
" BBOX[-43.7,96.76,-9.86,153.69]],\n"
6128+
" ID[\"EPSG\",5711]]]");
6129+
auto src = nn_dynamic_pointer_cast<CRS>(objSrc);
6130+
ASSERT_TRUE(src != nullptr);
6131+
6132+
// GDA2020 / MGA zone 55 + AHD height
6133+
auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");
6134+
auto dstObj = createFromUserInput("EPSG:7855+5711", dbContext, false);
6135+
auto dst = nn_dynamic_pointer_cast<CRS>(dstObj);
6136+
ASSERT_TRUE(dst != nullptr);
6137+
6138+
auto ctxt = CoordinateOperationContext::create(
6139+
AuthorityFactory::create(dbContext, std::string()), nullptr, 0.0);
6140+
ctxt->setSpatialCriterion(
6141+
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
6142+
auto list = CoordinateOperationFactory::create()->createOperations(
6143+
NN_NO_CHECK(src), NN_NO_CHECK(dst), ctxt);
6144+
ASSERT_GE(list.size(), 1U);
6145+
6146+
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
6147+
"+proj=pipeline "
6148+
"+step +inv +proj=utm +zone=55 +south +ellps=aust_SA "
6149+
"+step +proj=push +v_3 "
6150+
"+step +proj=cart +ellps=aust_SA "
6151+
"+step +proj=helmert +x=-153.501992013804 +y=91.8979603422544 "
6152+
"+z=-76.7203604254192 +rx=-11.9539464931016 +ry=13.9315768664084 "
6153+
"+rz=-3.45993996519333 +s=-34.0628770629792 "
6154+
"+convention=position_vector "
6155+
"+step +inv +proj=cart +ellps=GRS80 "
6156+
"+step +proj=pop +v_3 "
6157+
"+step +proj=utm +zone=55 +south +ellps=GRS80");
6158+
}
6159+
6160+
// ---------------------------------------------------------------------------
6161+
60406162
TEST(
60416163
operation,
60426164
compoundCRS_to_compoundCRS_concatenated_operation_with_two_vert_transformation) {

0 commit comments

Comments
 (0)