/src/PROJ/src/fwd.cpp

Source
/******************************************************************************
 * Project:  PROJ.4
 * Purpose:  Forward operation invocation
 * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2018-01-02
 *           Based on material from Gerald Evenden (original pj_fwd)
 *           and Piyush Agram (original pj_fwd3d)
 *
 ******************************************************************************
 * Copyright (c) 2000, Frank Warmerdam
 * Copyright (c) 2018, Thomas Knudsen / SDFE
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#include <errno.h>
#include <math.h>

#include "proj_internal.h"
#include <math.h>

#define INPUT_UNITS P->left
#define OUTPUT_UNITS P->right

static void fwd_prepare(PJ *P, PJ_COORD &coo) {

    /* Check validity of angular input coordinates */
    if (INPUT_UNITS == PJ_IO_UNITS_RADIANS) {

        /* check for latitude or longitude over-range */
        if (std::fabs(coo.lp.phi) > M_HALFPI) {
            if (HUGE_VAL == coo.lp.lam || HUGE_VAL == coo.lp.phi) {
                coo = proj_coord_error();
                return;
            }

            if (coo.lp.phi > 0) {
                if (coo.lp.phi - M_HALFPI > PJ_EPS_LAT) {
                    proj_log_error(P, _("Invalid latitude"));
                    proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
                    coo = proj_coord_error();
                    return;
                }
                coo.lp.phi = M_HALFPI;
            } else {
                if (coo.lp.phi - M_HALFPI < -PJ_EPS_LAT) {
                    proj_log_error(P, _("Invalid latitude"));
                    proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
                    coo = proj_coord_error();
                    return;
                }
                coo.lp.phi = -M_HALFPI;
            }
        }

        // Longitude check
        if (std::fabs(coo.lp.lam) > M_PI) {
            if (std::fabs(coo.lp.lam) > 10) {
                proj_log_error(P, _("Invalid longitude"));
                proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
                coo = proj_coord_error();
                return;
            }

            /* Ensure longitude is in the -pi:pi range */
            if (0 == P->over)
                coo.lp.lam = adjlon(coo.lp.lam);
        }

        if (HUGE_VAL == coo.v[2]) {
            coo = proj_coord_error();
            return;
        }

        /* If input latitude is geocentrical, convert to geographical */
        if (P->geoc)
            coo = pj_geocentric_latitude(P, PJ_INV, coo);

        if (P->hgridshift) {
            coo = proj_trans(P->hgridshift, PJ_INV, coo);
            if (coo.lp.lam == HUGE_VAL)
                return;
        } else if (P->helmert ||
                   (P->cart_wgs84 != nullptr && P->cart != nullptr)) {
            coo = proj_trans(P->cart_wgs84, PJ_FWD,
                             coo); /* Go cartesian in WGS84 frame */
            if (P->helmert)
                coo = proj_trans(P->helmert, PJ_INV,
                                 coo); /* Step into local frame */
            coo = proj_trans(P->cart, PJ_INV,
                             coo); /* Go back to angular using local ellps */
            if (coo.lp.lam == HUGE_VAL)
                return;
        }

        if (P->vgridshift)
            coo = proj_trans(P->vgridshift, PJ_FWD,
                             coo); /* Go orthometric from geometric */

        /* Distance from central meridian, taking system zero meridian into
         * account
         */
        coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;

        /* Ensure longitude is in the -pi:pi range */
        if (0 == P->over)
            coo.lp.lam = adjlon(coo.lp.lam);

        return;
    }

    if (HUGE_VAL == coo.v[0] || HUGE_VAL == coo.v[1] || HUGE_VAL == coo.v[2]) {
        coo = proj_coord_error();
        return;
    }

    /* We do not support gridshifts on cartesian input */
    if (INPUT_UNITS == PJ_IO_UNITS_CARTESIAN && P->helmert)
        coo = proj_trans(P->helmert, PJ_INV, coo);
    return;
}

