Skip to content

Commit

Permalink
PRIMEM WKT handling: fixes on import for 'sexagesimal DMS' or from WK…
Browse files Browse the repository at this point in the history 
…T1:GDAL/ESRI when GEOGCS UNIT != Degree; morph to ESRI the PRIMEM name on export
  • Loading branch information
@rouault
rouault committed Nov 27, 2020
1 parent 4be79fd commit a8ad205
Show file tree
Hide file tree
Showing 4 changed files with 159 additions and 4 deletions.
17 changes: 17 additions & 0 deletions src/iso19111/datum.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -352,6 +352,23 @@ void PrimeMeridian::_exportToWKT(
if (!(isWKT2 && formatter->primeMeridianOmittedIfGreenwich() &&
l_name == "Greenwich")) {
formatter->startNode(io::WKTConstants::PRIMEM, !identifiers().empty());

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

formatter->addQuotedString(l_name);
const auto &l_long = longitude();
if (formatter->primeMeridianInDegree()) {
Expand Down
76 changes: 72 additions & 4 deletions src/iso19111/io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1979,14 +1979,80 @@ PrimeMeridianNNPtr WKTParser::Private::buildPrimeMeridian(
try {
double angleValue = asDouble(children[1]);

// Correct for GDAL WKT1 departure
// Correct for GDAL WKT1 and WKT1-ESRI departure
if (name == "Paris" && std::fabs(angleValue - 2.33722917) < 1e-8 &&
unit == UnitOfMeasure::GRAD) {
unit._isEquivalentTo(UnitOfMeasure::GRAD,
util::IComparable::Criterion::EQUIVALENT)) {
angleValue = 2.5969213;
} else {
static const struct {
const char *name;
int deg;
int min;
double sec;
} primeMeridiansDMS[] = {
{"Lisbon", -9, 7, 54.862}, {"Bogota", -74, 4, 51.3},
{"Madrid", -3, 41, 14.55}, {"Rome", 12, 27, 8.4},
{"Bern", 7, 26, 22.5}, {"Jakarta", 106, 48, 27.79},
{"Ferro", -17, 40, 0}, {"Brussels", 4, 22, 4.71},
{"Stockholm", 18, 3, 29.8}, {"Athens", 23, 42, 58.815},
{"Oslo", 10, 43, 22.5}, {"Paris RGS", 2, 20, 13.95},
{"Paris_RGS", 2, 20, 13.95}};

// Current epsg.org output may use the EPSG:9110 "sexagesimal DMS"
// unit and a DD.MMSSsss value, but this will likely be changed to
// use decimal degree.
// Or WKT1 may for example use the Paris RGS decimal degree value
// but with a GEOGCS with UNIT["Grad"]
for (const auto &pmDef : primeMeridiansDMS) {
if (name == pmDef.name) {
double dmsAsDecimalValue =
(pmDef.deg >= 0 ? 1 : -1) *
(std::abs(pmDef.deg) + pmDef.min / 100. +
pmDef.sec / 10000.);
double dmsAsDecimalDegreeValue =
(pmDef.deg >= 0 ? 1 : -1) *
(std::abs(pmDef.deg) + pmDef.min / 60. +
pmDef.sec / 3600.);
if (std::fabs(angleValue - dmsAsDecimalValue) < 1e-8 ||
std::fabs(angleValue - dmsAsDecimalDegreeValue) <
1e-8) {
angleValue = dmsAsDecimalDegreeValue;
unit = UnitOfMeasure::DEGREE;
}
break;
}
}
}

auto &properties = buildProperties(node);
if (dbContext_ && esriStyle_) {
std::string outTableName;
std::string codeFromAlias;
std::string authNameFromAlias;
auto authFactory = AuthorityFactory::create(NN_NO_CHECK(dbContext_),
std::string());
auto officialName = authFactory->getOfficialNameFromAlias(
name, "prime_meridian", "ESRI", false, outTableName,
authNameFromAlias, codeFromAlias);
if (!officialName.empty()) {
properties.set(IdentifiedObject::NAME_KEY, officialName);
if (!authNameFromAlias.empty()) {
auto identifiers = ArrayOfBaseObject::create();
identifiers->add(Identifier::create(
codeFromAlias,
PropertyMap()
.set(Identifier::CODESPACE_KEY, authNameFromAlias)
.set(Identifier::AUTHORITY_KEY,
authNameFromAlias)));
properties.set(IdentifiedObject::IDENTIFIERS_KEY,
identifiers);
}
}
}

Angle angle(angleValue, unit);
return PrimeMeridian::create(buildProperties(node), angle);
return PrimeMeridian::create(properties, angle);
} catch (const std::exception &e) {
throw buildRethrow(__FUNCTION__, e);
}
Expand Down Expand Up @@ -2734,6 +2800,9 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
throw ParsingException("Missing DATUM or ENSEMBLE node");
}

// Do that now so that esriStyle_ can be set before buildPrimeMeridian()
auto props = buildProperties(node);

auto &dynamicNode = nodeP->lookForChild(WKTConstants::DYNAMIC);

auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
Expand Down Expand Up @@ -2771,7 +2840,6 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) {
angularUnit = primeMeridian->longitude().unit();
}

auto props = buildProperties(node);
addExtensionProj4ToProp(nodeP, props);

// No explicit AXIS node ? (WKT1)
Expand Down
14 changes: 14 additions & 0 deletions test/unit/test_crs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -414,6 +414,20 @@ TEST(crs, EPSG_4326_as_WKT1_ESRI_with_database) {

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

TEST(crs, EPSG_4901_as_WKT1_ESRI_with_PRIMEM_unit_name_morphing) {
auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto crs = factory->createCoordinateReferenceSystem("4901");
WKTFormatterNNPtr f(WKTFormatter::create(
WKTFormatter::Convention::WKT1_ESRI, DatabaseContext::create()));
EXPECT_EQ(crs->exportToWKT(f.get()),
"GEOGCS[\"GCS_ATF_Paris\",DATUM[\"D_ATF\","
"SPHEROID[\"Plessis_1817\",6376523.0,308.64]],"
"PRIMEM[\"Paris_RGS\",2.33720833333333],"
"UNIT[\"Grad\",0.0157079632679489]]");
}

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

TEST(crs, EPSG_4326_as_WKT1_ESRI_without_database) {
auto crs = GeographicCRS::EPSG_4326;
WKTFormatterNNPtr f(
Expand Down
56 changes: 56 additions & 0 deletions test/unit/test_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -488,6 +488,62 @@ TEST(wkt_parse, wkt1_EPSG_4807_grad_mess) {

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

TEST(wkt_parse, wkt1_esri_EPSG_4901_grad) {
auto obj =
WKTParser()
.attachDatabaseContext(DatabaseContext::create())
.createFromWKT("GEOGCS[\"GCS_ATF_Paris\",DATUM[\"D_ATF\","
"SPHEROID[\"Plessis_1817\",6376523.0,308.64]],"
"PRIMEM[\"Paris_RGS\",2.33720833333333],"
"UNIT[\"Grad\",0.0157079632679489]]");
auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
ASSERT_TRUE(crs != nullptr);

auto datum = crs->datum();
auto primem = datum->primeMeridian();
EXPECT_EQ(primem->nameStr(), "Paris RGS");
// The PRIMEM is really in degree
EXPECT_EQ(primem->longitude().unit(), UnitOfMeasure::DEGREE);
EXPECT_NEAR(primem->longitude().value(), 2.33720833333333, 1e-14);
}

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

TEST(wkt_parse, wkt2_epsg_org_EPSG_4901_PRIMEM_weird_sexagesimal_DMS) {
// Current epsg.org output may use the EPSG:9110 "sexagesimal DMS"
// unit and a DD.MMSSsss value, but this will likely be changed to
// use decimal degree.
auto obj = WKTParser().createFromWKT(
"GEOGCRS[\"ATF (Paris)\","
" DATUM[\"Ancienne Triangulation Francaise (Paris)\","
" ELLIPSOID[\"Plessis 1817\",6376523,308.64,"
" LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]],"
" ID[\"EPSG\",7027]],"
" ID[\"EPSG\",6901]],"
" PRIMEM[\"Paris RGS\",2.201395,"
" ANGLEUNIT[\"sexagesimal DMS\",1,ID[\"EPSG\",9110]],"
" ID[\"EPSG\",8914]],"
" CS[ellipsoidal,2,"
" ID[\"EPSG\",6403]],"
" AXIS[\"Geodetic latitude (Lat)\",north,"
" ORDER[1]],"
" AXIS[\"Geodetic longitude (Lon)\",east,"
" ORDER[2]],"
" ANGLEUNIT[\"grad\",0.015707963267949,ID[\"EPSG\",9105]],"
" USAGE[SCOPE[\"Geodesy.\"],AREA[\"France - mainland onshore.\"],"
" BBOX[42.33,-4.87,51.14,8.23]],"
"ID[\"EPSG\",4901]]");
auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj);
ASSERT_TRUE(crs != nullptr);

auto datum = crs->datum();
auto primem = datum->primeMeridian();
EXPECT_EQ(primem->longitude().unit(), UnitOfMeasure::DEGREE);
EXPECT_NEAR(primem->longitude().value(), 2.33720833333333, 1e-14);
}

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

TEST(wkt_parse, wkt1_geographic_old_datum_name_from_EPSG_code) {
auto wkt =
"GEOGCS[\"S-JTSK (Ferro)\",\n"
Expand Down

0 comments on commit a8ad205

Please sign in to comment.