diff --git a/src/iso19111/operation/coordinateoperationfactory.cpp b/src/iso19111/operation/coordinateoperationfactory.cpp index 0b21343ca9..98d4a7e32b 100644 --- a/src/iso19111/operation/coordinateoperationfactory.cpp +++ b/src/iso19111/operation/coordinateoperationfactory.cpp @@ -1671,6 +1671,10 @@ void CoordinateOperationFactory::Private::buildCRSIds( std::vector allowedObjects; auto geogCRS = dynamic_cast(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 @@ -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() + diff --git a/test/unit/test_operationfactory.cpp b/test/unit/test_operationfactory.cpp index 08ce7bb341..ad5f684e15 100644 --- a/test/unit/test_operationfactory.cpp +++ b/test/unit/test_operationfactory.cpp @@ -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(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(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");