Skip to content

Commit

Permalink
Merge pull request #4149 from rouault/fix_4148
Browse files Browse the repository at this point in the history 
Improve detection of compatible/incompatible celestial bodies
  • Loading branch information
@rouault
rouault authored May 24, 2024
2 parents ea00f18 + ed11397 commit 99a51d2
Show file tree
Hide file tree
Showing 8 changed files with 215 additions and 39 deletions.
3 changes: 2 additions & 1 deletion include/proj/datum.hpp
View file
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
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
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
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
Original file line number Diff line number Diff line change
Expand Up @@ -2105,15 +2105,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 @@ -7146,15 +7148,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 @@ -10745,7 +10750,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
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
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
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

0 comments on commit 99a51d2

Please sign in to comment.