Skip to content

Commit

Permalink
DerivedProjectedCRS WKT import: workaround lack of explicit CS in BAS…
Browse files Browse the repository at this point in the history 
…EPROJCRS

Current WKT2 grammar doesn't allow an explicit CS (Coordinate System) in the BASEPROJCRS of
a DERIVEDPROJCRS. This leads to ambiguities when instanciating the
BASEPROJCRS. Cf OSGeo/gdal#9732 (comment)
for a longer analysis

When importing such a DerivedProjectedCRS, we apply the following
strategy to try to reconstrut the CS of the baseProjCRS:
- if the baseProjCRS has an ID[], do a datebase lookup for this object,
  and if found, use its CS
- otherwise, do a database lookup from the name of the baseProjCRS, and
  if found, use its CS
- otherwise, browse through the PARAMETERs of the CONVERSION of the
  BASEPROJCRS, and if there is at least one such PARAMETER that has a
  linear unit, and if all PARAMETERs with a linear unit have the same
  one, then use it to reconstruct a Easting/Northing CartesianCS.
- otherwise, fallback to current behaviour, that is return a
  Easting/Northing CartesianCS with metre unit.

Fixes OSGeo/gdal#9732
  • Loading branch information
@rouault
rouault committed Apr 23, 2024
1 parent 08ee492 commit 1d6203e
Show file tree
Hide file tree
Showing 2 changed files with 248 additions and 2 deletions.
74 changes: 72 additions & 2 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4487,10 +4487,80 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
!ci_equal(nodeValue, WKTConstants::BASEPROJCRS)) {
ThrowMissing(WKTConstants::CS_);
}
auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);

const std::string projCRSName = stripQuotes(nodeP->children()[0]);

auto cs = [this, &projCRSName, &nodeP, &props, &csNode, &node, &nodeValue,
&conversionNode]() -> CoordinateSystemNNPtr {
if (isNull(csNode) && ci_equal(nodeValue, WKTConstants::BASEPROJCRS) &&
!isNull(conversionNode)) {
// A BASEPROJCRS (as of WKT2 18-010r11) normally lacks an explicit
// CS[] which cause issues to properly instanciate it. So we first
// start by trying to identify the BASEPROJCRS by its id or name.
// And fallback to exploring the conversion parameters to infer the
// CS AXIS unit from the linear parameter unit... Not fully bullet
// proof.
if (dbContext_) {
// Get official name from database if ID is present
auto &idNode = nodeP->lookForChild(WKTConstants::ID);
if (!isNull(idNode)) {
try {
auto id = buildId(idNode, false, false);
auto authFactory = AuthorityFactory::create(
NN_NO_CHECK(dbContext_), *id->codeSpace());
auto projCRS =
authFactory->createProjectedCRS(id->code());
return projCRS->coordinateSystem();
} catch (const std::exception &) {
}
}

auto authFactory = AuthorityFactory::create(
NN_NO_CHECK(dbContext_), std::string());
auto res = authFactory->createObjectsFromName(
projCRSName, {AuthorityFactory::ObjectType::PROJECTED_CRS},
false, 2);
if (res.size() == 1) {
auto projCRS =
dynamic_cast<const ProjectedCRS *>(res.front().get());
if (projCRS) {
return projCRS->coordinateSystem();
}
}
}

auto conv = buildConversion(conversionNode, UnitOfMeasure::METRE,
UnitOfMeasure::DEGREE);
UnitOfMeasure linearUOM = UnitOfMeasure::NONE;
for (const auto &genOpParamvalue : conv->parameterValues()) {
auto opParamvalue =
dynamic_cast<const operation::OperationParameterValue *>(
genOpParamvalue.get());
if (opParamvalue) {
const auto &parameterValue = opParamvalue->parameterValue();
if (parameterValue->type() ==
operation::ParameterValue::Type::MEASURE) {
const auto &measure = parameterValue->value();
const auto &unit = measure.unit();
if (unit.type() == UnitOfMeasure::Type::LINEAR) {
if (linearUOM == UnitOfMeasure::NONE) {
linearUOM = unit;
} else if (linearUOM != unit) {
linearUOM = UnitOfMeasure::NONE;
break;
}
}
}
}
}
if (linearUOM != UnitOfMeasure::NONE) {
return CartesianCS::createEastingNorthing(linearUOM);
}
}
return buildCS(csNode, node, UnitOfMeasure::NONE);
}();
auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);

if (esriStyle_ && dbContext_) {
if (cartesianCS) {
std::string outTableName;
Expand Down
176 changes: 176 additions & 0 deletions test/unit/test_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5643,6 +5643,182 @@ TEST(wkt_parse, DerivedProjectedCRS) {

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

TEST(wkt_parse, DerivedProjectedCRS_base_crs_cs_non_metre_from_conversion) {
auto wkt =
"DERIVEDPROJCRS[\"Ground for NAD83(2011) / Idaho West (ftUS)\",\n"
" BASEPROJCRS[\"foo\",\n"
" BASEGEOGCRS[\"NAD83(2011)\",\n"
" DATUM[\"NAD83 (National Spatial Reference System "
"2011)\",\n"
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
" CONVERSION[\"SPCS83 Idaho West zone (US Survey feet)\",\n"
" METHOD[\"Transverse Mercator\",\n"
" ID[\"EPSG\",9807]],\n"
" PARAMETER[\"Latitude of natural "
"origin\",41.6666666666667,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural origin\",-115.75,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"Scale factor at natural "
"origin\",0.999933333,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",8805]],\n"
" PARAMETER[\"False easting\",2624666.667,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",0,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8807]]]],\n"
" DERIVINGCONVERSION[\"Grid to ground\",\n"
" METHOD[\"Similarity transformation\",\n"
" ID[\"EPSG\",9621]],\n"
" PARAMETER[\"Ordinate 1 of evaluation point in target "
"CRS\",1000,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8621]],\n"
" PARAMETER[\"Ordinate 2 of evaluation point in target "
"CRS\",0,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8622]],\n"
" PARAMETER[\"Scale factor for source CRS axes\",1,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",1061]],\n"
" PARAMETER[\"Rotation angle of source CRS axes\",0,\n"
" ANGLEUNIT[\"degree\",0],\n"
" ID[\"EPSG\",8614]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"easting (X)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n"
" AXIS[\"northing (Y)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]]]";

