/src/PROJ/src/projections/imoll.cpp

Source


#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: 2026-05-16 07:15

Line	Count	Source
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	408	static PJ_XY imoll_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
62	408	using namespace pj_imoll_ns;
63
64	408	PJ_XY xy;
65	408	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(P->opaque);
66	408	int z;
67
68	408	if (lp.phi >= 0) { /* 1\|2 */
69	102	z = (lp.lam <= -d40 ? 1 : 2);
70	306	} else { /* 3\|4\|5\|6 */
71	306	if (lp.lam <= -d100)
72	51	z = 3; /* 3 */
73	255	else if (lp.lam <= -d20)
74	102	z = 4; /* 4 */
75	153	else if (lp.lam <= d80)
76	102	z = 5; /* 5 */
77	51	else
78	51	z = 6; /* 6 */
79	306	}
80
81	408	lp.lam -= Q->pj[z - 1]->lam0;
82	408	xy = Q->pj[z - 1]->fwd(lp, Q->pj[z - 1]);
83	408	xy.x += Q->pj[z - 1]->x0;
84	408	xy.y += Q->pj[z - 1]->y0;
85
86	408	return xy;
87	408	}
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	51	static PJ pj_imoll_destructor(PJ P, int errlev) {
158	51	using namespace pj_imoll_ns;
159
160	51	int i;
161	51	if (nullptr == P)
162	0	return nullptr;
163
164	51	if (nullptr == P->opaque)
165	0	return pj_default_destructor(P, errlev);
166
167	51	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(P->opaque);
168
169	357	for (i = 0; i < 6; ++i) {
170	306	if (Q->pj[i])
171	306	Q->pj[i]->destructor(Q->pj[i], errlev);
172	306	}
173
174	51	return pj_default_destructor(P, errlev);
175	51	}
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	306	double lon_0) {
200	306	if (!(Q->pj[n - 1] = proj_ptr(nullptr)))
201	0	return false;
202	306	if (!(Q->pj[n - 1] = proj_ptr(Q->pj[n - 1])))
203	0	return false;
204	306	Q->pj[n - 1]->ctx = P->ctx;
205	306	Q->pj[n - 1]->x0 = x_0;
206	306	Q->pj[n - 1]->y0 = y_0;
207	306	Q->pj[n - 1]->lam0 = lon_0;
208	306	return true;
209	306	}
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	255	double phi2) {
214	255	PJ_LP lp1, lp2;
215	255	PJ_XY xy1, xy2;
216
217	255	lp1.lam = lam - (Q->pj[zone1 - 1]->lam0);
218	255	lp1.phi = phi1;
219	255	lp2.lam = lam - (Q->pj[zone2 - 1]->lam0);
220	255	lp2.phi = phi2;
221	255	xy1 = Q->pj[zone1 - 1]->fwd(lp1, Q->pj[zone1 - 1]);
222	255	xy2 = Q->pj[zone2 - 1]->fwd(lp2, Q->pj[zone2 - 1]);
223	255	return (xy2.x + Q->pj[zone2 - 1]->x0) - (xy1.x + Q->pj[zone1 - 1]->x0);
224	255	}
225
226	204	static double compute_zone_x_boundary(PJ *P, double lam, double phi) {
227	204	PJ_LP lp1, lp2;
228	204	PJ_XY xy1, xy2;
229
230	204	lp1.lam = lam - pj_imoll_ns::EPSLN;
231	204	lp1.phi = phi;
232	204	lp2.lam = lam + pj_imoll_ns::EPSLN;
233	204	lp2.phi = phi;
234	204	xy1 = imoll_s_forward(lp1, P);
235	204	xy2 = imoll_s_forward(lp2, P);
236	204	return (xy1.x + xy2.x) / 2.;
237	204	}
238
239	51	PJ *PJ_PROJECTION(imoll) {
240	51	using namespace pj_imoll_ns;
241
242	51	struct pj_imoll_data Q = static_cast<struct pj_imoll_data >(
243	51	calloc(1, sizeof(struct pj_imoll_data)));
244	51	if (nullptr == Q)
245	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
246	51	P->opaque = Q;
247
248		/* Setup zones */
249	51	if (!setup_zone(P, Q, 1, pj_moll, -d100, 0, -d100) \|\|
250	51	!setup_zone(P, Q, 2, pj_moll, d30, 0, d30) \|\|
251	51	!setup_zone(P, Q, 3, pj_moll, -d160, 0, -d160) \|\|
252	51	!setup_zone(P, Q, 4, pj_moll, -d60, 0, -d60) \|\|
253	51	!setup_zone(P, Q, 5, pj_moll, d20, 0, d20) \|\|
254	51	!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	51	Q->pj[2]->x0 +=
262	51	compute_zone_offset(Q, 3, 1, -d160, 0.0 - EPSLN, 0.0 + EPSLN);
263
264		/* Match 2 (north-east) to 1 (north-west) */
265	51	Q->pj[1]->x0 +=
266	51	compute_zone_offset(Q, 2, 1, -d40, 0.0 + EPSLN, 0.0 + EPSLN);
267
268		/* Match 4 (south) to 1 (north) */
269	51	Q->pj[3]->x0 +=
270	51	compute_zone_offset(Q, 4, 1, -d100, 0.0 - EPSLN, 0.0 + EPSLN);
271
272		/* Match 5 (south) to 2 (north) */
273	51	Q->pj[4]->x0 +=
274	51	compute_zone_offset(Q, 5, 2, -d20, 0.0 - EPSLN, 0.0 + EPSLN);
275
276		/* Match 6 (south) to 2 (north) */
277	51	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	51	Q->boundary12 = compute_zone_x_boundary(P, -d40, 0.0 + EPSLN);
285	51	Q->boundary34 = compute_zone_x_boundary(P, -d100, 0.0 - EPSLN);
286	51	Q->boundary45 = compute_zone_x_boundary(P, -d20, 0.0 - EPSLN);
287	51	Q->boundary56 = compute_zone_x_boundary(P, d80, 0.0 - EPSLN);
288
289	51	P->inv = imoll_s_inverse;
290	51	P->fwd = imoll_s_forward;
291	51	P->destructor = pj_imoll_destructor;
292	51	P->es = 0.;
293
294	51	return P;
295	51	}