[Proj] confused from krovak

Tue Mar 21 10:16:33 PST 2006

Hallo,
I tested the whole problem one more time and I would say now, that I did
not compile  PROJ on MS Windows with the patch corectly. I setuped the
same configuration on Linux and everything worked well.

Here is the part of the mapfile:

 EXTENT  -603055.584603 -1163534.631291 -578334.193874 -1138794.167715
 PROJECTION
	"proj=krovak"
	"ellps=bessel"
	"a=6377397.155"
	"es=0.0066746722"
	"f=299.1528128"
	"towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56"
 END

And the resulting GetCapabilities from Linux server with applayed patch
on PJ_krovak.c from Radim Blazek

[...]
 <LatLonBoundingBox minx="16.5098" miny="49.1661" maxx="16.8839" maxy="49.411" />
[...]

Which is what I would expect.

The same input but different output I get from Proj and MapServer
installed on MS Windows server:
[...]
 <LatLonBoundingBox minx="37.9432" miny="69.5973" maxx="38.6756" maxy="69.8515" />
[...]

The way, how I try to apply the patch on MS Windows is following:

cygwin: > cd /cygdrive/c/
cygwin: > rm -rf proj
cygwin: > wget ftp://ftp.remotesensing.org/proj/proj-4.4.9.zip
cygwin: > unzip proj-4.4.9.zip
cygwin: > mv proj-4.4.9 proj
cygwin: > cd proj/src/
cygwin: > cat PJ_krovak.c |grep "\(xy\.y\)\|\(xy\.x\)"

        xy.y = ro * cos(eps) / a;
        xy.x = ro * sin(eps) / a;
        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);

cygwin: > wget http://mpa.itc.it/radim/jtsk/PJ_krovak.c.patch
cygwin: > patch PJ_krovak.c PJ_krovak.c.patch
cygwin: > cat PJ_krovak.c |grep "\(xy\.y\)\|\(xy\.x\)"

        xy.y = -ro * cos(eps) / a;
        xy.x = -ro * sin(eps) / a;
        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);

Right, now in Windows command line

wincmd: > cd c:\proj\src
wincmd: > vcvars32.bat
wincmd: > nmake /f makefile.vc all

Now the test:

wget -O - "http://indica.mendelu.cz/cgi-bin/mapserv.exe?map=c:\mapservdata\data\data_slp\krtiny_krovak.map&version=1.1.1&service=WMS&request=GetCapabilities" | grep LatLonBoundingBox

<LatLonBoundingBox minx="37.9432" miny="69.5973" maxx="38.6756" maxy="69.8515" />

and the same with linux server

wget -O - "http://mapserver.mendelu.cz/wms/?&service=WMS&request=GetCapabilities" | grep LatLonBoundingBox

<LatLonBoundingBox minx="16.5098" miny="49.1661" maxx="16.8839" maxy="49.411" />

Could somebody tell me, what I am doing wrong?

Thank you

Jachym

On Fri, Mar 17, 2006 at 11:00:28AM +0100, Oscar van Vlijmen wrote:
> > From: Jachym Cepicky
> > Subject: [Proj] confused from krovak
> 
> > My data are stored in Krovak projection (known as S-JTSK). All the data
> > are in "GIS version of Krovak", which means for the axes, that X=-y and
> > Y=-x. I used patch of PJ_krovak.c from Radim Blazek's site.
> > Output of cs2cs gives me expected result:
> > echo "-603056 -1.16353e+006" | cs2cs +proj=krovak +ellps=bessel
> > +a=6377397.1550000003 +es=0.0066743722 +f=299.1528128000
> > +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to +proj=latlong
> > +datum=WGS84 +ellps=WGS84 +no_defs
> > 16d32'43.318"E  49d9'55.923"N 44.706
> 
> I can confirm this:
> With my krovak function (not derived from (lib)proj) and the input
> x=1163530; y=603056; (note the +)
> I arrive at lat, lon wrt Greenwich:
> lat=49.1661377; lon=16.5466812 decimal degrees
> which after datum transform with your parameters gives:
> lat=49d 09m 55.9232s; lon=16d 32m 43.3178s; h=44.7065 m
> 
> > ....
> > To use EPSG code seems not to be working -- Krovak is defined with
> > <2065> +proj=krovak +lat_0=49.5 +lon_0=60.16666666666667
> > +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro
> > +units=m +no_defs
> 
> In the mean time EPSG has updated the parameters (alpha for 0.1
> milliarcsec), but this won't make any difference.
> 
> > and using this parameters with cs2cs gives
> > echo "-603056 -1.16353e+006" | cs2cs +proj=krovak +lat_0=49.5
> > +lon_0=60.16666666666667 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0
> > +ellps=bessel +pm=ferro +units=m +no_defs +to +proj=latlong +datum=WGS84
> > +ellps=WGS84 +no_defs
> > 34d12'48.052"E  49d10'0.227"N -701.834
> > which is definitely bad result
> 
> Remains to be seen. With my 'krovak' set to Ferro as prime meridian I
> obtain:
> lat=49d 09m 58.0958s; lon=34d 12m 48.0524s
> This is still on the Bessel ellipsoid!!!
> For your excessive height (depth) I have no explanation. Your latitude
> differs by 2" with mine, which is not very significant.
> Probably this information will help you to cure the cs2cs command line?
> 
> 
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj

-- 
Jachym Cepicky                                  
e-mail: jachym.cepicky at centrum.cz           
URL: http://les-ejk.cz                      
GPG: http://les-ejk.cz/gnupg_public_key/    
-----------------------------------------   
OFFICE:                                     
Department of Geoinformation Technologies   
LDF MZLU v Brně                             
Zemědělská 3                                
613 00 Brno                                 
e-mail: xcepicky at node.mendelu.cz            
URL:    http://mapserver.mendelu.cz
Tel.:   +420 545 134 514