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