1 /******************************************************************************
2 * Project: PROJ.4
3 * Purpose: Forward operation invocation
4 * Author: Thomas Knudsen, thokn@sdfe.dk, 2018-01-02
5 * Based on material from Gerald Evenden (original pj_fwd)
6 * and Piyush Agram (original pj_fwd3d)
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
31 #include <errno.h>
32 #include <math.h>
33
34 #include "proj_internal.h"
35 #include <math.h>
36
37 #define INPUT_UNITS P->left
38 #define OUTPUT_UNITS P->right
39
40
fwd_prepare(PJ * P,PJ_COORD coo)41 static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
42 if (HUGE_VAL==coo.v[0] || HUGE_VAL==coo.v[1] || HUGE_VAL==coo.v[2])
43 return proj_coord_error ();
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 /* Check validity of angular input coordinates */
50 if (INPUT_UNITS==PJ_IO_UNITS_RADIANS) {
51 double t;
52
53 /* check for latitude or longitude over-range */
54 t = (coo.lp.phi < 0 ? -coo.lp.phi : coo.lp.phi) - M_HALFPI;
55 if (t > PJ_EPS_LAT || coo.lp.lam > 10 || coo.lp.lam < -10) {
56 proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
57 return proj_coord_error ();
58 }
59
60 /* Clamp latitude to -90..90 degree range */
61 if (coo.lp.phi > M_HALFPI)
62 coo.lp.phi = M_HALFPI;
63 if (coo.lp.phi < -M_HALFPI)
64 coo.lp.phi = -M_HALFPI;
65
66 /* If input latitude is geocentrical, convert to geographical */
67 if (P->geoc)
68 coo = pj_geocentric_latitude (P, PJ_INV, coo);
69
70 /* Ensure longitude is in the -pi:pi range */
71 if (0==P->over)
72 coo.lp.lam = adjlon(coo.lp.lam);
73
74 if (P->hgridshift)
75 coo = proj_trans (P->hgridshift, PJ_INV, coo);
76 else if (P->helmert || (P->cart_wgs84 != nullptr && P->cart != nullptr)) {
77 coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
78 if( P->helmert )
79 coo = proj_trans (P->helmert, PJ_INV, coo); /* Step into local frame */
80 coo = proj_trans (P->cart, PJ_INV, coo); /* Go back to angular using local ellps */
81 }
82 if (coo.lp.lam==HUGE_VAL)
83 return coo;
84 if (P->vgridshift)
85 coo = proj_trans (P->vgridshift, PJ_FWD, coo); /* Go orthometric from geometric */
86
87 /* Distance from central meridian, taking system zero meridian into account */
88 coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
89
90 /* Ensure longitude is in the -pi:pi range */
91 if (0==P->over)
92 coo.lp.lam = adjlon(coo.lp.lam);
93
94 return coo;
95 }
96
97
98 /* We do not support gridshifts on cartesian input */
99 if (INPUT_UNITS==PJ_IO_UNITS_CARTESIAN && P->helmert)
100 return proj_trans (P->helmert, PJ_INV, coo);
101 return coo;
102 }
103
104
fwd_finalize(PJ * P,PJ_COORD coo)105 static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) {
106
107 switch (OUTPUT_UNITS) {
108
109 /* Handle false eastings/northings and non-metric linear units */
110 case PJ_IO_UNITS_CARTESIAN:
111
112 if (P->is_geocent) {
113 coo = proj_trans (P->cart, PJ_FWD, coo);
114 }
115 coo.xyz.x *= P->fr_meter;
116 coo.xyz.y *= P->fr_meter;
117 coo.xyz.z *= P->fr_meter;
118
119 break;
120
121 /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
122 case PJ_IO_UNITS_CLASSIC:
123 coo.xy.x *= P->a;
124 coo.xy.y *= P->a;
125
126 /* Falls through */ /* (<-- GCC warning silencer) */
127 /* to continue processing in common with PJ_IO_UNITS_PROJECTED */
128 case PJ_IO_UNITS_PROJECTED:
129 coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
130 coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
131 coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
132 break;
133
134 case PJ_IO_UNITS_WHATEVER:
135 break;
136
137 case PJ_IO_UNITS_DEGREES:
138 break;
139
140 case PJ_IO_UNITS_RADIANS:
141 coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);
142
143 if( P->is_long_wrap_set ) {
144 if( coo.lpz.lam != HUGE_VAL ) {
145 coo.lpz.lam = P->long_wrap_center +
146 adjlon(coo.lpz.lam - P->long_wrap_center);
147 }
148 }
149
150 break;
151 }
152
153 if (P->axisswap)
154 coo = proj_trans (P->axisswap, PJ_FWD, coo);
155
156 return coo;
157 }
158
159
error_or_coord(PJ * P,PJ_COORD coord,int last_errno)160 static PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
161 if (proj_errno(P))
162 return proj_coord_error();
163
164 proj_errno_restore(P, last_errno);
165 return coord;
166 }
167
168
pj_fwd(PJ_LP lp,PJ * P)169 PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
170 int last_errno;
171 PJ_COORD coo = {{0,0,0,0}};
172 coo.lp = lp;
173
174 last_errno = proj_errno_reset(P);
175
176 if (!P->skip_fwd_prepare)
177 coo = fwd_prepare (P, coo);
178 if (HUGE_VAL==coo.v[0] || HUGE_VAL==coo.v[1])
179 return proj_coord_error ().xy;
180
181 /* Do the transformation, using the lowest dimensional transformer available */
182 if (P->fwd)
183 coo.xy = P->fwd(coo.lp, P);
184 else if (P->fwd3d)
185 coo.xyz = P->fwd3d (coo.lpz, P);
186 else if (P->fwd4d)
187 coo = P->fwd4d (coo, P);
188 else {
189 proj_errno_set (P, EINVAL);
190 return proj_coord_error ().xy;
191 }
192 if (HUGE_VAL==coo.v[0])
193 return proj_coord_error ().xy;
194
195 if (!P->skip_fwd_finalize)
196 coo = fwd_finalize (P, coo);
197
198 return error_or_coord(P, coo, last_errno).xy;
199 }
200
201
202
pj_fwd3d(PJ_LPZ lpz,PJ * P)203 PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
204 int last_errno;
205 PJ_COORD coo = {{0,0,0,0}};
206 coo.lpz = lpz;
207
208 last_errno = proj_errno_reset(P);
209
210 if (!P->skip_fwd_prepare)
211 coo = fwd_prepare (P, coo);
212 if (HUGE_VAL==coo.v[0])
213 return proj_coord_error ().xyz;
214
215 /* Do the transformation, using the lowest dimensional transformer feasible */
216 if (P->fwd3d)
217 coo.xyz = P->fwd3d(coo.lpz, P);
218 else if (P->fwd4d)
219 coo = P->fwd4d (coo, P);
220 else if (P->fwd)
221 coo.xy = P->fwd (coo.lp, P);
222 else {
223 proj_errno_set (P, EINVAL);
224 return proj_coord_error ().xyz;
225 }
226 if (HUGE_VAL==coo.v[0])
227 return proj_coord_error ().xyz;
228
229 if (!P->skip_fwd_finalize)
230 coo = fwd_finalize (P, coo);
231
232 return error_or_coord(P, coo, last_errno).xyz;
233 }
234
235
236
pj_fwd4d(PJ_COORD coo,PJ * P)237 PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
238 int last_errno = proj_errno_reset(P);
239
240 if (!P->skip_fwd_prepare)
241 coo = fwd_prepare (P, coo);
242 if (HUGE_VAL==coo.v[0])
243 return proj_coord_error ();
244
245 /* Call the highest dimensional converter available */
246 if (P->fwd4d)
247 coo = P->fwd4d (coo, P);
248 else if (P->fwd3d)
249 coo.xyz = P->fwd3d (coo.lpz, P);
250 else if (P->fwd)
251 coo.xy = P->fwd (coo.lp, P);
252 else {
253 proj_errno_set (P, EINVAL);
254 return proj_coord_error ();
255 }
256 if (HUGE_VAL==coo.v[0])
257 return proj_coord_error ();
258
259 if (!P->skip_fwd_finalize)
260 coo = fwd_finalize (P, coo);
261
262 return error_or_coord(P, coo, last_errno);
263 }
264