/src/proj/src/projections/geos.cpp

Source (jump to first uncovered line)
/*
** libproj -- library of cartographic projections
**
** Copyright (c) 2004   Gerald I. Evenden
** Copyright (c) 2012   Martin Raspaud
**
** See also (section 4.4.3.2):
**   https://www.cgms-info.org/documents/pdf_cgms_03.pdf
**
** 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 <stddef.h>

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

namespace { // anonymous namespace
struct pj_geos_data {
    double h;
    double radius_p;
    double radius_p2;
    double radius_p_inv2;
    double radius_g;
    double radius_g_1;
    double C;
    int flip_axis;
};
} // anonymous namespace

PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";

static PJ_XY geos_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_geos_data *Q = static_cast<struct pj_geos_data *>(P->opaque);
    double Vx, Vy, Vz, tmp;

    /* Calculation of the three components of the vector from satellite to
    ** position on earth surface (long,lat).*/
    tmp = cos(lp.phi);
    Vx = cos(lp.lam) * tmp;
    Vy = sin(lp.lam) * tmp;
    Vz = sin(lp.phi);

    /* Check visibility*/

    /* Calculation based on view angles from satellite.*/
    tmp = Q->radius_g - Vx;

    if (Q->flip_axis) {
        xy.x = Q->radius_g_1 * atan(Vy / hypot(Vz, tmp));
        xy.y = Q->radius_g_1 * atan(Vz / tmp);
    } else {
        xy.x = Q->radius_g_1 * atan(Vy / tmp);
        xy.y = Q->radius_g_1 * atan(Vz / hypot(Vy, tmp));
    }

    return xy;
}

static PJ_XY geos_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_geos_data *Q = static_cast<struct pj_geos_data *>(P->opaque);
    double r, Vx, Vy, Vz, tmp;

    /* Calculation of geocentric latitude. */
    lp.phi = atan(Q->radius_p2 * tan(lp.phi));

    /* Calculation of the three components of the vector from satellite to
    ** position on earth surface (long,lat).*/
    r = (Q->radius_p) / hypot(Q->radius_p * cos(lp.phi), sin(lp.phi));
    Vx = r * cos(lp.lam) * cos(lp.phi);
    Vy = r * sin(lp.lam) * cos(lp.phi);
    Vz = r * sin(lp.phi);

    /* Check visibility. */
    if (((Q->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * Q->radius_p_inv2) < 0.) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return xy;
    }

    /* Calculation based on view angles from satellite. */
    tmp = Q->radius_g - Vx;

    if (Q->flip_axis) {
        xy.x = Q->radius_g_1 * atan(Vy / hypot(Vz, tmp));
        xy.y = Q->radius_g_1 * atan(Vz / tmp);
    } else {
        xy.x = Q->radius_g_1 * atan(Vy / tmp);
        xy.y = Q->radius_g_1 * atan(Vz / hypot(Vy, tmp));
    }

    return xy;
}

static PJ_LP geos_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_geos_data *Q = static_cast<struct pj_geos_data *>(P->opaque);
    double Vx, Vy, Vz, a, b, k;

    /* Setting three components of vector from satellite to position.*/
    Vx = -1.0;
    if (Q->flip_axis) {
        Vz = tan(xy.y / Q->radius_g_1);
        Vy = tan(xy.x / Q->radius_g_1) * sqrt(1.0 + Vz * Vz);
    } else {
        Vy = tan(xy.x / Q->radius_g_1);
        Vz = tan(xy.y / Q->radius_g_1) * sqrt(1.0 + Vy * Vy);
    }

    /* Calculation of terms in cubic equation and determinant.*/
    a = Vy * Vy + Vz * Vz + Vx * Vx;
    b = 2 * Q->radius_g * Vx;
    const double det = (b * b) - 4 * a * Q->C;
    if (det < 0.) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return lp;
    }

    /* Calculation of three components of vector from satellite to position.*/
    k = (-b - sqrt(det)) / (2 * a);
    Vx = Q->radius_g + k * Vx;
    Vy *= k;
    Vz *= k;

    /* Calculation of longitude and latitude.*/
    lp.lam = atan2(Vy, Vx);
    // Initial formula:
    // lp.phi = atan(Vz * cos(lp.lam) / Vx);
    // Given that cos(atan2(y,x) = x / hypot(x,y)
    lp.phi = atan(Vz / hypot(Vx, Vy));

    return lp;
}

