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

[Backport 9.4] Improve detection of compatible/incompatible celestial bodies #4150

Merged
merged 2 commits into from
May 24, 2024
Merged

[Backport 9.4] Improve detection of compatible/incompatible celestial bodies #4150

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
35a343f
WKT/PROJJSON importer: improve heuristics to guess the celestial body…
rouault May 24, 2024
38f218d
createOperations(): improve detection of compatible/incompatible cele…
rouault May 24, 2024
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
3 changes: 2 additions & 1 deletion include/proj/datum.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -338,7 +338,8 @@ class PROJ_GCC_DLL Ellipsoid final : public common::IdentifiedObject,
//! @endcond

PROJ_INTERNAL static std::string
guessBodyName(const io::DatabaseContextPtr &dbContext, double a);
guessBodyName(const io::DatabaseContextPtr &dbContext, double a,
const std::string &ellpsName = std::string());

PROJ_INTERNAL bool lookForProjWellKnownEllps(std::string &projEllpsName,
std::string &ellpsName) const;
Expand Down
7 changes: 7 additions & 0 deletions include/proj/internal/datum_internal.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,4 +35,11 @@

#define UNKNOWN_ENGINEERING_DATUM "Unknown engineering datum"

constexpr const char *NON_EARTH_BODY = "Non-Earth body";

// Mars (2015) - Sphere uses R=3396190
// and Mars polar radius (as used by HIRISE JPEG2000) is 3376200m
// which is a 0.59% relative difference.
constexpr double REL_ERROR_FOR_SAME_CELESTIAL_BODY = 0.007;

