/src/PROJ/src/projections/igh_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(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 `igh_o_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 | | constexpr double igh_o_phi_boundary = |
37 | | (40 + 44 / 60. + 11.8 / 3600.) * DEG_TO_RAD; |
38 | | |
39 | | namespace pj_igh_o_ns { |
40 | | struct pj_igh_o_data { |
41 | | struct PJconsts *pj[12]; |
42 | | double dy0; |
43 | | }; |
44 | | |
45 | | constexpr double d10 = 10 * DEG_TO_RAD; |
46 | | constexpr double d20 = 20 * DEG_TO_RAD; |
47 | | constexpr double d40 = 40 * DEG_TO_RAD; |
48 | | constexpr double d50 = 50 * DEG_TO_RAD; |
49 | | constexpr double d60 = 60 * DEG_TO_RAD; |
50 | | constexpr double d90 = 90 * DEG_TO_RAD; |
51 | | constexpr double d100 = 100 * DEG_TO_RAD; |
52 | | constexpr double d110 = 110 * DEG_TO_RAD; |
53 | | constexpr double d140 = 140 * DEG_TO_RAD; |
54 | | constexpr double d150 = 150 * DEG_TO_RAD; |
55 | | constexpr double d160 = 160 * DEG_TO_RAD; |
56 | | constexpr double d130 = 130 * DEG_TO_RAD; |
57 | | constexpr double d180 = 180 * DEG_TO_RAD; |
58 | | |
59 | | constexpr double EPSLN = |
60 | | 1.e-10; /* allow a little 'slack' on zone edge positions */ |
61 | | |
62 | | } // namespace pj_igh_o_ns |
63 | | |
64 | | /* |
65 | | Assign an integer index representing each of the 12 |
66 | | sub-projection zones based on latitude (phi) and |
67 | | longitude (lam) ranges. |
68 | | */ |
69 | | |
70 | 0 | static PJ_XY igh_o_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ |
71 | 0 | using namespace pj_igh_o_ns; |
72 | |
|
73 | 0 | PJ_XY xy; |
74 | 0 | struct pj_igh_o_data *Q = static_cast<struct pj_igh_o_data *>(P->opaque); |
75 | 0 | int z; |
76 | |
|
77 | 0 | if (lp.phi >= igh_o_phi_boundary) { |
78 | 0 | if (lp.lam <= -d90) |
79 | 0 | z = 1; |
80 | 0 | else if (lp.lam >= d60) |
81 | 0 | z = 3; |
82 | 0 | else |
83 | 0 | z = 2; |
84 | 0 | } else if (lp.phi >= 0) { |
85 | 0 | if (lp.lam <= -d90) |
86 | 0 | z = 4; |
87 | 0 | else if (lp.lam >= d60) |
88 | 0 | z = 6; |
89 | 0 | else |
90 | 0 | z = 5; |
91 | 0 | } else if (lp.phi >= -igh_o_phi_boundary) { |
92 | 0 | if (lp.lam <= -d60) |
93 | 0 | z = 7; |
94 | 0 | else if (lp.lam >= d90) |
95 | 0 | z = 9; |
96 | 0 | else |
97 | 0 | z = 8; |
98 | 0 | } else { |
99 | 0 | if (lp.lam <= -d60) |
100 | 0 | z = 10; |
101 | 0 | else if (lp.lam >= d90) |
102 | 0 | z = 12; |
103 | 0 | else |
104 | 0 | z = 11; |
105 | 0 | } |
106 | |
|
107 | 0 | lp.lam -= Q->pj[z - 1]->lam0; |
108 | 0 | xy = Q->pj[z - 1]->fwd(lp, Q->pj[z - 1]); |
109 | 0 | xy.x += Q->pj[z - 1]->x0; |
110 | 0 | xy.y += Q->pj[z - 1]->y0; |
111 | |
|
112 | 0 | return xy; |
113 | 0 | } |
114 | | |
115 | 0 | static PJ_LP igh_o_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ |
116 | 0 | using namespace pj_igh_o_ns; |
117 | |
|
118 | 0 | PJ_LP lp = {0.0, 0.0}; |
119 | 0 | struct pj_igh_o_data *Q = static_cast<struct pj_igh_o_data *>(P->opaque); |
120 | 0 | const double y90 = |
121 | 0 | Q->dy0 + sqrt(2.0); /* lt=90 corresponds to y=y0+sqrt(2) */ |
122 | |
|
123 | 0 | int z = 0; |
124 | 0 | if (xy.y > y90 + EPSLN || xy.y < -y90 + EPSLN) /* 0 */ |
125 | 0 | z = 0; |
126 | 0 | else if (xy.y >= igh_o_phi_boundary) |
127 | 0 | if (xy.x <= -d90) |
128 | 0 | z = 1; |
129 | 0 | else if (xy.x >= d60) |
130 | 0 | z = 3; |
131 | 0 | else |
132 | 0 | z = 2; |
133 | 0 | else if (xy.y >= 0) |
134 | 0 | if (xy.x <= -d90) |
135 | 0 | z = 4; |
136 | 0 | else if (xy.x >= d60) |
137 | 0 | z = 6; |
138 | 0 | else |
139 | 0 | z = 5; |
140 | 0 | else if (xy.y >= -igh_o_phi_boundary) { |
141 | 0 | if (xy.x <= -d60) |
142 | 0 | z = 7; |
143 | 0 | else if (xy.x >= d90) |
144 | 0 | z = 9; |
145 | 0 | else |
146 | 0 | z = 8; |
147 | 0 | } else { |
148 | 0 | if (xy.x <= -d60) |
149 | 0 | z = 10; |
150 | 0 | else if (xy.x >= d90) |
151 | 0 | z = 12; |
152 | 0 | else |
153 | 0 | z = 11; |
154 | 0 | } |
155 | |
|
156 | 0 | if (z) { |
157 | 0 | bool ok = false; |
158 | |
|
159 | 0 | xy.x -= Q->pj[z - 1]->x0; |
160 | 0 | xy.y -= Q->pj[z - 1]->y0; |
161 | 0 | lp = Q->pj[z - 1]->inv(xy, Q->pj[z - 1]); |
162 | 0 | lp.lam += Q->pj[z - 1]->lam0; |
163 | |
|
164 | 0 | switch (z) { |
165 | | /* Plot projectable ranges with exetension lobes in zones 1 & 3 */ |
166 | 0 | case 1: |
167 | 0 | ok = (lp.lam >= -d180 - EPSLN && lp.lam <= -d90 + EPSLN) || |
168 | 0 | ((lp.lam >= d160 - EPSLN && lp.lam <= d180 + EPSLN) && |
169 | 0 | (lp.phi >= d50 - EPSLN && lp.phi <= d90 + EPSLN)); |
170 | 0 | break; |
171 | 0 | case 2: |
172 | 0 | ok = (lp.lam >= -d90 - EPSLN && lp.lam <= d60 + EPSLN); |
173 | 0 | break; |
174 | 0 | case 3: |
175 | 0 | ok = (lp.lam >= d60 - EPSLN && lp.lam <= d180 + EPSLN) || |
176 | 0 | ((lp.lam >= -d180 - EPSLN && lp.lam <= -d160 + EPSLN) && |
177 | 0 | (lp.phi >= d50 - EPSLN && lp.phi <= d90 + EPSLN)); |
178 | 0 | break; |
179 | 0 | case 4: |
180 | 0 | ok = (lp.lam >= -d180 - EPSLN && lp.lam <= -d90 + EPSLN); |
181 | 0 | break; |
182 | 0 | case 5: |
183 | 0 | ok = (lp.lam >= -d90 - EPSLN && lp.lam <= d60 + EPSLN); |
184 | 0 | break; |
185 | 0 | case 6: |
186 | 0 | ok = (lp.lam >= d60 - EPSLN && lp.lam <= d180 + EPSLN); |
187 | 0 | break; |
188 | 0 | case 7: |
189 | 0 | ok = (lp.lam >= -d180 - EPSLN && lp.lam <= -d60 + EPSLN); |
190 | 0 | break; |
191 | 0 | case 8: |
192 | 0 | ok = (lp.lam >= -d60 - EPSLN && lp.lam <= d90 + EPSLN); |
193 | 0 | break; |
194 | 0 | case 9: |
195 | 0 | ok = (lp.lam >= d90 - EPSLN && lp.lam <= d180 + EPSLN); |
196 | 0 | break; |
197 | 0 | case 10: |
198 | 0 | ok = (lp.lam >= -d180 - EPSLN && lp.lam <= -d60 + EPSLN); |
199 | 0 | break; |
200 | 0 | case 11: |
201 | 0 | ok = (lp.lam >= -d60 - EPSLN && lp.lam <= d90 + EPSLN) || |
202 | 0 | ((lp.lam >= d90 - EPSLN && lp.lam <= d100 + EPSLN) && |
203 | 0 | (lp.phi >= -d90 - EPSLN && lp.phi <= -d40 + EPSLN)); |
204 | 0 | break; |
205 | 0 | case 12: |
206 | 0 | ok = (lp.lam >= d90 - EPSLN && lp.lam <= d180 + EPSLN); |
207 | 0 | break; |
208 | 0 | } |
209 | 0 | z = (!ok ? 0 : z); /* projectable? */ |
210 | 0 | } |
211 | | |
212 | 0 | if (!z) |
213 | 0 | lp.lam = HUGE_VAL; |
214 | 0 | if (!z) |
215 | 0 | lp.phi = HUGE_VAL; |
216 | |
|
217 | 0 | return lp; |
218 | 0 | } |
219 | | |
220 | 106 | static PJ *pj_igh_o_destructor(PJ *P, int errlev) { |
221 | 106 | using namespace pj_igh_o_ns; |
222 | | |
223 | 106 | int i; |
224 | 106 | if (nullptr == P) |
225 | 0 | return nullptr; |
226 | | |
227 | 106 | if (nullptr == P->opaque) |
228 | 0 | return pj_default_destructor(P, errlev); |
229 | | |
230 | 106 | struct pj_igh_o_data *Q = static_cast<struct pj_igh_o_data *>(P->opaque); |
231 | | |
232 | 1.37k | for (i = 0; i < 12; ++i) { |
233 | 1.27k | if (Q->pj[i]) |
234 | 1.27k | Q->pj[i]->destructor(Q->pj[i], errlev); |
235 | 1.27k | } |
236 | | |
237 | 106 | return pj_default_destructor(P, errlev); |
238 | 106 | } |
239 | | |
240 | | /* |
241 | | Zones: |
242 | | |
243 | | -180 -90 60 180 |
244 | | +---------+----------------+-------------+ Zones 1,2,3,10,11 & 12: |
245 | | |1 |2 |3 | Mollweide projection |
246 | | | | | | |
247 | | +---------+----------------+-------------+ Zones 4,5,6,7,8 & 9: |
248 | | |4 |5 |6 | Sinusoidal projection |
249 | | | | | | |
250 | | 0 +---------+--+-------------+--+----------+ |
251 | | |7 |8 |9 | |
252 | | | | | | |
253 | | +------------+----------------+----------+ |
254 | | |10 |11 |12 | |
255 | | | | | | |
256 | | +------------+----------------+----------+ |
257 | | -180 -60 90 180 |
258 | | */ |
259 | | |
260 | | static bool pj_igh_o_setup_zone(PJ *P, struct pj_igh_o_ns::pj_igh_o_data *Q, |
261 | | int n, PJ *(*proj_ptr)(PJ *), double x_0, |
262 | 1.27k | double y_0, double lon_0) { |
263 | 1.27k | if (!(Q->pj[n - 1] = proj_ptr(nullptr))) |
264 | 0 | return false; |
265 | 1.27k | if (!(Q->pj[n - 1] = proj_ptr(Q->pj[n - 1]))) |
266 | 0 | return false; |
267 | 1.27k | Q->pj[n - 1]->ctx = P->ctx; |
268 | 1.27k | Q->pj[n - 1]->x0 = x_0; |
269 | 1.27k | Q->pj[n - 1]->y0 = y_0; |
270 | 1.27k | Q->pj[n - 1]->lam0 = lon_0; |
271 | 1.27k | return true; |
272 | 1.27k | } |
273 | | |
274 | 106 | PJ *PJ_PROJECTION(igh_o) { |
275 | 106 | using namespace pj_igh_o_ns; |
276 | | |
277 | 106 | PJ_XY xy1, xy4; |
278 | 106 | PJ_LP lp = {0, igh_o_phi_boundary}; |
279 | 106 | struct pj_igh_o_data *Q = static_cast<struct pj_igh_o_data *>( |
280 | 106 | calloc(1, sizeof(struct pj_igh_o_data))); |
281 | 106 | if (nullptr == Q) |
282 | 0 | return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
283 | 106 | P->opaque = Q; |
284 | | |
285 | | /* sinusoidal zones */ |
286 | 106 | if (!pj_igh_o_setup_zone(P, Q, 4, pj_sinu, -d140, 0, -d140) || |
287 | 106 | !pj_igh_o_setup_zone(P, Q, 5, pj_sinu, -d10, 0, -d10) || |
288 | 106 | !pj_igh_o_setup_zone(P, Q, 6, pj_sinu, d130, 0, d130) || |
289 | 106 | !pj_igh_o_setup_zone(P, Q, 7, pj_sinu, -d110, 0, -d110) || |
290 | 106 | !pj_igh_o_setup_zone(P, Q, 8, pj_sinu, d20, 0, d20) || |
291 | 106 | !pj_igh_o_setup_zone(P, Q, 9, pj_sinu, d150, 0, d150)) { |
292 | 0 | return pj_igh_o_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
293 | 0 | } |
294 | | |
295 | | /* mollweide zones */ |
296 | 106 | if (!pj_igh_o_setup_zone(P, Q, 1, pj_moll, -d140, 0, -d140)) { |
297 | 0 | return pj_igh_o_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
298 | 0 | } |
299 | | |
300 | | /* y0 ? */ |
301 | 106 | xy1 = Q->pj[0]->fwd(lp, Q->pj[0]); /* zone 1 */ |
302 | 106 | xy4 = Q->pj[3]->fwd(lp, Q->pj[3]); /* zone 4 */ |
303 | | /* y0 + xy1.y = xy4.y for lt = 40d44'11.8" */ |
304 | 106 | Q->dy0 = xy4.y - xy1.y; |
305 | | |
306 | 106 | Q->pj[0]->y0 = Q->dy0; |
307 | | |
308 | | /* mollweide zones (cont'd) */ |
309 | 106 | if (!pj_igh_o_setup_zone(P, Q, 2, pj_moll, -d10, Q->dy0, -d10) || |
310 | 106 | !pj_igh_o_setup_zone(P, Q, 3, pj_moll, d130, Q->dy0, d130) || |
311 | 106 | !pj_igh_o_setup_zone(P, Q, 10, pj_moll, -d110, -Q->dy0, -d110) || |
312 | 106 | !pj_igh_o_setup_zone(P, Q, 11, pj_moll, d20, -Q->dy0, d20) || |
313 | 106 | !pj_igh_o_setup_zone(P, Q, 12, pj_moll, d150, -Q->dy0, d150)) { |
314 | 0 | return pj_igh_o_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
315 | 0 | } |
316 | | |
317 | 106 | P->inv = igh_o_s_inverse; |
318 | 106 | P->fwd = igh_o_s_forward; |
319 | 106 | P->destructor = pj_igh_o_destructor; |
320 | 106 | P->es = 0.; |
321 | | |
322 | 106 | return P; |
323 | 106 | } |