Skip to content

Commit

Permalink
Fix imoll and imoll_o zone calculations to correct inverse transforma…
Browse files Browse the repository at this point in the history 
…tions near the "seams". (#4159)

Fixes #4155

* Fix x zone boundaries when doing inverse imoll transform.

* Fix x zone boundaries when doing inverse imoll_o transform.

* Add roundtrip tests for imoll near the "seams".

* Add roundtrip tests for imoll_o near the "seams".
  • Loading branch information
@erykoff
erykoff authored Jun 3, 2024
1 parent 82af828 commit 929f946
Show file tree
Hide file tree
Showing 3 changed files with 113 additions and 8 deletions.
36 changes: 32 additions & 4 deletions src/projections/imoll.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,11 @@ C_NAMESPACE PJ *pj_moll(PJ *);
namespace pj_imoll_ns {
struct pj_imoll_data {
struct PJconsts *pj[6];
// We need to know the inverse boundary locations of the "seams".
double boundary12;
double boundary34;
double boundary45;
double boundary56;
};

constexpr double d20 = 20 * DEG_TO_RAD;
Expand Down Expand Up @@ -92,13 +97,13 @@ static PJ_LP imoll_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
if (xy.y > y90 + EPSLN || xy.y < -y90 + EPSLN) /* 0 */
z = 0;
else if (xy.y >= 0)
z = (xy.x <= -d40 ? 1 : 2); /* 1|2 */
z = (xy.x <= (Q->boundary12) ? 1 : 2); /* 1|2 */
else {
if (xy.x <= -d100)
if (xy.x <= Q->boundary34)
z = 3; /* 3 */
else if (xy.x <= -d20)
else if (xy.x <= Q->boundary45)
z = 4; /* 4 */
else if (xy.x <= d80)
else if (xy.x <= Q->boundary56)
z = 5; /* 5 */
else
z = 6; /* 6 */
Expand Down Expand Up @@ -218,6 +223,19 @@ static double compute_zone_offset(struct pj_imoll_ns::pj_imoll_data *Q,
return (xy2.x + Q->pj[zone2 - 1]->x0) - (xy1.x + Q->pj[zone1 - 1]->x0);
}

static double compute_zone_x_boundary(PJ *P, double lam, double phi) {
PJ_LP lp1, lp2;
PJ_XY xy1, xy2;

lp1.lam = lam - pj_imoll_ns::EPSLN;
lp1.phi = phi;
lp2.lam = lam + pj_imoll_ns::EPSLN;
lp2.phi = phi;
xy1 = imoll_s_forward(lp1, P);
xy2 = imoll_s_forward(lp2, P);
return (xy1.x + xy2.x)/2.;
}

PJ *PJ_PROJECTION(imoll) {
using namespace pj_imoll_ns;

Expand Down Expand Up @@ -258,6 +276,16 @@ PJ *PJ_PROJECTION(imoll) {
/* Match 6 (south) to 2 (north) */
Q->pj[5]->x0 += compute_zone_offset(Q, 6, 2, d80, 0.0 - EPSLN, 0.0 + EPSLN);

/*
The most straightforward way of computing the x locations of the "seams"
in the interrupted projection is to compute the forward transform at the
seams and record these values.
*/
Q->boundary12 = compute_zone_x_boundary(P, -d40, 0.0 + EPSLN);
Q->boundary34 = compute_zone_x_boundary(P, -d100, 0.0 - EPSLN);
Q->boundary45 = compute_zone_x_boundary(P, -d20, 0.0 - EPSLN);
Q->boundary56 = compute_zone_x_boundary(P, d80, 0.0 - EPSLN);

P->inv = imoll_s_inverse;
P->fwd = imoll_s_forward;
P->destructor = pj_imoll_destructor;
Expand Down
37 changes: 33 additions & 4 deletions src/projections/imoll_o.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,11 @@ C_NAMESPACE PJ *pj_moll(PJ *);
namespace pj_imoll_o_ns {
struct pj_imoll_o_data {
struct PJconsts *pj[6];
// We need to know the inverse boundary locations of the "seams".
double boundary12;
double boundary23;
double boundary45;
double boundary56;
};

/* SIMPLIFY THIS */
Expand Down Expand Up @@ -100,16 +105,16 @@ static PJ_LP imoll_o_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
if (xy.y > y90 + EPSLN || xy.y < -y90 + EPSLN) /* 0 */
z = 0;
else if (xy.y >= 0) {
if (xy.x <= -d90)
if (xy.x <= Q->boundary12)
z = 1;
else if (xy.x >= d60)
else if (xy.x >= Q->boundary23)
z = 3;
else
z = 2;
} else {
if (xy.x <= -d60)
if (xy.x <= Q->boundary45)
z = 4;
else if (xy.x >= d90)
else if (xy.x >= Q->boundary56)
z = 6;
else
z = 5;
Expand Down Expand Up @@ -232,6 +237,20 @@ pj_imoll_o_compute_zone_offset(struct pj_imoll_o_ns::pj_imoll_o_data *Q,
return (xy2.x + Q->pj[zone2 - 1]->x0) - (xy1.x + Q->pj[zone1 - 1]->x0);
}

static double
pj_imoll_o_compute_zone_x_boundary(PJ *P, double lam, double phi) {
PJ_LP lp1, lp2;
PJ_XY xy1, xy2;

lp1.lam = lam - pj_imoll_o_ns::EPSLN;
lp1.phi = phi;
lp2.lam = lam + pj_imoll_o_ns::EPSLN;
lp2.phi = phi;
xy1 = imoll_o_s_forward(lp1, P);
xy2 = imoll_o_s_forward(lp2, P);
return (xy1.x + xy2.x)/2.;
}

PJ *PJ_PROJECTION(imoll_o) {
using namespace pj_imoll_o_ns;

Expand Down Expand Up @@ -273,6 +292,16 @@ PJ *PJ_PROJECTION(imoll_o) {
Q->pj[5]->x0 +=
pj_imoll_o_compute_zone_offset(Q, 6, 3, d90, 0.0 - EPSLN, 0.0 + EPSLN);

/*
The most straightforward way of computing the x locations of the "seams"
in the interrupted projection is to compute the forward transform at the
seams and record these values.
*/
Q->boundary12 = pj_imoll_o_compute_zone_x_boundary(P, -d90, 0.0 + EPSLN);
Q->boundary23 = pj_imoll_o_compute_zone_x_boundary(P, d60, 0.0 + EPSLN);
Q->boundary45 = pj_imoll_o_compute_zone_x_boundary(P, -d60, 0.0 - EPSLN);
Q->boundary56 = pj_imoll_o_compute_zone_x_boundary(P, d90, 0.0 - EPSLN);

P->inv = imoll_o_s_inverse;
P->fwd = imoll_o_s_forward;
P->destructor = pj_imoll_o_destructor;
Expand Down
48 changes: 48 additions & 0 deletions test/gie/builtins.gie
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2575,6 +2575,30 @@ roundtrip 1
accept -2 -1
expect -1314402.165573241888 -124066.283433859542
roundtrip 1
accept -39.99 0.1
expect -5135117.0707450127 12406.8672748194
roundtrip 1
accept -40.01 0.1
expect -5137140.6776947584 12406.8672748194
roundtrip 1
accept -99.99 -0.1
expect -11169097.7713819221 -12406.8672748194
roundtrip 1
accept -100.01 -0.1
expect -11171118.5438199658 -12406.8672748194
roundtrip 1
accept -19.99 -0.1
expect -3123793.9498816459 -12406.8672748194
roundtrip 1
accept -20.01 -0.1
expect -3125812.8326452221 -12406.8672748194
roundtrip 1
accept 79.99 -0.1
expect 6930815.0545556545 -12406.8672748194
roundtrip 1
accept 80.01 -0.1
expect 6932837.7166681662 -12406.8672748194
roundtrip 1

accept -100.0 22.0
expect -11170107.212763708085 2703699.326638640370
Expand Down Expand Up @@ -2627,6 +2651,30 @@ roundtrip 1
accept -2 -1
expect -1759793.139928588411 -124066.283433859542
roundtrip 1
accept -89.99 0.1
expect -10608821.9887007959 12406.8672748194
roundtrip 1
accept -90.01 0.1
expect -10610845.5956505425 12406.8672748194
roundtrip 1
accept 59.99 0.1
expect 4474097.1799880061 12406.8672748194
roundtrip 1
accept 60.01 0.1
expect 4476121.7317749839 12406.8672748194
roundtrip 1
accept -59.99 -0.1
expect -7591833.0556381932 -12406.8672748194
roundtrip 1
accept -60.01 -0.1
expect -7593856.6625879407 -12406.8672748194
roundtrip 1
accept 89.99 -0.1
expect 7491086.1130506080 -12406.8672748194
roundtrip 1
accept 90.01 -0.1
expect 7493109.7200003546 -12406.8672748194
roundtrip 1

accept -140.0 22.0
expect -15638150.097869191319 2703699.326638640370
Expand Down

0 comments on commit 929f946

Please sign in to comment.