diff --git a/src/projections/imoll.cpp b/src/projections/imoll.cpp index 3a1e89fb55..c312e8af54 100644 --- a/src/projections/imoll.cpp +++ b/src/projections/imoll.cpp @@ -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; @@ -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 */ @@ -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; @@ -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; diff --git a/src/projections/imoll_o.cpp b/src/projections/imoll_o.cpp index 5b9d27f7dc..a3adfb35f6 100644 --- a/src/projections/imoll_o.cpp +++ b/src/projections/imoll_o.cpp @@ -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 */ @@ -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; @@ -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; @@ -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; diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 0c29235325..98ea577449 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -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 @@ -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