[OSRS-PROJ] PROJ4 and Krovak projection (Czech Republic)?
Thomas Flemming
tf at ttqv.com
Thu Apr 4 03:32:40 EST 2002
Hi,
here it is.
The code was translated from a fortran77 source provided from the
czechic Research Institute of Geodesy, Topography and Cartography to c.
The main problem for me was how to extend the proj library with a new
projection.
I don't think, that I solved this perfectly, but it works...
Hope, this will help.
Tom
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* pj_qv_krovak.c
*
* Krovak Projection for Proj.4
*
* last change 16.1.2002
*
* Copyright (c) 2001, Thomas Flemming, tf at ttqv.com
*
* 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.
*
*/
#ifndef lint
static const char SCCSID[]="@(#)pj_qv_krovak.c 4.1 94/02/15 GIE REL";
#endif
#define PROJ_PARMS__ \
double C_x;
#define PJ_LIB__
#include <projects.h>
//
// insert this line to pj_list.h
//
PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Sph.";
FORWARD(s_forward); /* spheroid */
//xy von lat/lon ausrechnen
char errmess[255];
char tmp[16];
// Constants, identisch wie in der Umkehrfunktion
double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a,
s0, n;
double gfi, u, lon17, lon42, lamdd, deltav, s, d, eps, ro;
s45 = 0.785398163397448; // 45°
s90 = 2 * s45;
fi0 = 0.863937979737193; //Latitude of projection centre 49° 30//
/* Ellipsoid wird als Übergabeparameter in for.c und inv.c
versrbeitet
daher muß a hier auf 1 gesetzt werden
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.15281,
e2=0.006674372230614;
*/
a = 1; //6377397.155;
e2 = P->es; //0.006674372230614;
e = sqrt(e2);
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
uq = 1.04216856380474; //DU(2, 59, 42, 42.69689)
u0 = asin(sin(fi0) / alfa);
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) ,
alfa * e / 2. );
k = tan( u0 / 2. + s45) /
pow (tan(fi0 / 2. + s45) , alfa) *
g;
k1 = 0.9999;
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
s0 = 1.37008346281555; //Latitude of pseudo standard parallel
78° 30'00" N
n = sin(s0);
ro0 = k1 * n0 / tan(s0);
ad = s90 - uq;
// Transformation
gfi =pow ( ((1. + e * sin(lp.phi)) /
(1. - e * sin(lp.phi))) , (alfa * e / 2.));
u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
lon17 = 0.308341501185665; // Longitude
of Ferro is 17° 40'00" West of Greenwich
lon42 = 0.74176493209759; // Longitude
of Origin 42° 30'00" East of Ferro
lamdd = lp.lam + lon17;
deltav = alfa * (lon42 - lamdd);
s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
d = asin(cos(u) * sin(deltav) / cos(s));
eps = n * d;
ro = ro0 * pow(tan(s0 / 2. + s45) , n) /
pow(tan(s / 2. + s45) , n) ;
// x und y vertauscht!
xy.y = ro * cos(eps);
xy.x = ro * sin(eps);
/* //debug
strcpy(errmess,"a: ");
strcpy(tmp," ");
ltoa((long)(a*1000000000),tmp,10);
strcat(errmess,tmp);
strcat(errmess,"e2: ");
strcpy(tmp," ");
ltoa((long)(e2*1000000000),tmp,10);
strcat(errmess,tmp);
MessageBox(NULL, errmess, NULL, 0);
*/
return (xy);
}
INVERSE(s_inverse); /* spheroid */
// lat/lon von xy ausrechnen
// Constants, identisch wie in der Umkehrfunktion
double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a,
s0, n;
double u, l24, lamdd, deltav, s, d, eps, ro, fi1, xy0;
int ok;
s45 = 0.785398163397448; // 45°
s90 = 2 * s45;
fi0 = 0.863937979737193; //Latitude of projection centre 49° 30//
/* Ellipsoid wird als Übergabeparameter in for.c und inv.c
versrbeitet
daher muß a hier auf 1 gesetzt werden
Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.15281,
e2=0.006674372230614;
*/
a = 1; //6377397.155;
e2 = P->es; //0.006674372230614;
e = sqrt(e2);
alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
uq = 1.04216856380474; //DU(2, 59, 42, 42.69689)
u0 = asin(sin(fi0) / alfa);
g = pow( (1. + e * sin(fi0)) / (1. - e * sin(fi0)) ,
alfa * e / 2. );
k = tan( u0 / 2. + s45) /
pow (tan(fi0 / 2. + s45) , alfa) * g;
k1 = 0.9999;
n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
s0 = 1.37008346281555; //Latitude of pseudo standard parallel
78° 30'00" N
n = sin(s0);
ro0 = k1 * n0 / tan(s0);
ad = s90 - uq;
// Transformation
//vertauschen
xy0=xy.x;
xy.x=xy.y;
xy.y=xy0;
ro = sqrt(xy.x * xy.x + xy.y * xy.y);
eps = atan2(xy.y, xy.x);
d = eps / sin(s0);
s = 2. * (atan( pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
deltav = asin(cos(s) * sin(d) / cos(u));
l24 = 0.433423430911925; //DU(2, 24, 50, 0.)
lp.lam = l24 - deltav / alfa;
// ITERATION FOR lp.phi
fi1 = u;
ok = 0;
do
{
lp.phi = 2. * ( atan( pow( k, -1. / alfa) *
pow( tan(u / 2. + s45) , 1. / alfa) *
pow( (1. + e * sin(fi1)) / (1. - e *
sin(fi1)) , e / 2.)
) - s45);
if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1;
fi1 = lp.phi;
}
while (ok==0);
return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
ENTRY0(krovak)
double ts;
// irgendwelche Parameter auslesen,
// hier Latitude Truescale
ts = pj_param(P->params, "rlat_ts").f;
P->C_x = ts;
// immer gleich
//P->es = 0.;
P->inv = s_inverse; P->fwd = s_forward;
ENDENTRY(P)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
Markus Neteler schrieb:
>
> On Wed, Apr 03, 2002 at 07:35:59PM +0200, Thomas Flemming wrote:
> > Hi Markus,
> >
> > last year I asked the same question, because we need this projection in
> > our GPS-Software.
> > I've got no answer, so I figured it out myself. It works, but I'm not
> > sure, if it is really correct integrated in proj.4.
> > For me, it is still very hard, to understand, how proj works. Anyway,
> > this is not my job.
>
> Hi Thomas,
>
> also we have a hard time currently as we want to integrate the ISIN
> projection for MODIS/TERRA stellite data. That's far from trivial.
>
> > If you like, I send you the c-sourcefile for the Krovak-Projection for
> > Proj.4.
>
> Oh great - I'll appreciate that. You efforts will also help the GRASS
> GIS software where our Czech users may need it.
>
> Thanks in advance,
>
> Markus
>
> > Markus Neteler schrieb:
> > >
> > > Hi,
> > >
> > > I am curious if it is possible to define the Krovak projection
> > > in PROJ4 (finally I want to do that in GRASS).
> > >
> > > The definition is:
> > > http://www.ihsenergy.com/epsg/guid7.html#1.4.3
> > >
> > > Parameters:
> > > Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.15281
> > > then e = 0.081696831 e2 = 0.006674372
> > > Latitude of projection centre (j_c) 49° 30'00" N = 0.604756586 rad
> > > Longitude of Origin 42° 30'00" East of Ferro
> > > Longitude of Ferro is 17° 40'00" West of Greenwich
> > > Longitude of Origin (l_c) 24° 50'00" East of Greenwich=0.433423431 rad
> > > Latitude of pseudo standard parallel (j_1) 78° 30'00" N
> > > Azimuth of centre line (a_c) 30° 17' 17.303"
> > > Scale factor on pseudo Standard Parallel (ko) 0.99990
> > > Easting at projection centre (Ec) 0.00 m
> > > Northing at projection centre (Nc) 0.00 m
> > >
> > > Projection constants:
> > > B = 1.000597498
> > > A = 6380703.61
> > > g 0 = 0.863239103
> > > t0 = 1.003419164
> > > n = 0.979924705
> > > r0 = 1298039.005
> > >
> > > It might be also called "S-JTSK Ferro Krovak".
> > >
> > > As being only experienced a bit with Transverse Mercator, I am a bit
> > > lost with this Krovak projection.
> > >
> > > Thanks in advance for any hint (or a proj command line...)
> > >
> > > Markus Neteler
> > >
> > > PS: There is also another description with all history:
> > > http://www.asprs.org/asprs/resources/grids/index.html
> > > -> The Czech Republic
> > > ----------------------------------------
> > > PROJ.4 Discussion List
> > > See http://www.remotesensing.org/proj for subscription, unsubscription
> > > and other information.
> >
> > --
> >
> > ****************************************
> > ** Dipl.-Ing. Thomas Flemming
> > ** Software Development
> > ** thomas.flemming at touratech.de
> > ** tf at ttqv.com
> > ** http://www.ttqv.com
> > **
> > ** ... und immer dem Pfeil nach!
> > ***************************************
> > ----------------------------------------
> > PROJ.4 Discussion List
> > See http://www.remotesensing.org/proj for subscription, unsubscription
> > and other information.
>
> --
> Markus Neteler
>
> ITC-irst, Istituto per la Ricerca Scientifica e Tecnologica
> Project on Predictive Models for the Environment
> Via Sommarive, 18 - 38050 Povo (Trento), Italia
> tel +39 0461 314 -520 (fax -591) http://mpa.itc.it
> ----------------------------------------
> PROJ.4 Discussion List
> See http://www.remotesensing.org/proj for subscription, unsubscription
> and other information.
--
****************************************
** Dipl.-Ing. Thomas Flemming
** Software Development
** thomas.flemming at touratech.de
** tf at ttqv.com
** http://www.ttqv.com
**
** ... und immer dem Pfeil nach!
***************************************
----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.
More information about the Proj
mailing list