Skip to content

Commit

Permalink
Merge 2bebb6b into 561722c
Browse files Browse the repository at this point in the history
  • Loading branch information
@rouault
rouault committed Apr 4, 2023
2 parents 561722c + 2bebb6b commit 7cbc6a0
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 57 deletions.
2 changes: 2 additions & 0 deletions data/sql/grid_alternatives.sql
View file
Expand Up @@ -207,9 +207,11 @@ VALUES

-- nc_dittt - Gouvernement de Nouvelle Calédonie - DITTT
('Ranc08_Circe.mnt','nc_dittt_Ranc08_Circe.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/nc_dittt_Ranc08_Circe.tif',1,1,NULL),
('RANC15.tac','nc_dittt_RANC15.tif',NULL,'GTiff','geoid_like',0,NULL,'https://cdn.proj.org/nc_dittt_RANC15.tif',1,1,NULL),
('gr3dnc01b.mnt','nc_dittt_gr3dnc01b.tif',NULL,'GTiff','geocentricoffset',0,NULL,'https://cdn.proj.org/nc_dittt_gr3dnc01b.tif',1,1,NULL),
('gr3dnc02b.mnt','nc_dittt_gr3dnc02b.tif',NULL,'GTiff','geocentricoffset',0,NULL,'https://cdn.proj.org/nc_dittt_gr3dnc02b.tif',1,1,NULL),
('gr3dnc03a.mnt','nc_dittt_gr3dnc03a.tif',NULL,'GTiff','geocentricoffset',0,NULL,'https://cdn.proj.org/nc_dittt_gr3dnc03a.tif',1,1,NULL),
('gr3dncl08.tac','nc_dittt_gr3dncI08.tif',NULL,'GTiff','geocentricoffset',0,NULL,'https://cdn.proj.org/nc_dittt_gr3dncI08.tif',1,1,NULL),

-- Netherlands / RDNAP (non-free grids). See https://salsa.debian.org/debian-gis-team/proj-rdnap/raw/master/debian/copyright
-- Netherlands / RDNAP 2018
Expand Down
149 changes: 92 additions & 57 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Expand Up @@ -2978,9 +2978,14 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
const crs::CRSNNPtr &candidateDstGeod,
const CoordinateOperationNNPtr &opFirst,
bool isNullFirst) {
bool isNullFirst,
bool useOnlyDirectRegistryOp) {
bool resNonEmptyBeforeFiltering;
const auto opsSecond =
createOperations(candidateSrcGeod, candidateDstGeod, context);
useOnlyDirectRegistryOp
? findOpsInRegistryDirect(candidateSrcGeod, candidateDstGeod,
context, resNonEmptyBeforeFiltering)
: createOperations(candidateSrcGeod, candidateDstGeod, context);
const auto opsThird = createOperations(
sourceAndTargetAre3D
? candidateDstGeod->promoteTo3D(std::string(), dbContext)
Expand Down Expand Up @@ -3099,9 +3104,13 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
}
};

// The below logic is thus quite fragile, and attempts at changing it
// result in degraded results for other use cases...
//
// Start in priority with candidates that have exactly the same name as
// the sourcCRS and targetCRS. Typically for the case of init=IGNF:XXXX

// the sourceCRS and targetCRS (Typically for the case of init=IGNF:XXXX);
// and then attempt candidate geodetic CRS with different names
//
// Transformation from IGNF:NTFP to IGNF:RGF93G,
// using
// NTF geographiques Paris (gr) vers NTF GEOGRAPHIQUES GREENWICH (DMS) +
Expand All @@ -3110,72 +3119,98 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
// of IGNF:RGF93G being returned before IGNF:RGF93GEO in candidatesDstGeod.
// If RGF93GEO is returned before then we go through WGS84 and use
// instead a Helmert transformation.
// The below logic is thus quite fragile, and attempts at changing it
// result in degraded results for other use cases...

