/src/proj/src/projections/imoll.cpp

Source (jump to first uncovered line)


#include <errno.h>
#include <math.h>

#include "proj.h"
#include "proj_internal.h"

PROJ_HEAD(imoll, "Interrupted Mollweide") "\n\tPCyl, Sph";

/*
This projection is a compilation of 6 separate sub-projections.
It is related to the Interrupted Goode Homolosine projection,
but uses Mollweide at all latitudes.

Each sub-projection is assigned an integer label
numbered 1 through 6. Most of this code contains logic to assign
the labels based on latitude (phi) and longitude (lam) regions.
The code is adapted from igh.cpp.

Original Reference:
J. Paul Goode (1919) STUDIES IN PROJECTIONS: ADAPTING THE
    HOMOLOGRAPHIC PROJECTION TO THE PORTRAYAL OF THE EARTH'S
    ENTIRE SURFACE, Bul. Geog. SOC.Phila., Vol. XWIJNO.3.
    July, 1919, pp. 103-113.
*/

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;
constexpr double d30 = 30 * DEG_TO_RAD;
constexpr double d40 = 40 * DEG_TO_RAD;
constexpr double d60 = 60 * DEG_TO_RAD;
constexpr double d80 = 80 * DEG_TO_RAD;
constexpr double d100 = 100 * DEG_TO_RAD;
constexpr double d140 = 140 * DEG_TO_RAD;
constexpr double d160 = 160 * DEG_TO_RAD;
constexpr double d180 = 180 * DEG_TO_RAD;

constexpr double EPSLN =
    1.e-10; /* allow a little 'slack' on zone edge positions */

} // namespace pj_imoll_ns

/*
Assign an integer index representing each of the 6
sub-projection zones based on latitude (phi) and
longitude (lam) ranges.
*/

static PJ_XY imoll_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
    using namespace pj_imoll_ns;

    PJ_XY xy;
    struct pj_imoll_data *Q = static_cast<struct pj_imoll_data *>(P->opaque);
    int z;

    if (lp.phi >= 0) { /* 1|2 */
        z = (lp.lam <= -d40 ? 1 : 2);
    } else { /* 3|4|5|6 */
        if (lp.lam <= -d100)
            z = 3; /* 3 */
        else if (lp.lam <= -d20)
            z = 4; /* 4 */
        else if (lp.lam <= d80)
            z = 5; /* 5 */
        else
            z = 6; /* 6 */
    }

    lp.lam -= Q->pj[z - 1]->lam0;
    xy = Q->pj[z - 1]->fwd(lp, Q->pj[z - 1]);
    xy.x += Q->pj[z - 1]->x0;
    xy.y += Q->pj[z - 1]->y0;

    return xy;
}

static PJ_LP imoll_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    using namespace pj_imoll_ns;

    PJ_LP lp = {0.0, 0.0};
    struct pj_imoll_data *Q = static_cast<struct pj_imoll_data *>(P->opaque);
    const double y90 = sqrt(2.0); /* lt=90 corresponds to y=sqrt(2) */

    int z = 0;
    if (xy.y > y90 + EPSLN || xy.y < -y90 + EPSLN) /* 0 */
        z = 0;
    else if (xy.y >= 0)
        z = (xy.x <= (Q->boundary12) ? 1 : 2); /* 1|2 */
    else {
        if (xy.x <= Q->boundary34)
            z = 3; /* 3 */
        else if (xy.x <= Q->boundary45)
            z = 4; /* 4 */
        else if (xy.x <= Q->boundary56)
            z = 5; /* 5 */
        else
            z = 6; /* 6 */
    }

    if (z) {
        bool ok = false;

        xy.x -= Q->pj[z - 1]->x0;
        xy.y -= Q->pj[z - 1]->y0;
        lp = Q->pj[z - 1]->inv(xy, Q->pj[z - 1]);
        lp.lam += Q->pj[z - 1]->lam0;

        switch (z) {
        case 1:
            ok = ((lp.lam >= -d180 - EPSLN && lp.lam <= -d40 + EPSLN) &&
                  (lp.phi >= 0.0 - EPSLN));
            break;
        case 2:
            ok = ((lp.lam >= -d40 - EPSLN && lp.lam <= d180 + EPSLN) &&
                  (lp.phi >= 0.0 - EPSLN));
            break;
        case 3:
            ok = ((lp.lam >= -d180 - EPSLN && lp.lam <= -d100 + EPSLN) &&
                  (lp.phi <= 0.0 + EPSLN));
            break;
        case 4:
            ok = ((lp.lam >= -d100 - EPSLN && lp.lam <= -d20 + EPSLN) &&
                  (lp.phi <= 0.0 + EPSLN));
            break;
        case 5:
            ok = ((lp.lam >= -d20 - EPSLN && lp.lam <= d80 + EPSLN) &&
                  (lp.phi <= 0.0 + EPSLN));
            break;
        case 6:
            ok = ((lp.lam >= d80 - EPSLN && lp.lam <= d180 + EPSLN) &&
                  (lp.phi <= 0.0 + EPSLN));
            break;
        }
        z = (!ok ? 0 : z); /* projectable? */
    }

    if (!z)
        lp.lam = HUGE_VAL;
    if (!z)
        lp.phi = HUGE_VAL;

    return lp;
}

