/src/PROJ/src/fwd.cpp

Source (jump to first uncovered line)
/******************************************************************************
 * 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) {
    if (HUGE_VAL == coo.v[0] || HUGE_VAL == coo.v[1] || HUGE_VAL == coo.v[2]) {
        coo = proj_coord_error();
        return;
    }

    /* The helmert datum shift will choke unless it gets a sensible 4D
     * coordinate
     */
    if (HUGE_VAL == coo.v[2] && P->helmert)
        coo.v[2] = 0.0;
    if (HUGE_VAL == coo.v[3] && P->helmert)
        coo.v[3] = 0.0;

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

        /* check for latitude or longitude over-range */
        t = (coo.lp.phi < 0 ? -coo.lp.phi : coo.lp.phi) - M_HALFPI;
        if (t > PJ_EPS_LAT) {
            proj_log_error(P, _("Invalid latitude"));
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_INVALID_COORD);
            coo = proj_coord_error();
            return;
        }
        if (coo.lp.lam > 10 || 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;
        }

        /* Clamp latitude to -90..90 degree range */
        if (coo.lp.phi > M_HALFPI)
            coo.lp.phi = M_HALFPI;
        if (coo.lp.phi < -M_HALFPI)
            coo.lp.phi = -M_HALFPI;

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

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

        if (P->hgridshift)
            coo = proj_trans(P->hgridshift, PJ_INV, coo);
        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;
    }

    /* 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: 2025-07-18 07:18

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