Coverage Report

Created: 2026-01-09 07:03

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