Skip to content

Commit

Permalink
Merge pull request #3876 from rouault/fix_3854
Browse files Browse the repository at this point in the history 
CoordinateOperationFactory: deal with CompoundToCompound with a horizontal similarity transformation and a ballpark vertical (fixes #3854)
  • Loading branch information
@rouault
rouault authored Sep 5, 2023
2 parents 8b4fa7d + 89943df commit d750d68
Show file tree
Hide file tree
Showing 6 changed files with 168 additions and 40 deletions.
44 changes: 23 additions & 21 deletions src/4D_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,22 +322,14 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
// onshore and offshore Tunisia area of uses, but is slightly
// onshore. So in a general way, prefer a onshore area to a
// offshore one.
if (iBest < 0 || (alt.accuracy >= 0 &&
(alt.accuracy < bestAccuracy ||
// If two operations have the same accuracy, use
// the one that is contained within a larger one
(alt.accuracy == bestAccuracy &&
alt.minxSrc >= opList[iBest].minxSrc &&
alt.minySrc >= opList[iBest].minySrc &&
alt.maxxSrc <= opList[iBest].maxxSrc &&
alt.maxySrc <= opList[iBest].maxySrc &&
// check that this is not equality
!(alt.minxSrc == opList[iBest].minxSrc &&
alt.minySrc == opList[iBest].minySrc &&
alt.maxxSrc == opList[iBest].maxxSrc &&
alt.maxySrc == opList[iBest].maxySrc) &&
!opList[iBest].isPriorityOp)) &&
!alt.isOffshore)) {
if (iBest < 0 ||
(((alt.accuracy >= 0 && alt.accuracy < bestAccuracy) ||
// If two operations have the same accuracy, use
// the one that has the smallest area
(alt.accuracy == bestAccuracy &&
alt.pseudoArea < opList[iBest].pseudoArea &&
!opList[iBest].isPriorityOp)) &&
!alt.isOffshore)) {

if (skipNonInstantiable && !alt.isInstantiable()) {
continue;
Expand Down Expand Up @@ -1713,6 +1705,16 @@ static PJ *add_coord_op_to_list(
double maxxDst;
double maxyDst;

double w = west_lon / 180 * M_PI;
double s = south_lat / 180 * M_PI;
double e = east_lon / 180 * M_PI;
double n = north_lat / 180 * M_PI;
if (w > e) {
e += 2 * M_PI;
}
// Integrate cos(lat) between south_lat and north_lat
const double pseudoArea = (e - w) * (std::sin(n) - std::sin(s));

if (pjSrcGeocentricToLonLat) {
minxSrc = west_lon;
minySrc = south_lat;
Expand Down Expand Up @@ -1740,8 +1742,8 @@ static PJ *add_coord_op_to_list(
const double accuracy = proj_coordoperation_get_accuracy(op->ctx, op);
altCoordOps.emplace_back(
idxInOriginalList, minxSrc, minySrc, maxxSrc, maxySrc, minxDst,
minyDst, maxxDst, maxyDst, op, name, accuracy, isOffshore,
pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
minyDst, maxxDst, maxyDst, op, name, accuracy, pseudoArea,
isOffshore, pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
op = nullptr;
}
return op;
Expand Down Expand Up @@ -2987,13 +2989,13 @@ PJCoordOperation::PJCoordOperation(
int idxInOriginalListIn, double minxSrcIn, double minySrcIn,
double maxxSrcIn, double maxySrcIn, double minxDstIn, double minyDstIn,
double maxxDstIn, double maxyDstIn, PJ *pjIn, const std::string &nameIn,
double accuracyIn, bool isOffshoreIn, const PJ *pjSrcGeocentricToLonLatIn,
const PJ *pjDstGeocentricToLonLatIn)
double accuracyIn, double pseudoAreaIn, bool isOffshoreIn,
const PJ *pjSrcGeocentricToLonLatIn, const PJ *pjDstGeocentricToLonLatIn)
: idxInOriginalList(idxInOriginalListIn), minxSrc(minxSrcIn),
minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn),
minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn),
maxyDst(maxyDstIn), pj(pjIn), name(nameIn), accuracy(accuracyIn),
isOffshore(isOffshoreIn),
pseudoArea(pseudoAreaIn), isOffshore(isOffshoreIn),
isPriorityOp(isSpecialCaseForNAD83_to_NAD83HARN(name) ||
isSpecialCaseForGDA94_to_WGS84(name) ||
isSpecialCaseForWGS84_to_GDA2020(name)),
Expand Down
3 changes: 2 additions & 1 deletion src/iso19111/c_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9110,7 +9110,8 @@ PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) {
alt.idxInOriginalList, minxSrc, minySrc, maxxSrc,
maxySrc, minxDst, minyDst, maxxDst, maxyDst,
pjNormalized, co->nameStr(), alt.accuracy,
alt.isOffshore, alt.pjSrcGeocentricToLonLat,
alt.pseudoArea, alt.isOffshore,
alt.pjSrcGeocentricToLonLat,
alt.pjDstGeocentricToLonLat);
}
}
Expand Down
135 changes: 121 additions & 14 deletions src/iso19111/operation/coordinateoperationfactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1384,16 +1384,23 @@ struct FilterResults {
}

#if 0
std::cerr << op->nameStr() << " ";
std::cerr << area << " ";
std::cerr << getAccuracy(op) << " ";
std::cerr << isPROJExportable << " ";
std::cerr << hasGrids << " ";
std::cerr << gridsAvailable << " ";
std::cerr << gridsKnown << " ";
std::cerr << stepCount << " ";
std::cerr << op->hasBallparkTransformation() << " ";
std::cerr << isNullTransformation(op->nameStr()) << " ";
std::cerr << "name=" << op->nameStr() << " ";
std::cerr << "area=" << area << " ";
std::cerr << "accuracy=" << getAccuracy(op) << " ";
std::cerr << "isPROJExportable=" << isPROJExportable << " ";
std::cerr << "hasGrids=" << hasGrids << " ";
std::cerr << "gridsAvailable=" << gridsAvailable << " ";
std::cerr << "gridsKnown=" << gridsKnown << " ";
std::cerr << "stepCount=" << stepCount << " ";
std::cerr << "projStepCount=" << projStepCount << " ";
std::cerr << "ballpark=" << op->hasBallparkTransformation() << " ";
std::cerr << "vertBallpark="
<< (op->nameStr().find(
BALLPARK_VERTICAL_TRANSFORMATION) !=
std::string::npos)
<< " ";
std::cerr << "isNull=" << isNullTransformation(op->nameStr())
<< " ";
std::cerr << std::endl;
#endif
map[op.get()] = PrecomputedOpCharacteristics(
Expand Down Expand Up @@ -2378,6 +2385,29 @@ struct MyPROJStringExportableHorizVerticalHorizPROJBased final
MyPROJStringExportableHorizVerticalHorizPROJBased::
~MyPROJStringExportableHorizVerticalHorizPROJBased() = default;

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

struct MyPROJStringExportableHorizNullVertical final
: public io::IPROJStringExportable {
CoordinateOperationPtr horizTransform{};

MyPROJStringExportableHorizNullVertical(
const CoordinateOperationPtr &horizTransformIn)
: horizTransform(horizTransformIn) {}

~MyPROJStringExportableHorizNullVertical() override;

void
// cppcheck-suppress functionStatic
_exportToPROJString(io::PROJStringFormatter *formatter) const override {

horizTransform->_exportToPROJString(formatter);
}
};

MyPROJStringExportableHorizNullVertical::
~MyPROJStringExportableHorizNullVertical() = default;

//! @endcond

} // namespace operation
Expand All @@ -2388,6 +2418,7 @@ namespace dropbox{ namespace oxygen {
template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableGeodToGeod>>::~nn() = default;
template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableHorizVertical>>::~nn() = default;
template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableHorizVerticalHorizPROJBased>>::~nn() = default;
template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableHorizNullVertical>>::~nn() = default;
}}
#endif

Expand Down Expand Up @@ -2662,6 +2693,40 @@ static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
nullptr, accuracies, hasBallparkTransformation);
}

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

