1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
|
/*****************************************************************************
Lambert Conformal Conic Alternative
-----------------------------------
This is Gerald Evenden's 2003 implementation of an alternative
"almost" LCC, which has been in use historically, but which
should NOT be used for new projects - i.e: use this implementation
if you need interoperability with old data represented in this
projection, but not in any other case.
The code was originally discussed on the PROJ.4 mailing list in
a thread archived over at
http://lists.maptools.org/pipermail/proj/2003-March/000644.html
It was discussed again in the thread starting at
http://lists.maptools.org/pipermail/proj/2017-October/007828.html
and continuing at
http://lists.maptools.org/pipermail/proj/2017-November/007831.html
which prompted Clifford J. Mugnier to add these clarifying notes:
The French Army Truncated Cubic Lambert (partially conformal) Conic
projection is the Legal system for the projection in France between
the late 1800s and 1948 when the French Legislature changed the law
to recognize the fully conformal version.
It was (might still be in one or two North African prior French
Colonies) used in North Africa in Algeria, Tunisia, & Morocco, as
well as in Syria during the Levant.
Last time I have seen it used was about 30+ years ago in
Algeria when it was used to define Lease Block boundaries for
Petroleum Exploration & Production.
(signed)
Clifford J. Mugnier, c.p., c.m.s.
Chief of Geodesy
LSU Center for GeoInformatics
Dept. of Civil Engineering
LOUISIANA STATE UNIVERSITY
*****************************************************************************/
#define PJ_LIB__
#include <errno.h>
#include <math.h>
#include "proj.h"
#include "proj_internal.h"
PROJ_HEAD(lcca, "Lambert Conformal Conic Alternative")
"\n\tConic, Sph&Ell\n\tlat_0=";
#define MAX_ITER 10
#define DEL_TOL 1e-12
namespace { // anonymous namespace
struct pj_opaque {
double *en;
double r0, l, M0;
double C;
};
} // anonymous namespace
static double fS(double S, double C) { /* func to compute dr */
return S * ( 1. + S * S * C);
}
static double fSp(double S, double C) { /* deriv of fs */
return 1. + 3.* S * S * C;
}
static PJ_XY lcca_e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
PJ_XY xy = {0.0,0.0};
struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
double S, r, dr;
S = pj_mlfn(lp.phi, sin(lp.phi), cos(lp.phi), Q->en) - Q->M0;
dr = fS(S, Q->C);
r = Q->r0 - dr;
const double lam_mul_l = lp.lam * Q->l;
xy.x = P->k0 * (r * sin(lam_mul_l));
xy.y = P->k0 * (Q->r0 - r * cos(lam_mul_l) );
return xy;
}
static PJ_LP lcca_e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
PJ_LP lp = {0.0,0.0};
struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
double theta, dr, S, dif;
int i;
xy.x /= P->k0;
xy.y /= P->k0;
theta = atan2(xy.x , Q->r0 - xy.y);
dr = xy.y - xy.x * tan(0.5 * theta);
lp.lam = theta / Q->l;
S = dr;
for (i = MAX_ITER; i ; --i) {
S -= (dif = (fS(S, Q->C) - dr) / fSp(S, Q->C));
if (fabs(dif) < DEL_TOL) break;
}
if (!i) {
proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
return lp;
}
lp.phi = pj_inv_mlfn(P->ctx, S + Q->M0, P->es, Q->en);
return lp;
}
static PJ *destructor (PJ *P, int errlev) {
if (nullptr==P)
return nullptr;
if (nullptr==P->opaque)
return pj_default_destructor (P, errlev);
free (static_cast<struct pj_opaque*>(P->opaque)->en);
return pj_default_destructor (P, errlev);
}
PJ *PROJECTION(lcca) {
double s2p0, N0, R0, tan0;
struct pj_opaque *Q = static_cast<struct pj_opaque*>(calloc (1, sizeof (struct pj_opaque)));
if (nullptr==Q)
return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
P->opaque = Q;
(Q->en = pj_enfn(P->es));
if (!Q->en)
return pj_default_destructor (P, PROJ_ERR_OTHER /*ENOMEM*/);
if (P->phi0 == 0.) {
proj_log_error(P, _("Invalid value for lat_0: it should be different from 0."));
return destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
}
Q->l = sin(P->phi0);
Q->M0 = pj_mlfn(P->phi0, Q->l, cos(P->phi0), Q->en);
s2p0 = Q->l * Q->l;
R0 = 1. / (1. - P->es * s2p0);
N0 = sqrt(R0);
R0 *= P->one_es * N0;
tan0 = tan(P->phi0);
Q->r0 = N0 / tan0;
Q->C = 1. / (6. * R0 * N0);
P->inv = lcca_e_inverse;
P->fwd = lcca_e_forward;
P->destructor = destructor;
return P;
}
|
