/src/proj/src/projections/aea.cpp

Source (jump to first uncovered line)
/******************************************************************************
 * Project:  PROJ.4
 * Purpose:  Implementation of the aea (Albers Equal Area) projection.
 *           and the leac (Lambert Equal Area Conic) projection
 * Author:   Gerald Evenden (1995)
 *           Thomas Knudsen (2016) - revise/add regression tests
 *
 ******************************************************************************
 * Copyright (c) 1995, Gerald Evenden
 *
 * 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 "proj.h"
#include "proj_internal.h"
#include <errno.h>
#include <math.h>

#define EPS10 1.e-10
#define TOL7 1.e-7

PROJ_HEAD(aea, "Albers Equal Area") "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
PROJ_HEAD(leac, "Lambert Equal Area Conic")
"\n\tConic, Sph&Ell\n\tlat_1= south";

namespace { // anonymous namespace
struct pj_aea {
    double ec;
    double n;
    double c;
    double dd;
    double n2;
    double rho0;
    double rho;
    double phi1;
    double phi2;
    int ellips;
    double *apa;
    double qp;
};
} // anonymous namespace

static PJ *pj_aea_destructor(PJ *P, int errlev) { /* Destructor */
    if (nullptr == P)
        return nullptr;

    if (nullptr == P->opaque)
        return pj_default_destructor(P, errlev);

    free(static_cast<struct pj_aea *>(P->opaque)->apa);

    return pj_default_destructor(P, errlev);
}

static PJ_XY aea_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoid/spheroid, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_aea *Q = static_cast<struct pj_aea *>(P->opaque);
    Q->rho = Q->c - (Q->ellips ? Q->n * pj_authalic_lat_q(sin(lp.phi), P)
                               : Q->n2 * sin(lp.phi));
    if (Q->rho < 0.) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return xy;
    }
    Q->rho = Q->dd * sqrt(Q->rho);
    lp.lam *= Q->n;
    xy.x = Q->rho * sin(lp.lam);
    xy.y = Q->rho0 - Q->rho * cos(lp.lam);
    return xy;
}