auto obj = WKTParser().createFromWKT(wkt);
auto crs = nn_dynamic_pointer_cast<DerivedProjectedCRS>(obj);
ASSERT_TRUE(crs != nullptr);

const auto &axisList = crs->baseCRS()->coordinateSystem()->axisList();
ASSERT_EQ(axisList.size(), 2U);
EXPECT_EQ(axisList[0]->unit(), UnitOfMeasure::US_FOOT);
}

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

TEST(
wkt_parse,
DerivedProjectedCRS_base_crs_cs_non_metre_from_conversion_context_from_baseprojcrs_name) {
auto wkt =
"DERIVEDPROJCRS[\"Ground for NAD83(2011) / Idaho West (ftUS)\",\n"
" BASEPROJCRS[\"NAD83(2011) / Idaho West (ftUS)\",\n"
" BASEGEOGCRS[\"NAD83(2011)\",\n"
" DATUM[\"NAD83 (National Spatial Reference System "
"2011)\",\n"
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
" CONVERSION[\"incomplete conversion\",\n"
" METHOD[\"Transverse Mercator\",\n"
" ID[\"EPSG\",9807]]]],\n"
" DERIVINGCONVERSION[\"Grid to ground\",\n"
" METHOD[\"Similarity transformation\",\n"
" ID[\"EPSG\",9621]],\n"
" PARAMETER[\"Ordinate 1 of evaluation point in target "
"CRS\",1000,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8621]],\n"
" PARAMETER[\"Ordinate 2 of evaluation point in target "
"CRS\",0,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8622]],\n"
" PARAMETER[\"Scale factor for source CRS axes\",1,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",1061]],\n"
" PARAMETER[\"Rotation angle of source CRS axes\",0,\n"
" ANGLEUNIT[\"degree\",0],\n"
" ID[\"EPSG\",8614]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"easting (X)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n"
" AXIS[\"northing (Y)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]]]";

auto obj = WKTParser()
.attachDatabaseContext(DatabaseContext::create())
.createFromWKT(wkt);
auto crs = nn_dynamic_pointer_cast<DerivedProjectedCRS>(obj);
ASSERT_TRUE(crs != nullptr);

const auto &axisList = crs->baseCRS()->coordinateSystem()->axisList();
ASSERT_EQ(axisList.size(), 2U);
EXPECT_EQ(axisList[0]->unit(), UnitOfMeasure::US_FOOT);
}

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

TEST(
wkt_parse,
DerivedProjectedCRS_base_crs_cs_non_metre_from_conversion_context_from_baseprojcrs_id) {
auto wkt =
"DERIVEDPROJCRS[\"Ground for NAD83(2011) / Idaho West (ftUS)\",\n"
" BASEPROJCRS[\"foo\",\n"
" BASEGEOGCRS[\"NAD83(2011)\",\n"
" DATUM[\"NAD83 (National Spatial Reference System "
"2011)\",\n"
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
" CONVERSION[\"incomplete conversion\",\n"
" METHOD[\"Transverse Mercator\",\n"
" ID[\"EPSG\",9807]]],\n"
" ID[\"EPSG\",6453]],\n"
" DERIVINGCONVERSION[\"Grid to ground\",\n"
" METHOD[\"Similarity transformation\",\n"
" ID[\"EPSG\",9621]],\n"
" PARAMETER[\"Ordinate 1 of evaluation point in target "
"CRS\",1000,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8621]],\n"
" PARAMETER[\"Ordinate 2 of evaluation point in target "
"CRS\",0,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8622]],\n"
" PARAMETER[\"Scale factor for source CRS axes\",1,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",1061]],\n"
" PARAMETER[\"Rotation angle of source CRS axes\",0,\n"
" ANGLEUNIT[\"degree\",0],\n"
" ID[\"EPSG\",8614]]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"easting (X)\",east,\n"
" ORDER[1],\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n"
" AXIS[\"northing (Y)\",north,\n"
" ORDER[2],\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]]]";

auto obj = WKTParser()
.attachDatabaseContext(DatabaseContext::create())
.createFromWKT(wkt);
auto crs = nn_dynamic_pointer_cast<DerivedProjectedCRS>(obj);
ASSERT_TRUE(crs != nullptr);

const auto &axisList = crs->baseCRS()->coordinateSystem()->axisList();
ASSERT_EQ(axisList.size(), 2U);
EXPECT_EQ(axisList[0]->unit(), UnitOfMeasure::US_FOOT);
}

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

TEST(wkt_parse, DerivedProjectedCRS_ordinal) {
auto wkt = "DERIVEDPROJCRS[\"derived projectedCRS\",\n"
" BASEPROJCRS[\"BASEPROJCRS\",\n"
Expand Down

0 comments on commit 1d6203e

Please sign in to comment.