static void fwd_finalize(PJ *P, PJ_COORD &coo) {

    switch (OUTPUT_UNITS) {

    /* Handle false eastings/northings and non-metric linear units */
    case PJ_IO_UNITS_CARTESIAN:

        if (P->is_geocent) {
            coo = proj_trans(P->cart, PJ_FWD, coo);
        }
        coo.xyz.x *= P->fr_meter;
        coo.xyz.y *= P->fr_meter;
        coo.xyz.z *= P->fr_meter;

        break;

    /* Classic proj.4 functions return plane coordinates in units of the
     * semimajor axis */
    case PJ_IO_UNITS_CLASSIC:
        coo.xy.x *= P->a;
        coo.xy.y *= P->a;
        PROJ_FALLTHROUGH;

    /* to continue processing in common with PJ_IO_UNITS_PROJECTED */
    case PJ_IO_UNITS_PROJECTED:
        coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
        coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
        coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
        break;

    case PJ_IO_UNITS_WHATEVER:
        break;

    case PJ_IO_UNITS_DEGREES:
        break;

    case PJ_IO_UNITS_RADIANS:
        coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);

        if (P->is_long_wrap_set) {
            if (coo.lpz.lam != HUGE_VAL) {
                coo.lpz.lam = P->long_wrap_center +
                              adjlon(coo.lpz.lam - P->long_wrap_center);
            }
        }

        break;
    }

    if (P->axisswap)
        coo = proj_trans(P->axisswap, PJ_FWD, coo);
}

static inline PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
    if (P->ctx->last_errno)
        return proj_coord_error();

    P->ctx->last_errno = last_errno;

    return coord;
}

PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
    PJ_COORD coo = {{0, 0, 0, 0}};
    coo.lp = lp;

    const int last_errno = P->ctx->last_errno;
    P->ctx->last_errno = 0;

    if (!P->skip_fwd_prepare)
        fwd_prepare(P, coo);
    if (HUGE_VAL == coo.v[0] || HUGE_VAL == coo.v[1])
        return proj_coord_error().xy;

    /* Do the transformation, using the lowest dimensional transformer available
     */
    if (P->fwd) {
        const auto xy = P->fwd(coo.lp, P);
        coo.xy = xy;
    } else if (P->fwd3d) {
        const auto xyz = P->fwd3d(coo.lpz, P);
        coo.xyz = xyz;
    } else if (P->fwd4d)
        P->fwd4d(coo, P);
    else {
        proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
        return proj_coord_error().xy;
    }
    if (HUGE_VAL == coo.v[0])
        return proj_coord_error().xy;

    if (!P->skip_fwd_finalize)
        fwd_finalize(P, coo);

    return error_or_coord(P, coo, last_errno).xy;
}

PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
    PJ_COORD coo = {{0, 0, 0, 0}};
    coo.lpz = lpz;

    const int last_errno = P->ctx->last_errno;
    P->ctx->last_errno = 0;

    if (!P->skip_fwd_prepare)
        fwd_prepare(P, coo);
    if (HUGE_VAL == coo.v[0])
        return proj_coord_error().xyz;

    /* Do the transformation, using the lowest dimensional transformer feasible
     */
    if (P->fwd3d) {
        const auto xyz = P->fwd3d(coo.lpz, P);
        coo.xyz = xyz;
    } else if (P->fwd4d)
        P->fwd4d(coo, P);
    else if (P->fwd) {
        const auto xy = P->fwd(coo.lp, P);
        coo.xy = xy;
    } else {
        proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
        return proj_coord_error().xyz;
    }
    if (HUGE_VAL == coo.v[0])
        return proj_coord_error().xyz;

    if (!P->skip_fwd_finalize)
        fwd_finalize(P, coo);

    return error_or_coord(P, coo, last_errno).xyz;
}