static PJ_LP aea_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoid/spheroid, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_aea *Q = static_cast<struct pj_aea *>(P->opaque);
    xy.y = Q->rho0 - xy.y;
    Q->rho = hypot(xy.x, xy.y);
    if (Q->rho != 0.0) {
        if (Q->n < 0.) {
            Q->rho = -Q->rho;
            xy.x = -xy.x;
            xy.y = -xy.y;
        }
        lp.phi = Q->rho / Q->dd;
        if (Q->ellips) {
            const double qs = (Q->c - lp.phi * lp.phi) / Q->n;
            if (fabs(Q->ec - fabs(qs)) > TOL7) {
                if (fabs(qs) > 2) {
                    proj_errno_set(
                        P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                    return lp;
                }
                lp.phi = pj_authalic_lat_inverse(asin(qs / Q->qp), Q->apa,
                                                 P, Q->qp);
                if (lp.phi == HUGE_VAL) {
                    proj_errno_set(
                        P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
                    return lp;
                }
            } else
                lp.phi = qs < 0. ? -M_HALFPI : M_HALFPI;
        } else {
            const double qs_div_2 = (Q->c - lp.phi * lp.phi) / Q->n2;
            if (fabs(qs_div_2) <= 1.)
                lp.phi = asin(qs_div_2);
            else
                lp.phi = qs_div_2 < 0. ? -M_HALFPI : M_HALFPI;
        }
        lp.lam = atan2(xy.x, xy.y) / Q->n;
    } else {
        lp.lam = 0.;
        lp.phi = Q->n > 0. ? M_HALFPI : -M_HALFPI;
    }
    return lp;
}

static PJ *setup(PJ *P) {
    struct pj_aea *Q = static_cast<struct pj_aea *>(P->opaque);

    P->inv = aea_e_inverse;
    P->fwd = aea_e_forward;

    if (fabs(Q->phi1) > M_HALFPI) {
        proj_log_error(P,
                       _("Invalid value for lat_1: |lat_1| should be <= 90°"));
        return pj_aea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    if (fabs(Q->phi2) > M_HALFPI) {
        proj_log_error(P,
                       _("Invalid value for lat_2: |lat_2| should be <= 90°"));
        return pj_aea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    if (fabs(Q->phi1 + Q->phi2) < EPS10) {
        proj_log_error(P, _("Invalid value for lat_1 and lat_2: |lat_1 + "
                            "lat_2| should be > 0"));
        return pj_aea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    double sinphi = sin(Q->phi1);
    Q->n = sinphi;
    double cosphi = cos(Q->phi1);
    const int secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
    Q->ellips = (P->es > 0.);
    if (Q->ellips) {
        double ml1, m1;

        Q->apa = pj_authalic_lat_compute_coeffs(P->n);
        if (Q->apa == nullptr)
            return pj_aea_destructor(P, 0);
        Q->qp = pj_authalic_lat_q(1.0, P);
        m1 = pj_msfn(sinphi, cosphi, P->es);
        ml1 = pj_authalic_lat_q(sinphi, P);
        if (secant) { /* secant cone */
            double ml2, m2;

            sinphi = sin(Q->phi2);
            cosphi = cos(Q->phi2);
            m2 = pj_msfn(sinphi, cosphi, P->es);
            ml2 = pj_authalic_lat_q(sinphi, P);
            if (ml2 == ml1)
                return pj_aea_destructor(P, 0);

            Q->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
            if (Q->n == 0) {
                // Not quite, but es is very close to 1...
                proj_log_error(P, _("Invalid value for eccentricity"));
                return pj_aea_destructor(P,
                                         PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
            }
        }
        Q->ec = 1. - .5 * P->one_es * log((1. - P->e) / (1. + P->e)) / P->e;
        Q->c = m1 * m1 + Q->n * ml1;
        Q->dd = 1. / Q->n;
        Q->rho0 = Q->dd *
                  sqrt(Q->c - Q->n * pj_authalic_lat_q(sin(P->phi0), P));
    } else {
        if (secant)
            Q->n = .5 * (Q->n + sin(Q->phi2));
        Q->n2 = Q->n + Q->n;
        Q->c = cosphi * cosphi + Q->n2 * sinphi;
        Q->dd = 1. / Q->n;
        Q->rho0 = Q->dd * sqrt(Q->c - Q->n2 * sin(P->phi0));
    }

    return P;
}

PJ *PJ_PROJECTION(aea) {
    struct pj_aea *Q =
        static_cast<struct pj_aea *>(calloc(1, sizeof(struct pj_aea)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_aea_destructor;

    Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
    Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
    return setup(P);
}

PJ *PJ_PROJECTION(leac) {
    struct pj_aea *Q =
        static_cast<struct pj_aea *>(calloc(1, sizeof(struct pj_aea)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_aea_destructor;

    Q->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
    Q->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? -M_HALFPI : M_HALFPI;
    return setup(P);
}

#undef EPS10
#undef TOL7

Coverage Report

Created: 2025-06-13 06:29

Line	Count	Source (jump to first uncovered line)
1		/******************************************************************************
2		* Project: PROJ.4
3		* Purpose: Implementation of the aea (Albers Equal Area) projection.
4		* and the leac (Lambert Equal Area Conic) projection
5		* Author: Gerald Evenden (1995)
6		* Thomas Knudsen (2016) - revise/add regression tests
7		*
8		******************************************************************************
9		* Copyright (c) 1995, Gerald Evenden
10		*
11		* Permission is hereby granted, free of charge, to any person obtaining a
12		* copy of this software and associated documentation files (the "Software"),
13		* to deal in the Software without restriction, including without limitation
14		* the rights to use, copy, modify, merge, publish, distribute, sublicense,
15		* and/or sell copies of the Software, and to permit persons to whom the
16		* Software is furnished to do so, subject to the following conditions:
17		*
18		* The above copyright notice and this permission notice shall be included
19		* in all copies or substantial portions of the Software.
20		*
21		* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22		* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23		* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24		* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25		* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26		* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27		* DEALINGS IN THE SOFTWARE.
28		*****************************************************************************/
29
30		#include "proj.h"
31		#include "proj_internal.h"
32		#include <errno.h>
33		#include <math.h>
34
35	1.80k	#define EPS10 1.e-10
36	0	#define TOL7 1.e-7
37
38		PROJ_HEAD(aea, "Albers Equal Area") "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
39		PROJ_HEAD(leac, "Lambert Equal Area Conic")
40		"\n\tConic, Sph&Ell\n\tlat_1= south";
41
42		namespace { // anonymous namespace
43		struct pj_aea {
44		double ec;
45		double n;
46		double c;
47		double dd;
48		double n2;
49		double rho0;
50		double rho;
51		double phi1;
52		double phi2;
53		int ellips;
54		double *apa;
55		double qp;
56		};
57		} // anonymous namespace
58
59	908	static PJ pj_aea_destructor(PJ P, int errlev) { /* Destructor */
60	908	if (nullptr == P)
61	0	return nullptr;
62
63	908	if (nullptr == P->opaque)
64	0	return pj_default_destructor(P, errlev);
65
66	908	free(static_cast<struct pj_aea *>(P->opaque)->apa);
67
68	908	return pj_default_destructor(P, errlev);
69	908	}
70
71	0	static PJ_XY aea_e_forward(PJ_LP lp, PJ P) { / Ellipsoid/spheroid, forward */
72	0	PJ_XY xy = {0.0, 0.0};
73	0	struct pj_aea Q = static_cast<struct pj_aea >(P->opaque);
74	0	Q->rho = Q->c - (Q->ellips ? Q->n * pj_authalic_lat_q(sin(lp.phi), P)
75	0	: Q->n2 * sin(lp.phi));
76	0	if (Q->rho < 0.) {
77	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
78	0	return xy;
79	0	}
80	0	Q->rho = Q->dd * sqrt(Q->rho);
81	0	lp.lam *= Q->n;
82	0	xy.x = Q->rho * sin(lp.lam);
83	0	xy.y = Q->rho0 - Q->rho * cos(lp.lam);
84	0	return xy;
85	0	}
86
87	0	static PJ_LP aea_e_inverse(PJ_XY xy, PJ P) { / Ellipsoid/spheroid, inverse */
88	0	PJ_LP lp = {0.0, 0.0};
89	0	struct pj_aea Q = static_cast<struct pj_aea >(P->opaque);
90	0	xy.y = Q->rho0 - xy.y;
91	0	Q->rho = hypot(xy.x, xy.y);
92	0	if (Q->rho != 0.0) {
93	0	if (Q->n < 0.) {
94	0	Q->rho = -Q->rho;
95	0	xy.x = -xy.x;
96	0	xy.y = -xy.y;
97	0	}
98	0	lp.phi = Q->rho / Q->dd;
99	0	if (Q->ellips) {
100	0	const double qs = (Q->c - lp.phi * lp.phi) / Q->n;
101	0	if (fabs(Q->ec - fabs(qs)) > TOL7) {
102	0	if (fabs(qs) > 2) {
103	0	proj_errno_set(
104	0	P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
105	0	return lp;
106	0	}
107	0	lp.phi = pj_authalic_lat_inverse(asin(qs / Q->qp), Q->apa,
108	0	P, Q->qp);
109	0	if (lp.phi == HUGE_VAL) {
110	0	proj_errno_set(
111	0	P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
112	0	return lp;
113	0	}
114	0	} else
115	0	lp.phi = qs < 0. ? -M_HALFPI : M_HALFPI;
116	0	} else {
117	0	const double qs_div_2 = (Q->c - lp.phi * lp.phi) / Q->n2;
118	0	if (fabs(qs_div_2) <= 1.)
119	0	lp.phi = asin(qs_div_2);
120	0	else
121	0	lp.phi = qs_div_2 < 0. ? -M_HALFPI : M_HALFPI;
122	0	}
123	0	lp.lam = atan2(xy.x, xy.y) / Q->n;
124	0	} else {
125	0	lp.lam = 0.;
126	0	lp.phi = Q->n > 0. ? M_HALFPI : -M_HALFPI;
127	0	}
128	0	return lp;
129	0	}
130
131	908	static PJ setup(PJ P) {
132	908	struct pj_aea Q = static_cast<struct pj_aea >(P->opaque);
133
134	908	P->inv = aea_e_inverse;
135	908	P->fwd = aea_e_forward;
136
137	908	if (fabs(Q->phi1) > M_HALFPI) {
138	0	proj_log_error(P,
139	0	_("Invalid value for lat_1: \|lat_1\| should be <= 90°"));
140	0	return pj_aea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
141	0	}
142	908	if (fabs(Q->phi2) > M_HALFPI) {
143	1	proj_log_error(P,
144	1	_("Invalid value for lat_2: \|lat_2\| should be <= 90°"));
145	1	return pj_aea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
146	1	}
147	907	if (fabs(Q->phi1 + Q->phi2) < EPS10) {
148	6	proj_log_error(P, _("Invalid value for lat_1 and lat_2: \|lat_1 + "
149	6	"lat_2\| should be > 0"));
150	6	return pj_aea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
151	6	}
152	901	double sinphi = sin(Q->phi1);
153	901	Q->n = sinphi;
154	901	double cosphi = cos(Q->phi1);
155	901	const int secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
156	901	Q->ellips = (P->es > 0.);
157	901	if (Q->ellips) {
158	390	double ml1, m1;
159
160	390	Q->apa = pj_authalic_lat_compute_coeffs(P->n);
161	390	if (Q->apa == nullptr)
162	0	return pj_aea_destructor(P, 0);
163	390	Q->qp = pj_authalic_lat_q(1.0, P);
164	390	m1 = pj_msfn(sinphi, cosphi, P->es);
165	390	ml1 = pj_authalic_lat_q(sinphi, P);
166	390	if (secant) { /* secant cone */
167	389	double ml2, m2;
168
169	389	sinphi = sin(Q->phi2);
170	389	cosphi = cos(Q->phi2);
171	389	m2 = pj_msfn(sinphi, cosphi, P->es);
172	389	ml2 = pj_authalic_lat_q(sinphi, P);
173	389	if (ml2 == ml1)
174	0	return pj_aea_destructor(P, 0);
175
176	389	Q->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
177	389	if (Q->n == 0) {
178		// Not quite, but es is very close to 1...
179	0	proj_log_error(P, _("Invalid value for eccentricity"));
180	0	return pj_aea_destructor(P,
181	0	PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
182	0	}
183	389	}
184	390	Q->ec = 1. - .5 * P->one_es * log((1. - P->e) / (1. + P->e)) / P->e;
185	390	Q->c = m1 * m1 + Q->n * ml1;
186	390	Q->dd = 1. / Q->n;
187	390	Q->rho0 = Q->dd *
188	390	sqrt(Q->c - Q->n * pj_authalic_lat_q(sin(P->phi0), P));
189	511	} else {
190	511	if (secant)
191	501	Q->n = .5 * (Q->n + sin(Q->phi2));
192	511	Q->n2 = Q->n + Q->n;
193	511	Q->c = cosphi * cosphi + Q->n2 * sinphi;
194	511	Q->dd = 1. / Q->n;
195	511	Q->rho0 = Q->dd * sqrt(Q->c - Q->n2 * sin(P->phi0));
196	511	}
197
198	901	return P;
199	901	}
200
201	6	PJ *PJ_PROJECTION(aea) {
202	6	struct pj_aea *Q =
203	6	static_cast<struct pj_aea *>(calloc(1, sizeof(struct pj_aea)));
204	6	if (nullptr == Q)
205	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
206	6	P->opaque = Q;
207	6	P->destructor = pj_aea_destructor;
208
209	6	Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
210	6	Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
211	6	return setup(P);
212	6	}
213
214	902	PJ *PJ_PROJECTION(leac) {
215	902	struct pj_aea *Q =
216	902	static_cast<struct pj_aea *>(calloc(1, sizeof(struct pj_aea)));
217	902	if (nullptr == Q)
218	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
219	902	P->opaque = Q;
220	902	P->destructor = pj_aea_destructor;
221
222	902	Q->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
223	902	Q->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? -M_HALFPI : M_HALFPI;
224	902	return setup(P);
225	902	}
226
227		#undef EPS10
228		#undef TOL7