/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-06-30 06:29

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.33M	#define INPUT_UNITS P->left
38	4.75M	#define OUTPUT_UNITS P->right
39
40	3.26M	static void fwd_prepare(PJ *P, PJ_COORD &coo) {
41
42		/* Check validity of angular input coordinates */
43	3.26M	if (INPUT_UNITS == PJ_IO_UNITS_RADIANS) {
44
45		/* check for latitude or longitude over-range */
46	2.19M	if (std::fabs(coo.lp.phi) > M_HALFPI) {
47	3.02k	if (HUGE_VAL == coo.lp.lam \|\| HUGE_VAL == coo.lp.phi) {
48	0	coo = proj_coord_error();
49	0	return;
50	0	}
51
52	3.02k	if (coo.lp.phi > 0) {
53	2.85k	if (coo.lp.phi - M_HALFPI > PJ_EPS_LAT) {
54	2.85k	proj_log_error(P, _("Invalid latitude"));
55	2.85k	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
56	2.85k	coo = proj_coord_error();
57	2.85k	return;
58	2.85k	}
59	0	coo.lp.phi = M_HALFPI;
60	168	} else {
61	168	if (coo.lp.phi - M_HALFPI < -PJ_EPS_LAT) {
62	168	proj_log_error(P, _("Invalid latitude"));
63	168	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
64	168	coo = proj_coord_error();
65	168	return;
66	168	}
67	0	coo.lp.phi = -M_HALFPI;
68	0	}
69	3.02k	}
70
71		// Longitude check
72	2.18M	if (std::fabs(coo.lp.lam) > M_PI) {
73	24	if (std::fabs(coo.lp.lam) > 10) {
74	1	proj_log_error(P, _("Invalid longitude"));
75	1	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
76	1	coo = proj_coord_error();
77	1	return;
78	1	}
79
80		/* Ensure longitude is in the -pi:pi range */
81	23	if (0 == P->over)
82	23	coo.lp.lam = adjlon(coo.lp.lam);
83	23	}
84
85	2.18M	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	2.18M	if (P->geoc)
92	0	coo = pj_geocentric_latitude(P, PJ_INV, coo);
93
94	2.18M	if (P->hgridshift) {
95	46.4k	coo = proj_trans(P->hgridshift, PJ_INV, coo);
96	46.4k	if (coo.lp.lam == HUGE_VAL)
97	0	return;
98	2.14M	} else if (P->helmert \|\|
99	1.81M	(P->cart_wgs84 != nullptr && P->cart != nullptr)) {
100	342k	coo = proj_trans(P->cart_wgs84, PJ_FWD,
101	342k	coo); /* Go cartesian in WGS84 frame */
102	342k	if (P->helmert)
103	333k	coo = proj_trans(P->helmert, PJ_INV,
104	333k	coo); /* Step into local frame */
105	342k	coo = proj_trans(P->cart, PJ_INV,
106	342k	coo); /* Go back to angular using local ellps */
107	342k	if (coo.lp.lam == HUGE_VAL)
108	0	return;
109	342k	}
110
111	2.18M	if (P->vgridshift)
112	7.39k	coo = proj_trans(P->vgridshift, PJ_FWD,
113	7.39k	coo); /* Go orthometric from geometric */
114
115		/* Distance from central meridian, taking system zero meridian into
116		* account
117		*/
118	2.18M	coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
119
120		/* Ensure longitude is in the -pi:pi range */
121	2.18M	if (0 == P->over)
122	2.18M	coo.lp.lam = adjlon(coo.lp.lam);
123
124	2.18M	return;
125	2.18M	}
126
127	1.06M	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.06M	if (INPUT_UNITS == PJ_IO_UNITS_CARTESIAN && P->helmert)
134	0	coo = proj_trans(P->helmert, PJ_INV, coo);
135	1.06M	return;
136	1.06M	}
137
138	4.75M	static void fwd_finalize(PJ *P, PJ_COORD &coo) {
139
140	4.75M	switch (OUTPUT_UNITS) {
141
142		/* Handle false eastings/northings and non-metric linear units */
143	55.3k	case PJ_IO_UNITS_CARTESIAN:
144
145	55.3k	if (P->is_geocent) {
146	0	coo = proj_trans(P->cart, PJ_FWD, coo);
147	0	}
148	55.3k	coo.xyz.x *= P->fr_meter;
149	55.3k	coo.xyz.y *= P->fr_meter;
150	55.3k	coo.xyz.z *= P->fr_meter;
151
152	55.3k	break;
153
154		/* Classic proj.4 functions return plane coordinates in units of the
155		* semimajor axis */
156	1.26M	case PJ_IO_UNITS_CLASSIC:
157	1.26M	coo.xy.x *= P->a;
158	1.26M	coo.xy.y *= P->a;
159	1.26M	PROJ_FALLTHROUGH;
160
161		/* to continue processing in common with PJ_IO_UNITS_PROJECTED */
162	1.31M	case PJ_IO_UNITS_PROJECTED:
163	1.31M	coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
164	1.31M	coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
165	1.31M	coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
166	1.31M	break;
167
168	708k	case PJ_IO_UNITS_WHATEVER:
169	708k	break;
170
171	367k	case PJ_IO_UNITS_DEGREES:
172	367k	break;
173
174	2.29M	case PJ_IO_UNITS_RADIANS:
175	2.29M	coo.lpz.z = P->vfr_meter * (coo.lpz.z + P->z0);
176
177	2.29M	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.29M	break;
185	4.75M	}
186
187	4.75M	if (P->axisswap)
188	15.2k	coo = proj_trans(P->axisswap, PJ_FWD, coo);
189	4.75M	}
190
191	23	static inline PJ_COORD error_or_coord(PJ *P, PJ_COORD coord, int last_errno) {
192	23	if (P->ctx->last_errno)
193	0	return proj_coord_error();
194
195	23	P->ctx->last_errno = last_errno;
196
197	23	return coord;
198	23	}
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	24	PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
236	24	PJ_COORD coo = {{0, 0, 0, 0}};
237	24	coo.lpz = lpz;
238
239	24	const int last_errno = P->ctx->last_errno;
240	24	P->ctx->last_errno = 0;
241
242	24	if (!P->skip_fwd_prepare)
243	24	fwd_prepare(P, coo);
244	24	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	23	if (P->fwd3d) {
250	23	const auto xyz = P->fwd3d(coo.lpz, P);
251	23	coo.xyz = xyz;
252	23	} 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	23	if (HUGE_VAL == coo.v[0])
262	0	return proj_coord_error().xyz;
263
264	23	if (!P->skip_fwd_finalize)
265	23	fwd_finalize(P, coo);
266
267	23	return error_or_coord(P, coo, last_errno).xyz;
268	23	}
269
270	8.30M	bool pj_fwd4d(PJ_COORD &coo, PJ *P) {
271
272	8.30M	const int last_errno = P->ctx->last_errno;
273	8.30M	P->ctx->last_errno = 0;
274
275	8.30M	if (!P->skip_fwd_prepare)
276	3.26M	fwd_prepare(P, coo);
277	8.30M	if (HUGE_VAL == coo.v[0]) {
278	3.02k	coo = proj_coord_error();
279	3.02k	return false;
280	3.02k	}
281
282		/* Call the highest dimensional converter available */
283	8.29M	if (P->fwd4d)
284	5.25M	P->fwd4d(coo, P);
285	3.04M	else if (P->fwd3d) {
286	557k	const auto xyz = P->fwd3d(coo.lpz, P);
287	557k	coo.xyz = xyz;
288	2.48M	} else if (P->fwd) {
289	1.74M	const auto xy = P->fwd(coo.lp, P);
290	1.74M	coo.xy = xy;
291	1.74M	} else {
292	739k	proj_errno_set(P, PROJ_ERR_OTHER_NO_INVERSE_OP);
293	739k	coo = proj_coord_error();
294	739k	return false;
295	739k	}
296	7.55M	if (HUGE_VAL == coo.v[0]) {
297	12.2k	coo = proj_coord_error();
298	12.2k	return false;
299	12.2k	}
300
301	7.54M	if (!P->skip_fwd_finalize)
302	4.75M	fwd_finalize(P, coo);
303
304	7.54M	if (P->ctx->last_errno) {
305	0	coo = proj_coord_error();
306	0	return false;
307	0	}
308
309	7.54M	P->ctx->last_errno = last_errno;
310	7.54M	return true;
311	7.54M	}