static PJ *pj_imoll_destructor(PJ *P, int errlev) {
    using namespace pj_imoll_ns;

    int i;
    if (nullptr == P)
        return nullptr;

    if (nullptr == P->opaque)
        return pj_default_destructor(P, errlev);

    struct pj_imoll_data *Q = static_cast<struct pj_imoll_data *>(P->opaque);

    for (i = 0; i < 6; ++i) {
        if (Q->pj[i])
            Q->pj[i]->destructor(Q->pj[i], errlev);
    }

    return pj_default_destructor(P, errlev);
}

/*
  Zones:

    -180            -40                       180
      +--------------+-------------------------+
      |1             |2                        |
      |              |                         |
      |              |                         |
      |              |                         |
      |              |                         |
    0 +-------+------+-+-----------+-----------+
      |3      |4       |5          |6          |
      |       |        |           |           |
      |       |        |           |           |
      |       |        |           |           |
      |       |        |           |           |
      +-------+--------+-----------+-----------+
    -180    -100      -20         80          180
*/

static bool setup_zone(PJ *P, struct pj_imoll_ns::pj_imoll_data *Q, int n,
                       PJ *(*proj_ptr)(PJ *), double x_0, double y_0,
                       double lon_0) {
    if (!(Q->pj[n - 1] = proj_ptr(nullptr)))
        return false;
    if (!(Q->pj[n - 1] = proj_ptr(Q->pj[n - 1])))
        return false;
    Q->pj[n - 1]->ctx = P->ctx;
    Q->pj[n - 1]->x0 = x_0;
    Q->pj[n - 1]->y0 = y_0;
    Q->pj[n - 1]->lam0 = lon_0;
    return true;
}

