1 /*
2 * Project: PROJ.4
3 * Purpose: Implementation of the krovak (Krovak) projection.
4 * Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
5 * Author: Thomas Flemming, tf@ttqv.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
9 *
10 * Permission is hereby granted, free of charge, to any person obtaining a
11 * copy of this software and associated documentation files (the "Software"),
12 * to deal in the Software without restriction, including without limitation
13 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 * and/or sell copies of the Software, and to permit persons to whom the
15 * Software is furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included
18 * in all copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
21 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
22 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
23 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
24 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
25 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
26 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
27 * SOFTWARE.
28 ******************************************************************************
29 * A description of the (forward) projection is found in:
30 *
31 * Bohuslav Veverka,
32 *
33 * KROVAK’S PROJECTION AND ITS USE FOR THE
34 * CZECH REPUBLIC AND THE SLOVAK REPUBLIC,
35 *
36 * 50 years of the Research Institute of
37 * and the Slovak Republic Geodesy, Topography and Cartography
38 *
39 * which can be found via the Wayback Machine:
40 *
41 * https://web.archive.org/web/20150216143806/https://www.vugtk.cz/odis/sborniky/sb2005/Sbornik_50_let_VUGTK/Part_1-Scientific_Contribution/16-Veverka.pdf
42 *
43 * Further info, including the inverse projection, is given by EPSG:
44 *
45 * Guidance Note 7 part 2
46 * Coordinate Conversions and Transformations including Formulas
47 *
48 * http://www.iogp.org/pubs/373-07-2.pdf
49 *
50 * Variable names in this file mostly follows what is used in the
51 * paper by Veverka.
52 *
53 * According to EPSG the full Krovak projection method should have
54 * the following parameters. Within PROJ.4 the azimuth, and pseudo
55 * standard parallel are hardcoded in the algorithm and can't be
56 * altered from outside. The others all have defaults to match the
57 * common usage with Krovak projection.
58 *
59 * lat_0 = latitude of centre of the projection
60 *
61 * lon_0 = longitude of centre of the projection
62 *
63 * ** = azimuth (true) of the centre line passing through the
64 * centre of the projection
65 *
66 * ** = latitude of pseudo standard parallel
67 *
68 * k = scale factor on the pseudo standard parallel
69 *
70 * x_0 = False Easting of the centre of the projection at the
71 * apex of the cone
72 *
73 * y_0 = False Northing of the centre of the projection at
74 * the apex of the cone
75 *
76 *****************************************************************************/
77
78
79 #define PJ_LIB__
80 #include <projects.h>
81
82 PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Ellps.";
83
84 #define EPS 1e-15
85 #define S45 0.785398163397448 /* 45 deg */
86 #define S90 1.570796326794896 /* 90 deg */
87 #define UQ 1.04216856380474 /* DU(2, 59, 42, 42.69689) */
88 #define S0 1.37008346281555 /* Latitude of pseudo standard parallel 78deg 30'00" N */
89
90 struct pj_opaque {
91 double alpha;
92 double k;
93 double n;
94 double rho0;
95 double ad;
96 int czech;
97 };
98
99
e_forward(LP lp,PJ * P)100 static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
101 struct pj_opaque *Q = P->opaque;
102 XY xy = {0.0,0.0};
103
104 double gfi, u, deltav, s, d, eps, rho;
105
106 gfi = pow ( (1. + P->e * sin(lp.phi)) / (1. - P->e * sin(lp.phi)), Q->alpha * P->e / 2.);
107
108 u = 2. * (atan(Q->k * pow( tan(lp.phi / 2. + S45), Q->alpha) / gfi)-S45);
109 deltav = -lp.lam * Q->alpha;
110
111 s = asin(cos(Q->ad) * sin(u) + sin(Q->ad) * cos(u) * cos(deltav));
112 d = asin(cos(u) * sin(deltav) / cos(s));
113
114 eps = Q->n * d;
115 rho = Q->rho0 * pow(tan(S0 / 2. + S45) , Q->n) / pow(tan(s / 2. + S45) , Q->n);
116
117 xy.y = rho * cos(eps);
118 xy.x = rho * sin(eps);
119
120 xy.y *= Q->czech;
121 xy.x *= Q->czech;
122
123 return xy;
124 }
125
126
e_inverse(XY xy,PJ * P)127 static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
128 struct pj_opaque *Q = P->opaque;
129 LP lp = {0.0,0.0};
130
131 double u, deltav, s, d, eps, rho, fi1, xy0;
132 int ok;
133
134 xy0 = xy.x;
135 xy.x = xy.y;
136 xy.y = xy0;
137
138 xy.x *= Q->czech;
139 xy.y *= Q->czech;
140
141 rho = sqrt(xy.x * xy.x + xy.y * xy.y);
142 eps = atan2(xy.y, xy.x);
143
144 d = eps / sin(S0);
145 s = 2. * (atan( pow(Q->rho0 / rho, 1. / Q->n) * tan(S0 / 2. + S45)) - S45);
146
147 u = asin(cos(Q->ad) * sin(s) - sin(Q->ad) * cos(s) * cos(d));
148 deltav = asin(cos(s) * sin(d) / cos(u));
149
150 lp.lam = P->lam0 - deltav / Q->alpha;
151
152 /* ITERATION FOR lp.phi */
153 fi1 = u;
154
155 ok = 0;
156 do {
157 lp.phi = 2. * ( atan( pow( Q->k, -1. / Q->alpha) *
158 pow( tan(u / 2. + S45) , 1. / Q->alpha) *
159 pow( (1. + P->e * sin(fi1)) / (1. - P->e * sin(fi1)) , P->e / 2.)
160 ) - S45);
161
162 if (fabs(fi1 - lp.phi) < EPS) ok=1;
163 fi1 = lp.phi;
164 } while (ok==0);
165
166 lp.lam -= P->lam0;
167
168 return lp;
169 }
170
171
freeup_new(PJ * P)172 static void *freeup_new (PJ *P) { /* Destructor */
173 if (0==P)
174 return 0;
175 if (0==P->opaque)
176 return pj_dealloc(P);
177
178 pj_dealloc(P->opaque);
179 return pj_dealloc(P);
180 }
181
freeup(PJ * P)182 static void freeup (PJ *P) {
183 freeup_new (P);
184 return;
185 }
186
187
PROJECTION(krovak)188 PJ *PROJECTION(krovak) {
189 double u0, n0, g;
190 struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
191 if (0==Q)
192 return freeup_new (P);
193 P->opaque = Q;
194
195 /* we want Bessel as fixed ellipsoid */
196 P->a = 6377397.155;
197 P->e = sqrt(P->es = 0.006674372230614);
198
199 /* if latitude of projection center is not set, use 49d30'N */
200 if (!pj_param(P->ctx, P->params, "tlat_0").i)
201 P->phi0 = 0.863937979737193;
202
203 /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
204 /* that will correspond to using longitudes relative to greenwich */
205 /* as input and output, instead of lat/long relative to Ferro */
206 if (!pj_param(P->ctx, P->params, "tlon_0").i)
207 P->lam0 = 0.7417649320975901 - 0.308341501185665;
208
209 /* if scale not set default to 0.9999 */
210 if (!pj_param(P->ctx, P->params, "tk").i)
211 P->k0 = 0.9999;
212
213 Q->czech = 1;
214 if( !pj_param(P->ctx, P->params, "tczech").i )
215 Q->czech = -1;
216
217 /* Set up shared parameters between forward and inverse */
218 Q->alpha = sqrt(1. + (P->es * pow(cos(P->phi0), 4)) / (1. - P->es));
219 u0 = asin(sin(P->phi0) / Q->alpha);
220 g = pow( (1. + P->e * sin(P->phi0)) / (1. - P->e * sin(P->phi0)) , Q->alpha * P->e / 2. );
221 Q->k = tan( u0 / 2. + S45) / pow (tan(P->phi0 / 2. + S45) , Q->alpha) * g;
222 n0 = sqrt(1. - P->es) / (1. - P->es * pow(sin(P->phi0), 2));
223 Q->n = sin(S0);
224 Q->rho0 = P->k0 * n0 / tan(S0);
225 Q->ad = S90 - UQ;
226
227 P->inv = e_inverse;
228 P->fwd = e_forward;
229
230 return P;
231 }
232
233
234 #ifndef PJ_SELFTEST
pj_krovak_selftest(void)235 int pj_krovak_selftest (void) {return 0;}
236 #else
237
pj_krovak_selftest(void)238 int pj_krovak_selftest (void) {
239 double tolerance_lp = 1e-10;
240 double tolerance_xy = 1e-7;
241
242 char e_args[] = {"+proj=krovak +ellps=GRS80"};
243
244 LP fwd_in[] = {
245 { 2, 1},
246 { 2,-1},
247 {-2, 1},
248 {-2,-1}
249 };
250
251 XY e_fwd_expect[] = {
252 {-3196535.2325636409, -6617878.8675514441},
253 {-3260035.4405521089, -6898873.6148780314},
254 {-3756305.3288691747, -6478142.5615715114},
255 {-3831703.6585019818, -6759107.1701553948},
256 };
257
258 XY inv_in[] = {
259 { 200, 100},
260 { 200,-100},
261 {-200, 100},
262 {-200,-100}
263 };
264
265 LP e_inv_expect[] = {
266 {24.836218918719162, 59.758403933233858},
267 {24.836315484509566, 59.756888425730189},
268 {24.830447747947495, 59.758403933233858},
269 {24.830351182157091, 59.756888425730189},
270 };
271
272 return pj_generic_selftest (e_args, 0, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, 0, inv_in, e_inv_expect, 0);
273 }
274
275
276 #endif
277