bool pj_fwd4d(PJ_COORD &coo, PJ *P) {

    const int last_errno = P->ctx->last_errno;
    P->ctx->last_errno = 0;

    if (!P->skip_fwd_prepare)
        fwd_prepare(P, coo);
    if (HUGE_VAL == coo.v[0]) {
        coo = proj_coord_error();
        return false;
    }

    /* Call the highest dimensional converter available */
    if (P->fwd4d)
        P->fwd4d(coo, P);
    else if (P->fwd3d) {
        const auto xyz = P->fwd3d(coo.lpz, P);
        coo.xyz = xyz;
    } else if (P->fwd) {
        const auto xy = P->fwd(coo.lp, P);
        coo.xy = xy;
    } else {
        proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
        coo = proj_coord_error();
        return false;
    }
    if (HUGE_VAL == coo.v[0]) {
        coo = proj_coord_error();
        return false;
    }

    if (!P->skip_fwd_finalize)
        fwd_finalize(P, coo);

    if (P->ctx->last_errno) {
        coo = proj_coord_error();
        return false;
    }

    P->ctx->last_errno = last_errno;
    return true;
}

Coverage Report

Created: 2026-09-28 07:23

Line	Count	Source
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	4.17M	#define INPUT_UNITS P->left
38	4.33M	#define OUTPUT_UNITS P->right
39
40	3.07M	static void fwd_prepare(PJ *P, PJ_COORD &coo) {
41
42		/* Check validity of angular input coordinates */
43	3.07M	if (INPUT_UNITS == PJ_IO_UNITS_RADIANS) {
44
45		/* check for latitude or longitude over-range */
46	1.97M	if (std::fabs(coo.lp.phi) > M_HALFPI) {
47	4.31k	if (HUGE_VAL == coo.lp.lam \|\| HUGE_VAL == coo.lp.phi) {
48	0	coo = proj_coord_error();
49	0	return;
50	0	}
51
52	4.31k	if (coo.lp.phi > 0) {
53	3.24k	if (coo.lp.phi - M_HALFPI > PJ_EPS_LAT) {
54	3.24k	proj_log_error(P, _("Invalid latitude"));
55	3.24k	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
56	3.24k	coo = proj_coord_error();
57	3.24k	return;
58	3.24k	}
59	0	coo.lp.phi = M_HALFPI;
60	1.06k	} else {
61	1.06k	if (coo.lp.phi - M_HALFPI < -PJ_EPS_LAT) {
62	1.06k	proj_log_error(P, _("Invalid latitude"));
63	1.06k	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
64	1.06k	coo = proj_coord_error();
65	1.06k	return;
66	1.06k	}
67	0	coo.lp.phi = -M_HALFPI;
68	0	}
69	4.31k	}
70
71		// Longitude check
72	1.97M	if (std::fabs(coo.lp.lam) > M_PI) {
73	318	if (std::fabs(coo.lp.lam) > 10) {
74	226	proj_log_error(P, _("Invalid longitude"));
75	226	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
76	226	coo = proj_coord_error();
77	226	return;
78	226	}
79
80		/* Ensure longitude is in the -pi:pi range */
81	92	if (0 == P->over)
82	92	coo.lp.lam = adjlon(coo.lp.lam);
83	92	}
84
85	1.97M	if (HUGE_VAL == coo.v[2]) {
86	0	coo = proj_coord_error();
87	0	return;
88	0	}
89
90		/* If input latitude is geocentrical, convert to geographical */
91	1.97M	if (P->geoc)
92	1.00k	coo = pj_geocentric_latitude(P, PJ_INV, coo);
93
94	1.97M	if (P->hgridshift) {
95	56.4k	coo = proj_trans(P->hgridshift, PJ_INV, coo);
96	56.4k	if (coo.lp.lam == HUGE_VAL)
97	0	return;
98	1.91M	} else if (P->helmert \|\|
99	1.63M	(P->cart_wgs84 != nullptr && P->cart != nullptr)) {
100	294k	coo = proj_trans(P->cart_wgs84, PJ_FWD,
101	294k	coo); /* Go cartesian in WGS84 frame */
102	294k	if (P->helmert)
103	281k	coo = proj_trans(P->helmert, PJ_INV,
104	281k	coo); /* Step into local frame */
105	294k	coo = proj_trans(P->cart, PJ_INV,
106	294k	coo); /* Go back to angular using local ellps */
107	294k	if (coo.lp.lam == HUGE_VAL)
108	0	return;
109	294k	}
110
111	1.97M	if (P->vgridshift)
112	12.7k	coo = proj_trans(P->vgridshift, PJ_FWD,
113	12.7k	coo); /* Go orthometric from geometric */
114
115		/* Distance from central meridian, taking system zero meridian into
116		* account
117		*/
118	1.97M	coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
119
120		/* Ensure longitude is in the -pi:pi range */
121	1.97M	if (0 == P->over)
122	1.97M	coo.lp.lam = adjlon(coo.lp.lam);
123
124	1.97M	return;
125	1.97M	}
126
127	1.09M	if (HUGE_VAL == coo.v[0] \|\| HUGE_VAL == coo.v[1] \|\| HUGE_VAL == coo.v[2]) {
128	0	coo = proj_coord_error();
129	0	return;
130	0	}
131
132		/* We do not support gridshifts on cartesian input */
133	1.09M	if (INPUT_UNITS == PJ_IO_UNITS_CARTESIAN && P->helmert)
134	0	coo = proj_trans(P->helmert, PJ_INV, coo);
135	1.09M	return;
136	1.09M	}
137
138	4.33M	static void fwd_finalize(PJ *P, PJ_COORD &coo) {
139
140	4.33M	switch (OUTPUT_UNITS) {
141
142		/* Handle false eastings/northings and non-metric linear units */
143	37.3k	case PJ_IO_UNITS_CARTESIAN:
144
145	37.3k	if (P->is_geocent) {
146	0	coo = proj_trans(P->cart, PJ_FWD, coo);
147	0	}
148	37.3k	coo.xyz.x *= P->fr_meter;
149	37.3k	coo.xyz.y *= P->fr_meter;
150	37.3k	coo.xyz.z *= P->fr_meter;
151
152	37.3k	break;
153
154		/* Classic proj.4 functions return plane coordinates in units of the
155		* semimajor axis */
156	1.12M	case PJ_IO_UNITS_CLASSIC:
157	1.12M	coo.xy.x *= P->a;
158	1.12M	coo.xy.y *= P->a;
159	1.12M	PROJ_FALLTHROUGH;
160
161		/* to continue processing in common with PJ_IO_UNITS_PROJECTED */
162	1.21M	case PJ_IO_UNITS_PROJECTED:
163	1.21M	coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
164	1.21M	coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
165	1.21M	coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
166	1.21M	break;
167
168	632k	case PJ_IO_UNITS_WHATEVER:
169	632k	break;
170
171	334k	case PJ_IO_UNITS_DEGREES:
172	334k	break;
173
174	2.11M	case PJ_IO_UNITS_RADIANS:
175	2.11M	coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);
176
177	2.11M	if (P->is_long_wrap_set) {
178	0	if (coo.lpz.lam != HUGE_VAL) {
179	0	coo.lpz.lam = P->long_wrap_center +
180	0	adjlon(coo.lpz.lam - P->long_wrap_center);
181	0	}
182	0	}
183
184	2.11M	break;
185	4.33M	}
186
187	4.33M	if (P->axisswap)
188	16.2k	coo = proj_trans(P->axisswap, PJ_FWD, coo);
189	4.33M	}
190
191	30	static inline PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
192	30	if (P->ctx->last_errno)
193	0	return proj_coord_error();
194
195	30	P->ctx->last_errno = last_errno;
196
197	30	return coord;
198	30	}
199
200	0	PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
201	0	PJ_COORD coo = {{0, 0, 0, 0}};
202	0	coo.lp = lp;
203
204	0	const int last_errno = P->ctx->last_errno;
205	0	P->ctx->last_errno = 0;
206
207	0	if (!P->skip_fwd_prepare)
208	0	fwd_prepare(P, coo);
209	0	if (HUGE_VAL == coo.v[0] \|\| HUGE_VAL == coo.v[1])
210	0	return proj_coord_error().xy;
211
212		/* Do the transformation, using the lowest dimensional transformer available
213		*/
214	0	if (P->fwd) {
215	0	const auto xy = P->fwd(coo.lp, P);
216	0	coo.xy = xy;
217	0	} else if (P->fwd3d) {
218	0	const auto xyz = P->fwd3d(coo.lpz, P);
219	0	coo.xyz = xyz;
220	0	} else if (P->fwd4d)
221	0	P->fwd4d(coo, P);
222	0	else {
223	0	proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
224	0	return proj_coord_error().xy;
225	0	}
226	0	if (HUGE_VAL == coo.v[0])
227	0	return proj_coord_error().xy;
228
229	0	if (!P->skip_fwd_finalize)
230	0	fwd_finalize(P, coo);
231
232	0	return error_or_coord(P, coo, last_errno).xy;
233	0	}
234
235	31	PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
236	31	PJ_COORD coo = {{0, 0, 0, 0}};
237	31	coo.lpz = lpz;
238
239	31	const int last_errno = P->ctx->last_errno;
240	31	P->ctx->last_errno = 0;
241
242	31	if (!P->skip_fwd_prepare)
243	31	fwd_prepare(P, coo);
244	31	if (HUGE_VAL == coo.v[0])
245	1	return proj_coord_error().xyz;
246
247		/* Do the transformation, using the lowest dimensional transformer feasible
248		*/
249	30	if (P->fwd3d) {
250	30	const auto xyz = P->fwd3d(coo.lpz, P);
251	30	coo.xyz = xyz;
252	30	} else if (P->fwd4d)
253	0	P->fwd4d(coo, P);
254	0	else if (P->fwd) {
255	0	const auto xy = P->fwd(coo.lp, P);
256	0	coo.xy = xy;
257	0	} else {
258	0	proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
259	0	return proj_coord_error().xyz;
260	0	}
261	30	if (HUGE_VAL == coo.v[0])
262	0	return proj_coord_error().xyz;
263
264	30	if (!P->skip_fwd_finalize)
265	30	fwd_finalize(P, coo);
266
267	30	return error_or_coord(P, coo, last_errno).xyz;
268	30	}
269
270	7.72M	bool pj_fwd4d(PJ_COORD &coo, PJ *P) {
271
272	7.72M	const int last_errno = P->ctx->last_errno;
273	7.72M	P->ctx->last_errno = 0;
274
275	7.72M	if (!P->skip_fwd_prepare)
276	3.07M	fwd_prepare(P, coo);
277	7.72M	if (HUGE_VAL == coo.v[0]) {
278	4.53k	coo = proj_coord_error();
279	4.53k	return false;
280	4.53k	}
281
282		/* Call the highest dimensional converter available */
283	7.72M	if (P->fwd4d)
284	4.88M	P->fwd4d(coo, P);
285	2.83M	else if (P->fwd3d) {
286	491k	const auto xyz = P->fwd3d(coo.lpz, P);
287	491k	coo.xyz = xyz;
288	2.34M	} else if (P->fwd) {
289	1.53M	const auto xy = P->fwd(coo.lp, P);
290	1.53M	coo.xy = xy;
291	1.53M	} else {
292	804k	proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
293	804k	coo = proj_coord_error();
294	804k	return false;
295	804k	}
296	6.91M	if (HUGE_VAL == coo.v[0]) {
297	13.0k	coo = proj_coord_error();
298	13.0k	return false;
299	13.0k	}
300
301	6.90M	if (!P->skip_fwd_finalize)
302	4.33M	fwd_finalize(P, coo);
303
304	6.90M	if (P->ctx->last_errno) {
305	0	coo = proj_coord_error();
306	0	return false;
307	0	}
308
309	6.90M	P->ctx->last_errno = last_errno;
310	6.90M	return true;
311	6.90M	}