/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: 2025-10-10 06:31

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