Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Backport 9.3] createOperations(): fix issue with a obscure case involving CompoundCRS of unknown horizontal datum + boundCRS of vertical #3934

Merged
merged 1 commit into from
Oct 30, 2023
Merged

[Backport 9.3] createOperations(): fix issue with a obscure case involving CompoundCRS of unknown horizontal datum + boundCRS of vertical #3934

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 12 additions & 6 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1671,6 +1671,10 @@ void CoordinateOperationFactory::Private::buildCRSIds(
std::vector<io::AuthorityFactory::ObjectType> allowedObjects;
auto geogCRS = dynamic_cast<const crs::GeographicCRS *>(crs.get());
if (geogCRS) {
if (geogCRS->datumNonNull(authFactory->databaseContext())
->nameStr() == "unknown") {
return;
}
allowedObjects.push_back(
geogCRS->coordinateSystem()->axisList().size() == 2
? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS
Expand Down Expand Up @@ -5765,14 +5769,16 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
}
}

const bool hasOnlyOneOp =
srcToInterpOps.size() == 1 && interpToTargetOps.size() == 1;
for (const auto &srcToInterp : srcToInterpOps) {
for (const auto &interpToTarget : interpToTargetOps) {

if ((srcAndTargetGeogAreSame &&
mapSetDatumsUsed[srcToInterp.get()] !=
mapSetDatumsUsed[interpToTarget.get()]) ||
!useCompatibleTransformationsForSameSourceTarget(
srcToInterp, interpToTarget)) {
if (!hasOnlyOneOp &&
((srcAndTargetGeogAreSame &&
mapSetDatumsUsed[srcToInterp.get()] !=
mapSetDatumsUsed[interpToTarget.get()]) ||
!useCompatibleTransformationsForSameSourceTarget(
srcToInterp, interpToTarget))) {
#ifdef TRACE_CREATE_OPERATIONS
logTrace(
"Considering that '" + srcToInterp->nameStr() +
Expand Down
176 changes: 176 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4690,6 +4690,182 @@ TEST(operation,

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

TEST(operation, compoundCRS_with_vert_bound_to_bound_geog3D) {
// Test case of https://github.com/OSGeo/PROJ/issues/3927

auto dbContext = DatabaseContext::create();

const char *wktSrc =
"COMPOUNDCRS[\"ENU (-77.410692720411:39.4145340892321) + EGM96 geoid "
"height\",\n"
" PROJCRS[\"ENU (-77.410692720411:39.4145340892321)\",\n"
" BASEGEOGCRS[\"WGS 84\",\n"
" DATUM[\"unknown\",\n"
" ELLIPSOID[\"WGS84\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n"
" CONVERSION[\"unnamed\",\n"
" METHOD[\"Orthographic\",\n"
" ID[\"EPSG\",9840]],\n"
" PARAMETER[\"Latitude of natural "
"origin\",39.4145340892321,\n"
" ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural "
"origin\",-77.410692720411,\n"
" ANGLEUNIT[\"Degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"False easting\",0,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",0,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8807]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"easting\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]],\n"
" AXIS[\"northing\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]]],\n"
" BOUNDCRS[\n"
" SOURCECRS[\n"
" VERTCRS[\"EGM96 geoid height\",\n"
" VDATUM[\"EGM96 geoid\"],\n"
" CS[vertical,1],\n"
" AXIS[\"up\",up,\n"
" LENGTHUNIT[\"m\",1]],\n"
" ID[\"EPSG\",5773]]],\n"
" TARGETCRS[\n"
" GEOGCRS[\"WGS 84\",\n"
" DATUM[\"World Geodetic System 1984\",\n"
" ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" CS[ellipsoidal,3],\n"
" AXIS[\"geodetic latitude (Lat)\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"geodetic longitude (Lon)\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"ellipsoidal height (h)\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",4979]]],\n"
" ABRIDGEDTRANSFORMATION[\"EGM96 geoid height to WGS 84 "
"ellipsoidal height\",\n"
" METHOD[\"GravityRelatedHeight to Geographic3D\"],\n"
" PARAMETERFILE[\"Geoid (height correction) model "
"file\",\"egm96_15.gtx\",\n"
" ID[\"EPSG\",8666]]]]]";
auto objSrc =
WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc);
auto srcCRS = nn_dynamic_pointer_cast<CompoundCRS>(objSrc);
ASSERT_TRUE(srcCRS != nullptr);

const char *wktDst =
"BOUNDCRS[\n"
" SOURCECRS[\n"
" GEOGCRS[\"WGS84 Coordinate System\",\n"
" DATUM[\"World Geodetic System 1984\",\n"
" ELLIPSOID[\"WGS 1984\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",6326]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" CS[ellipsoidal,3],\n"
" AXIS[\"geodetic latitude (Lat)\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"geodetic longitude (Lon)\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"ellipsoidal height (h)\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]],\n"
" REMARK[\"Promoted to 3D from EPSG:4326\"]]],\n"
" TARGETCRS[\n"
" GEOGCRS[\"WGS 84\",\n"
" ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n"
" MEMBER[\"World Geodetic System 1984 (Transit)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G730)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G873)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1150)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1674)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G1762)\"],\n"
" MEMBER[\"World Geodetic System 1984 (G2139)\"],\n"
" ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ENSEMBLEACCURACY[2.0]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" CS[ellipsoidal,3],\n"
" AXIS[\"geodetic latitude (Lat)\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"geodetic longitude (Lon)\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"ellipsoidal height (h)\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" USAGE[\n"
" SCOPE[\"Geodesy. Navigation and positioning using GPS "
"satellite system.\"],\n"
" AREA[\"World.\"],\n"
" BBOX[-90,-180,90,180]],\n"
" ID[\"EPSG\",4979]]],\n"
" ABRIDGEDTRANSFORMATION[\"Transformation from WGS84 Coordinate "
"System to WGS84\",\n"
" METHOD[\"Position Vector transformation (geog2D domain)\",\n"
" ID[\"EPSG\",9606]],\n"
" PARAMETER[\"X-axis translation\",0,\n"
" ID[\"EPSG\",8605]],\n"
" PARAMETER[\"Y-axis translation\",0,\n"
" ID[\"EPSG\",8606]],\n"
" PARAMETER[\"Z-axis translation\",0,\n"
" ID[\"EPSG\",8607]],\n"
" PARAMETER[\"X-axis rotation\",0,\n"
" ID[\"EPSG\",8608]],\n"
" PARAMETER[\"Y-axis rotation\",0,\n"
" ID[\"EPSG\",8609]],\n"
" PARAMETER[\"Z-axis rotation\",0,\n"
" ID[\"EPSG\",8610]],\n"
" PARAMETER[\"Scale difference\",1,\n"
" ID[\"EPSG\",8611]]]]";
auto objDst =
WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst);
auto dstCRS = nn_dynamic_pointer_cast<CRS>(objDst);
ASSERT_TRUE(dstCRS != nullptr);

auto authFactory = AuthorityFactory::create(dbContext, std::string());
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt);
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +inv +proj=ortho +lat_0=39.4145340892321 "
"+lon_0=-77.410692720411 +x_0=0 +y_0=0 +ellps=WGS84 "
"+step +proj=vgridshift +grids=egm96_15.gtx +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

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

TEST(operation, geocent_to_compoundCRS) {
auto objSrc = PROJStringParser().createFromPROJString(
"+proj=geocent +datum=WGS84 +units=m +type=crs");
Expand Down