static double compute_zone_offset(struct pj_imoll_ns::pj_imoll_data *Q,
                                  int zone1, int zone2, double lam, double phi1,
                                  double phi2) {
    PJ_LP lp1, lp2;
    PJ_XY xy1, xy2;

    lp1.lam = lam - (Q->pj[zone1 - 1]->lam0);
    lp1.phi = phi1;
    lp2.lam = lam - (Q->pj[zone2 - 1]->lam0);
    lp2.phi = phi2;
    xy1 = Q->pj[zone1 - 1]->fwd(lp1, Q->pj[zone1 - 1]);
    xy2 = Q->pj[zone2 - 1]->fwd(lp2, Q->pj[zone2 - 1]);
    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;

    struct pj_imoll_data *Q = static_cast<struct pj_imoll_data *>(
        calloc(1, sizeof(struct pj_imoll_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    /* Setup zones */
    if (!setup_zone(P, Q, 1, pj_moll, -d100, 0, -d100) ||
        !setup_zone(P, Q, 2, pj_moll, d30, 0, d30) ||
        !setup_zone(P, Q, 3, pj_moll, -d160, 0, -d160) ||
        !setup_zone(P, Q, 4, pj_moll, -d60, 0, -d60) ||
        !setup_zone(P, Q, 5, pj_moll, d20, 0, d20) ||
        !setup_zone(P, Q, 6, pj_moll, d140, 0, d140)) {
        return pj_imoll_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    }

    /* Adjust zones */

    /* Match 3 (south) to 1 (north) */
    Q->pj[2]->x0 +=
        compute_zone_offset(Q, 3, 1, -d160, 0.0 - EPSLN, 0.0 + EPSLN);

    /* Match 2 (north-east) to 1 (north-west) */
    Q->pj[1]->x0 +=
        compute_zone_offset(Q, 2, 1, -d40, 0.0 + EPSLN, 0.0 + EPSLN);

    /* Match 4 (south) to 1 (north) */
    Q->pj[3]->x0 +=
        compute_zone_offset(Q, 4, 1, -d100, 0.0 - EPSLN, 0.0 + EPSLN);

    /* Match 5 (south) to 2 (north) */
    Q->pj[4]->x0 +=
        compute_zone_offset(Q, 5, 2, -d20, 0.0 - EPSLN, 0.0 + EPSLN);

    /* 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;
    P->es = 0.;

    return P;
}

Coverage Report

Created: 2025-08-28 06:57

Line	Count	Source (jump to first uncovered line)
1
2
3		#include <errno.h>
4		#include <math.h>
5
6		#include "proj.h"
7		#include "proj_internal.h"
8
9		PROJ_HEAD(imoll, "Interrupted Mollweide") "\n\tPCyl, Sph";
10
11		/*
12		This projection is a compilation of 6 separate sub-projections.
13		It is related to the Interrupted Goode Homolosine projection,
14		but uses Mollweide at all latitudes.
15
16		Each sub-projection is assigned an integer label
17		numbered 1 through 6. Most of this code contains logic to assign
18		the labels based on latitude (phi) and longitude (lam) regions.
19		The code is adapted from igh.cpp.
20
21		Original Reference:
22		J. Paul Goode (1919) STUDIES IN PROJECTIONS: ADAPTING THE
23		HOMOLOGRAPHIC PROJECTION TO THE PORTRAYAL OF THE EARTH'S
24		ENTIRE SURFACE, Bul. Geog. SOC.Phila., Vol. XWIJNO.3.
25		July, 1919, pp. 103-113.
26		*/
27
28		C_NAMESPACE PJ pj_moll(PJ );
29
30		namespace pj_imoll_ns {
31		struct pj_imoll_data {
32		struct PJconsts *pj[6];
33		// We need to know the inverse boundary locations of the "seams".
34		double boundary12;
35		double boundary34;
36		double boundary45;
37		double boundary56;
38		};
39
40		constexpr double d20 = 20 * DEG_TO_RAD;
41		constexpr double d30 = 30 * DEG_TO_RAD;
42		constexpr double d40 = 40 * DEG_TO_RAD;
43		constexpr double d60 = 60 * DEG_TO_RAD;
44		constexpr double d80 = 80 * DEG_TO_RAD;
45		constexpr double d100 = 100 * DEG_TO_RAD;
46		constexpr double d140 = 140 * DEG_TO_RAD;
47		constexpr double d160 = 160 * DEG_TO_RAD;
48		constexpr double d180 = 180 * DEG_TO_RAD;
49
50		constexpr double EPSLN =
51		1.e-10; /* allow a little 'slack' on zone edge positions */
52
53		} // namespace pj_imoll_ns
54
55		/*
56		Assign an integer index representing each of the 6
57		sub-projection zones based on latitude (phi) and
58		longitude (lam) ranges.
59		*/
60
61	0	static PJ_XY imoll_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
62	0	using namespace pj_imoll_ns;
63
64	0	PJ_XY xy;
65	0	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(P->opaque);
66	0	int z;
67
68	0	if (lp.phi >= 0) { /* 1\|2 */
69	0	z = (lp.lam <= -d40 ? 1 : 2);
70	0	} else { /* 3\|4\|5\|6 */
71	0	if (lp.lam <= -d100)
72	0	z = 3; /* 3 */
73	0	else if (lp.lam <= -d20)
74	0	z = 4; /* 4 */
75	0	else if (lp.lam <= d80)
76	0	z = 5; /* 5 */
77	0	else
78	0	z = 6; /* 6 */
79	0	}
80
81	0	lp.lam -= Q->pj[z - 1]->lam0;
82	0	xy = Q->pj[z - 1]->fwd(lp, Q->pj[z - 1]);
83	0	xy.x += Q->pj[z - 1]->x0;
84	0	xy.y += Q->pj[z - 1]->y0;
85
86	0	return xy;
87	0	}
88
89	0	static PJ_LP imoll_s_inverse(PJ_XY xy, PJ P) { / Spheroidal, inverse */
90	0	using namespace pj_imoll_ns;
91
92	0	PJ_LP lp = {0.0, 0.0};
93	0	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(P->opaque);
94	0	const double y90 = sqrt(2.0); /* lt=90 corresponds to y=sqrt(2) */
95
96	0	int z = 0;
97	0	if (xy.y > y90 + EPSLN \|\| xy.y < -y90 + EPSLN) /* 0 */
98	0	z = 0;
99	0	else if (xy.y >= 0)
100	0	z = (xy.x <= (Q->boundary12) ? 1 : 2); /* 1\|2 */
101	0	else {
102	0	if (xy.x <= Q->boundary34)
103	0	z = 3; /* 3 */
104	0	else if (xy.x <= Q->boundary45)
105	0	z = 4; /* 4 */
106	0	else if (xy.x <= Q->boundary56)
107	0	z = 5; /* 5 */
108	0	else
109	0	z = 6; /* 6 */
110	0	}
111
112	0	if (z) {
113	0	bool ok = false;
114
115	0	xy.x -= Q->pj[z - 1]->x0;
116	0	xy.y -= Q->pj[z - 1]->y0;
117	0	lp = Q->pj[z - 1]->inv(xy, Q->pj[z - 1]);
118	0	lp.lam += Q->pj[z - 1]->lam0;
119
120	0	switch (z) {
121	0	case 1:
122	0	ok = ((lp.lam >= -d180 - EPSLN && lp.lam <= -d40 + EPSLN) &&
123	0	(lp.phi >= 0.0 - EPSLN));
124	0	break;
125	0	case 2:
126	0	ok = ((lp.lam >= -d40 - EPSLN && lp.lam <= d180 + EPSLN) &&
127	0	(lp.phi >= 0.0 - EPSLN));
128	0	break;
129	0	case 3:
130	0	ok = ((lp.lam >= -d180 - EPSLN && lp.lam <= -d100 + EPSLN) &&
131	0	(lp.phi <= 0.0 + EPSLN));
132	0	break;
133	0	case 4:
134	0	ok = ((lp.lam >= -d100 - EPSLN && lp.lam <= -d20 + EPSLN) &&
135	0	(lp.phi <= 0.0 + EPSLN));
136	0	break;
137	0	case 5:
138	0	ok = ((lp.lam >= -d20 - EPSLN && lp.lam <= d80 + EPSLN) &&
139	0	(lp.phi <= 0.0 + EPSLN));
140	0	break;
141	0	case 6:
142	0	ok = ((lp.lam >= d80 - EPSLN && lp.lam <= d180 + EPSLN) &&
143	0	(lp.phi <= 0.0 + EPSLN));
144	0	break;
145	0	}
146	0	z = (!ok ? 0 : z); /* projectable? */
147	0	}
148
149	0	if (!z)
150	0	lp.lam = HUGE_VAL;
151	0	if (!z)
152	0	lp.phi = HUGE_VAL;
153
154	0	return lp;
155	0	}
156
157	0	static PJ pj_imoll_destructor(PJ P, int errlev) {
158	0	using namespace pj_imoll_ns;
159
160	0	int i;
161	0	if (nullptr == P)
162	0	return nullptr;
163
164	0	if (nullptr == P->opaque)
165	0	return pj_default_destructor(P, errlev);
166
167	0	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(P->opaque);
168
169	0	for (i = 0; i < 6; ++i) {
170	0	if (Q->pj[i])
171	0	Q->pj[i]->destructor(Q->pj[i], errlev);
172	0	}
173
174	0	return pj_default_destructor(P, errlev);
175	0	}
176
177		/*
178		Zones:
179
180		-180 -40 180
181		+--------------+-------------------------+
182		\|1 \|2 \|
183		\| \| \|
184		\| \| \|
185		\| \| \|
186		\| \| \|
187		0 +-------+------+-+-----------+-----------+
188		\|3 \|4 \|5 \|6 \|
189		\| \| \| \| \|
190		\| \| \| \| \|
191		\| \| \| \| \|
192		\| \| \| \| \|
193		+-------+--------+-----------+-----------+
194		-180 -100 -20 80 180
195		*/
196
197		static bool setup_zone(PJ P, struct pj_imoll_ns::pj_imoll_data Q, int n,
198		PJ (proj_ptr)(PJ *), double x_0, double y_0,
199	0	double lon_0) {
200	0	if (!(Q->pj[n - 1] = proj_ptr(nullptr)))
201	0	return false;
202	0	if (!(Q->pj[n - 1] = proj_ptr(Q->pj[n - 1])))
203	0	return false;
204	0	Q->pj[n - 1]->ctx = P->ctx;
205	0	Q->pj[n - 1]->x0 = x_0;
206	0	Q->pj[n - 1]->y0 = y_0;
207	0	Q->pj[n - 1]->lam0 = lon_0;
208	0	return true;
209	0	}
210
211		static double compute_zone_offset(struct pj_imoll_ns::pj_imoll_data *Q,
212		int zone1, int zone2, double lam, double phi1,
213	0	double phi2) {
214	0	PJ_LP lp1, lp2;
215	0	PJ_XY xy1, xy2;
216
217	0	lp1.lam = lam - (Q->pj[zone1 - 1]->lam0);
218	0	lp1.phi = phi1;
219	0	lp2.lam = lam - (Q->pj[zone2 - 1]->lam0);
220	0	lp2.phi = phi2;
221	0	xy1 = Q->pj[zone1 - 1]->fwd(lp1, Q->pj[zone1 - 1]);
222	0	xy2 = Q->pj[zone2 - 1]->fwd(lp2, Q->pj[zone2 - 1]);
223	0	return (xy2.x + Q->pj[zone2 - 1]->x0) - (xy1.x + Q->pj[zone1 - 1]->x0);
224	0	}
225
226	0	static double compute_zone_x_boundary(PJ *P, double lam, double phi) {
227	0	PJ_LP lp1, lp2;
228	0	PJ_XY xy1, xy2;
229
230	0	lp1.lam = lam - pj_imoll_ns::EPSLN;
231	0	lp1.phi = phi;
232	0	lp2.lam = lam + pj_imoll_ns::EPSLN;
233	0	lp2.phi = phi;
234	0	xy1 = imoll_s_forward(lp1, P);
235	0	xy2 = imoll_s_forward(lp2, P);
236	0	return (xy1.x + xy2.x) / 2.;
237	0	}
238
239	0	PJ *PJ_PROJECTION(imoll) {
240	0	using namespace pj_imoll_ns;
241
242	0	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(
243	0	calloc(1, sizeof(struct pj_imoll_data)));
244	0	if (nullptr == Q)
245	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
246	0	P->opaque = Q;
247
248		/* Setup zones */
249	0	if (!setup_zone(P, Q, 1, pj_moll, -d100, 0, -d100) \|\|
250	0	!setup_zone(P, Q, 2, pj_moll, d30, 0, d30) \|\|
251	0	!setup_zone(P, Q, 3, pj_moll, -d160, 0, -d160) \|\|
252	0	!setup_zone(P, Q, 4, pj_moll, -d60, 0, -d60) \|\|
253	0	!setup_zone(P, Q, 5, pj_moll, d20, 0, d20) \|\|
254	0	!setup_zone(P, Q, 6, pj_moll, d140, 0, d140)) {
255	0	return pj_imoll_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
256	0	}
257
258		/* Adjust zones */
259
260		/* Match 3 (south) to 1 (north) */
261	0	Q->pj[2]->x0 +=
262	0	compute_zone_offset(Q, 3, 1, -d160, 0.0 - EPSLN, 0.0 + EPSLN);
263
264		/* Match 2 (north-east) to 1 (north-west) */
265	0	Q->pj[1]->x0 +=
266	0	compute_zone_offset(Q, 2, 1, -d40, 0.0 + EPSLN, 0.0 + EPSLN);
267
268		/* Match 4 (south) to 1 (north) */
269	0	Q->pj[3]->x0 +=
270	0	compute_zone_offset(Q, 4, 1, -d100, 0.0 - EPSLN, 0.0 + EPSLN);
271
272		/* Match 5 (south) to 2 (north) */
273	0	Q->pj[4]->x0 +=
274	0	compute_zone_offset(Q, 5, 2, -d20, 0.0 - EPSLN, 0.0 + EPSLN);
275
276		/* Match 6 (south) to 2 (north) */
277	0	Q->pj[5]->x0 += compute_zone_offset(Q, 6, 2, d80, 0.0 - EPSLN, 0.0 + EPSLN);
278
279		/*
280		The most straightforward way of computing the x locations of the "seams"
281		in the interrupted projection is to compute the forward transform at the
282		seams and record these values.
283		*/
284	0	Q->boundary12 = compute_zone_x_boundary(P, -d40, 0.0 + EPSLN);
285	0	Q->boundary34 = compute_zone_x_boundary(P, -d100, 0.0 - EPSLN);
286	0	Q->boundary45 = compute_zone_x_boundary(P, -d20, 0.0 - EPSLN);
287	0	Q->boundary56 = compute_zone_x_boundary(P, d80, 0.0 - EPSLN);
288
289	0	P->inv = imoll_s_inverse;
290	0	P->fwd = imoll_s_forward;
291	0	P->destructor = pj_imoll_destructor;
292	0	P->es = 0.;
293
294	0	return P;
295	0	}