[FreeGeoCZ] Krovak v GRASSu podporovan

["Jachym Cepicky" <jachym.cepicky AT gmail.com>]

Zdravím,
v Trentu právě proběhla jedna z mnoha hacking sessions a zdá se, že
souřadný systém S-JTSK je od nynějška kompletně podporovaný v GRASSu.

Obrázek v příloze, nějaké výpisy zde

GRASS6.3.cvs krovak_z_EPSG/jachym: grass6 > g.region -p
projection: 99 (Krovak)
zone:       0
datum:      hermannskogel
ellipsoid:  bessel
north:      1
south:      0
west:       0
east:       1
nsres:      1
ewres:      1
rows:       1
cols:       1
cells:      1
GRASS6.3.cvs krovak_z_EPSG/jachym: grass6 > g.proj -p
-PROJ_INFO-------------------------------------------------
name       : Krovak
proj       : krovak
datum      : hermannskogel
ellps      : bessel
lat_0      : 49.5
lon_0      : 42.5
alpha      : 30.28813972222222
k          : 0.9999
x_0        : 0
y_0        : 0
pm         : ferro
no_defs    : defined
towgs84    : 570.8,85.7,462.8,4.998,1.587,5.261,3.56
-PROJ_UNITS------------------------------------------------
unit       : metre
units      : metres
meters     : 1
GRASS6.3.cvs krovak_z_EPSG/jachym: grass6 >

zbývá rozjet r.in.wms  ;-)

Jáchym

---------- Forwarded message ----------
From: Markus Neteler <neteler AT itc.it>
Date: 6.3.2007 18:45
Subject: Re: GRASS: prime meridian recognition?
To: Martin Landa <landa AT itc.it>
Cc: Paul Kelly <paul-grass AT stjohnspoint.co.uk>, Jachym Cepicky 
<cepicky AT itc.it>


Hi,

I have added support for the Krovak datum to get queried
once the user creates a location from EPSG code (at least
is solved now).

If we reproject CZ from LatLong, the coordinates are right:
g.region vect=czech_republic  -p
projection: 99 (Krovak)
zone:       0
datum:      hermannskogel
ellipsoid:  bessel
north:      -934850.34795202
south:      -1227231.04178299
west:       -904253.20215554
east:       -431692.75129609

