[OSRSPROJ] Krovak and GIS usage
Markus Neteler
neteler at itc.it
Mon May 27 10:22:53 EDT 2002
Dear list,
while adding Krovak Projection to GRASS we came across a
problem. In general the Krovak projection is defined as
westing and southing (not easting and northing):
http://krovak.webpark.cz/index_kartograf.htm
See figure:
http://krovak.webpark.cz/obr/index_1.gif
y <+



x
The current PJ_krovak.c implementation supports positive coordinate values,
however, Czech and Slovak governmental data require support for negative
coordinate values with switched axis (here is the problem).
In an email exchange with Thomas Flemming who implemented PJ_krovak.c we
found out that the current implementation in PROJ4 uses positive values:
PJ_krovak.c is implemented as follows (read from bottom):
^
 user interface: x=westing y=southing
 
 mid layer: switch of axes
 
 Inner layer: y=westing, x=southing
+
Unfortunately westing/northing cannot be implemented at least in GRASS (in
any GIS?) as north must be up etc. For the GIS data we have, next is needed
(for GIS maps with negative coordinate values):
+y


21
x +> +x
34
GIS 
maps 

y
So the GIS maps are located in the 3rd quadrant.
A solution could be to optionally multiply/divide internally with 1:
cvs diff u PJ_krovak.c
Index: PJ_krovak.c
===================================================================
RCS file: /cvsroot/osrs/proj/src/PJ_krovak.c,v
retrieving revision 1.1
diff u r1.1 PJ_krovak.c
 PJ_krovak.c 16 Apr 2002 12:16:07 0000 1.1
+++ PJ_krovak.c 27 May 2002 13:59:53 0000
@@ 100,8 +100,9 @@
ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n)
;
/* x and y are reverted! */
 xy.y = ro * cos(eps) / a;
 xy.x = ro * sin(eps) / a;
+ /* when P>k0 == 1, then negative coordinates can be used */
+ xy.y = P>k0 * ro * cos(eps) / a;
+ xy.x = P>k0 * ro * sin(eps) / a;
#ifdef DEBUG
strcpy(errmess,"a: ");
@@ 162,9 +163,10 @@
/* Transformation */
/* revert y, x*/
+ /* when P>k0 == 1, then negative coordinates can be used */
xy0=xy.x;
 xy.x=xy.y;
 xy.y=xy0;
+ xy.x=xy.y / P>k0;
+ xy.y=xy0 / P>k0;
ro = sqrt(xy.x * xy.x + xy.y * xy.y);
eps = atan2(xy.y, xy.x);
For GRASS we have temporarily implemented "Kavork" (read backwards) which
is Krovak with switched coordinate axes. The new Kavork support negative
values, eg:
GRASS:~ > g.region pd
projection: 99 (Kavork)
zone: 0
datum: hermannskogel
ellipsoid: bessel
north: 1029600
south: 1040900
west: 629200
east: 618300
nsres: 10
ewres: 10
rows: 1130
cols: 1090
GRASS:~ > g.region l
long: 15.97793 lat: 50.33660 (north/west corner)
long: 16.13006 lat: 50.34784 (north/east corner)
long: 16.14807 lat: 50.24690 (south/east corner)
long: 15.99626 lat: 50.23569 (south/west corner)
rows: 1130
cols: 1090
Center Longitude: 16:03:47.083724E [16.06308]
Center latitude: 50:17:30.334822N [50.29176]
which are reasonable results (error because of missing datum
transformation). I also added datum "hermannskogel" to the PROJ4 version
in GRASS (although GRASS does not yet fully support datum transformations,
it wants to see the datum entry already).
So  what is Krovak? It must be well defined, either from east to west or
from west to east and exactly in this way it must be defined in PROJ4. We do
not have a solution now  should we implement two Krovak versions?
Alternate method might to utilize +k=1.0 to use above patch. But at time
+k=value must be value positive only.
Looking at other GIS:
MicroStation:
http://selectservices.bentley.com/technotes/faqs/6192.htm#8
ESRI:
http://arconline.esri.com/arconline/documentation/ims_/webhelp31/elements/pcs_list.htm
...they also implement several Krovak variations.
Hoping for a suggestion,
Markus Neteler

Markus Neteler
ITCirst, 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.
More information about the Proj
mailing list