Skip to content

Commit

Permalink
createOperations(): changes so that EPSG:9635 'Geog3D to Geog2D+Gravi…
Browse files Browse the repository at this point in the history 
…tyRelatedHeight (US .gtx)' method used by Slovakian geoids correctly deal with axis order and unit conversion, to be used as 'standalone'. Also improves when using directly 'Geographic3D to GravityRelatedHeight' method
  • Loading branch information
@rouault
rouault committed Jun 6, 2020
1 parent 325d8e1 commit a96e8f5
Show file tree
Hide file tree
Showing 5 changed files with 163 additions and 21 deletions.
4 changes: 4 additions & 0 deletions include/proj/io.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -439,6 +439,10 @@ class PROJ_GCC_DLL PROJStringFormatter {
PROJ_INTERNAL void popOmitZUnitConversion();
PROJ_INTERNAL bool omitZUnitConversion() const;

PROJ_INTERNAL void pushOmitHorizontalConversionInVertTransformation();
PROJ_INTERNAL void popOmitHorizontalConversionInVertTransformation();
PROJ_INTERNAL bool omitHorizontalConversionInVertTransformation() const;

PROJ_INTERNAL void setLegacyCRSToCRSContext(bool legacyContext);
PROJ_INTERNAL bool getLegacyCRSToCRSContext() const;

Expand Down
54 changes: 54 additions & 0 deletions src/iso19111/coordinateoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8662,6 +8662,7 @@ isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
"9663", // Geographic3D to GravityRelatedHeight (OSGM-GB)
"9664", // Geographic3D to GravityRelatedHeight (IGN1997)
"9665", // Geographic3D to GravityRelatedHeight (US .gtx)
"9635", // Geog3D to Geog2D+GravityRelatedHeight (US .gtx)
};