Also WKT is right:
g.proj -w
PROJCS["Krovak",
   GEOGCS["bessel",
       DATUM["Militar_Geographische_Institut",
           SPHEROID["Bessel_1841",6377397.155,299.1528128],
           TOWGS84[570.8,85.7,462.8,4.998,1.587,5.261,3.56]],
       PRIMEM["ferro",-17.666666666668],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Krovak"],
   PARAMETER["latitude_of_center",49.5],
   PARAMETER["longitude_of_center",24.833333333332],
   PARAMETER["azimuth",0],
   PARAMETER["pseudo_standard_parallel_1",0],
   PARAMETER["scale_factor",0.9999],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   UNIT["metre",1]]

We still observe the problem in 'g.region -l' but this
seems to be a bug in g.region itself (say, the underlying
G_proj_xx() functions).

Best,
Markus

On Tue, Mar 06, 2007 at 06:14:09PM +0100, Martin Landa wrote:
> Hi,
>
> On Tue, Mar 06, 2007 at 05:12:34PM +0100, Markus Neteler wrote:
>
> > On Tue, Mar 06, 2007 at 02:39:38PM +0000, Paul Kelly wrote:
> > > Hello Markus,
> > > This is an interesting problem. I vaguely remember something like this
> > > before not being resolved, and also that the prime meridian support in
> > > PROJ is a bit confused.
> > >
> > > The key parameter is the longitude of origin, which is given in the WKT
> > > file as 24.833 relative to the Ferro meridian (which is -17.667), which
> > > makes it +7.166 relative to the Greenwich meridian, right?
> >
> > Honestly, I don't know. Martin may have an idea.
>
> the longitude of origin used in Krovak projection is 42d30m relative
> to Ferro meridian. I have tried to create GRASS location based on EPSG
> 2065.
>
> g.proj -p
> -PROJ_INFO-------------------------------------------------
> name       : Krovak
> proj       : krovak
> ellps      : bessel
> lat_0      : 49.5
> lon_0      : 42.5                                                <-- OK 
> relative to Ferro !
> alpha      : 30.28813972222222
> k          : 0.9999
> x_0        : 0
> y_0        : 0
> pm         : ferro                                               <-- OK
> no_defs    : defined
> -PROJ_UNITS------------------------------------------------
> unit       : metre
> units      : metres
> meters     : 1
>
> g.proj -w
> PROJCS["Krovak",
>     GEOGCS["bessel",
>         DATUM["unknown",
>             SPHEROID["Bessel_1841",6377397.155,299.1528128]],
>         PRIMEM["ferro",-17.666666666668],                       <-- still 
> Ferro
>         UNIT["degree",0.0174532925199433]],
>     PROJECTION["Krovak"],
>     PARAMETER["latitude_of_center",49.5],
>     PARAMETER["longitude_of_center",24.833333333332],           <-- relative 
> to Greenwich
>     PARAMETER["azimuth",0],
>     PARAMETER["pseudo_standard_parallel_1",0],
>     PARAMETER["scale_factor",0.9999],
>     PARAMETER["false_easting",0],
>     PARAMETER["false_northing",0],
>     UNIT["metre",1]]
>
> Reprojecting data from/to the location (e.g. Krovak <-> WGS84)
> seems to work on my machine.
>
> Running g.region -l I get (wrong values)
>
> g.region -l
> north-west corner: long: 29:36:34.190968E lat: 50:49:22.537224N
> north-east corner: long: 36:18:19.791042E lat: 51:21:12.834112N
> south-east corner: long: 36:37:33.383801E lat: 48:44:13.412399N
> south-west corner: long: 30:16:25.239095E lat: 48:14:00.914015N
> center longitude:  33:12:13.151227E
> center latitude:   49:47:12.424438N
>
> instead of
>
> north-west corner: long: 11:56:36.977797E lat: 50:49:27.846078N
> north-east corner: long: 18:38:25.8542E lat: 51:21:18.085148N
> south-east corner: long: 18:57:39.307535E lat: 48:44:17.405995N
> south-west corner: long: 12:36:28.254847E lat: 48:14:04.96853N
> center longitude:  15:32:17.598595E
> center latitude:   49:47:17.076438N
>
> I hope it is less confusing then before;-)
>
> Best regards, Martin
>
> > >  if I run g.proj -p yeps=2065 I see the lon_0 value given as 42.5. This
> > > is actually what you would get if you subtracted (-17.667) from 24.833
> > > instead of adding them like I did above, but I can't see how that's
> > > correct. It's confusing me. +42.5 E is I suspect much closer to the
> > > Ukraine than to the Czech Republic.
> >
> > That's true and in fact confusing.
> >
> > > My hypothesis for the time being is that this value is being converted
> > > wrongly. If you manually edit the PROJ_INFO so that lon_0 is 7.16667
> > > rather than 42.5, does it reproject correctly? Can you try it with and
> > > without the
> > > pm: ferro
> > > line (I suspect it is ignored by PROJ but not very sure).
> >
> > I tried to no avail. Say, the 15 deg appear but then the CZ boundary
> > coordinates are looking wrong.
> >
> > > The results of those tests (if they make it clear what is the correct
> > > value) might help us to investigate a bit further where the conversion is
> > > going wrong.
> >
> > For playing around I have extracted CZ from VMAP0:
> > http://mpa.itc.it/markus/tmp/czech_republic_ll_wgs84.tar.gz
> >
> > Hopefully Martin picks it up and sends us a comment :-)
> > (he and Jachym are sitting 1.2meters away in our office)
> >
> > more soon,
> > Markus
> >
> > > Paul
> > >
> > > On Tue, 6 Mar 2007, Markus Neteler wrote:
> > >
> > > >Hello Paul,
> > > >
> > > >for a project I have to deal with Krovak projection:
> > > >
> > > >EPSG: 2065
> > > >
> > > >g.proj -w
> > > >PROJCS["Krovak",
> > > >   GEOGCS["bessel",
> > > >       DATUM["unknown",
> > > >           SPHEROID["Bessel_1841",6377397.155,299.1528128]],
> > > >       PRIMEM["ferro",-17.666666666668],
> > > >       UNIT["degree",0.0174532925199433]],
> > > >   PROJECTION["Krovak"],
> > > >   PARAMETER["latitude_of_center",49.5],
> > > >   PARAMETER["longitude_of_center",24.833333333332],
> > > >   PARAMETER["azimuth",0],
> > > >   PARAMETER["pseudo_standard_parallel_1",0],
> > > >   PARAMETER["scale_factor",0.9999],
> > > >   PARAMETER["false_easting",0],
> > > >   PARAMETER["false_northing",0],
> > > >   UNIT["metre",1]]
> > > >
> > > >GRASS seems to define things correctly. When I
> > > >reproject the Czech republic into the new location,
> > > >zoom to the map, I end up in Ukraine!
> > > >
> > > >g.region -lg
> > > >NV_long=29.60814214
> > > >No_lat=50.82631427
> > > >ENE_long=36.30134121
> > > >Be_lat=51.35684619
> > > >Ce_long=36.62241227
> > > >se_lat=48.73697405
> > > >Sc_long=30.27325184
> > > >sw_lat=48.23363642
> > > >center_long=33.20128687
> > > >center_lat=49.78844273
> > > >rows=474
> > > >cols=766
> > > >
> > > >But if I calculate
> > > >>33.20128687 - 17.666666666668
> > > >[1] 15.53462
> > > >
> > > >I am back in Czech republic. It seems that the prime
> > > >meridian is recorded but not applied.
> > > >
> > > >Do you have a suggestion how to make use of the PM?
> > > >
> > > >thanks,
> > > >Markus
> > > >
> > > >------------------
> > > >ITC -> hall'1 marzo 2007 Fondazione Bruno Kessler
> > > >ITC -> since 1 March 2007 Fondazione Bruno Kessler
> > > >------------------
> > > >
> > > >
> >
> > --
> > Markus Neteler  <neteler itc it>  http://mpa.itc.it/markus/
> > ITC-irst -  Centro per la Ricerca Scientifica e Tecnologica
> > MPBA - Predictive Models for Biol. & Environ. Data Analysis
> > Via Sommarive, 18        -       38050 Povo (Trento), Italy
> --
> Martin Landa  <landa itc it>
> ITC-irst -  Centro per la Ricerca Scientifica e Tecnologica
> MPBA - Predictive Models for Biol. & Environ. Data Analysis
> Via Sommarive, 18        -       38050 Povo (Trento), Italy
>
>  LocalWords:  reproject VMAP se sw ITC marzo Fondazione Kessler neteler

------------------
ITC -> dall'1 marzo 2007 Fondazione Bruno Kessler
ITC -> since 1 March 2007 Fondazione Bruno Kessler
------------------



--
Jachym Cepicky
e-mail: jachym.cepicky gmail com
URL: http://les-ejk.cz
GPG: http://www.les-ejk.cz/pgp/jachym_cepicky-gpg.pub

Attachment: krovak1.png
Description: PNG image