/src/proj/src/projections/cea.cpp
Line | Count | Source |
1 | | #include <errno.h> |
2 | | #include <math.h> |
3 | | |
4 | | #include "proj.h" |
5 | | #include "proj_internal.h" |
6 | | |
7 | | namespace { // anonymous namespace |
8 | | struct pj_cea_data { |
9 | | double qp; |
10 | | double *apa; |
11 | | }; |
12 | | } // anonymous namespace |
13 | | |
14 | | PROJ_HEAD(cea, "Equal Area Cylindrical") "\n\tCyl, Sph&Ell\n\tlat_ts="; |
15 | 0 | #define EPS 1e-10 |
16 | | |
17 | 0 | static PJ_XY cea_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ |
18 | 0 | PJ_XY xy = {0.0, 0.0}; |
19 | 0 | xy.x = P->k0 * lp.lam; |
20 | 0 | xy.y = 0.5 * pj_authalic_lat_q(sin(lp.phi), P) / P->k0; |
21 | 0 | return xy; |
22 | 0 | } |
23 | | |
24 | 0 | static PJ_XY cea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ |
25 | 0 | PJ_XY xy = {0.0, 0.0}; |
26 | 0 | xy.x = P->k0 * lp.lam; |
27 | 0 | xy.y = sin(lp.phi) / P->k0; |
28 | 0 | return xy; |
29 | 0 | } |
30 | | |
31 | 0 | static PJ_LP cea_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ |
32 | 0 | PJ_LP lp = {0.0, 0.0}; |
33 | 0 | const struct pj_cea_data *Q = |
34 | 0 | static_cast<const struct pj_cea_data *>(P->opaque); |
35 | 0 | lp.phi = pj_authalic_lat_inverse(asin(2. * xy.y * P->k0 / Q->qp), Q->apa, P, |
36 | 0 | Q->qp); |
37 | 0 | lp.lam = xy.x / P->k0; |
38 | 0 | return lp; |
39 | 0 | } |
40 | | |
41 | 0 | static PJ_LP cea_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ |
42 | 0 | PJ_LP lp = {0.0, 0.0}; |
43 | |
|
44 | 0 | xy.y *= P->k0; |
45 | 0 | const double t = fabs(xy.y); |
46 | 0 | if (t - EPS <= 1.) { |
47 | 0 | if (t >= 1.) |
48 | 0 | lp.phi = xy.y < 0. ? -M_HALFPI : M_HALFPI; |
49 | 0 | else |
50 | 0 | lp.phi = asin(xy.y); |
51 | 0 | lp.lam = xy.x / P->k0; |
52 | 0 | } else { |
53 | 0 | proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); |
54 | 0 | return lp; |
55 | 0 | } |
56 | 0 | return (lp); |
57 | 0 | } |
58 | | |
59 | 0 | static PJ *pj_cea_destructor(PJ *P, int errlev) { /* Destructor */ |
60 | 0 | if (nullptr == P) |
61 | 0 | return nullptr; |
62 | | |
63 | 0 | if (nullptr == P->opaque) |
64 | 0 | return pj_default_destructor(P, errlev); |
65 | | |
66 | 0 | free(static_cast<struct pj_cea_data *>(P->opaque)->apa); |
67 | 0 | return pj_default_destructor(P, errlev); |
68 | 0 | } |
69 | | |
70 | 0 | PJ *PJ_PROJECTION(cea) { |
71 | 0 | double t = 0.0; |
72 | 0 | struct pj_cea_data *Q = static_cast<struct pj_cea_data *>( |
73 | 0 | calloc(1, sizeof(struct pj_cea_data))); |
74 | 0 | if (nullptr == Q) |
75 | 0 | return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
76 | 0 | P->opaque = Q; |
77 | 0 | P->destructor = pj_cea_destructor; |
78 | |
|
79 | 0 | if (pj_param(P->ctx, P->params, "tlat_ts").i) { |
80 | 0 | t = pj_param(P->ctx, P->params, "rlat_ts").f; |
81 | 0 | P->k0 = cos(t); |
82 | 0 | if (P->k0 < 0.) { |
83 | 0 | proj_log_error( |
84 | 0 | P, _("Invalid value for lat_ts: |lat_ts| should be <= 90°")); |
85 | 0 | return pj_default_destructor(P, |
86 | 0 | PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE); |
87 | 0 | } |
88 | 0 | } |
89 | 0 | if (P->es != 0.0) { |
90 | 0 | t = sin(t); |
91 | 0 | P->k0 /= sqrt(1. - P->es * t * t); |
92 | 0 | P->e = sqrt(P->es); |
93 | 0 | Q->apa = pj_authalic_lat_compute_coeffs(P->n); |
94 | 0 | if (!(Q->apa)) |
95 | 0 | return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/); |
96 | | |
97 | 0 | Q->qp = pj_authalic_lat_q(1.0, P); |
98 | 0 | P->inv = cea_e_inverse; |
99 | 0 | P->fwd = cea_e_forward; |
100 | 0 | } else { |
101 | 0 | P->inv = cea_s_inverse; |
102 | 0 | P->fwd = cea_s_forward; |
103 | 0 | } |
104 | | |
105 | 0 | return P; |
106 | 0 | } |
107 | | |
108 | | #undef EPS |