[OSRS-PROJ] Krovak and GIS usage

[Markus Neteler]

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
               |
               |
              2|1
   -x ---------+------> +x
              3|4
        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

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.