#endif // DATUM_INTERNAL_HH_INCLUDED
27 changes: 19 additions & 8 deletions src/iso19111/datum.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1179,24 +1179,35 @@ bool Ellipsoid::_isEquivalentTo(const util::IComparable *other,
// ---------------------------------------------------------------------------

std::string Ellipsoid::guessBodyName(const io::DatabaseContextPtr &dbContext,
double a) {
// Mars (2015) - Sphere uses R=3396190
// and Mars polar radius (as used by HIRISE JPEG2000) is 3376200m
// which is a 0.59% relative difference.
constexpr double relError = 0.007;
double a, const std::string &ellpsName) {
constexpr double earthMeanRadius = 6375000.0;
if (std::fabs(a - earthMeanRadius) < relError * earthMeanRadius) {
if (std::fabs(a - earthMeanRadius) <
REL_ERROR_FOR_SAME_CELESTIAL_BODY * earthMeanRadius) {
return Ellipsoid::EARTH;
}
if (dbContext) {
try {
auto factory = io::AuthorityFactory::create(NN_NO_CHECK(dbContext),
std::string());
return factory->identifyBodyFromSemiMajorAxis(a, relError);
if (!ellpsName.empty()) {
auto matches = factory->createObjectsFromName(
ellpsName, {io::AuthorityFactory::ObjectType::ELLIPSOID},
true, 1);
if (!matches.empty()) {
auto ellps =
static_cast<const Ellipsoid *>(matches.front().get());
if (std::fabs(a - ellps->semiMajorAxis().getSIValue()) <
REL_ERROR_FOR_SAME_CELESTIAL_BODY * a) {
return ellps->celestialBody();
}
}
}
return factory->identifyBodyFromSemiMajorAxis(
a, REL_ERROR_FOR_SAME_CELESTIAL_BODY);
} catch (const std::exception &) {
}
}
return "Non-Earth body";
return NON_EARTH_BODY;
}

// ---------------------------------------------------------------------------
Expand Down
18 changes: 14 additions & 4 deletions src/iso19111/factory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4524,20 +4524,30 @@ std::string
AuthorityFactory::identifyBodyFromSemiMajorAxis(double semi_major_axis,
double tolerance) const {
auto res =
d->run("SELECT name, (ABS(semi_major_axis - ?) / semi_major_axis ) "
"AS rel_error FROM celestial_body WHERE rel_error <= ?",
d->run("SELECT DISTINCT name, "
"(ABS(semi_major_axis - ?) / semi_major_axis ) AS rel_error "
"FROM celestial_body WHERE rel_error <= ? "
"ORDER BY rel_error, name",
{semi_major_axis, tolerance});
if (res.empty()) {
throw FactoryException("no match found");
}
constexpr int IDX_NAME = 0;
if (res.size() > 1) {
constexpr int IDX_REL_ERROR = 1;
// If the first object has a relative error of 0 and the next one
// a non-zero error, then use the first one.
if (res.front()[IDX_REL_ERROR] == "0" &&
(*std::next(res.begin()))[IDX_REL_ERROR] != "0") {
return res.front()[IDX_NAME];
}
for (const auto &row : res) {
if (row[0] != res.front()[0]) {
if (row[IDX_NAME] != res.front()[IDX_NAME]) {
throw FactoryException("more than one match found");
}
}
}
return res.front()[0];
return res.front()[IDX_NAME];
}

// ---------------------------------------------------------------------------
Expand Down
25 changes: 15 additions & 10 deletions src/iso19111/io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2083,15 +2083,17 @@ EllipsoidNNPtr WKTParser::Private::buildEllipsoid(const WKTNodeNNPtr &node) {
Scale invFlattening(invFlatteningChild->GP()->value() == "\"inf\""
? 0
: asDouble(invFlatteningChild));
const auto celestialBody(
Ellipsoid::guessBodyName(dbContext_, semiMajorAxis.getSIValue()));
const auto ellpsProperties = buildProperties(node);
std::string ellpsName;
ellpsProperties.getStringValue(IdentifiedObject::NAME_KEY, ellpsName);
const auto celestialBody(Ellipsoid::guessBodyName(
dbContext_, semiMajorAxis.getSIValue(), ellpsName));
if (invFlattening.getSIValue() == 0) {
return Ellipsoid::createSphere(buildProperties(node), semiMajorAxis,
return Ellipsoid::createSphere(ellpsProperties, semiMajorAxis,
celestialBody);
} else {
return Ellipsoid::createFlattenedSphere(
buildProperties(node), semiMajorAxis, invFlattening,
celestialBody);
ellpsProperties, semiMajorAxis, invFlattening, celestialBody);
}
} catch (const std::exception &e) {
throw buildRethrow(__FUNCTION__, e);
Expand Down Expand Up @@ -7122,15 +7124,18 @@ PrimeMeridianNNPtr JSONParser::buildPrimeMeridian(const json &j) {
EllipsoidNNPtr JSONParser::buildEllipsoid(const json &j) {
if (j.contains("semi_major_axis")) {
auto semiMajorAxis = getLength(j, "semi_major_axis");
const auto celestialBody(
Ellipsoid::guessBodyName(dbContext_, semiMajorAxis.getSIValue()));
const auto ellpsProperties = buildProperties(j);
std::string ellpsName;
ellpsProperties.getStringValue(IdentifiedObject::NAME_KEY, ellpsName);
const auto celestialBody(Ellipsoid::guessBodyName(
dbContext_, semiMajorAxis.getSIValue(), ellpsName));
if (j.contains("semi_minor_axis")) {
return Ellipsoid::createTwoAxis(buildProperties(j), semiMajorAxis,
return Ellipsoid::createTwoAxis(ellpsProperties, semiMajorAxis,
getLength(j, "semi_minor_axis"),
celestialBody);
} else if (j.contains("inverse_flattening")) {
return Ellipsoid::createFlattenedSphere(
buildProperties(j), semiMajorAxis,
ellpsProperties, semiMajorAxis,
Scale(getNumber(j, "inverse_flattening")), celestialBody);
} else {
throw ParsingException(
Expand Down Expand Up @@ -10717,7 +10722,7 @@ PrimeMeridianNNPtr PROJStringParser::Private::buildPrimeMeridian(Step &step) {
std::string PROJStringParser::Private::guessBodyName(double a) {

auto ret = Ellipsoid::guessBodyName(dbContext_, a);
if (ret == "Non-Earth body" && dbContext_ == nullptr && ctx_ != nullptr) {
if (ret == NON_EARTH_BODY && dbContext_ == nullptr && ctx_ != nullptr) {
dbContext_ =
ctx_->get_cpp_context()->getDatabaseContext().as_nullable();
if (dbContext_) {
Expand Down
39 changes: 23 additions & 16 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
#include "proj/metadata.hpp"
#include "proj/util.hpp"

#include "proj/internal/datum_internal.hpp"
#include "proj/internal/internal.hpp"
#include "proj/internal/io_internal.hpp"
#include "proj/internal/tracing.hpp"
Expand Down Expand Up @@ -4442,28 +4443,34 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod(

ENTER_FUNCTION();

if (geodSrc->ellipsoid()->celestialBody() !=
geodDst->ellipsoid()->celestialBody()) {
const auto &srcEllps = geodSrc->ellipsoid();
const auto &dstEllps = geodDst->ellipsoid();
if (srcEllps->celestialBody() == dstEllps->celestialBody() &&
srcEllps->celestialBody() != NON_EARTH_BODY) {
// Same celestial body, with a known name (that is not the generic
// NON_EARTH_BODY used when we can't guess it) ==> compatible
} else if ((srcEllps->celestialBody() != dstEllps->celestialBody() &&
srcEllps->celestialBody() != NON_EARTH_BODY &&
dstEllps->celestialBody() != NON_EARTH_BODY) ||
std::fabs(srcEllps->semiMajorAxis().getSIValue() -
dstEllps->semiMajorAxis().getSIValue()) >
REL_ERROR_FOR_SAME_CELESTIAL_BODY *
dstEllps->semiMajorAxis().getSIValue()) {
const char *envVarVal = getenv("PROJ_IGNORE_CELESTIAL_BODY");
if (envVarVal == nullptr) {
if (envVarVal == nullptr || ci_equal(envVarVal, "NO") ||
ci_equal(envVarVal, "FALSE") || ci_equal(envVarVal, "OFF")) {
std::string osMsg(
"Source and target ellipsoid do not belong to the same "
"celestial body (");
osMsg += geodSrc->ellipsoid()->celestialBody();
osMsg += srcEllps->celestialBody();
osMsg += " vs ";
osMsg += geodDst->ellipsoid()->celestialBody();
osMsg += "). You may override this check by setting the "
"PROJ_IGNORE_CELESTIAL_BODY environment variable to YES.";
throw util::UnsupportedOperationException(osMsg.c_str());
} else if (ci_equal(envVarVal, "NO") || ci_equal(envVarVal, "FALSE") ||
ci_equal(envVarVal, "OFF")) {
std::string osMsg(
"Source and target ellipsoid do not belong to the same "
"celestial body (");
osMsg += geodSrc->ellipsoid()->celestialBody();
osMsg += " vs ";
osMsg += geodDst->ellipsoid()->celestialBody();
osMsg += dstEllps->celestialBody();
osMsg += ").";
if (envVarVal == nullptr) {
osMsg += " You may override this check by setting the "
"PROJ_IGNORE_CELESTIAL_BODY environment variable "
"to YES.";
}
throw util::UnsupportedOperationException(osMsg.c_str());
}
}
Expand Down
26 changes: 26 additions & 0 deletions test/unit/test_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,32 @@ TEST(wkt_parse, datum_no_pm_not_earth) {

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

TEST(wkt_parse, guess_celestial_body_from_ellipsoid_name) {
auto obj = WKTParser()
.attachDatabaseContext(DatabaseContext::create())
.createFromWKT("DATUM[\"unnamed\",\n"
" ELLIPSOID[\"Ananke\",10000,0,\n"
" LENGTHUNIT[\"metre\",1]]]");
auto datum = nn_dynamic_pointer_cast<GeodeticReferenceFrame>(obj);
ASSERT_TRUE(datum != nullptr);
EXPECT_EQ(datum->ellipsoid()->celestialBody(), "Ananke");
}

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

TEST(wkt_parse, guess_celestial_body_from_ellipsoid_name_false_positive) {
auto obj = WKTParser()
.attachDatabaseContext(DatabaseContext::create())
.createFromWKT("DATUM[\"unnamed\",\n"
" ELLIPSOID[\"Ananke\",999999,0,\n"
" LENGTHUNIT[\"metre\",1]]]");
auto datum = nn_dynamic_pointer_cast<GeodeticReferenceFrame>(obj);
ASSERT_TRUE(datum != nullptr);
EXPECT_EQ(datum->ellipsoid()->celestialBody(), "Non-Earth body");
}

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

TEST(wkt_parse, dynamic_geodetic_reference_frame) {
auto obj = WKTParser().createFromWKT(
"GEOGCRS[\"WGS 84 (G1762)\","
Expand Down
109 changes: 109 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -637,6 +637,115 @@ TEST(operation, geogCRS_to_geogCRS_context_datum_ensemble) {

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

TEST(operation, geogCRS_to_geogCRS_context_incompatible_celestial_body) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), std::string());
auto authFactoryEPSG =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto authFactoryIAU_2015 =
AuthorityFactory::create(DatabaseContext::create(), "IAU_2015");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
EXPECT_THROW(
CoordinateOperationFactory::create()->createOperations(
authFactoryEPSG->createCoordinateReferenceSystem("4326"), // WGS 84
authFactoryIAU_2015->createCoordinateReferenceSystem(
"51200"), // Ananke
ctxt),
UnsupportedOperationException);
}

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

TEST(operation,
geogCRS_to_geogCRS_context_incompatible_celestial_body_but_same_radius) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), std::string());
auto authFactoryESRI =
AuthorityFactory::create(DatabaseContext::create(), "ESRI");
auto authFactoryIAU_2015 =
AuthorityFactory::create(DatabaseContext::create(), "IAU_2015");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
EXPECT_THROW(CoordinateOperationFactory::create()->createOperations(
authFactoryESRI->createCoordinateReferenceSystem(
"104936"), // GCS_Pan_2000
authFactoryIAU_2015->createCoordinateReferenceSystem(
"51200"), // Ananke
ctxt),
UnsupportedOperationException);
}

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

TEST(operation, geogCRS_to_geogCRS_compatible_unknown_celestial_body) {
auto dbContext = DatabaseContext::create();
auto authFactory = AuthorityFactory::create(dbContext, std::string());
auto objSrc =
createFromUserInput("+proj=longlat +R=10000 +type=crs", dbContext);
auto srcCRS = nn_dynamic_pointer_cast<CRS>(objSrc);
auto objDst =
createFromUserInput("+proj=longlat +R=10001 +type=crs", dbContext);
auto dstCRS = nn_dynamic_pointer_cast<CRS>(objDst);
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
auto list = CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt);
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=noop");
}

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

TEST(operation, geogCRS_to_geogCRS_incompatible_unknown_celestial_body) {
auto dbContext = DatabaseContext::create();
auto authFactory = AuthorityFactory::create(dbContext, std::string());
auto objSrc =
createFromUserInput("+proj=longlat +R=10000 +type=crs", dbContext);
auto srcCRS = nn_dynamic_pointer_cast<CRS>(objSrc);
auto objDst =
createFromUserInput("+proj=longlat +R=99999 +type=crs", dbContext);
auto dstCRS = nn_dynamic_pointer_cast<CRS>(objDst);
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
EXPECT_THROW(CoordinateOperationFactory::create()->createOperations(
NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt),
UnsupportedOperationException);
}

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

TEST(operation,
geogCRS_to_geogCRS_compatible_celestial_body_through_semi_major_axis) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), std::string());
auto authFactoryIAU_2015 =
AuthorityFactory::create(DatabaseContext::create(), "IAU_2015");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);

auto dst_wkt = "GEOGCRS[\"unknown\",\n"
" DATUM[\"unknown\",\n"
" ELLIPSOID[\"unknown\",9999,0,\n" // Ananke is 10000
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Reference_Meridian\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" CS[ellipsoidal,2],\n"
" AXIS[\"geodetic latitude (Lat)\",north,\n"
" ORDER[1],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
" AXIS[\"geodetic longitude (Lon)\",east,\n"
" ORDER[2],\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]]]";
auto dstObj = WKTParser().createFromWKT(dst_wkt);
auto dstCRS = nn_dynamic_pointer_cast<CRS>(dstObj);

auto list = CoordinateOperationFactory::create()->createOperations(
authFactoryIAU_2015->createCoordinateReferenceSystem("51200"), // Ananke
NN_NO_CHECK(dstCRS), ctxt);
ASSERT_EQ(list.size(), 1U);
EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=noop");
}

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

TEST(operation, geogCRS_to_derived_geogCRS_3D) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
Expand Down
Loading