1
/******************************************************************************
2
* $Id: PJ_krovak.c,v 1.4 2002/12/15 22:31:04 warmerda Exp $
5
* Purpose: Implementation of the krovak (Krovak) projection.
6
* Definition: http://www.ihsenergy.com/epsg/guid7.html#1.4.3
7
* Author: Thomas Flemming, tf@ttqv.com
9
******************************************************************************
10
* Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
12
* Permission is hereby granted, free of charge, to any person obtaining a
13
* copy of this software and associated documentation files (the "Software"),
14
* to deal in the Software without restriction, including without limitation
15
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
16
* and/or sell copies of the Software, and to permit persons to whom the
17
* Software is furnished to do so, subject to the following conditions:
19
* The above copyright notice and this permission notice shall be included
20
* in all copies or substantial portions of the Software.
22
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
25
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
26
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
27
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
28
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
30
******************************************************************************
32
* $Log: PJ_krovak.c,v $
33
* Revision 1.4 2002/12/15 22:31:04 warmerda
34
* handle lon_0, k, and prime meridian properly
36
* Revision 1.3 2002/12/15 00:13:30 warmerda
37
* lat_0 may now be set by user, but still defaults to 49d30N
39
* Revision 1.2 2002/12/14 19:35:21 warmerda
44
#define PROJ_PARMS__ \
52
PJ_CVSID("$Id: PJ_krovak.c,v 1.4 2002/12/15 22:31:04 warmerda Exp $");
54
PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Sph.";
57
NOTES: According to EPSG the full Krovak projection method should have
58
the following parameters. Within PROJ.4 the azimuth, and pseudo
59
standard parallel are hardcoded in the algorithm and can't be
60
altered from outside. The others all have defaults to match the
61
common usage with Krovak projection.
63
lat_0 = latitude of centre of the projection
65
lon_0 = longitude of centre of the projection
67
** = azimuth (true) of the centre line passing through the centre of the projection
69
** = latitude of pseudo standard parallel
71
k = scale factor on the pseudo standard parallel
73
x_0 = False Easting of the centre of the projection at the apex of the cone
75
y_0 = False Northing of the centre of the projection at the apex of the cone
81
FORWARD(s_forward); /* spheroid */
82
/* calculate xy from lat/lon */
87
/* Constants, identical to inverse transform function */
88
double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
89
double gfi, u, fi0, lon17, lamdd, deltav, s, d, eps, ro;
92
s45 = 0.785398163397448; /* 45� */
94
fi0 = P->phi0; /* Latitude of projection centre 49� 30' */
96
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
98
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
101
a = 1; /* 6377397.155; */
102
/* e2 = P->es;*/ /* 0.006674372230614; */
103
e2 = 0.006674372230614;
106
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
108
uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */
109
u0 = asin(sin(fi0) / alfa);
110
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. );
112
k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g;
115
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
116
s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78� 30'00" N */
118
ro0 = k1 * n0 / tan(s0);
123
gfi =pow ( ((1. + e * sin(lp.phi)) /
124
(1. - e * sin(lp.phi))) , (alfa * e / 2.));
126
u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
128
deltav = - lp.lam * alfa;
130
s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
131
d = asin(cos(u) * sin(deltav) / cos(s));
133
ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n) ;
135
/* x and y are reverted! */
136
xy.y = ro * cos(eps) / a;
137
xy.x = ro * sin(eps) / a;
140
strcpy(errmess,"a: ");
142
ltoa((long)(a*1000000000),tmp,10);
144
strcat(errmess,"e2: ");
146
ltoa((long)(e2*1000000000),tmp,10);
149
MessageBox(NULL, errmess, NULL, 0);
158
INVERSE(s_inverse); /* spheroid */
159
/* calculate lat/lon from xy */
161
/* Constants, identisch wie in der Umkehrfunktion */
162
double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
163
double u, l24, lamdd, deltav, s, d, eps, ro, fi1, xy0, lon17;
166
s45 = 0.785398163397448; /* 45� */
168
fi0 = P->phi0; /* Latitude of projection centre 49� 30' */
171
/* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
173
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
174
e2=0.006674372230614;
176
a = 1; /* 6377397.155; */
177
/* e2 = P->es; */ /* 0.006674372230614; */
178
e2 = 0.006674372230614;
181
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
182
uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */
183
u0 = asin(sin(fi0) / alfa);
184
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2. );
186
k = tan( u0 / 2. + s45) / pow (tan(fi0 / 2. + s45) , alfa) * g;
189
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
190
s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78� 30'00" N */
192
ro0 = k1 * n0 / tan(s0);
202
ro = sqrt(xy.x * xy.x + xy.y * xy.y);
203
eps = atan2(xy.y, xy.x);
205
s = 2. * (atan( pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
207
u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
208
deltav = asin(cos(s) * sin(d) / cos(u));
210
lp.lam = P->lam0 - deltav / alfa;
212
/* ITERATION FOR lp.phi */
218
lp.phi = 2. * ( atan( pow( k, -1. / alfa) *
219
pow( tan(u / 2. + s45) , 1. / alfa) *
220
pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
223
if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1;
234
FREEUP; if (P) pj_dalloc(P); }
238
/* read some Parameters,
239
* here Latitude Truescale */
241
ts = pj_param(P->params, "rlat_ts").f;
244
/* we want Bessel as fixed ellipsoid */
246
P->e = sqrt(P->es = 0.006674372230614);
248
/* if latitude of projection center is not set, use 49d30'N */
249
if (!pj_param(P->params, "tlat_0").i)
250
P->phi0 = 0.863937979737193;
252
/* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
253
/* that will correspond to using longitudes relative to greenwich */
254
/* as input and output, instead of lat/long relative to Ferro */
255
if (!pj_param(P->params, "tlon_0").i)
256
P->lam0 = 0.7417649320975901 - 0.308341501185665;
259
/* if scale not set default to 0.9999 */
260
if (!pj_param(P->params, "tk").i)
263
/* always the same */
added
removed
src/PJ_krovak.c
