Skip to content

Commit

Permalink
Merge pull request #2389 from rouault/wkt_parse_implict_compound_crs_…
Browse files Browse the repository at this point in the history 
…esri

WKT1_ESRI: fix import and export of CompoundCRS
  • Loading branch information
@rouault
rouault authored Oct 23, 2020
2 parents 4fa0c25 + 6351422 commit aaa1592
Show file tree
Hide file tree
Showing 11 changed files with 1,328 additions and 432 deletions.
857 changes: 685 additions & 172 deletions data/sql/esri.sql
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion data/sql/proj_db_table_defs.sql
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ FOR EACH ROW BEGIN
SELECT RAISE(ABORT, 'insert on usage violates constraint: extent must not be deprecated when object is not deprecated')
WHERE EXISTS (
SELECT 1 FROM extent JOIN object_view o WHERE
NOT (o.table_name IN ('projected_crs', 'conversion') AND o.auth_name = 'ESRI') AND
NOT (o.table_name IN ('projected_crs', 'vertical_crs', 'vertical_datum', 'conversion') AND o.auth_name = 'ESRI') AND
o.table_name = NEW.object_table_name AND
o.auth_name = NEW.object_auth_name AND
o.code = NEW.object_code AND
Expand Down
113 changes: 82 additions & 31 deletions scripts/build_db_from_esri.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -436,6 +436,16 @@ def import_datum():
sql = """INSERT INTO alias_name VALUES('geodetic_datum','EPSG','%s','%s','ESRI');""" % (
code, escape_literal(esri_name))
all_sql.append(sql)

description = row[idx_description]
deprecated = 1 if row[idx_deprecated] == 'yes' else 0

map_datum_esri_to_parameters[code] = {
'esri_name': esri_name,
'description': description,
'deprecated': deprecated
}

else:
assert authority.upper() == 'ESRI', row

Expand Down Expand Up @@ -1187,7 +1197,14 @@ def import_vertcs():
idx_rlon = header.index('rlon')
assert idx_rlon >= 0

datum_written = set()
vdatum_written = set()

sql = """-- vertical coordinate system for ellipsoidal height. Not really ISO 19111 valid..."""
all_sql.append(sql)
sql = """INSERT INTO "coordinate_system" VALUES('ESRI','ELLPS_HEIGHT_METRE','vertical',1);"""
all_sql.append(sql)
sql = """INSERT INTO "axis" VALUES('ESRI','ELLPS_HEIGHT_METRE','Ellipsoidal height','h','up','ESRI','ELLPS_HEIGHT_METRE',1,'EPSG','9001');"""
all_sql.append(sql)

while True:
try:
Expand Down Expand Up @@ -1221,51 +1238,81 @@ def import_vertcs():

wkt = row[idx_wkt]

if ',DATUM[' in wkt:
# FIXME ??
print('Skipping %s. Should be a CompoundCRS' % (esri_name))
sql = """-- Skipping %s. Should be a CompoundCRS""" % (
esri_name)
all_sql.append(sql)
continue

pos = wkt.find('VDATUM["')
assert pos >= 0
pos += len('VDATUM["')
end_pos = wkt[pos:].find('"')
assert end_pos >= 0
end_pos += pos
datum_name = wkt[pos:end_pos]
is_vdatum = True
if pos > 0:
pos += len('VDATUM["')
end_pos = wkt[pos:].find('"')
assert end_pos >= 0
end_pos += pos
datum_name = wkt[pos:end_pos]
if datum_name not in map_vdatum_esri_name_to_auth_code:
print('Skipping vertcs %s. Cannot find vertical datum %s' % (
str(row), datum_name))
sql = """-- Skipping vertcs %s. Cannot find vertical datum %s""" % (
esri_name, datum_name)
all_sql.append(sql)
continue
datum_auth_name, datum_code = map_vdatum_esri_name_to_auth_code[datum_name]
else:
pos = wkt.find(',DATUM[')
assert pos > 0
pos += len(',DATUM["')
is_vdatum = False
end_pos = wkt[pos:].find('"')
assert end_pos >= 0
end_pos += pos
datum_name = wkt[pos:end_pos]
if datum_name not in map_datum_esri_name_to_auth_code:
print('Skipping vertcs %s. Cannot find geodetic datum %s' % (
str(row), datum_name))
sql = """-- Skipping vertcs %s. Cannot find geodetic datum %s""" % (
esri_name, datum_name)
all_sql.append(sql)
continue
datum_auth_name, datum_code = map_datum_esri_name_to_auth_code[datum_name]

assert 'PARAMETER["Vertical_Shift",0.0]' in wkt, row

if datum_name not in map_vdatum_esri_name_to_auth_code:
print('Skipping vertcs %s. it expresses an ellipsoidal height regarding a horizontal datum' % (
str(row)))
sql = """-- Skipping vertcs %s. it expresses an ellipsoidal height regarding a horizontal datum""" % (
esri_name, datum_name)
all_sql.append(sql)
continue

datum_auth_name, datum_code = map_vdatum_esri_name_to_auth_code[datum_name]

deprecated = 1 if row[idx_deprecated] == 'yes' else 0

extent_auth_name, extent_code = find_extent(
row[idx_areaname], row[idx_slat], row[idx_nlat], row[idx_llon], row[idx_rlon])

if datum_auth_name == 'ESRI':
if datum_code not in datum_written:
datum_written.add(datum_code)
if not is_vdatum:
new_datum_code = 'from_geogdatum_' + datum_auth_name + '_' + datum_code
if new_datum_code not in vdatum_written:
vdatum_written.add(new_datum_code)

p = map_datum_esri_to_parameters[datum_code]

datum_code = new_datum_code

p = map_vdatum_esri_to_parameters[datum_code]
sql = """INSERT INTO "vertical_datum" VALUES('ESRI','%s','%s',NULL,NULL,NULL,NULL,%d);""" % (
datum_code, p['esri_name'], p['deprecated'])
all_sql.append(sql)
sql = """INSERT INTO "usage" VALUES('ESRI', '%s_USAGE','vertical_datum','ESRI','%s','%s','%s','%s','%s');""" % (datum_code, datum_code, extent_auth_name, extent_code, 'EPSG', '1024')
all_sql.append(sql)
else:
datum_code = new_datum_code

datum_auth_name = 'ESRI'

elif datum_auth_name == 'ESRI':
assert datum_code not in vdatum_written

vdatum_written.add(datum_code)

p = map_vdatum_esri_to_parameters[datum_code]
sql = """INSERT INTO "vertical_datum" VALUES('ESRI','%s','%s',NULL,NULL,NULL,NULL,%d);""" % (
datum_code, p['esri_name'], p['deprecated'])
all_sql.append(sql)
sql = """INSERT INTO "usage" VALUES('ESRI', '%s_USAGE','vertical_datum','ESRI','%s','%s','%s','%s','%s');""" % (datum_code, datum_code, extent_auth_name, extent_code, 'EPSG', '1024')
all_sql.append(sql)

map_vertcs_esri_name_to_auth_code[esri_name] = ['ESRI', code]

cs_auth = 'EPSG'
if 'PARAMETER["Direction",1.0]' in wkt and 'UNIT["Meter"' in wkt:
cs_code = 6499
elif 'PARAMETER["Direction",1.0]' in wkt and 'UNIT["Foot"' in wkt:
Expand All @@ -1275,9 +1322,13 @@ def import_vertcs():
else:
assert False, ('unknown coordinate system for %s' %
str(row))
if not is_vdatum:
assert cs_code == 6499
cs_auth = 'ESRI'
cs_code = 'ELLPS_HEIGHT_METRE'

sql = """INSERT INTO "vertical_crs" VALUES('ESRI','%s','%s',NULL,'EPSG','%s','%s','%s',%d);""" % (
code, esri_name, cs_code, datum_auth_name, datum_code, deprecated)
sql = """INSERT INTO "vertical_crs" VALUES('ESRI','%s','%s',NULL,'%s','%s','%s','%s',%d);""" % (
code, esri_name, cs_auth, cs_code, datum_auth_name, datum_code, deprecated)
all_sql.append(sql)
sql = """INSERT INTO "usage" VALUES('ESRI', '%s_USAGE','vertical_crs','ESRI','%s','%s','%s','%s','%s');""" % (code, code, extent_auth_name, extent_code, 'EPSG', '1024')
all_sql.append(sql)
Expand Down
122 changes: 109 additions & 13 deletions src/iso19111/crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1577,6 +1577,49 @@ GeodeticCRS::create(const util::PropertyMap &properties,

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

//! @cond Doxygen_Suppress

// Try to format a Geographic/ProjectedCRS 3D CRS as a
// GEOGCS[]/PROJCS[],VERTCS[...,DATUM[],...] if we find corresponding objects
static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
const CRS *self, const GeodeticCRS *geodCRS, io::WKTFormatter *formatter) {
const auto &dbContext = formatter->databaseContext();
if (!dbContext) {
return false;
}
const auto l_datum = geodCRS->datumNonNull(formatter->databaseContext());
auto l_alias = dbContext->getAliasFromOfficialName(
l_datum->nameStr(), "geodetic_datum", "ESRI");
if (l_alias.empty()) {
return false;
}
auto authFactory =
io::AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string());
auto list = authFactory->createObjectsFromName(
l_alias, {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
false /* approximate=false*/);
if (list.empty()) {
return false;
}
auto gdatum = util::nn_dynamic_pointer_cast<datum::Datum>(list.front());
if (gdatum == nullptr || gdatum->identifiers().empty()) {
return false;
}
const auto &gdatum_ids = gdatum->identifiers();
auto vertCRSList = authFactory->createVerticalCRSFromDatum(
"ESRI", "from_geogdatum_" + *gdatum_ids[0]->codeSpace() + '_' +
gdatum_ids[0]->code());
if (vertCRSList.size() != 1) {
return false;
}
self->demoteTo2D(std::string(), dbContext)->_exportToWKT(formatter);
vertCRSList.front()->_exportToWKT(formatter);
return true;
}
//! @endcond

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

//! @cond Doxygen_Suppress
void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
Expand All @@ -1589,11 +1632,21 @@ void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
auto l_name = nameStr();
const auto &dbContext = formatter->databaseContext();

if (formatter->useESRIDialect()) {
if (axisList.size() != 2) {
if (!isWKT2 && formatter->useESRIDialect() && axisList.size() == 3) {
if (!isGeographic) {
io::FormattingException::Throw(
"Only export of Geographic 2D CRS is supported in WKT1_ESRI");
"Geocentric CRS not supported in WKT1_ESRI");
}
// Try to format the Geographic 3D CRS as a GEOGCS[],VERTCS[...,DATUM[]]
// if we find corresponding objects
if (dbContext) {
if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(this, this,
formatter)) {
return;
}
}
io::FormattingException::Throw(
"Cannot export this Geographic 3D CRS in WKT1_ESRI");
}

if (!isWKT2 && formatter->isStrict() && isGeographic &&
Expand Down Expand Up @@ -2853,9 +2906,9 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
!identifiers().empty());

auto l_name = nameStr();
const auto &dbContext = formatter->databaseContext();
if (formatter->useESRIDialect()) {
bool aliasFound = false;
const auto &dbContext = formatter->databaseContext();
if (dbContext) {
auto l_alias = dbContext->getAliasFromOfficialName(
l_name, "vertical_crs", "ESRI");
Expand All @@ -2870,7 +2923,34 @@ void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
}

formatter->addQuotedString(l_name);
exportDatumOrDatumEnsembleToWkt(formatter);

const auto l_datum = datum();
if (formatter->useESRIDialect() && l_datum &&
l_datum->getWKT1DatumType() == "2002") {
bool foundMatch = false;
if (dbContext) {
auto authFactory = io::AuthorityFactory::create(
NN_NO_CHECK(dbContext), std::string());
auto list = authFactory->createObjectsFromName(
l_datum->nameStr(),
{io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
false /* approximate=false*/);
if (!list.empty()) {
auto gdatum =
util::nn_dynamic_pointer_cast<datum::Datum>(list.front());
if (gdatum) {
gdatum->_exportToWKT(formatter);
foundMatch = true;
}
}
}
if (!foundMatch) {
// We should export a geodetic datum, but we cannot really do better
l_datum->_exportToWKT(formatter);
}
} else {
exportDatumOrDatumEnsembleToWkt(formatter);
}
const auto &cs = SingleCRS::getPrivate()->coordinateSystem;
const auto &axisList = cs->axisList();

Expand Down Expand Up @@ -3497,6 +3577,16 @@ void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
}
}

if (formatter->useESRIDialect() && dbContext) {
// Try to format the ProjecteD 3D CRS as a
// PROJCS[],VERTCS[...,DATUM[]]
// if we find corresponding objects
if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
this, baseCRS().as_nullable().get(), formatter)) {
return;
}
}

