1 #define PJ_LIB__
2
3 #include <errno.h>
4 #include <math.h>
5
6 #include "proj.h"
7 #include "proj_internal.h"
8
9 PROJ_HEAD(igh_o, "Interrupted Goode Homolosine Oceanic View") "\n\tPCyl, Sph";
10
11 /*
12 This projection is a variant of the Interrupted Goode Homolosine
13 projection that emphasizes ocean areas. The projection is a
14 compilation of 12 separate sub-projections. Sinusoidal projections
15 are found near the equator and Mollweide projections are found at
16 higher latitudes. The transition between the two occurs at 40 degrees
17 latitude and is represented by `phi_boundary`.
18
19 Each sub-projection is assigned an integer label
20 numbered 1 through 12. Most of this code contains logic to assign
21 the labels based on latitude (phi) and longitude (lam) regions.
22
23 Original Reference:
24 J. Paul Goode (1925) THE HOMOLOSINE PROJECTION: A NEW DEVICE FOR
25 PORTRAYING THE EARTH'S SURFACE ENTIRE, Annals of the Association of
26 American Geographers, 15:3, 119-125, DOI: 10.1080/00045602509356949
27 */
28
29 C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *);
30
31 /*
32 Transition from sinusoidal to Mollweide projection
33 Latitude (phi): 40deg 44' 11.8"
34 */
35
36 static const double phi_boundary = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD;
37
38 static const double d10 = 10 * DEG_TO_RAD;
39 static const double d20 = 20 * DEG_TO_RAD;
40 static const double d40 = 40 * DEG_TO_RAD;
41 static const double d50 = 50 * DEG_TO_RAD;
42 static const double d60 = 60 * DEG_TO_RAD;
43 static const double d90 = 90 * DEG_TO_RAD;
44 static const double d100 = 100 * DEG_TO_RAD;
45 static const double d110 = 110 * DEG_TO_RAD;
46 static const double d140 = 140 * DEG_TO_RAD;
47 static const double d150 = 150 * DEG_TO_RAD;
48 static const double d160 = 160 * DEG_TO_RAD;
49 static const double d130 = 130 * DEG_TO_RAD;
50 static const double d180 = 180 * DEG_TO_RAD;
51
52 static const double EPSLN = 1.e-10; /* allow a little 'slack' on zone edge positions */
53
54 namespace { // anonymous namespace
55 struct pj_opaque {
56 struct PJconsts* pj[12]; \
57 double dy0;
58 };
59 } // anonymous namespace
60
61
62 /*
63 Assign an integer index representing each of the 12
64 sub-projection zones based on latitude (phi) and
65 longitude (lam) ranges.
66 */
67
igh_o_s_forward(PJ_LP lp,PJ * P)68 static PJ_XY igh_o_s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */
69 PJ_XY xy;
70 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
71 int z;
72
73 if (lp.phi >= phi_boundary) {
74 if (lp.lam <= -d90) z = 1;
75 else if (lp.lam >= d60) z = 3;
76 else z = 2;
77 }
78 else if (lp.phi >= 0) {
79 if (lp.lam <= -d90) z = 4;
80 else if (lp.lam >= d60) z = 6;
81 else z = 5;
82 }
83 else if (lp.phi >= -phi_boundary) {
84 if (lp.lam <= -d60) z = 7;
85 else if (lp.lam >= d90) z = 9;
86 else z = 8;
87 }
88 else {
89 if (lp.lam <= -d60) z = 10;
90 else if (lp.lam >= d90) z = 12;
91 else z = 11;
92 }
93
94 lp.lam -= Q->pj[z-1]->lam0;
95 xy = Q->pj[z-1]->fwd(lp, Q->pj[z-1]);
96 xy.x += Q->pj[z-1]->x0;
97 xy.y += Q->pj[z-1]->y0;
98
99 return xy;
100 }
101
102
igh_o_s_inverse(PJ_XY xy,PJ * P)103 static PJ_LP igh_o_s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
104 PJ_LP lp = {0.0,0.0};
105 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
106 const double y90 = Q->dy0 + sqrt(2.0); /* lt=90 corresponds to y=y0+sqrt(2) */
107
108 int z = 0;
109 if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) /* 0 */
110 z = 0;
111 else if (xy.y >= phi_boundary)
112 if (xy.x <= -d90) z = 1;
113 else if (xy.x >= d60) z = 3;
114 else z = 2;
115 else if (xy.y >= 0)
116 if (xy.x <= -d90) z = 4;
117 else if (xy.x >= d60) z = 6;
118 else z = 5;
119 else if (xy.y >= -phi_boundary) {
120 if (xy.x <= -d60) z = 7;
121 else if (xy.x >= d90) z = 9;
122 else z = 8;
123 }
124 else {
125 if (xy.x <= -d60) z = 10;
126 else if (xy.x >= d90) z = 12;
127 else z = 11;
128 }
129
130 if (z) {
131 bool ok = false;
132
133 xy.x -= Q->pj[z-1]->x0;
134 xy.y -= Q->pj[z-1]->y0;
135 lp = Q->pj[z-1]->inv(xy, Q->pj[z-1]);
136 lp.lam += Q->pj[z-1]->lam0;
137
138 switch (z) {
139 /* Plot projectable ranges with exetension lobes in zones 1 & 3 */
140 case 1: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d90+EPSLN) ||
141 ((lp.lam >= d160-EPSLN && lp.lam <= d180+EPSLN) &&
142 (lp.phi >= d50-EPSLN && lp.phi <= d90+EPSLN)); break;
143 case 2: ok = (lp.lam >= -d90-EPSLN && lp.lam <= d60+EPSLN); break;
144 case 3: ok = (lp.lam >= d60-EPSLN && lp.lam <= d180+EPSLN) ||
145 ((lp.lam >= -d180-EPSLN && lp.lam <= -d160+EPSLN) &&
146 (lp.phi >= d50-EPSLN && lp.phi <= d90+EPSLN)); break;
147 case 4: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d90+EPSLN); break;
148 case 5: ok = (lp.lam >= -d90-EPSLN && lp.lam <= d60+EPSLN); break;
149 case 6: ok = (lp.lam >= d60-EPSLN && lp.lam <= d180+EPSLN); break;
150 case 7: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d60+EPSLN); break;
151 case 8: ok = (lp.lam >= -d60-EPSLN && lp.lam <= d90+EPSLN); break;
152 case 9: ok = (lp.lam >= d90-EPSLN && lp.lam <= d180+EPSLN); break;
153 case 10: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d60+EPSLN); break;
154 case 11: ok = (lp.lam >= -d60-EPSLN && lp.lam <= d90+EPSLN) ||
155 ((lp.lam >= d90-EPSLN && lp.lam <= d100+EPSLN) &&
156 (lp.phi >= -d90-EPSLN && lp.phi <= -d40+EPSLN)); break;
157 case 12: ok = (lp.lam >= d90-EPSLN && lp.lam <= d180+EPSLN); break;
158
159 }
160 z = (!ok? 0: z); /* projectable? */
161 }
162
163 if (!z) lp.lam = HUGE_VAL;
164 if (!z) lp.phi = HUGE_VAL;
165
166 return lp;
167 }
168
169
destructor(PJ * P,int errlev)170 static PJ *destructor (PJ *P, int errlev) {
171 int i;
172 if (nullptr==P)
173 return nullptr;
174
175 if (nullptr==P->opaque)
176 return pj_default_destructor (P, errlev);
177
178 struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
179
180 for (i = 0; i < 12; ++i) {
181 if (Q->pj[i])
182 Q->pj[i]->destructor(Q->pj[i], errlev);
183 }
184
185 return pj_default_destructor(P, errlev);
186 }
187
188
189
190 /*
191 Zones:
192
193 -180 -90 60 180
194 +---------+----------------+-------------+ Zones 1,2,3,10,11 & 12:
195 |1 |2 |3 | Mollweide projection
196 | | | |
197 +---------+----------------+-------------+ Zones 4,5,6,7,8 & 9:
198 |4 |5 |6 | Sinusoidal projection
199 | | | |
200 0 +---------+--+-------------+--+----------+
201 |7 |8 |9 |
202 | | | |
203 +------------+----------------+----------+
204 |10 |11 |12 |
205 | | | |
206 +------------+----------------+----------+
207 -180 -60 90 180
208 */
209
setup_zone(PJ * P,struct pj_opaque * Q,int n,PJ * (* proj_ptr)(PJ *),double x_0,double y_0,double lon_0)210 static bool setup_zone(PJ *P, struct pj_opaque *Q, int n,
211 PJ*(*proj_ptr)(PJ*), double x_0,
212 double y_0, double lon_0) {
213 if (!(Q->pj[n-1] = proj_ptr(nullptr))) return false;
214 if (!(Q->pj[n-1] = proj_ptr(Q->pj[n-1]))) return false;
215 Q->pj[n-1]->ctx = P->ctx;
216 Q->pj[n-1]->x0 = x_0;
217 Q->pj[n-1]->y0 = y_0;
218 Q->pj[n-1]->lam0 = lon_0;
219 return true;
220 }
221
PROJECTION(igh_o)222 PJ *PROJECTION(igh_o) {
223 PJ_XY xy1, xy4;
224 PJ_LP lp = { 0, phi_boundary };
225 struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
226 if (nullptr==Q)
227 return pj_default_destructor (P, ENOMEM);
228 P->opaque = Q;
229
230
231 /* sinusoidal zones */
232 if (!setup_zone(P, Q, 4, pj_sinu, -d140, 0, -d140) ||
233 !setup_zone(P, Q, 5, pj_sinu, -d10, 0, -d10) ||
234 !setup_zone(P, Q, 6, pj_sinu, d130, 0, d130) ||
235 !setup_zone(P, Q, 7, pj_sinu, -d110, 0, -d110) ||
236 !setup_zone(P, Q, 8, pj_sinu, d20, 0, d20) ||
237 !setup_zone(P, Q, 9, pj_sinu, d150, 0, d150))
238 {
239 return destructor(P, ENOMEM);
240 }
241
242
243 /* mollweide zones */
244 if (!setup_zone(P, Q, 1, pj_moll, -d140, 0, -d140)) {
245 return destructor(P, ENOMEM);
246 }
247
248 /* y0 ? */
249 xy1 = Q->pj[0]->fwd(lp, Q->pj[0]); /* zone 1 */
250 xy4 = Q->pj[3]->fwd(lp, Q->pj[3]); /* zone 4 */
251 /* y0 + xy1.y = xy4.y for lt = 40d44'11.8" */
252 Q->dy0 = xy4.y - xy1.y;
253
254 Q->pj[0]->y0 = Q->dy0;
255
256 /* mollweide zones (cont'd) */
257 if (!setup_zone(P, Q, 2, pj_moll, -d10, Q->dy0, -d10) ||
258 !setup_zone(P, Q, 3, pj_moll, d130, Q->dy0, d130) ||
259 !setup_zone(P, Q, 10, pj_moll, -d110, -Q->dy0, -d110) ||
260 !setup_zone(P, Q, 11, pj_moll, d20, -Q->dy0, d20) ||
261 !setup_zone(P, Q, 12, pj_moll, d150, -Q->dy0, d150))
262 {
263 return destructor(P, ENOMEM);
264 }
265
266 P->inv = igh_o_s_inverse;
267 P->fwd = igh_o_s_forward;
268 P->destructor = destructor;
269 P->es = 0.;
270
271 return P;
272 }
273