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.

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

createOperations(): take into account axis unit and inversion of target DerivedProjectedCRS #3281

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
9 changes: 8 additions & 1 deletion include/proj/crs.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -661,6 +661,10 @@ class PROJ_GCC_DLL ProjectedCRS final : public DerivedCRS,
addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,
bool axisSpecFound) const;

PROJ_INTERNAL static void addUnitConvertAndAxisSwap(
const std::vector<cs::CoordinateSystemAxisNNPtr> &axisListIn,
io::PROJStringFormatter *formatter, bool axisSpecFound);

PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter)
const override; // throw(io::FormattingException)

Expand Down Expand Up @@ -1272,7 +1276,10 @@ class PROJ_GCC_DLL DerivedProjectedCRS final : public DerivedCRS {
//! @cond Doxygen_Suppress
PROJ_INTERNAL void _exportToWKT(io::WKTFormatter *formatter)
const override; // throw(io::FormattingException)
//! @endcond

PROJ_INTERNAL void
addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter) const;
//! @endcond

protected:
PROJ_INTERNAL
Expand Down
33 changes: 25 additions & 8 deletions src/iso19111/crs.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4474,9 +4474,16 @@ ProjectedCRS::alterParametersLinearUnit(const common::UnitOfMeasure &unit,
//! @cond Doxygen_Suppress
void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,
bool axisSpecFound) const {
const auto &axisList = d->coordinateSystem()->axisList();
const auto &unit = axisList[0]->unit();
const auto *zUnit = axisList.size() == 3 ? &(axisList[2]->unit()) : nullptr;
ProjectedCRS::addUnitConvertAndAxisSwap(d->coordinateSystem()->axisList(),
formatter, axisSpecFound);
}

