1 /******************************************************************************
2 * Project: PROJ.4
3 * Purpose: Inverse operation invocation
4 * Author: Thomas Knudsen, thokn@sdfe.dk, 2018-01-02
5 * Based on material from Gerald Evenden (original pj_inv)
6 * and Piyush Agram (original pj_inv3d)
7 *
8 ******************************************************************************
9 * Copyright (c) 2000, Frank Warmerdam
10 * Copyright (c) 2018, Thomas Knudsen / SDFE
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 *****************************************************************************/
30 #include <errno.h>
31 #include <math.h>
32
33 #include "proj_internal.h"
34 #include <math.h>
35
36 #define INPUT_UNITS P->right
37 #define OUTPUT_UNITS P->left
38
inv_prepare(PJ * P,PJ_COORD coo)39 static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
40 if (coo.v[0] == HUGE_VAL || coo.v[1] == HUGE_VAL || coo.v[2] == HUGE_VAL) {
41 proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
42 return proj_coord_error ();
43 }
44
45 /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
46 if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
47 if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
48
49 if (P->axisswap)
50 coo = proj_trans (P->axisswap, PJ_INV, coo);
51
52 /* Handle remaining possible input types */
53 switch (INPUT_UNITS) {
54 case PJ_IO_UNITS_WHATEVER:
55 return coo;
56
57 case PJ_IO_UNITS_DEGREES:
58 return coo;
59
60 /* de-scale and de-offset */
61 case PJ_IO_UNITS_CARTESIAN:
62 coo.xyz.x *= P->to_meter;
63 coo.xyz.y *= P->to_meter;
64 coo.xyz.z *= P->to_meter;
65 if (P->is_geocent) {
66 coo = proj_trans (P->cart, PJ_INV, coo);
67 }
68
69 return coo;
70
71 case PJ_IO_UNITS_PROJECTED:
72 case PJ_IO_UNITS_CLASSIC:
73 coo.xyz.x = P->to_meter * coo.xyz.x - P->x0;
74 coo.xyz.y = P->to_meter * coo.xyz.y - P->y0;
75 coo.xyz.z = P->vto_meter * coo.xyz.z - P->z0;
76 if (INPUT_UNITS==PJ_IO_UNITS_PROJECTED)
77 return coo;
78
79 /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
80 /* Multiplying by ra, rather than dividing by a because the CalCOFI projection */
81 /* stomps on a and hence (apparently) depends on this to roundtrip correctly */
82 /* (CalCOFI avoids further scaling by stomping - but a better solution is possible) */
83 coo.xyz.x *= P->ra;
84 coo.xyz.y *= P->ra;
85 return coo;
86
87 case PJ_IO_UNITS_RADIANS:
88 coo.lpz.z = P->vto_meter * coo.lpz.z - P->z0;
89 break;
90 }
91
92 /* Should not happen, so we could return pj_coord_err here */
93 return coo;
94 }
95
96
97
inv_finalize(PJ * P,PJ_COORD coo)98 static PJ_COORD inv_finalize (PJ *P, PJ_COORD coo) {
99 if (coo.xyz.x == HUGE_VAL) {
100 proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
101 return proj_coord_error ();
102 }
103
104 if (OUTPUT_UNITS==PJ_IO_UNITS_RADIANS) {
105
106 /* Distance from central meridian, taking system zero meridian into account */
107 coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
108
109 /* adjust longitude to central meridian */
110 if (0==P->over)
111 coo.lpz.lam = adjlon(coo.lpz.lam);
112
113 if (P->vgridshift)
114 coo = proj_trans (P->vgridshift, PJ_INV, coo); /* Go geometric from orthometric */
115 if (coo.lp.lam==HUGE_VAL)
116 return coo;
117 if (P->hgridshift)
118 coo = proj_trans (P->hgridshift, PJ_FWD, coo);
119 else if (P->helmert || (P->cart_wgs84 != nullptr && P->cart != nullptr)) {
120 coo = proj_trans (P->cart, PJ_FWD, coo); /* Go cartesian in local frame */
121 if( P->helmert )
122 coo = proj_trans (P->helmert, PJ_FWD, coo); /* Step into WGS84 */
123 coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
124 }
125 if (coo.lp.lam==HUGE_VAL)
126 return coo;
127
128 /* If input latitude was geocentrical, convert back to geocentrical */
129 if (P->geoc)
130 coo = pj_geocentric_latitude (P, PJ_FWD, coo);
131 }
132
133 return coo;
134 }
135
136
error_or_coord(PJ * P,PJ_COORD coord,int last_errno)137 static PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
138 if (proj_errno(P))
139 return proj_coord_error();
140
141 proj_errno_restore(P, last_errno);
142 return coord;
143 }
144
145
pj_inv(PJ_XY xy,PJ * P)146 PJ_LP pj_inv(PJ_XY xy, PJ *P) {
147 int last_errno;
148 PJ_COORD coo = {{0,0,0,0}};
149 coo.xy = xy;
150
151 last_errno = proj_errno_reset(P);
152
153 if (!P->skip_inv_prepare)
154 coo = inv_prepare (P, coo);
155 if (HUGE_VAL==coo.v[0])
156 return proj_coord_error ().lp;
157
158 /* Do the transformation, using the lowest dimensional transformer available */
159 if (P->inv)
160 coo.lp = P->inv(coo.xy, P);
161 else if (P->inv3d)
162 coo.lpz = P->inv3d (coo.xyz, P);
163 else if (P->inv4d)
164 coo = P->inv4d (coo, P);
165 else {
166 proj_errno_set (P, EINVAL);
167 return proj_coord_error ().lp;
168 }
169 if (HUGE_VAL==coo.v[0])
170 return proj_coord_error ().lp;
171
172 if (!P->skip_inv_finalize)
173 coo = inv_finalize (P, coo);
174
175 return error_or_coord(P, coo, last_errno).lp;
176 }
177
178
179
pj_inv3d(PJ_XYZ xyz,PJ * P)180 PJ_LPZ pj_inv3d (PJ_XYZ xyz, PJ *P) {
181 int last_errno;
182 PJ_COORD coo = {{0,0,0,0}};
183 coo.xyz = xyz;
184
185 last_errno = proj_errno_reset(P);
186
187 if (!P->skip_inv_prepare)
188 coo = inv_prepare (P, coo);
189 if (HUGE_VAL==coo.v[0])
190 return proj_coord_error ().lpz;
191
192 /* Do the transformation, using the lowest dimensional transformer feasible */
193 if (P->inv3d)
194 coo.lpz = P->inv3d (coo.xyz, P);
195 else if (P->inv4d)
196 coo = P->inv4d (coo, P);
197 else if (P->inv)
198 coo.lp = P->inv (coo.xy, P);
199 else {
200 proj_errno_set (P, EINVAL);
201 return proj_coord_error ().lpz;
202 }
203 if (HUGE_VAL==coo.v[0])
204 return proj_coord_error ().lpz;
205
206 if (!P->skip_inv_finalize)
207 coo = inv_finalize (P, coo);
208
209 return error_or_coord(P, coo, last_errno).lpz;
210 }
211
212
213
pj_inv4d(PJ_COORD coo,PJ * P)214 PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
215 int last_errno = proj_errno_reset(P);
216
217 if (!P->skip_inv_prepare)
218 coo = inv_prepare (P, coo);
219 if (HUGE_VAL==coo.v[0])
220 return proj_coord_error ();
221
222 /* Call the highest dimensional converter available */
223 if (P->inv4d)
224 coo = P->inv4d (coo, P);
225 else if (P->inv3d)
226 coo.lpz = P->inv3d (coo.xyz, P);
227 else if (P->inv)
228 coo.lp = P->inv (coo.xy, P);
229 else {
230 proj_errno_set (P, EINVAL);
231 return proj_coord_error ();
232 }
233 if (HUGE_VAL==coo.v[0])
234 return proj_coord_error ();
235
236 if (!P->skip_inv_finalize)
237 coo = inv_finalize (P, coo);
238
239 return error_or_coord(P, coo, last_errno);
240 }
241