for (const auto &candidateSrcGeod : candidatesSrcGeod) {
if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
auto sourceSrcGeodModified(
sourceAndTargetAre3D
? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
: candidateSrcGeod);
for (const auto &candidateDstGeod : candidatesDstGeod) {
if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
//
// Actually, in the general case, we do the lookup in 2 passes with the 2
// above steps in each pass:
// - one first pass where we only consider direct transformations (no
// other intermediate CRS)
// - a second pass where we allow transformation through another
// intermediate CRS.
// ... but when transforming between 2 IGNF CRS, we do just one single pass
// by allowing directly all transformation. There is no strong reason for
// that particular case, except that otherwise we'd get different results
// for thest test/cli/testIGNF script when transforming a point outside
// the area of validity... Not totally sure the behaviour we try to preserve
// here with the particular case is fundamentally better than the general
// case. The general case is needed typically for the RGNC91-93 -> RGNC15
// transformation where we we need to actually use a transformation between
// RGNC91-93 (lon-lat) -> RGNC15 (lon-lat), and not chain 2 no-op
// transformations RGNC91-93 -> WGS 84 and RGNC15 -> WGS84.

const auto isIGNF = [](const crs::CRSNNPtr &crs) {
const auto &ids = crs->identifiers();
return !ids.empty() && *(ids.front()->codeSpace()) == "IGNF";
};
const int nIters = (isIGNF(sourceCRS) && isIGNF(targetCRS)) ? 1 : 2;
for (int iter = 0; iter < nIters; ++iter) {
const bool useOnlyDirectRegistryOp = (iter == 0 && nIters == 2);
for (const auto &candidateSrcGeod : candidatesSrcGeod) {
if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
auto sourceSrcGeodModified(sourceAndTargetAre3D
? candidateSrcGeod->promoteTo3D(
std::string(), dbContext)
: candidateSrcGeod);
for (const auto &candidateDstGeod : candidatesDstGeod) {
if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
#ifdef TRACE_CREATE_OPERATIONS
ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) + "->" +
objectAsStr(candidateSrcGeod.get()) + "->" +
objectAsStr(candidateDstGeod.get()) + "->" +
objectAsStr(targetCRS.get()) + ")");
ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) +
"->" + objectAsStr(candidateSrcGeod.get()) +
"->" + objectAsStr(candidateDstGeod.get()) +
"->" + objectAsStr(targetCRS.get()) + ")");
#endif
const auto opsFirst = createOperations(
sourceCRS, sourceSrcGeodModified, context);
assert(!opsFirst.empty());
const bool isNullFirst =
isNullTransformation(opsFirst[0]->nameStr());
createTransformations(candidateSrcGeod, candidateDstGeod,
opsFirst[0], isNullFirst);
if (!res.empty()) {
if (hasResultSetOnlyResultsWithPROJStep(res)) {
continue;
const auto opsFirst = createOperations(
sourceCRS, sourceSrcGeodModified, context);
assert(!opsFirst.empty());
const bool isNullFirst =
isNullTransformation(opsFirst[0]->nameStr());
createTransformations(
candidateSrcGeod, candidateDstGeod, opsFirst[0],
isNullFirst, useOnlyDirectRegistryOp);
if (!res.empty()) {
if (hasResultSetOnlyResultsWithPROJStep(res)) {
continue;
}
return;
}
return;
}
}
}
}
}

for (const auto &candidateSrcGeod : candidatesSrcGeod) {
const bool bSameSrcName =
candidateSrcGeod->nameStr() == sourceCRS->nameStr();
for (const auto &candidateSrcGeod : candidatesSrcGeod) {
const bool bSameSrcName =
candidateSrcGeod->nameStr() == sourceCRS->nameStr();
#ifdef TRACE_CREATE_OPERATIONS
ENTER_BLOCK("");
ENTER_BLOCK("");
#endif
auto sourceSrcGeodModified(
sourceAndTargetAre3D
? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
: candidateSrcGeod);
const auto opsFirst =
createOperations(sourceCRS, sourceSrcGeodModified, context);
assert(!opsFirst.empty());
const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());

for (const auto &candidateDstGeod : candidatesDstGeod) {
if (bSameSrcName &&
candidateDstGeod->nameStr() == targetCRS->nameStr()) {
continue;
}
auto sourceSrcGeodModified(
sourceAndTargetAre3D
? candidateSrcGeod->promoteTo3D(std::string(), dbContext)
: candidateSrcGeod);
const auto opsFirst =
createOperations(sourceCRS, sourceSrcGeodModified, context);
assert(!opsFirst.empty());
const bool isNullFirst =
isNullTransformation(opsFirst[0]->nameStr());

for (const auto &candidateDstGeod : candidatesDstGeod) {
if (bSameSrcName &&
candidateDstGeod->nameStr() == targetCRS->nameStr()) {
continue;
}

#ifdef TRACE_CREATE_OPERATIONS
ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) + "->" +
objectAsStr(candidateSrcGeod.get()) + "->" +
objectAsStr(candidateDstGeod.get()) + "->" +
objectAsStr(targetCRS.get()) + ")");
ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) + "->" +
objectAsStr(candidateSrcGeod.get()) + "->" +
objectAsStr(candidateDstGeod.get()) + "->" +
objectAsStr(targetCRS.get()) + ")");
#endif
createTransformations(candidateSrcGeod, candidateDstGeod,
opsFirst[0], isNullFirst);
if (!res.empty() && !hasResultSetOnlyResultsWithPROJStep(res)) {
return;
createTransformations(candidateSrcGeod, candidateDstGeod,
opsFirst[0], isNullFirst,
useOnlyDirectRegistryOp);
if (!res.empty() && !hasResultSetOnlyResultsWithPROJStep(res)) {
return;
}
}
}
}
Expand Down
25 changes: 25 additions & 0 deletions test/unit/test_operationfactory.cpp
View file
Expand Up @@ -1046,6 +1046,31 @@ TEST(operation, geogCRS_to_geogCRS_longitude_rotation_context) {

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

TEST(operation, geogCRS_to_geogCRS_context_lonlat_vs_latlon_crs) {
auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0);
ctxt->setGridAvailabilityUse(
CoordinateOperationContext::GridAvailabilityUse::
IGNORE_GRID_AVAILABILITY);
ctxt->setSpatialCriterion(
CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
auto list = CoordinateOperationFactory::create()->createOperations(
authFactory->createCoordinateReferenceSystem("4749"), // RGNC91-93
authFactory->createCoordinateReferenceSystem("10310"), // RGNC15
ctxt);
ASSERT_EQ(list.size(), 3U);
// Check that we get direct transformation, and not through WGS 84
// The difficulty here is that the transformation is registered between
// "RGNC91-93 (lon-lat)" et "RGNC15 (lon-lat)"
EXPECT_EQ(list[0]->nameStr(), "axis order change (2D) + RGNC91-93 to "
"RGNC15 (2) + axis order change (2D)");
EXPECT_EQ(list[1]->nameStr(), "axis order change (2D) + RGNC91-93 to "
"RGNC15 (1) + axis order change (2D)");
}

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

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

0 comments on commit 7cbc6a0

Please sign in to comment.