Skip to content

Commit

Permalink
Merge pull request #2386 from rouault/improve_compoundcrs_identification
Browse files Browse the repository at this point in the history 
Improve CompoundCRS identification and name morphing in VerticalCRS with ESRI WKT1
  • Loading branch information
@rouault
rouault authored Oct 21, 2020
2 parents 6b8ef54 + 31f2727 commit b61a1d0
Show file tree
Hide file tree
Showing 7 changed files with 247 additions and 32 deletions.
1 change: 1 addition & 0 deletions include/proj/crs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,6 +154,7 @@ class PROJ_GCC_DLL CRS : public common::ObjectUsage,
const cs::CoordinateSystemAxisNNPtr &verticalAxisIfNotAlreadyPresent)
const;

PROJ_INTERNAL bool hasImplicitCS() const;
//! @endcond

protected:
Expand Down
58 changes: 40 additions & 18 deletions src/apps/projinfo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1198,26 +1198,48 @@ int main(int argc, char **argv) {
std::cout << *ids[0]->codeSpace() << ":"
<< ids[0]->code() << ": "
<< pair.second << " %" << std::endl;
} else {
auto boundCRS = dynamic_cast<BoundCRS *>(
identifiedCRS.get());
if (boundCRS &&
!boundCRS->baseCRS()
->identifiers()
.empty()) {
const auto &idsBase =
boundCRS->baseCRS()->identifiers();
std::cout << "BoundCRS of "
<< *idsBase[0]->codeSpace() << ":"
<< idsBase[0]->code() << ": "
<< pair.second << " %"
<< std::endl;
} else {
std::cout
<< "un-identifier CRS: " << pair.second
<< " %" << std::endl;
continue;
}

auto boundCRS =
dynamic_cast<BoundCRS *>(identifiedCRS.get());
if (boundCRS &&
!boundCRS->baseCRS()->identifiers().empty()) {
const auto &idsBase =
boundCRS->baseCRS()->identifiers();
std::cout << "BoundCRS of "
<< *idsBase[0]->codeSpace() << ":"
<< idsBase[0]->code() << ": "
<< pair.second << " %" << std::endl;
continue;
}

auto compoundCRS = dynamic_cast<CompoundCRS *>(
identifiedCRS.get());
if (compoundCRS) {
const auto &components =
compoundCRS->componentReferenceSystems();
if (components.size() == 2 &&
!components[0]->identifiers().empty() &&
!components[1]->identifiers().empty()) {
const auto &idH =
components[0]->identifiers().front();
const auto &idV =
components[1]->identifiers().front();
if (*idH->codeSpace() ==
*idV->codeSpace()) {
std::cout << *idH->codeSpace() << ":"
<< idH->code() << '+'
<< idV->code() << ": "
<< pair.second << " %"
<< std::endl;
continue;
}
}
}

std::cout << "un-identified CRS: " << pair.second
<< " %" << std::endl;
}
} catch (const std::exception &e) {
std::cerr << "Identification failed: " << e.what()
Expand Down
82 changes: 70 additions & 12 deletions src/iso19111/crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,16 @@ CRS::~CRS() = default;

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

//! @cond Doxygen_Suppress

/** \brief Return whether the CRS has an implicit coordinate system
* (e.g from ESRI WKT) */
bool CRS::hasImplicitCS() const { return d->implicitCS_; }

//! @endcond

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

/** \brief Return the BoundCRS potentially attached to this CRS.
*
* In the case this method is called on a object returned by
Expand Down Expand Up @@ -1971,7 +1981,7 @@ GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
io::DatabaseContextPtr dbContext =
authorityFactory ? authorityFactory->databaseContext().as_nullable()
: nullptr;
const bool l_implicitCS = CRS::getPrivate()->implicitCS_;
const bool l_implicitCS = hasImplicitCS();
const auto crsCriterion =
l_implicitCS
? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS
Expand Down Expand Up @@ -2841,7 +2851,25 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
? io::WKTConstants::VERTCS
: io::WKTConstants::VERT_CS,
!identifiers().empty());
formatter->addQuotedString(nameStr());

auto l_name = nameStr();
if (formatter->useESRIDialect()) {
bool aliasFound = false;
const auto &dbContext = formatter->databaseContext();
if (dbContext) {
auto l_alias = dbContext->getAliasFromOfficialName(
l_name, "vertical_crs", "ESRI");
if (!l_alias.empty()) {
l_name = l_alias;
aliasFound = true;
}
}
if (!aliasFound) {
l_name = io::WKTFormatter::morphNameToESRI(l_name);
}
}

formatter->addQuotedString(l_name);
exportDatumOrDatumEnsembleToWkt(formatter);
const auto &cs = SingleCRS::getPrivate()->coordinateSystem;
const auto &axisList = cs->axisList();
Expand Down Expand Up @@ -4020,7 +4048,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
}
}

const bool l_implicitCS = CRS::getPrivate()->implicitCS_;
const bool l_implicitCS = hasImplicitCS();
const auto addCRS = [&](const ProjectedCRSNNPtr &crs, const bool eqName) {
const auto &l_unit = cs->axisList()[0]->unit();
if (_isEquivalentTo(crs.get(), util::IComparable::Criterion::
Expand Down Expand Up @@ -4617,6 +4645,13 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {

const auto &thisName(nameStr());

const auto &components = componentReferenceSystems();
const bool l_implicitCS = components[0]->hasImplicitCS();
const auto crsCriterion =
l_implicitCS
? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS
: util::IComparable::Criterion::EQUIVALENT;

if (authorityFactory) {
const io::DatabaseContextNNPtr &dbContext =
authorityFactory->databaseContext();
Expand All @@ -4635,9 +4670,8 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
auto crs = io::AuthorityFactory::create(
dbContext, *id->codeSpace())
->createCompoundCRS(id->code());
bool match = _isEquivalentTo(
crs.get(), util::IComparable::Criterion::EQUIVALENT,
dbContext);
bool match =
_isEquivalentTo(crs.get(), crsCriterion, dbContext);
res.emplace_back(crs, match ? 100 : 25);
return res;
} catch (const std::exception &) {
Expand All @@ -4657,16 +4691,15 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
const bool eqName = metadata::Identifier::isEquivalentName(
thisName.c_str(), crs->nameStr().c_str());
foundEquivalentName |= eqName;
if (_isEquivalentTo(
crs.get(), util::IComparable::Criterion::EQUIVALENT,
dbContext)) {
if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) {
if (crs->nameStr() == thisName) {
res.clear();
res.emplace_back(crsNN, 100);
return res;
}
res.emplace_back(crsNN, eqName ? 90 : 70);
} else {

res.emplace_back(crsNN, 25);
}
}
Expand Down Expand Up @@ -4725,9 +4758,7 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
continue;
}

if (_isEquivalentTo(crs.get(),
util::IComparable::Criterion::EQUIVALENT,
dbContext)) {
if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) {
res.emplace_back(crs, unsignificantName ? 90 : 70);
} else {
res.emplace_back(crs, 25);
Expand All @@ -4737,6 +4768,33 @@ CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
res.sort(lambdaSort);
}

// If we didn't find a match for the CompoundCRS, check if the
// horizontal and vertical parts are not themselves well known.
if (identifiers().empty() && res.empty() && components.size() == 2) {
auto candidatesHorizCRS = components[0]->identify(authorityFactory);
auto candidatesVertCRS = components[1]->identify(authorityFactory);
if (candidatesHorizCRS.size() == 1 &&
candidatesVertCRS.size() == 1 &&
candidatesHorizCRS.front().second >= 70 &&
candidatesVertCRS.front().second >= 70) {
auto newCRS = CompoundCRS::create(
util::PropertyMap().set(
common::IdentifiedObject::NAME_KEY,
candidatesHorizCRS.front().first->nameStr() + " + " +
candidatesVertCRS.front().first->nameStr()),
{candidatesHorizCRS.front().first,
candidatesVertCRS.front().first});
const bool eqName = metadata::Identifier::isEquivalentName(
thisName.c_str(), newCRS->nameStr().c_str());
res.emplace_back(
newCRS,
std::min(thisName == newCRS->nameStr() ? 100 : eqName ? 90
: 70,
std::min(candidatesHorizCRS.front().second,
candidatesVertCRS.front().second)));
}
}

// Keep only results of the highest confidence
if (res.size() >= 2) {
const auto highestConfidence = res.front().second;
Expand Down
17 changes: 16 additions & 1 deletion src/iso19111/datum.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1933,8 +1933,23 @@ void VerticalReferenceFrame::_exportToWKT(
? io::WKTConstants::VDATUM
: io::WKTConstants::VERT_DATUM,
!identifiers().empty());
const auto &l_name = nameStr();
auto l_name = nameStr();
if (!l_name.empty()) {
if (!isWKT2 && formatter->useESRIDialect()) {
bool aliasFound = false;
const auto &dbContext = formatter->databaseContext();
if (dbContext) {
auto l_alias = dbContext->getAliasFromOfficialName(
l_name, "vertical_datum", "ESRI");
if (!l_alias.empty()) {
l_name = l_alias;
aliasFound = true;
}
}
if (!aliasFound) {
l_name = io::WKTFormatter::morphNameToESRI(l_name);
}
}
formatter->addQuotedString(l_name);
} else {
formatter->addQuotedString("unnamed");
Expand Down
41 changes: 41 additions & 0 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4006,6 +4006,22 @@ VerticalReferenceFrameNNPtr WKTParser::Private::buildVerticalReferenceFrame(
const auto *nodeP = node->GP();
const std::string &name(nodeP->value());
auto &props = buildProperties(node);

if (esriStyle_ && dbContext_) {
std::string outTableName;
std::string authNameFromAlias;
std::string codeFromAlias;
auto authFactory =
AuthorityFactory::create(NN_NO_CHECK(dbContext_), std::string());
const std::string datumName = stripQuotes(nodeP->children()[0]);
auto officialName = authFactory->getOfficialNameFromAlias(
datumName, "vertical_datum", "ESRI", false, outTableName,
authNameFromAlias, codeFromAlias);
if (!officialName.empty()) {
props.set(IdentifiedObject::NAME_KEY, officialName);
}
}

if (ci_equal(name, WKTConstants::VERT_DATUM)) {
const auto &children = nodeP->children();
if (children.size() >= 2) {
Expand Down Expand Up @@ -4137,6 +4153,16 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {
throw ParsingException("Missing VDATUM or ENSEMBLE node");
}

for (const auto &childNode : nodeP->children()) {
const auto &childNodeChildren = childNode->GP()->children();
if (childNodeChildren.size() == 2 &&
ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER) &&
childNodeChildren[0]->GP()->value() == "\"Vertical_Shift\"") {
esriStyle_ = true;
break;
}
}

auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC);
auto datum =
!isNull(datumNode)
Expand All @@ -4162,6 +4188,21 @@ CRSNNPtr WKTParser::Private::buildVerticalCRS(const WKTNodeNNPtr &node) {

auto &props = buildProperties(node);

if (esriStyle_ && dbContext_) {
std::string outTableName;
std::string authNameFromAlias;
std::string codeFromAlias;
auto authFactory =
AuthorityFactory::create(NN_NO_CHECK(dbContext_), std::string());
const std::string vertCRSName = stripQuotes(nodeP->children()[0]);
auto officialName = authFactory->getOfficialNameFromAlias(
vertCRSName, "vertical_crs", "ESRI", false, outTableName,
authNameFromAlias, codeFromAlias);
if (!officialName.empty()) {
props.set(IdentifiedObject::NAME_KEY, officialName);
}
}

// Deal with Lidar WKT1 VertCRS that embeds geoid model in CRS name,
// following conventions from
// https://pubs.usgs.gov/tm/11b4/pdf/tm11-B4.pdf
Expand Down
55 changes: 54 additions & 1 deletion test/unit/test_crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3593,7 +3593,7 @@ TEST(crs, verticalCRS_as_WKT1_GDAL) {

TEST(crs, verticalCRS_as_WKT1_ESRI) {
auto crs = createVerticalCRS();
auto expected = "VERTCS[\"ODN height\",VDATUM[\"Ordnance Datum Newlyn\"],"
auto expected = "VERTCS[\"ODN_height\",VDATUM[\"Ordnance_Datum_Newlyn\"],"
"PARAMETER[\"Vertical_Shift\",0.0],"
"PARAMETER[\"Direction\",1.0],"
"UNIT[\"Meter\",1.0]]";
Expand All @@ -3603,6 +3603,23 @@ TEST(crs, verticalCRS_as_WKT1_ESRI) {
WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI).get()),
expected);
}

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

TEST(crs, verticalCRS_as_WKT1_ESRI_context) {
auto crs = createVerticalCRS();
auto expected = "VERTCS[\"Newlyn\",VDATUM[\"Ordnance_Datum_Newlyn\"],"
"PARAMETER[\"Vertical_Shift\",0.0],"
"PARAMETER[\"Direction\",1.0],"
"UNIT[\"Meter\",1.0]]";

EXPECT_EQ(crs->exportToWKT(
WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI,
DatabaseContext::create())
.get()),
expected);
}

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

TEST(crs, verticalCRS_down_as_WKT1_ESRI) {
Expand Down Expand Up @@ -4058,6 +4075,42 @@ TEST(crs, compoundCRS_identify_db) {
// Just check we don't get an exception
crs->identify(factory);
}
// Identify from ESRI WKT
{
auto sourceCRS = factory->createCompoundCRS("7405");
auto wkt = sourceCRS->exportToWKT(
WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI, dbContext)
.get());
auto obj =
WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj);
ASSERT_TRUE(crs != nullptr);
auto res = crs->identify(factory);
ASSERT_EQ(res.size(), 1U);
EXPECT_EQ(res.front().first->getEPSGCode(), 7405);
EXPECT_EQ(res.front().second, 100);
}
// Identify a CompoundCRS whose horizontal and vertical parts are known
// but not the composition.
{
auto obj = createFromUserInput("EPSG:4326+3855", dbContext);
auto sourceCRS = nn_dynamic_pointer_cast<CompoundCRS>(obj);
ASSERT_TRUE(sourceCRS != nullptr);
auto wkt = sourceCRS->exportToWKT(
WKTFormatter::create(WKTFormatter::Convention::WKT1_ESRI, dbContext)
.get());
auto obj2 =
WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt);
auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj2);
ASSERT_TRUE(crs != nullptr);
auto res = crs->identify(factory);
ASSERT_EQ(res.size(), 1U);
EXPECT_EQ(res.front().first->getEPSGCode(), 0);
EXPECT_EQ(res.front().second, 100);
const auto &components = res.front().first->componentReferenceSystems();
EXPECT_EQ(components[0]->getEPSGCode(), 4326);
EXPECT_EQ(components[1]->getEPSGCode(), 3855);
}
}

// ---------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit b61a1d0

Please sign in to comment.