if (ci_find(methodName, "Geographic3D to GravityRelatedHeight") == 0) {
Expand Down Expand Up @@ -9697,6 +9698,21 @@ void Transformation::_exportToPROJString(

const auto &heightFilename = _getHeightToGeographic3DFilename(this, true);
if (!heightFilename.empty()) {
auto targetCRSGeog =
extractGeographicCRSIfGeographicCRSOrEquivalent(targetCRS());
if (!targetCRSGeog) {
throw io::FormattingException(
concat("Can apply ", methodName, " only to GeographicCRS"));
}

if (!formatter->omitHorizontalConversionInVertTransformation()) {
formatter->startInversion();
formatter->pushOmitZUnitConversion();
targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
formatter->popOmitZUnitConversion();
formatter->stopInversion();
}

if (isMethodInverseOf) {
formatter->startInversion();
}
Expand All @@ -9706,6 +9722,13 @@ void Transformation::_exportToPROJString(
if (isMethodInverseOf) {
formatter->stopInversion();
}

if (!formatter->omitHorizontalConversionInVertTransformation()) {
formatter->pushOmitZUnitConversion();
targetCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
formatter->popOmitZUnitConversion();
}

return;
}

Expand All @@ -9716,6 +9739,22 @@ void Transformation::_exportToPROJString(
if (fileParameter &&
fileParameter->type() == ParameterValue::Type::FILENAME) {
auto filename = fileParameter->valueFile();

auto sourceCRSGeog =
extractGeographicCRSIfGeographicCRSOrEquivalent(sourceCRS());
if (!sourceCRSGeog) {
throw io::FormattingException(
concat("Can apply ", methodName, " only to GeographicCRS"));
}

if (!formatter->omitHorizontalConversionInVertTransformation()) {
formatter->startInversion();
formatter->pushOmitZUnitConversion();
sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
formatter->popOmitZUnitConversion();
formatter->stopInversion();
}

bool doInversion = isMethodInverseOf;
// The EPSG Geog3DToHeight is the reverse convention of PROJ !
doInversion = !doInversion;
Expand All @@ -9728,6 +9767,13 @@ void Transformation::_exportToPROJString(
if (doInversion) {
formatter->stopInversion();
}

if (!formatter->omitHorizontalConversionInVertTransformation()) {
formatter->pushOmitZUnitConversion();
sourceCRSGeog->addAngularUnitConvertAndAxisSwap(formatter);
formatter->popOmitZUnitConversion();
}

return;
}
}
Expand Down Expand Up @@ -12286,6 +12332,8 @@ createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,

//! @cond Doxygen_Suppress

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

struct MyPROJStringExportableGeodToGeod final
: public io::IPROJStringExportable {
crs::GeodeticCRSPtr geodSrc{};
Expand Down Expand Up @@ -12332,14 +12380,18 @@ struct MyPROJStringExportableHorizVertical final
_exportToPROJString(io::PROJStringFormatter *formatter) const override {

formatter->pushOmitZUnitConversion();

horizTransform->_exportToPROJString(formatter);

formatter->startInversion();
geogDst->addAngularUnitConvertAndAxisSwap(formatter);
formatter->stopInversion();

formatter->popOmitZUnitConversion();

formatter->pushOmitHorizontalConversionInVertTransformation();
verticalTransform->_exportToPROJString(formatter);
formatter->popOmitHorizontalConversionInVertTransformation();

formatter->pushOmitZUnitConversion();
geogDst->addAngularUnitConvertAndAxisSwap(formatter);
Expand Down Expand Up @@ -12385,7 +12437,9 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final

formatter->popOmitZUnitConversion();

formatter->pushOmitHorizontalConversionInVertTransformation();
verticalTransform->_exportToPROJString(formatter);
formatter->popOmitHorizontalConversionInVertTransformation();

formatter->pushOmitZUnitConversion();

Expand Down
20 changes: 20 additions & 0 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6635,6 +6635,7 @@ struct PROJStringFormatter::Private {
std::vector<InversionStackElt> inversionStack_{InversionStackElt()};
bool omitProjLongLatIfPossible_ = false;
std::vector<bool> omitZUnitConversion_{false};
std::vector<bool> omitHorizontalConversionInVertTransformation_{false};
DatabaseContextPtr dbContext_{};
bool useApproxTMerc_ = false;
bool addNoDefs_ = true;
Expand Down Expand Up @@ -7739,6 +7740,25 @@ bool PROJStringFormatter::omitZUnitConversion() const {

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

void PROJStringFormatter::pushOmitHorizontalConversionInVertTransformation() {
d->omitHorizontalConversionInVertTransformation_.push_back(true);
}

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

void PROJStringFormatter::popOmitHorizontalConversionInVertTransformation() {
assert(d->omitHorizontalConversionInVertTransformation_.size() > 1);
d->omitHorizontalConversionInVertTransformation_.pop_back();
}

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

bool PROJStringFormatter::omitHorizontalConversionInVertTransformation() const {
return d->omitHorizontalConversionInVertTransformation_.back();
}

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

void PROJStringFormatter::setLegacyCRSToCRSContext(bool legacyContext) {
d->legacyCRSToCRSContext_ = legacyContext;
}
Expand Down
10 changes: 5 additions & 5 deletions test/cli/testprojinfo_out.dist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -903,7 +903,7 @@ Operation No. 1:
INVERSE(DERIVED_FROM(PROJ)):EPSG_4977_TO_EPSG_5613, Inverse of SWEREF99 to RH2000 height, unknown accuracy, Sweden - onshore

PROJ string:
+proj=vgridshift +grids=se_lantmateriet_SWEN17_RH2000.tif +multiplier=1
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=se_lantmateriet_SWEN17_RH2000.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

WKT2:2019 string:
COORDINATEOPERATION["Inverse of SWEREF99 to RH2000 height",
Expand Down Expand Up @@ -967,15 +967,15 @@ Operation No. 1:
INVERSE(DERIVED_FROM(EPSG)):8885, Inverse of RGF93 to NGF-IGN69 height (3), 0.01 m, France - mainland onshore

PROJ string:
+proj=vgridshift +grids=fr_ign_RAF18.tif +multiplier=1
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=fr_ign_RAF18.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 2:

INVERSE(EPSG):10000, Inverse of RGF93 to NGF-IGN69 height (1), 0.5 m, France - mainland onshore

PROJ string:
+proj=vgridshift +grids=ggf97a.txt +multiplier=1
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=ggf97a.txt +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

Testing EPSG:32631 --3d
PROJ.4 string:
Expand Down Expand Up @@ -1081,7 +1081,7 @@ Operation No. 1:
DERIVED_FROM(EPSG):5656, GDA94 to AHD height (49), 0.03 m, Australia - mainland

PROJ string:
+proj=pipeline +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid09_V1.01.tif +multiplier=1
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid09_V1.01.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

Testing -s "GDA2020" -t "AHD height" --grid-check none -o PROJ --spatial-test intersects
Candidate operations found: 1
Expand All @@ -1091,7 +1091,7 @@ Operation No. 1:
DERIVED_FROM(EPSG):8451, GDA2020 to AHD height (1), 0.03 m, Australia Christmas and Cocos - onshore

PROJ string:
+proj=pipeline +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid2020_20180201.tif +multiplier=1
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=au_ga_AUSGeoid2020_20180201.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

Testing -k ellipsoid WGS84
PROJ string:
Expand Down
96 changes: 80 additions & 16 deletions test/unit/test_operation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4819,12 +4819,18 @@ TEST(operation, vertCRS_to_geogCRS_context) {
authFactory->createCoordinateReferenceSystem("4979"), // WGS 84
ctxt);
ASSERT_EQ(list.size(), 2U);
EXPECT_EQ(list[1]->exportToPROJString(
PROJStringFormatter::create(
PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1");
EXPECT_EQ(
list[1]->exportToPROJString(
PROJStringFormatter::create(
PROJStringFormatter::Convention::PROJ_5,
authFactory->databaseContext())
.get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
{
auto ctxt =
Expand All @@ -4837,7 +4843,12 @@ TEST(operation, vertCRS_to_geogCRS_context) {
ASSERT_EQ(list.size(), 2U);
EXPECT_EQ(
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1");
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
{
auto ctxt =
Expand All @@ -4850,9 +4861,13 @@ TEST(operation, vertCRS_to_geogCRS_context) {
ASSERT_EQ(list.size(), 2U);
EXPECT_EQ(
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline +step +inv +proj=vgridshift "
"+grids=us_nga_egm08_25.tif "
"+multiplier=1");
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +inv +proj=vgridshift +grids=us_nga_egm08_25.tif "
"+multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
{
auto ctxt =
Expand All @@ -4877,7 +4892,41 @@ TEST(operation, vertCRS_to_geogCRS_context) {
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=vgridshift +grids=nz_linz_nzgeoid2016.tif +multiplier=1");
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=nz_linz_nzgeoid2016.tif "
"+multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
}

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

TEST(operation, geog3DCRS_to_geog2DCRS_plus_vertCRS_context) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
{
auto ctxt =
CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
auto list = CoordinateOperationFactory::create()->createOperations(
// ETRS89 (3D)
authFactory->createCoordinateReferenceSystem("4937"),
// ETRS89 + Baltic 1957 height
authFactory->createCoordinateReferenceSystem("8360"), ctxt);
ASSERT_GE(list.size(), 1U);
EXPECT_EQ(
list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +inv +proj=vgridshift "
"+grids=sk_gku_Slovakia_ETRS89h_to_Baltic1957.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}
}

Expand Down Expand Up @@ -6774,7 +6823,13 @@ static BoundCRSNNPtr createBoundVerticalCRS() {
TEST(operation, transformation_height_to_PROJ_string) {
auto transf = createBoundVerticalCRS()->transformation();
EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1");
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +proj=vgridshift +grids=us_nga_egm08_25.tif "
"+multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");

auto grids = transf->gridsNeeded(DatabaseContext::create(), false);
ASSERT_EQ(grids.size(), 1U);
Expand Down Expand Up @@ -6843,8 +6898,13 @@ TEST(operation, transformation_Geographic3D_to_GravityRelatedHeight_gtx) {
// Check that we correctly inverse files in the case of
// "Geographic3D to GravityRelatedHeight (US .gtx)"
EXPECT_EQ(transf->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline +step +inv +proj=vgridshift "
"+grids=naptrans2008.gtx +multiplier=1");
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +inv +proj=vgridshift "
"+grids=naptrans2008.gtx +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -10710,6 +10770,10 @@ TEST(operation,
// Test that even if the .gtx file is unknown, we export in the correct
// direction
EXPECT_EQ(crs->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline +step +inv +proj=vgridshift +grids=foo.gtx "
"+multiplier=1");
"+proj=pipeline "
"+step +proj=axisswap +order=2,1 "
"+step +proj=unitconvert +xy_in=deg +xy_out=rad "
"+step +inv +proj=vgridshift +grids=foo.gtx +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

0 comments on commit a96e8f5

Please sign in to comment.