void ProjectedCRS::addUnitConvertAndAxisSwap(
const std::vector<cs::CoordinateSystemAxisNNPtr> &axisListIn,
io::PROJStringFormatter *formatter, bool axisSpecFound) {
const auto &unit = axisListIn[0]->unit();
const auto *zUnit =
axisListIn.size() == 3 ? &(axisListIn[2]->unit()) : nullptr;
if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE,
util::IComparable::Criterion::EQUIVALENT) ||
(zUnit &&
Expand Down Expand Up @@ -4518,8 +4525,8 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,

if (!axisSpecFound &&
(!formatter->getCRSExport() || formatter->getLegacyCRSToCRSContext())) {
const auto &dir0 = axisList[0]->direction();
const auto &dir1 = axisList[1]->direction();
const auto &dir0 = axisListIn[0]->direction();
const auto &dir1 = axisListIn[1]->direction();
if (!(&dir0 == &cs::AxisDirection::EAST &&
&dir1 == &cs::AxisDirection::NORTH) &&
// For polar projections, that have south+south direction,
Expand All @@ -4528,7 +4535,7 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,

const char *order[2] = {nullptr, nullptr};
for (int i = 0; i < 2; i++) {
const auto &dir = axisList[i]->direction();
const auto &dir = axisListIn[i]->direction();
if (&dir == &cs::AxisDirection::WEST)
order[i] = "-1";
else if (&dir == &cs::AxisDirection::EAST)
Expand All @@ -4546,8 +4553,8 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,
formatter->addParam("order", orderStr);
}
} else {
const auto &name0 = axisList[0]->nameStr();
const auto &name1 = axisList[1]->nameStr();
const auto &name0 = axisListIn[0]->nameStr();
const auto &name1 = axisListIn[1]->nameStr();
const bool northingEasting = ci_starts_with(name0, "northing") &&
ci_starts_with(name1, "easting");
// case of EPSG:32661 ["WGS 84 / UPS North (N,E)]"
Expand Down Expand Up @@ -6548,6 +6555,16 @@ bool DerivedProjectedCRS::_isEquivalentTo(

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

//! @cond Doxygen_Suppress
void DerivedProjectedCRS::addUnitConvertAndAxisSwap(
io::PROJStringFormatter *formatter) const {
ProjectedCRS::addUnitConvertAndAxisSwap(coordinateSystem()->axisList(),
formatter, false);
}
//! @endcond

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

//! @cond Doxygen_Suppress
struct TemporalCRS::Private {};
//! @endcond
Expand Down
12 changes: 11 additions & 1 deletion src/iso19111/operation/conversion.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4086,8 +4086,12 @@ void Conversion::_exportToPROJString(
}
}

auto derivedProjCRS =
dynamic_cast<const crs::DerivedProjectedCRS *>(horiz);

// horiz != nullptr: only to make clang static analyzer happy
if (!bEllipsoidParametersDone && horiz != nullptr) {
if (!bEllipsoidParametersDone && horiz != nullptr &&
derivedProjCRS == nullptr) {
auto targetGeodCRS = horiz->extractGeodeticCRS();
auto targetGeogCRS =
std::dynamic_pointer_cast<crs::GeographicCRS>(targetGeodCRS);
Expand Down Expand Up @@ -4119,6 +4123,12 @@ void Conversion::_exportToPROJString(
if (projCRS->hasOver()) {
formatter->addParam("over");
}
} else {
if (derivedProjCRS) {
formatter->pushOmitZUnitConversion();
derivedProjCRS->addUnitConvertAndAxisSwap(formatter);
formatter->popOmitZUnitConversion();
}
}

auto derivedGeographicCRS =
Expand Down
62 changes: 62 additions & 0 deletions test/unit/test_operationfactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8411,3 +8411,65 @@ TEST(operation, createOperation_ossfuzz_47873_simplified_if_i_might_say) {
} catch (const std::exception &) {
}
}

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

TEST(operation, createOperation_derived_projected_crs) {

auto authFactory =
AuthorityFactory::create(DatabaseContext::create(), "EPSG");
auto src = authFactory->createCoordinateReferenceSystem("6507");

auto wkt =
"DERIVEDPROJCRS[\"Custom Site Calibrated CRS\",\n"
" BASEPROJCRS[\"NAD83(2011) / Mississippi East (ftUS)\",\n"
" BASEGEOGCRS[\"NAD83(2011)\",\n"
" DATUM[\"NAD83 (National Spatial Reference System "
"2011)\",\n"
" ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n"
" LENGTHUNIT[\"metre\",1]]],\n"
" PRIMEM[\"Greenwich\",0,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433]]],\n"
" CONVERSION[\"SPCS83 Mississippi East zone (US Survey "
"feet)\",\n"
" METHOD[\"Transverse Mercator\",\n"
" ID[\"EPSG\",9807]],\n"
" PARAMETER[\"Latitude of natural origin\",29.5,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8801]],\n"
" PARAMETER[\"Longitude of natural "
"origin\",-88.8333333333333,\n"
" ANGLEUNIT[\"degree\",0.0174532925199433],\n"
" ID[\"EPSG\",8802]],\n"
" PARAMETER[\"Scale factor at natural origin\",0.99995,\n"
" SCALEUNIT[\"unity\",1],\n"
" ID[\"EPSG\",8805]],\n"
" PARAMETER[\"False easting\",984250,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8806]],\n"
" PARAMETER[\"False northing\",0,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219],\n"
" ID[\"EPSG\",8807]]]],\n"
" DERIVINGCONVERSION[\"Affine transformation as PROJ-based\",\n"
" METHOD[\"PROJ-based operation method: "
"+proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=us-ft "
"+step +proj=affine +xoff=20 "
"+step +proj=unitconvert +xy_in=us-ft +xy_out=m\"]],\n"
" CS[Cartesian,2],\n"
" AXIS[\"northing (Y)\",north,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n"
" AXIS[\"easting (X)\",east,\n"
" LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n"
" REMARK[\"EPSG:6507 with 20 feet offset and axis inversion\"]]";
auto objDst = WKTParser().createFromWKT(wkt);
auto dst = nn_dynamic_pointer_cast<CRS>(objDst);
ASSERT_TRUE(dst != nullptr);

auto op = CoordinateOperationFactory::create()->createOperation(
src, NN_NO_CHECK(dst));
ASSERT_TRUE(op != nullptr);
EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()),
"+proj=pipeline "
"+step +proj=affine +xoff=20 "
"+step +proj=axisswap +order=2,1");
}