Skip to content

Commit

Permalink
pj_get_suggested_operation(): handle longitudes outside of [-180,180]…
Browse files Browse the repository at this point in the history 
… for coordinate operation selection (fixes #3594)
  • Loading branch information
@rouault
rouault committed Jan 26, 2023
1 parent f7e8883 commit 32b04ef
Show file tree
Hide file tree
Showing 3 changed files with 176 additions and 8 deletions.
96 changes: 96 additions & 0 deletions src/4D_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,19 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
PJ_COORD coord)
/**************************************************************************************/
{
const auto normalizeLongitude = [](double x) {
if (x > 180.0) {
x -= 360.0;
if (x > 180.0)
x = fmod(x + 180.0, 360.0) - 180.0;
} else if (x < -180.0) {
x += 360.0;
if (x < -180.0)
x = fmod(x + 180.0, 360.0) - 180.0;
}
return x;
};

// Select the operations that match the area of use
// and has the best accuracy.
int iBest = -1;
Expand All @@ -228,11 +241,39 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
if (coord.xyzt.x >= alt.minxSrc && coord.xyzt.y >= alt.minySrc &&
coord.xyzt.x <= alt.maxxSrc && coord.xyzt.y <= alt.maxySrc) {
spatialCriterionOK = true;
} else if (alt.srcIsLonLatDegree && coord.xyzt.y >= alt.minySrc &&
coord.xyzt.y <= alt.maxySrc) {
const double normalizedLon = normalizeLongitude(coord.xyzt.x);
if (normalizedLon >= alt.minxSrc &&
normalizedLon <= alt.maxxSrc) {
spatialCriterionOK = true;
}
} else if (alt.srcIsLatLonDegree && coord.xyzt.x >= alt.minxSrc &&
coord.xyzt.x <= alt.maxxSrc) {
const double normalizedLon = normalizeLongitude(coord.xyzt.y);
if (normalizedLon >= alt.minySrc &&
normalizedLon <= alt.maxySrc) {
spatialCriterionOK = true;
}
}
} else {
if (coord.xyzt.x >= alt.minxDst && coord.xyzt.y >= alt.minyDst &&
coord.xyzt.x <= alt.maxxDst && coord.xyzt.y <= alt.maxyDst) {
spatialCriterionOK = true;
} else if (alt.dstIsLonLatDegree && coord.xyzt.y >= alt.minyDst &&
coord.xyzt.y <= alt.maxyDst) {
const double normalizedLon = normalizeLongitude(coord.xyzt.x);
if (normalizedLon >= alt.minxDst &&
normalizedLon <= alt.maxxDst) {
spatialCriterionOK = true;
}
} else if (alt.dstIsLatLonDegree && coord.xyzt.x >= alt.minxDst &&
coord.xyzt.x <= alt.maxxDst) {
const double normalizedLon = normalizeLongitude(coord.xyzt.y);
if (normalizedLon >= alt.minyDst &&
normalizedLon <= alt.maxyDst) {
spatialCriterionOK = true;
}
}
}

Expand Down Expand Up @@ -2589,3 +2630,58 @@ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {

return factors;
}

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

//! @cond Doxygen_Suppress
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)
: 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) {

const auto IsLonLatOrLatLon = [](const PJ *crs, bool &isLonLatDegreeOut,
bool &isLatLonDegreeOut) {
const auto eType = proj_get_type(crs);
if (eType == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
eType == PJ_TYPE_GEOGRAPHIC_3D_CRS) {
const auto cs = proj_crs_get_coordinate_system(crs->ctx, crs);
const char *direction = "";
double conv_factor = 0;
constexpr double EPS = 1e-14;
if (proj_cs_get_axis_info(crs->ctx, cs, 0, nullptr, nullptr,
&direction, &conv_factor, nullptr,
nullptr, nullptr) &&
ci_equal(direction, "East")) {
isLonLatDegreeOut = fabs(conv_factor - M_PI / 180) < EPS;
} else if (proj_cs_get_axis_info(crs->ctx, cs, 1, nullptr, nullptr,
&direction, &conv_factor, nullptr,
nullptr, nullptr) &&
ci_equal(direction, "East")) {
isLatLonDegreeOut = fabs(conv_factor - M_PI / 180) < EPS;
}
proj_destroy(cs);
}
};

const auto source = proj_get_source_crs(pj->ctx, pj);
if (source) {
IsLonLatOrLatLon(source, srcIsLonLatDegree, srcIsLatLonDegree);
proj_destroy(source);
}

const auto target = proj_get_target_crs(pj->ctx, pj);
if (target) {
IsLonLatOrLatLon(target, dstIsLonLatDegree, dstIsLatLonDegree);
proj_destroy(target);
}
}

