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

Fix issue when transforming from/into a WKT2 Bound VerticalCRS with a 'Geographic3D to GravityRelatedHeight' method (fixes #3354) #3355

Merged
merged 1 commit into from
Sep 30, 2022
Merged

Fix issue when transforming from/into a WKT2 Bound VerticalCRS with a 'Geographic3D to GravityRelatedHeight' method (fixes #3354) #3355

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
3 changes: 3 additions & 0 deletions include/proj/coordinateoperation.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1617,6 +1617,9 @@ class PROJ_GCC_DLL Transformation : public SingleOperation {
demoteTo2D(const std::string &newName,
const io::DatabaseContextPtr &dbContext) const;

PROJ_INTERNAL static bool
isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
bool allowInverse);
//! @endcond

protected:
Expand Down
91 changes: 65 additions & 26 deletions src/iso19111/io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1472,10 +1472,6 @@ struct WKTParser::Private {

CRSPtr buildCRS(const WKTNodeNNPtr &node);

CRSPtr dealWithEPSGCodeForInterpolationCRSParameter(
std::vector<OperationParameterNNPtr> &parameters,
std::vector<ParameterValueNNPtr> &values);

TransformationNNPtr buildCoordinateOperation(const WKTNodeNNPtr &node);

ConcatenatedOperationNNPtr
Expand Down Expand Up @@ -3337,6 +3333,11 @@ void WKTParser::Private::consumeParameters(

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

static CRSPtr dealWithEPSGCodeForInterpolationCRSParameter(
DatabaseContextPtr &dbContext,
std::vector<OperationParameterNNPtr> &parameters,
std::vector<ParameterValueNNPtr> &values);

ConversionNNPtr
WKTParser::Private::buildConversion(const WKTNodeNNPtr &node,
const UnitOfMeasure &defaultLinearUnit,
Expand Down Expand Up @@ -3371,21 +3372,22 @@ WKTParser::Private::buildConversion(const WKTNodeNNPtr &node,
->inverse()));
}
auto conv = Conversion::create(convProps, methodProps, parameters, values);
auto interpolationCRS =
dealWithEPSGCodeForInterpolationCRSParameter(parameters, values);
auto interpolationCRS = dealWithEPSGCodeForInterpolationCRSParameter(
dbContext_, parameters, values);
if (interpolationCRS)
conv->setInterpolationCRS(interpolationCRS);
return conv;
}

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

CRSPtr WKTParser::Private::dealWithEPSGCodeForInterpolationCRSParameter(
static CRSPtr dealWithEPSGCodeForInterpolationCRSParameter(
DatabaseContextPtr &dbContext,
std::vector<OperationParameterNNPtr> &parameters,
std::vector<ParameterValueNNPtr> &values) {
// Transform EPSG hacky PARAMETER["EPSG code for Interpolation CRS",
// crs_epsg_code] into proper interpolation CRS
if (dbContext_ != nullptr) {
if (dbContext != nullptr) {
for (size_t i = 0; i < parameters.size(); ++i) {
const auto &l_name = parameters[i]->nameStr();
const auto epsgCode = parameters[i]->getEPSGCode();
Expand All @@ -3398,7 +3400,7 @@ CRSPtr WKTParser::Private::dealWithEPSGCodeForInterpolationCRSParameter(
static_cast<int>(values[i]->value().getSIValue());
try {
auto authFactory = AuthorityFactory::create(
NN_NO_CHECK(dbContext_), Identifier::EPSG);
NN_NO_CHECK(dbContext), Identifier::EPSG);
auto interpolationCRS =
authFactory
->createGeographicCRS(internal::toString(code))
Expand Down Expand Up @@ -3461,8 +3463,8 @@ WKTParser::Private::buildCoordinateOperation(const WKTNodeNNPtr &node) {
defaultAngularUnit);

if (interpolationCRS == nullptr)
interpolationCRS =
dealWithEPSGCodeForInterpolationCRSParameter(parameters, values);
interpolationCRS = dealWithEPSGCodeForInterpolationCRSParameter(
dbContext_, parameters, values);

std::vector<PositionalAccuracyNNPtr> accuracies;
auto &accuracyNode = nodeP->lookForChild(WKTConstants::OPERATIONACCURACY);
Expand Down Expand Up @@ -4883,6 +4885,50 @@ CRSNNPtr WKTParser::Private::buildCompoundCRS(const WKTNodeNNPtr &node) {

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

static TransformationNNPtr buildTransformationForBoundCRS(
DatabaseContextPtr &dbContext,
const util::PropertyMap &abridgedNodeProperties,
const util::PropertyMap &methodNodeProperties, const CRSNNPtr &sourceCRS,
const CRSNNPtr &targetCRS, std::vector<OperationParameterNNPtr> &parameters,
std::vector<ParameterValueNNPtr> &values) {

auto interpolationCRS = dealWithEPSGCodeForInterpolationCRSParameter(
dbContext, parameters, values);

const auto sourceTransformationCRS(
createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS));
auto transformation = Transformation::create(
abridgedNodeProperties, sourceTransformationCRS, targetCRS,
interpolationCRS, methodNodeProperties, parameters, values,
std::vector<PositionalAccuracyNNPtr>());

// If the transformation is a "Geographic3D to GravityRelatedHeight" one,
// then the sourceCRS is expected to be a GeographicCRS and the target a
// VerticalCRS. Due to how things work in a BoundCRS, we have the opposite,
// so use our "GravityRelatedHeight to Geographic3D" method instead.
if (Transformation::isGeographic3DToGravityRelatedHeight(
transformation->method(), true) &&
dynamic_cast<VerticalCRS *>(sourceTransformationCRS.get()) &&
dynamic_cast<GeographicCRS *>(targetCRS.get())) {
auto fileParameter = transformation->parameterValue(
EPSG_NAME_PARAMETER_GEOID_CORRECTION_FILENAME,
EPSG_CODE_PARAMETER_GEOID_CORRECTION_FILENAME);
if (fileParameter &&
fileParameter->type() == ParameterValue::Type::FILENAME) {
auto filename = fileParameter->valueFile();

transformation =
Transformation::createGravityRelatedHeightToGeographic3D(
abridgedNodeProperties, sourceTransformationCRS, targetCRS,
interpolationCRS, filename,
std::vector<PositionalAccuracyNNPtr>());
}
}
return transformation;
}

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

BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) {
const auto *nodeP = node->GP();
auto &abridgedNode =
Expand Down Expand Up @@ -4926,15 +4972,11 @@ BoundCRSNNPtr WKTParser::Private::buildBoundCRS(const WKTNodeNNPtr &node) {
consumeParameters(abridgedNode, true, parameters, values, defaultLinearUnit,
defaultAngularUnit);

auto interpolationCRS =
dealWithEPSGCodeForInterpolationCRSParameter(parameters, values);

const auto sourceTransformationCRS(
createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS));
auto transformation = Transformation::create(
buildProperties(abridgedNode), sourceTransformationCRS,
NN_NO_CHECK(targetCRS), interpolationCRS, buildProperties(methodNode),
parameters, values, std::vector<PositionalAccuracyNNPtr>());
const auto nnSourceCRS = NN_NO_CHECK(sourceCRS);
const auto nnTargetCRS = NN_NO_CHECK(targetCRS);
const auto transformation = buildTransformationForBoundCRS(
dbContext_, buildProperties(abridgedNode), buildProperties(methodNode),
nnSourceCRS, nnTargetCRS, parameters, values);

return BoundCRS::create(buildProperties(node, false, false),
NN_NO_CHECK(sourceCRS), NN_NO_CHECK(targetCRS),
Expand Down Expand Up @@ -6236,12 +6278,9 @@ BoundCRSNNPtr JSONParser::buildBoundCRS(const json &j) {
values.emplace_back(ParameterValue::create(getMeasure(param)));
}

const auto sourceTransformationCRS(
createBoundCRSSourceTransformationCRS(sourceCRS, targetCRS));
auto transformation = Transformation::create(
buildProperties(transformationJ), sourceTransformationCRS, targetCRS,
nullptr, buildProperties(methodJ), parameters, values,
std::vector<PositionalAccuracyNNPtr>());
const auto transformation = buildTransformationForBoundCRS(
dbContext_, buildProperties(transformationJ), buildProperties(methodJ),
sourceCRS, targetCRS, parameters, values);

return BoundCRS::create(buildProperties(j,
/* removeInverseOf= */ false,
Expand Down
9 changes: 4 additions & 5 deletions src/iso19111/operation/singleoperation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1842,9 +1842,8 @@ _getHeightToGeographic3DFilename(const SingleOperation *op, bool allowInverse) {
// ---------------------------------------------------------------------------

//! @cond Doxygen_Suppress
static bool
isGeographic3DToGravityRelatedHeight(const OperationMethodNNPtr &method,
bool allowInverse) {
bool Transformation::isGeographic3DToGravityRelatedHeight(
const OperationMethodNNPtr &method, bool allowInverse) {
const auto &methodName = method->nameStr();
static const char *const methodCodes[] = {
"1025", // Geographic3D to GravityRelatedHeight (EGM2008)
Expand Down Expand Up @@ -2155,7 +2154,7 @@ TransformationNNPtr SingleOperation::substitutePROJAlternativeGridNames(
}
}

if (isGeographic3DToGravityRelatedHeight(method(), false)) {
if (Transformation::isGeographic3DToGravityRelatedHeight(method(), false)) {
const auto &fileParameter =
parameterValue(EPSG_NAME_PARAMETER_GEOID_CORRECTION_FILENAME,
EPSG_CODE_PARAMETER_GEOID_CORRECTION_FILENAME);
Expand Down Expand Up @@ -3648,7 +3647,7 @@ bool SingleOperation::exportToPROJStringGeneric(
return true;
}

if (isGeographic3DToGravityRelatedHeight(method(), true)) {
if (Transformation::isGeographic3DToGravityRelatedHeight(method(), true)) {
auto fileParameter =
parameterValue(EPSG_NAME_PARAMETER_GEOID_CORRECTION_FILENAME,
EPSG_CODE_PARAMETER_GEOID_CORRECTION_FILENAME);
Expand Down
145 changes: 145 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 @@ -6646,6 +6646,151 @@ TEST(operation,

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

TEST(operation, compoundCRS_of_bound_horizCRS_and_bound_vertCRS_to_geogCRS) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
auto wkt =
"COMPOUNDCRS[\"CH1903 / LV03 + EGM96 height\",\n"
" BOUNDCRS[\n"
" SOURCECRS[\n"
" PROJCRS[\"CH1903 / LV03\",\n"
" BASEGEOGCRS[\"CH1903\",\n"
" DATUM[\"CH1903\",\n"
" ELLIPSOID[\"Bessel "
"1841\",6377397.155,299.1528128,\n"
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" ID[\"EPSG\",4149]],\n"
" CONVERSION[\"Swiss Oblique Mercator 1903M\",\n"
" METHOD[\"Hotine Oblique Mercator (variant B)\",\n"
" ID[\"EPSG\",9815]],\n"
" PARAMETER[\"Latitude of projection "
"centre\",46.9524055555556,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8811]],\n"
" PARAMETER[\"Longitude of projection "
"centre\",7.43958333333333,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8812]],\n"
" PARAMETER[\"Azimuth of initial line\",90,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8813]],\n"
" PARAMETER[\"Angle from Rectified to Skew "
"Grid\",90,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8814]],\n"
" PARAMETER[\"Scale factor on initial line\",1,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",8815]],\n"
" PARAMETER[\"Easting at projection "
"centre\",600000,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8816]],\n"
" PARAMETER[\"Northing at projection "
"centre\",200000,\n"
" LENGTHUNIT[\"metre\",1],\n"
" ID[\"EPSG\",8817]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"easting (Y)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" AXIS[\"northing (X)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"metre\",1]]]],\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,2],\n"
" AXIS[\"latitude\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"longitude\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" ID[\"EPSG\",4326]]],\n"
" ABRIDGEDTRANSFORMATION[\"CH1903 to WGS 84 (2)\",\n"
" VERSION[\"BfL-CH 2\"],\n"
" METHOD[\"Geocentric translations (geog2D domain)\",\n"
" ID[\"EPSG\",9603]],\n"
" PARAMETER[\"X-axis translation\",674.374,\n"
" ID[\"EPSG\",8605]],\n"
" PARAMETER[\"Y-axis translation\",15.056,\n"
" ID[\"EPSG\",8606]],\n"
" PARAMETER[\"Z-axis translation\",405.346,\n"
" ID[\"EPSG\",8607]]]],\n"
" BOUNDCRS[\n"
" SOURCECRS[\n"
" VERTCRS[\"EGM96 height\",\n"
" VDATUM[\"EGM96 geoid\"],\n"
" CS[vertical,1],\n"
" AXIS[\"gravity-related height (H)\",up,\n"
" LENGTHUNIT[\"metre\",1]],\n"
" USAGE[\n"
" SCOPE[\"Geodesy.\"],\n"
" AREA[\"World.\"],\n"
" BBOX[-90,-180,90,180]],\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[\"latitude\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"longitude\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"ellipsoidal height\",up,\n"
" ORDER[3],\n"
" LENGTHUNIT[\"metre\",1]],\n"
" ID[\"EPSG\",4979]]],\n"
" ABRIDGEDTRANSFORMATION[\"WGS 84 to EGM96 height (1)\",\n"
" METHOD[\"Geographic3D to GravityRelatedHeight (EGM)\",\n"
" ID[\"EPSG\",9661]],\n"
" PARAMETERFILE[\"Geoid (height correction) model "
"file\",\"us_nga_egm96_15.tif\"]]]]";
auto srcObj =
createFromUserInput(wkt, authFactory->databaseContext(), false);
auto src = nn_dynamic_pointer_cast<CRS>(srcObj);
ASSERT_TRUE(src != nullptr);
// WGS 84 3D
auto dst = authFactory->createCoordinateReferenceSystem("4979");

auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(src), dst, ctxt);
ASSERT_EQ(list.size(), 1U);
auto op_proj =
list[0]->exportToPROJString(PROJStringFormatter::create().get());
EXPECT_EQ(op_proj,
"+proj=pipeline "
"+step +inv +proj=somerc +lat_0=46.9524055555556 "
"+lon_0=7.43958333333333 +k_0=1 "
"+x_0=600000 +y_0=200000 +ellps=bessel "
"+step +proj=push +v_3 "
"+step +proj=cart +ellps=bessel "
"+step +proj=helmert +x=674.374 +y=15.056 +z=405.346 "
"+step +inv +proj=cart +ellps=WGS84 "
"+step +proj=pop +v_3 "
"+step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1 "
"+step +proj=unitconvert +xy_in=rad +xy_out=deg "
"+step +proj=axisswap +order=2,1");
}

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

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