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

createOperations(): fix issues when transforming between Geog3D and DerivedGeog3D CRS with Geographic3D offsets method #3411

Merged
merged 1 commit into from
Oct 23, 2022
Merged

createOperations(): fix issues when transforming between Geog3D and DerivedGeog3D CRS with Geographic3D offsets method #3411

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
8 changes: 8 additions & 0 deletions src/iso19111/io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3201,6 +3201,14 @@ CRSNNPtr WKTParser::Private::buildDerivedGeodeticCRS(const WKTNodeNNPtr &node) {

auto ellipsoidalCS = nn_dynamic_pointer_cast<EllipsoidalCS>(cs);
if (ellipsoidalCS) {

if (ellipsoidalCS->axisList().size() == 3 &&
baseGeodCRS->coordinateSystem()->axisList().size() == 2) {
baseGeodCRS =
NN_NO_CHECK(util::nn_dynamic_pointer_cast<GeodeticCRS>(
baseGeodCRS->promoteTo3D(std::string(), dbContext_)));
}

return DerivedGeographicCRS::create(buildProperties(node), baseGeodCRS,
derivingConversion,
NN_NO_CHECK(ellipsoidalCS));
Expand Down
6 changes: 5 additions & 1 deletion src/iso19111/operation/conversion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3465,12 +3465,16 @@ void Conversion::_exportToPROJString(
methodEPSGCode == EPSG_CODE_METHOD_AFFINE_PARAMETRIC_TRANSFORMATION;
const bool isGeographicGeocentric =
methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_GEOCENTRIC;
const bool isGeographicOffsets =
methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC2D_OFFSETS ||
methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC3D_OFFSETS ||
methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC2D_WITH_HEIGHT_OFFSETS;
const bool isHeightDepthReversal =
methodEPSGCode == EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL;
const bool applySourceCRSModifiers =
!isZUnitConversion && !isAffineParametric &&
!isAxisOrderReversal(methodEPSGCode) && !isGeographicGeocentric &&
!isHeightDepthReversal;
!isGeographicOffsets && !isHeightDepthReversal;
bool applyTargetCRSModifiers = applySourceCRSModifiers;

if (formatter->getCRSExport()) {
Expand Down
67 changes: 67 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 @@ -635,6 +635,73 @@ TEST(operation, geogCRS_to_geogCRS_context_datum_ensemble) {

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

TEST(operation, geogCRS_to_derived_geogCRS_3D) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);

auto dst_wkt =
"GEOGCRS[\"CH1903+ with 10m offset on ellipsoidal height\",\n"
" BASEGEOGCRS[\"CH1903+\",\n"
" DATUM[\"CH1903+\",\n"
" ELLIPSOID[\"Bessel 1841\",6377397.155,299.1528128,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",6150]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8901]]],\n"
" DERIVINGCONVERSION[\"Offset on ellipsoidal height\",\n"
" METHOD[\"Geographic3D offsets\",\n"
" ID[\"EPSG\",9660]],\n"
" PARAMETER[\"Latitude offset\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8601]],\n"
" PARAMETER[\"Longitude offset\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8602]],\n"
" PARAMETER[\"Vertical Offset\",10,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8603]]],\n"
" CS[ellipsoidal,3],\n"
" AXIS[\"geodetic latitude (Lat)\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433,\n"
" ID[\"EPSG\",9122]]],\n"
" AXIS[\"geodetic longitude (Lon)\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433,\n"
" ID[\"EPSG\",9122]]],\n"
" AXIS[\"ellipsoidal height (h)\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1,\n"
" ID[\"EPSG\",9001]]]]";
auto dstObj = WKTParser().createFromWKT(dst_wkt);
auto dstCRS = nn_dynamic_pointer_cast<CRS>(dstObj);
ASSERT_TRUE(dstCRS != nullptr);

auto list = CoordinateOperationFactory::create()->createOperations(
authFactory->createCoordinateReferenceSystem("4979"), // WGS 84 3D
NN_NO_CHECK(dstCRS), ctxt);
ASSERT_GE(list.size(), 1U);
EXPECT_EQ(
list[0]->nameStr(),
"Inverse of CH1903+ to WGS 84 (1) + Offset on ellipsoidal height");
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m "
"+step +proj=cart +ellps=WGS84 "
"+step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346 "
"+step +inv +proj=cart +ellps=bessel "
"+step +proj=geogoffset +dlat=0 +dlon=0 +dh=10 "
"+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m "
"+step +proj=axisswap +order=2,1");
}

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

TEST(operation, vertCRS_to_geogCRS_context) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
Expand Down