//! @endcond
23 changes: 15 additions & 8 deletions src/proj_internal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -312,17 +312,16 @@ struct PJCoordOperation {
std::string name{};
double accuracy = -1.0;
bool isOffshore = false;
bool srcIsLonLatDegree = false;
bool srcIsLatLonDegree = false;
bool dstIsLonLatDegree = false;
bool dstIsLatLonDegree = false;

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)
: 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) {}
double accuracyIn, bool isOffshoreIn);

PJCoordOperation(const PJCoordOperation &) = delete;

Expand All @@ -333,15 +332,23 @@ 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) {}
isOffshore(other.isOffshore),
srcIsLonLatDegree(other.srcIsLonLatDegree),
srcIsLatLonDegree(other.srcIsLatLonDegree),
dstIsLonLatDegree(other.dstIsLonLatDegree),
dstIsLatLonDegree(other.dstIsLatLonDegree) {}

PJCoordOperation(PJCoordOperation &&other)
: idxInOriginalList(other.idxInOriginalList), minxSrc(other.minxSrc),
minySrc(other.minySrc), maxxSrc(other.maxxSrc),
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) {
accuracy(other.accuracy), isOffshore(other.isOffshore),
srcIsLonLatDegree(other.srcIsLonLatDegree),
srcIsLatLonDegree(other.srcIsLatLonDegree),
dstIsLonLatDegree(other.dstIsLonLatDegree),
dstIsLatLonDegree(other.dstIsLatLonDegree) {
pj = other.pj;
other.pj = nullptr;
}
Expand Down
65 changes: 65 additions & 0 deletions test/unit/gie_self_tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1041,6 +1041,71 @@ TEST(gie, proj_create_crs_to_crs_with_area_large) {

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

TEST(gie, proj_create_crs_to_crs_with_longitude_outside_minus_180_180) {

// Test bugfix for /~https://github.com/OSGeo/PROJ/issues/3594
auto P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, "EPSG:4277", "EPSG:4326",
nullptr);
ASSERT_TRUE(P != nullptr);
PJ_COORD c;

c.xyzt.x = 50; // Lat in deg
c.xyzt.y = -2 + 360; // Long in deg
c.xyzt.z = 0;
c.xyzt.t = HUGE_VAL;
c = proj_trans(P, PJ_FWD, c);
EXPECT_NEAR(c.xy.x, 50.00065628, 1e-8);
EXPECT_NEAR(c.xy.y, -2.00133989, 1e-8);

c.xyzt.x = 50; // Lat in deg
c.xyzt.y = -2 - 360; // Long in deg
c.xyzt.z = 0;
c.xyzt.t = HUGE_VAL;
c = proj_trans(P, PJ_FWD, c);
EXPECT_NEAR(c.xy.x, 50.00065628, 1e-8);
EXPECT_NEAR(c.xy.y, -2.00133989, 1e-8);

c.xyzt.x = 50.00065628; // Lat in deg
c.xyzt.y = -2.00133989 + 360; // Long in deg
c.xyzt.z = 0;
c.xyzt.t = HUGE_VAL;
c = proj_trans(P, PJ_INV, c);
EXPECT_NEAR(c.xy.x, 50, 1e-8);
EXPECT_NEAR(c.xy.y, -2, 1e-8);

c.xyzt.x = 50.00065628; // Lat in deg
c.xyzt.y = -2.00133989 - 360; // Long in deg
c.xyzt.z = 0;
c.xyzt.t = HUGE_VAL;
c = proj_trans(P, PJ_INV, c);
EXPECT_NEAR(c.xy.x, 50, 1e-8);
EXPECT_NEAR(c.xy.y, -2, 1e-8);

auto Pnormalized = proj_normalize_for_visualization(PJ_DEFAULT_CTX, P);

c.xyzt.x = -2 + 360; // Long in deg
c.xyzt.y = 50; // Lat in deg
c.xyzt.z = 0;
c.xyzt.t = HUGE_VAL;
c = proj_trans(Pnormalized, PJ_FWD, c);
EXPECT_NEAR(c.xy.x, -2.00133989, 1e-8);
EXPECT_NEAR(c.xy.y, 50.00065628, 1e-8);

c.xyzt.x = -2.00133989 + 360; // Long in deg
c.xyzt.y = 50.00065628; // Lat in deg
c.xyzt.z = 0;
c.xyzt.t = HUGE_VAL;
c = proj_trans(Pnormalized, PJ_INV, c);
EXPECT_NEAR(c.xy.x, -2, 1e-8);
EXPECT_NEAR(c.xy.y, 50, 1e-8);

proj_destroy(Pnormalized);

proj_destroy(P);
}

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

TEST(gie, proj_trans_generic) {
// GDA2020 to WGS84 (G1762)
auto P = proj_create(
Expand Down

0 comments on commit 32b04ef

Please sign in to comment.