static PJ_LP geos_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_geos_data *Q = static_cast<struct pj_geos_data *>(P->opaque);
    double Vx, Vy, Vz, a, b, k;

    /* Setting three components of vector from satellite to position.*/
    Vx = -1.0;

    if (Q->flip_axis) {
        Vz = tan(xy.y / Q->radius_g_1);
        Vy = tan(xy.x / Q->radius_g_1) * hypot(1.0, Vz);
    } else {
        Vy = tan(xy.x / Q->radius_g_1);
        Vz = tan(xy.y / Q->radius_g_1) * hypot(1.0, Vy);
    }

    /* Calculation of terms in cubic equation and determinant.*/
    a = Vz / Q->radius_p;
    a = Vy * Vy + a * a + Vx * Vx;
    b = 2 * Q->radius_g * Vx;
    const double det = (b * b) - 4 * a * Q->C;
    if (det < 0.) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return lp;
    }

    /* Calculation of three components of vector from satellite to position.*/
    k = (-b - sqrt(det)) / (2. * a);
    Vx = Q->radius_g + k * Vx;
    Vy *= k;
    Vz *= k;

    /* Calculation of longitude and latitude.*/
    lp.lam = atan2(Vy, Vx);
    // Initial formula:
    // lp.phi = atan(Q->radius_p_inv2 * Vz * cos(lp.lam) / Vx);
    // Given that cos(atan2(y,x) = x / hypot(x,y)
    lp.phi = atan(Q->radius_p_inv2 * Vz / hypot(Vx, Vy));

    return lp;
}