if (!formatter->useESRIDialect() &&
CRS::getPrivate()->allowNonConformantWKT1Export_) {
formatter->startNode(io::WKTConstants::COMPD_CS, false);
Expand Down Expand Up @@ -4528,15 +4618,21 @@ CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties,
//! @cond Doxygen_Suppress
void CompoundCRS::_exportToWKT(io::WKTFormatter *formatter) const {
const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS
: io::WKTConstants::COMPD_CS,
!identifiers().empty());
formatter->addQuotedString(nameStr());
for (const auto &crs : componentReferenceSystems()) {
crs->_exportToWKT(formatter);
const auto &l_components = componentReferenceSystems();
if (!isWKT2 && formatter->useESRIDialect() && l_components.size() == 2) {
l_components[0]->_exportToWKT(formatter);
l_components[1]->_exportToWKT(formatter);
} else {
formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS
: io::WKTConstants::COMPD_CS,
!identifiers().empty());
formatter->addQuotedString(nameStr());
for (const auto &crs : l_components) {
crs->_exportToWKT(formatter);
}
ObjectUsage::baseExportToWKT(formatter);
formatter->endNode();
}
ObjectUsage::baseExportToWKT(formatter);
formatter->endNode();
}
//! @endcond

Expand Down
19 changes: 15 additions & 4 deletions src/iso19111/factory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2100,6 +2100,9 @@ AuthorityFactory::createVerticalDatum(const std::string &code) const {
if (!publication_date.empty()) {
props.set("PUBLICATION_DATE", publication_date);
}
if (d->authority() == "ESRI" && starts_with(code, "from_geogdatum_")) {
props.set("VERT_DATUM_TYPE", "2002");
}
auto anchor = util::optional<std::string>();
if (frame_reference_epoch.empty()) {
return datum::VerticalReferenceFrame::create(props, anchor);
Expand Down Expand Up @@ -5719,17 +5722,23 @@ AuthorityFactory::createObjectsFromNameEx(
"SELECT table_name, auth_name, code, name, deprecated, is_alias "
"FROM (");

const auto getTableAndTypeConstraints = [&allowedObjectTypes]() {
const auto getTableAndTypeConstraints = [&allowedObjectTypes,
&searchedName]() {
typedef std::pair<std::string, std::string> TableType;
std::list<TableType> res;
// Hide ESRI D_ vertical datums
const bool startsWithDUnderscore = starts_with(searchedName, "D_");
if (allowedObjectTypes.empty()) {
for (const auto &tableName :
{"prime_meridian", "ellipsoid", "geodetic_datum",
"vertical_datum", "geodetic_crs", "projected_crs",
"vertical_crs", "compound_crs", "conversion",
"helmert_transformation", "grid_transformation",
"other_transformation", "concatenated_operation"}) {
res.emplace_back(TableType(tableName, std::string()));
if (!(startsWithDUnderscore &&
strcmp(tableName, "vertical_datum") == 0)) {
res.emplace_back(TableType(tableName, std::string()));
}
}
} else {
for (const auto type : allowedObjectTypes) {
Expand All @@ -5744,8 +5753,10 @@ AuthorityFactory::createObjectsFromNameEx(
case ObjectType::DATUM:
res.emplace_back(
TableType("geodetic_datum", std::string()));
res.emplace_back(
TableType("vertical_datum", std::string()));
if (!startsWithDUnderscore) {
res.emplace_back(
TableType("vertical_datum", std::string()));
}
break;
case ObjectType::GEODETIC_REFERENCE_FRAME:
res.emplace_back(
Expand Down
Loading

0 comments on commit aaa1592

Please sign in to comment.