static CoordinateOperationNNPtr createHorizNullVerticalPROJBased(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
const operation::CoordinateOperationNNPtr &horizTransform,
const operation::CoordinateOperationNNPtr &verticalTransform) {

auto exportable =
util::nn_make_shared<MyPROJStringExportableHorizNullVertical>(
horizTransform);

std::vector<CoordinateOperationNNPtr> ops = {horizTransform,
verticalTransform};
const std::string opName = computeConcatenatedName(ops);

auto properties = util::PropertyMap();
properties.set(common::IdentifiedObject::NAME_KEY, opName);

bool emptyIntersection = false;
auto extent = getExtent(ops, false, emptyIntersection);
if (extent) {
properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
NN_NO_CHECK(extent));
}

const auto remarks = getRemarks(ops);
if (!remarks.empty()) {
properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
}

return createPROJBased(properties, exportable, sourceCRS, targetCRS,
nullptr, {}, /*hasBallparkTransformation=*/true);
}

//! @endcond

// ---------------------------------------------------------------------------
Expand Down Expand Up @@ -4829,6 +4894,18 @@ void CoordinateOperationFactory::Private::createOperationsBoundToVert(

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

static std::string
getBallparkTransformationVertToVert(const crs::CRSNNPtr &sourceCRS,
const crs::CRSNNPtr &targetCRS) {
auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
name += " (";
name += BALLPARK_VERTICAL_TRANSFORMATION;
name += ')';
return name;
}

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

void CoordinateOperationFactory::Private::createOperationsVertToVert(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &context, const crs::VerticalCRS *vertSrc,
Expand Down Expand Up @@ -4874,10 +4951,8 @@ void CoordinateOperationFactory::Private::createOperationsVertToVert(
: metadata::Extent::WORLD);

if (!equivalentVDatum) {
auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
name += " (";
name += BALLPARK_VERTICAL_TRANSFORMATION;
name += ')';
const auto name =
getBallparkTransformationVertToVert(sourceCRS, targetCRS);
auto conv = Transformation::createChangeVerticalUnit(
map.set(common::IdentifiedObject::NAME_KEY, name), sourceCRS,
targetCRS,
Expand Down Expand Up @@ -6096,6 +6171,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
}

for (const auto &verticalTransform : verticalTransforms) {

auto interpolationGeogCRS = NN_NO_CHECK(srcGeog);
auto interpTransformCRS = verticalTransform->interpolationCRS();
if (interpTransformCRS) {
Expand Down Expand Up @@ -6138,6 +6214,37 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
}
}
}

if (verticalTransforms.size() == 1U &&
verticalTransform->hasBallparkTransformation() &&
context.context->getAuthorityFactory() &&
dynamic_cast<crs::ProjectedCRS *>(componentsSrc[0].get()) &&
dynamic_cast<crs::ProjectedCRS *>(componentsDst[0].get()) &&
verticalTransform->nameStr() ==
getBallparkTransformationVertToVert(componentsSrc[1],
componentsDst[1])) {
// e.g EPSG:3912+EPSG:5779 to EPSG:3794+EPSG:8690
// "MGI 1901 / Slovene National Grid + SVS2000 height" to
// "Slovenia 1996 / Slovene National Grid + SVS2010 height"
// using the "D48/GK to D96/TM (xx)" family of horizontal
// transformatoins between the projected CRS
// Cf
// https://github.com/OSGeo/PROJ/issues/3854#issuecomment-1689964773
// We restrict to a ballpark vertical transformation for now,
// but ideally we should deal with a regular vertical transformation
// but that would involve doing a transformation to its
// interpolation CRS and we don't have such cases for now.

std::vector<CoordinateOperationNNPtr> opsHoriz;
createOperationsFromDatabase(
componentsSrc[0], componentsDst[0], context, srcGeog.get(),
dstGeog.get(), srcGeog.get(), dstGeog.get(),
/*vertSrc=*/nullptr, /*vertDst=*/nullptr, opsHoriz);
for (const auto &opHoriz : opsHoriz) {
res.emplace_back(createHorizNullVerticalPROJBased(
sourceCRS, targetCRS, opHoriz, verticalTransform));
}
}
}

if (verticalTransforms.empty()) {
Expand Down
10 changes: 6 additions & 4 deletions src/proj_internal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,6 +342,7 @@ struct PJCoordOperation {
PJ *pj = nullptr;
std::string name{};
double accuracy = -1.0;
double pseudoArea = 0.0;
bool isOffshore = false;
bool isPriorityOp = false;
bool srcIsLonLatDegree = false;
Expand All @@ -363,7 +364,7 @@ struct PJCoordOperation {
double minySrcIn, double maxxSrcIn, double maxySrcIn,
double minxDstIn, double minyDstIn, double maxxDstIn,
double maxyDstIn, PJ *pjIn, const std::string &nameIn,
double accuracyIn, bool isOffshoreIn,
double accuracyIn, double pseudoAreaIn, bool isOffshoreIn,
const PJ *pjSrcGeocentricToLonLatIn,
const PJ *pjDstGeocentricToLonLatIn);

Expand All @@ -376,7 +377,8 @@ struct PJCoordOperation {
minyDst(other.minyDst), maxxDst(other.maxxDst),
maxyDst(other.maxyDst), pj(proj_clone(ctx, other.pj)),
name(std::move(other.name)), accuracy(other.accuracy),
isOffshore(other.isOffshore), isPriorityOp(other.isPriorityOp),
pseudoArea(other.pseudoArea), isOffshore(other.isOffshore),
isPriorityOp(other.isPriorityOp),
srcIsLonLatDegree(other.srcIsLonLatDegree),
srcIsLatLonDegree(other.srcIsLatLonDegree),
dstIsLonLatDegree(other.dstIsLonLatDegree),
Expand All @@ -396,8 +398,8 @@ struct PJCoordOperation {
maxySrc(other.maxySrc), minxDst(other.minxDst),
minyDst(other.minyDst), maxxDst(other.maxxDst),
maxyDst(other.maxyDst), name(std::move(other.name)),
accuracy(other.accuracy), isOffshore(other.isOffshore),
isPriorityOp(other.isPriorityOp),
accuracy(other.accuracy), pseudoArea(other.pseudoArea),
isOffshore(other.isOffshore), isPriorityOp(other.isPriorityOp),
srcIsLonLatDegree(other.srcIsLonLatDegree),
srcIsLatLonDegree(other.srcIsLatLonDegree),
dstIsLonLatDegree(other.dstIsLonLatDegree),
Expand Down
12 changes: 12 additions & 0 deletions test/cli/testvarious
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1320,6 +1320,18 @@ $EXE -d 3 EPSG:25831 EPSG:23031 -E >>${OUT} 2>&1 <<EOF
299905.060 4499796.515
EOF

echo "##############################################################" >> ${OUT}
echo "Test Similarity Transformation through CompoundCRS" >> ${OUT}
# Cf https://github.com/OSGeo/PROJ/issues/3854#issuecomment-1689964773
#
$EXE -d 3 EPSG:3912 EPSG:3794 -E >>${OUT} 2>&1 <<EOF
477134.28 95134.21
EOF

$EXE -d 3 EPSG:3912+EPSG:5779 EPSG:3794+EPSG:8690 -E >>${OUT} 2>&1 <<EOF
477134.28 95134.21 5
EOF


# Done!
# do 'diff' with distribution results
Expand Down
4 changes: 4 additions & 0 deletions test/cli/tv_out.dist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -633,3 +633,7 @@ Test Similarity Transformation (example from EPSG Guidance Note 7.2)
##############################################################
Test inverse of Similarity Transformation (example from EPSG Guidance Note 7.2)
299905.060 4499796.515 300000.000 4500000.000 0.000
##############################################################
Test Similarity Transformation through CompoundCRS
477134.28 95134.21 476763.303 95620.222 0.000
477134.28 95134.21 5 476763.303 95620.222 5.000

0 comments on commit d750d68

Please sign in to comment.