Skip to content

Commit

Permalink
createOperations(): fix issue with a obscure case involving CompoundC…
Browse files Browse the repository at this point in the history 
…RS of unknown horizontal datum + boundCRS of vertical (fixes #3927)
  • Loading branch information
@rouault @github-actions
rouault authored and github-actions[bot] committed Oct 27, 2023
1 parent 2245779 commit 44b6fc3
Show file tree
Hide file tree
Showing 2 changed files with 188 additions and 6 deletions.
18 changes: 12 additions & 6 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
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
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

0 comments on commit 44b6fc3

Please sign in to comment.