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

WKT1_ESRI: fix import and export of CompoundCRS #2389

Merged
merged 6 commits into from
Oct 23, 2020
Merged

WKT1_ESRI: fix import and export of CompoundCRS #2389

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
698e51f
WKT parser: accept implicit compoundCRS from ESRI WKT, like "PROJCS[.…
rouault Oct 22, 2020
f2cd2b3
WKT parser: accept ESRI VERTCS[...,DATUM[...,SPHEROID[]] syntax to ex…
rouault Oct 22, 2020
82cd218
WKT1_ESRI export: export CompoundCRS as PROJCS[...],VERTCS[...] or GE…
rouault Oct 22, 2020
9a11293
WKT1_ESRI export: generate VERTCS[...,DATUM[...,SPHEROID[]] syntax wh…
rouault Oct 22, 2020
4969076
Database: import ESRI VERTCS that uses a (geodetic) datum to express …
rouault Oct 22, 2020
6351422
WKT1_ESRI export: try to export Geographic3D and Projected3D CRS when…
rouault Oct 22, 2020
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
857 changes: 685 additions & 172 deletions data/sql/esri.sql
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion data/sql/proj_db_table_defs.sql
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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