[GRASSLIST:5954] Re: Krovak / PROJ4 (+axis option?)

Radim Blazek wrote:

Hi Frank,

as you probably know, Krovak defines X=-y, Y=-x, i.e. X is vertical axis, oriented down and Y is horizontal oriented to the left. But all
GIS/GPS/... systems must use usual coordinate system where coordinates are in 3. quadrant.
For example:
14.39862613611E 50.465671925N LL
X=1001120.05 Y=738666.75 Krovak
x=-738666.75 y=-1001120.05 In usual coor system and all GIS

Because:
1) Evenden's implementation of Krovak (PJ_kocc.c) returns by default
   values in 3. quadrant (-738666.75 -1001120.05) and Krovak X,Y
   is returned only with +czech option
2) All GIS applications require coor in 3. quadrant and there is no reasonable way to switch axis outside PROJ4 only for Krovak
   (For GRASS for example, I used code PJ_krovak.c from PROJ4
    but modified the code to get negative signs.)
3) Current behaviour in PROJ4 is wrong anyway in my opinion, because
   returns 738666.76,1001120.13 (switched signs but not axis)

I suggest to change PJ_krovak.c in PROJ.4 as you can see below.
That means to get values that most of applications require.

Radim

136,137c139,140
< xy.y = ro * cos(eps) / a;
< xy.x = ro * sin(eps) / a;
---

     xy.y = -ro * cos(eps) / a;
     xy.x = -ro * sin(eps) / a;

199,200c202,203
< xy.x=xy.y;
< xy.y=xy0;
---

     xy.x = -xy.y;
     xy.y = -xy0;

Radim,

Rather than have a +czech switch specifically for the Krovak projections, and
instead of altering the sign conventions from the "usual", I would like to
introduce a generic mechanism at the pj_transform() level for alternate axis
orientations.

Something along the line of a +axis switch, where "+axis=en" means that
x axis = 'e'ast, and y axis = 'n'orth. Then a coordinate system string
could specify alternate axis orientations as +axis=sw indicating that
X is a "southing" (or negative northing) and Y is a westing (or negative
easting). I think this would be the right organization for Krovak right?

My thinking is that this isn't a "projection" issue, but rather at a higher
level a "coordinate system" issue which is why I think it should be handled
by pj_transform() and not effect the low level projections routines.

I would add that Krovak is not the only coordinate system with unusual
axis conventions but this is an area I have been "sweeping under the carpet"
in most of my work in PROJ.4, and OGRSpatialReference, so far.

I would also note that I intend to adopt Gerald Krovak code at some point.

I have no schedule for implementing generic axis support or migrating to
use Geralds Krovak (or simplified libproj). I will not have time to even
look into this for several weeks. Thus, for your own purposes you should hack
your local GRASS built in a manner suitable to your needs in the meantime.

I have taken the liberty of cc:ing this reply to the PROJ list where I hope
any discussion of a +axis switch can take place. GRASS folks might also be
interested as it would presumably be the way that alternate axis orientations
are selected in GRASS.

Best regards,

--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

> as you probably know, Krovak defines X=-y, Y=-x, i.e. X is vertical
> axis, oriented down and Y is horizontal oriented to the left. But
> all GIS/GPS/... systems must use usual coordinate system where
> coordinates are in 3. quadrant.

[...]

>
> Radim
>

Something along the line of a +axis switch, where "+axis=en" means
that x axis = 'e'ast, and y axis = 'n'orth. Then a coordinate system
string could specify alternate axis orientations as +axis=sw
indicating that X is a "southing" (or negative northing) and Y is a
westing (or negative easting). I think this would be the right
organization for Krovak right?

Hi - just to point out how Matlab handles this issue, which always made
good sense to me when I have needed to do this sort of thing on
occasion:

    AXIS IJ puts MATLAB into its "matrix" axes mode. The coordinate
       system origin is at the upper left corner. The i axis is
       vertical and is numbered from top to bottom. The j axis is
       horizontal and is numbered from left to right.
    AXIS XY puts MATLAB into its default "Cartesian" axes mode. The
       coordinate system origin is at the lower left corner. The x
       axis is horizontal and is numbered from left to right. The y
       axis is vertical and is numbered from bottom to top.
[(c) Mathworks 199something I suppose]

I find thinking of it this way requires less mental gymnastics than
trying to visualize rotated planes, which in turn makes debugging a lot
more straight forward.

regards,
Hamish