PJ *PJ_PROJECTION(geos) {
    char *sweep_axis;
    struct pj_geos_data *Q = static_cast<struct pj_geos_data *>(
        calloc(1, sizeof(struct pj_geos_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;

    Q->h = pj_param(P->ctx, P->params, "dh").f;

    sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
    if (sweep_axis == nullptr)
        Q->flip_axis = 0;
    else {
        if ((sweep_axis[0] != 'x' && sweep_axis[0] != 'y') ||
            sweep_axis[1] != '\0') {
            proj_log_error(
                P, _("Invalid value for sweep: it should be equal to x or y."));
            return pj_default_destructor(P,
                                         PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }

        if (sweep_axis[0] == 'x')
            Q->flip_axis = 1;
        else
            Q->flip_axis = 0;
    }

    Q->radius_g_1 = Q->h / P->a;
    if (Q->radius_g_1 <= 0 || Q->radius_g_1 > 1e10) {
        proj_log_error(P, _("Invalid value for h."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    Q->radius_g = 1. + Q->radius_g_1;
    Q->C = Q->radius_g * Q->radius_g - 1.0;
    if (P->es != 0.0) {
        Q->radius_p = sqrt(P->one_es);
        Q->radius_p2 = P->one_es;
        Q->radius_p_inv2 = P->rone_es;
        P->inv = geos_e_inverse;
        P->fwd = geos_e_forward;
    } else {
        Q->radius_p = Q->radius_p2 = Q->radius_p_inv2 = 1.0;
        P->inv = geos_s_inverse;
        P->fwd = geos_s_forward;
    }

    return P;
}

Coverage Report

Created: 2025-06-22 06:59

Line	Count	Source (jump to first uncovered line)
1		/*
2		** libproj -- library of cartographic projections
3		**
4		** Copyright (c) 2004 Gerald I. Evenden
5		** Copyright (c) 2012 Martin Raspaud
6		**
7		** See also (section 4.4.3.2):
8		** https://www.cgms-info.org/documents/pdf_cgms_03.pdf
9		**
10		** Permission is hereby granted, free of charge, to any person obtaining
11		** a copy of this software and associated documentation files (the
12		** "Software"), to deal in the Software without restriction, including
13		** without limitation the rights to use, copy, modify, merge, publish,
14		** distribute, sublicense, and/or sell copies of the Software, and to
15		** permit persons to whom the Software is furnished to do so, subject to
16		** the following conditions:
17		**
18		** The above copyright notice and this permission notice shall be
19		** included in all copies or substantial portions of the Software.
20		**
21		** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
22		** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
23		** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
24		** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
25		** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
26		** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
27		** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
28		*/
29
30		#include <errno.h>
31		#include <math.h>
32		#include <stddef.h>
33
34		#include "proj.h"
35		#include "proj_internal.h"
36		#include <math.h>
37
38		namespace { // anonymous namespace
39		struct pj_geos_data {
40		double h;
41		double radius_p;
42		double radius_p2;
43		double radius_p_inv2;
44		double radius_g;
45		double radius_g_1;
46		double C;
47		int flip_axis;
48		};
49		} // anonymous namespace
50
51		PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";
52
53	0	static PJ_XY geos_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
54	0	PJ_XY xy = {0.0, 0.0};
55	0	struct pj_geos_data Q = static_cast<struct pj_geos_data >(P->opaque);
56	0	double Vx, Vy, Vz, tmp;
57
58		/* Calculation of the three components of the vector from satellite to
59		** position on earth surface (long,lat).*/
60	0	tmp = cos(lp.phi);
61	0	Vx = cos(lp.lam) * tmp;
62	0	Vy = sin(lp.lam) * tmp;
63	0	Vz = sin(lp.phi);
64
65		/* Check visibility*/
66
67		/* Calculation based on view angles from satellite.*/
68	0	tmp = Q->radius_g - Vx;
69
70	0	if (Q->flip_axis) {
71	0	xy.x = Q->radius_g_1 * atan(Vy / hypot(Vz, tmp));
72	0	xy.y = Q->radius_g_1 * atan(Vz / tmp);
73	0	} else {
74	0	xy.x = Q->radius_g_1 * atan(Vy / tmp);
75	0	xy.y = Q->radius_g_1 * atan(Vz / hypot(Vy, tmp));
76	0	}
77
78	0	return xy;
79	0	}
80
81	0	static PJ_XY geos_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
82	0	PJ_XY xy = {0.0, 0.0};
83	0	struct pj_geos_data Q = static_cast<struct pj_geos_data >(P->opaque);
84	0	double r, Vx, Vy, Vz, tmp;
85
86		/* Calculation of geocentric latitude. */
87	0	lp.phi = atan(Q->radius_p2 * tan(lp.phi));
88
89		/* Calculation of the three components of the vector from satellite to
90		** position on earth surface (long,lat).*/
91	0	r = (Q->radius_p) / hypot(Q->radius_p * cos(lp.phi), sin(lp.phi));
92	0	Vx = r * cos(lp.lam) * cos(lp.phi);
93	0	Vy = r * sin(lp.lam) * cos(lp.phi);
94	0	Vz = r * sin(lp.phi);
95
96		/* Check visibility. */
97	0	if (((Q->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * Q->radius_p_inv2) < 0.) {
98	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
99	0	return xy;
100	0	}
101
102		/* Calculation based on view angles from satellite. */
103	0	tmp = Q->radius_g - Vx;
104
105	0	if (Q->flip_axis) {
106	0	xy.x = Q->radius_g_1 * atan(Vy / hypot(Vz, tmp));
107	0	xy.y = Q->radius_g_1 * atan(Vz / tmp);
108	0	} else {
109	0	xy.x = Q->radius_g_1 * atan(Vy / tmp);
110	0	xy.y = Q->radius_g_1 * atan(Vz / hypot(Vy, tmp));
111	0	}
112
113	0	return xy;
114	0	}
115
116	0	static PJ_LP geos_s_inverse(PJ_XY xy, PJ P) { / Spheroidal, inverse */
117	0	PJ_LP lp = {0.0, 0.0};
118	0	struct pj_geos_data Q = static_cast<struct pj_geos_data >(P->opaque);
119	0	double Vx, Vy, Vz, a, b, k;
120
121		/* Setting three components of vector from satellite to position.*/
122	0	Vx = -1.0;
123	0	if (Q->flip_axis) {
124	0	Vz = tan(xy.y / Q->radius_g_1);
125	0	Vy = tan(xy.x / Q->radius_g_1) * sqrt(1.0 + Vz * Vz);
126	0	} else {
127	0	Vy = tan(xy.x / Q->radius_g_1);
128	0	Vz = tan(xy.y / Q->radius_g_1) * sqrt(1.0 + Vy * Vy);
129	0	}
130
131		/* Calculation of terms in cubic equation and determinant.*/
132	0	a = Vy * Vy + Vz * Vz + Vx * Vx;
133	0	b = 2 * Q->radius_g * Vx;
134	0	const double det = (b * b) - 4 * a * Q->C;
135	0	if (det < 0.) {
136	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
137	0	return lp;
138	0	}
139
140		/* Calculation of three components of vector from satellite to position.*/
141	0	k = (-b - sqrt(det)) / (2 * a);
142	0	Vx = Q->radius_g + k * Vx;
143	0	Vy *= k;
144	0	Vz *= k;
145
146		/* Calculation of longitude and latitude.*/
147	0	lp.lam = atan2(Vy, Vx);
148		// Initial formula:
149		// lp.phi = atan(Vz * cos(lp.lam) / Vx);
150		// Given that cos(atan2(y,x) = x / hypot(x,y)
151	0	lp.phi = atan(Vz / hypot(Vx, Vy));
152
153	0	return lp;
154	0	}
155
156	0	static PJ_LP geos_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
157	0	PJ_LP lp = {0.0, 0.0};
158	0	struct pj_geos_data Q = static_cast<struct pj_geos_data >(P->opaque);
159	0	double Vx, Vy, Vz, a, b, k;
160
161		/* Setting three components of vector from satellite to position.*/
162	0	Vx = -1.0;
163
164	0	if (Q->flip_axis) {
165	0	Vz = tan(xy.y / Q->radius_g_1);
166	0	Vy = tan(xy.x / Q->radius_g_1) * hypot(1.0, Vz);
167	0	} else {
168	0	Vy = tan(xy.x / Q->radius_g_1);
169	0	Vz = tan(xy.y / Q->radius_g_1) * hypot(1.0, Vy);
170	0	}
171
172		/* Calculation of terms in cubic equation and determinant.*/
173	0	a = Vz / Q->radius_p;
174	0	a = Vy * Vy + a * a + Vx * Vx;
175	0	b = 2 * Q->radius_g * Vx;
176	0	const double det = (b * b) - 4 * a * Q->C;
177	0	if (det < 0.) {
178	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
179	0	return lp;
180	0	}
181
182		/* Calculation of three components of vector from satellite to position.*/
183	0	k = (-b - sqrt(det)) / (2. * a);
184	0	Vx = Q->radius_g + k * Vx;
185	0	Vy *= k;
186	0	Vz *= k;
187
188		/* Calculation of longitude and latitude.*/
189	0	lp.lam = atan2(Vy, Vx);
190		// Initial formula:
191		// lp.phi = atan(Q->radius_p_inv2 * Vz * cos(lp.lam) / Vx);
192		// Given that cos(atan2(y,x) = x / hypot(x,y)
193	0	lp.phi = atan(Q->radius_p_inv2 * Vz / hypot(Vx, Vy));
194
195	0	return lp;
196	0	}
197
198	181	PJ *PJ_PROJECTION(geos) {
199	181	char *sweep_axis;
200	181	struct pj_geos_data Q = static_cast<struct pj_geos_data >(
201	181	calloc(1, sizeof(struct pj_geos_data)));
202	181	if (nullptr == Q)
203	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
204	181	P->opaque = Q;
205
206	181	Q->h = pj_param(P->ctx, P->params, "dh").f;
207
208	181	sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
209	181	if (sweep_axis == nullptr)
210	105	Q->flip_axis = 0;
211	76	else {
212	76	if ((sweep_axis[0] != 'x' && sweep_axis[0] != 'y') \|\|
213	76	sweep_axis[1] != '\0') {
214	5	proj_log_error(
215	5	P, _("Invalid value for sweep: it should be equal to x or y."));
216	5	return pj_default_destructor(P,
217	5	PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
218	5	}
219
220	71	if (sweep_axis[0] == 'x')
221	15	Q->flip_axis = 1;
222	56	else
223	56	Q->flip_axis = 0;
224	71	}
225
226	176	Q->radius_g_1 = Q->h / P->a;
227	176	if (Q->radius_g_1 <= 0 \|\| Q->radius_g_1 > 1e10) {
228	6	proj_log_error(P, _("Invalid value for h."));
229	6	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
230	6	}
231	170	Q->radius_g = 1. + Q->radius_g_1;
232	170	Q->C = Q->radius_g * Q->radius_g - 1.0;
233	170	if (P->es != 0.0) {
234	99	Q->radius_p = sqrt(P->one_es);
235	99	Q->radius_p2 = P->one_es;
236	99	Q->radius_p_inv2 = P->rone_es;
237	99	P->inv = geos_e_inverse;
238	99	P->fwd = geos_e_forward;
239	99	} else {
240	71	Q->radius_p = Q->radius_p2 = Q->radius_p_inv2 = 1.0;
241	71	P->inv = geos_s_inverse;
242	71	P->fwd = geos_s_forward;
243	71	}
244
245	170	return P;
246	176	}