S-JTSK (zastaralý článek)

Zdroje:

Radim Blažek: GIS GRASS a JTSK (Krovak)
Jan Ježek: Vývoj programového modulu pro převod souřadnic mezi kartografickými zobrazeními], diplomová práce, 2003, ČVUT v Praze, Fakulta stavební

!!! Nova transformace v PROJ.4 pomoci GRIDu je k dispozici na: S-JTSK-Grid

Obsah

1 Křovákovo zobrazení
2 Knihovna PROJ.4
3 Gdalwarp a ogr2ogr
4 GRASS lokace S-JTSK
- 4.1 Založení nové GRASS lokace
  - 4.1.1 wxGUI - Location wizard
  - 4.1.2 Textové rozhraní
5 QGIS a projekt v S-JTSK
6 Související články

Křovákovo zobrazení

Souřadnicový systém S-JTSK je mimojiné definován vlastním kartografickým zobrazením – tzv. Křovákovým zobrazením.

Křovákovo zobrazení je konformní kuželové zobrazení v obecné poloze, které v roce 1922 navrhl Ing. Josef Křovák. Transformace souřadnic <math>\varphi, \lambda</math> na pravoúhlé <math>X, Y\,</math> se provádí v několika krocích. Nejprve je provedeno Gaussovo konformní zobrazení Besselova elipsoidu na kouli a poté konformní zobrazení na kuželovou plochu obecně položenou.

Schéma zobrazení:

<math>(\varphi, \lambda) \rightarrow (U, V) \rightarrow (\check{S}, D) \rightarrow (\rho, \epsilon) \rightarrow (X, Y)</math>

V první fázi je zobrazen Besselův elipsoid na kouli. Vzhledem k tomu, že tento krok předchází zobrazení na kuželovou plochu, nazývá se Křovákovo zobrazení dvojité. Podmínkou zobrazení je minimální délkové zkreslení kolem základní rovnoběžky, která byla zvolena <math>\varphi_{o}=49^{o}30'</math> (Pro <math>\varphi_{o}</math> je hodnota délkového zkreslení <math>m_{o}=1\,</math>).

Křovák zvolil za zobrazovací plochu kužel v obecné poloze. Empiricky zjistil, že ČSR lze ohraničit horizontálními kružnicemi do pásu o šířce <math>2^{o}31'\,</math>, přičemž zkreslení na okrajích dosahovalo hodnot +24cm/km. V případě normální polohy kužele by byl pás široký <math>3^{o}20'\,</math> a zkreslení v maximální vzdálenosti od střední rovnoběžky až +42 cm/km. Při definitivní úpravě byla zvolena jako základní rovnoběžka kartografická rovnoběžka <math> \check{S}_{o}=78^{o}30'</math> a okrajové <math>\check{S}_{1}=77^{o}13'</math> a <math>\check{S}_{2}=79^{o}44'</math>. Základní rovnoběžka je kolmá na zeměpisný poledník <math>\lambda=42^{o}30'\,</math> východně od Ferra a jejich průsečík <math>A</math> má šířku <math>\varphi=48^{o}15'</math>.

Referenční koule je dále konformně zobrazena na kužel v obecné poloze.

Křovák umístil osu X do obrazu základního poledníku <math>\lambda=42^{o}30'\,</math> v.F. a počátek souřadnic do vrcholu kužele Q. Tím byla celá ČSR umístěna do jediného kvadrantu (viz obrázek níže).

Křovákovo zobrazení

Knihovna PROJ.4

Práci s kartografickými zobrazeními má (nejen) v GRASSu na starost knihovna PROJ.4. Křovákovo zobrazení v PROJ.4 odpovídá jeho definici, nicméně v této podobě není použitelné v GIS produktech (x, y). Aby bylo možné použít Křovákovo zobrazení z PROJ.4 v GRASSu (tj. pracovat ve 3. kvadrantu) je třeba aplikovat patch na soubor PJ_krovak.c.

Na Radimův požadavek aplikovat tento patch na PROJ.4 odpověděl Frank Warmerdam takto:

... With all respect to the Czech nation, I don't feel the projection is widely
used enough to cause a lot of people problems.

Na Veliký pátek 2005, požádal Radim znovu Franka Warmerdama o začlenění patche do oficialní distribuce PROJ.4, v jeho odpovědi věnujte pozornost zejména poslednímu odstavci a pokud potřebujete "KrovakGIS", přidávejte poznámky do bugzilly.

Velkopáteční odpověď Franka Warmerdama

On Fri, 25 Mar 2005 09:03:57 +0100, Radim Blazek <blazek@itc.it> wrote:

>> Last time you refused the patch with
>> 
>>    "... With all respect to the Czech nation, I don't feel
>>     the projection is widely used enough to cause alot of
>>     people problems."

Radim, 

I am surprised the above was my primary objection.  I would
have thought my primary objection would have been that I 
intend to do a "proper" approach to CRS'es with axis orientations.
However, I have clearly *not* done that yet.  I would have expected
the above objection would instead have been an excuse for why
I thought I could leave Krovak "out in the cold" while I worked on
such a solution. 

Looking at your patch, it just changes the sign of the x and y by 
the looks of it. I had thought the x and y axes were also swapped
in Krovak.  I guess I am mistaken on that.  Or perhaps that is already
taken into account? 

>> Do we have to organize a petition to convince you? I thought that open
>> software development respects reasonable user needs. Am I wrong?

A petition is a questionable tactic since it is hard for all the
folks who might start running into problems when I apply an
ill-considered patch to get organized and counter-petition in
advance of the change actually being made. 

As for whether open software development respects reasonable
user needs, that is not at all clear.  What is "open software 
development"?  OSS?  FOSS?  I certainly provide public access
to PROJ.4's development status, and provide mechanisms (bugzilla
and the mailing list) for users to provide feedback.  But it is not a
democratic environment.  One, or a few oligarchs control the what
actually happens.  

I like to think I am open to feedback and it is rare that I don't
apply a patch unless I am concerned that it will have negative
consequences.  In this case I am concerned, though my concern
is at best ill-formed.  It is a vague sense of ill-ease that I am going
to be sorry for this "hard coded" solution when I actually try to
implement real axis orientation support.  Also, that I will have
broken things for a bunch of people already using Krovak support
successfully.  Or that I am breaking compatibility with ESRI or
existing EPSG formulations or something. 

Has Gerald patched libproj4 according to your needs?  One 
reasonably compelling argument would be that I might as well
come into line with what he does since I will eventually abandon
all the low level projection code in favor of that in libproj4. 

Feel free to add notes to the bugzilla bug, raise it's priority and
harrass me every week or two.  I will try and resolve this before
the next public release of PROJ.4. 

Best regards,

On Tue, 29 Mar 2005 10:50:42 +0200, Radim Blazek <blazek@itc.it> wrote:

>> Bunch of people?! Nobody can use it successfully in GIS.
>> In theory it is possible to use it in application which is aware of
>> Krovak axis orientation and labels, but it is not because current
>> implementation is wrong (labels only switched).

Radim,

Well, it must work for someone.  Someone submitted the code, and
claimed it worked.  

>>> > Has Gerald patched libproj4 according to your needs?  One
>>> > reasonably compelling argument would be that I might as well
>>> > come into line with what he does since I will eventually abandon
>>> > all the low level projection code in favor of that in libproj4.
>
>> 
>> I don't know how Krovak is implemented in his version, I have never got
>> it compiled.

I would suggest you submit the patch to him.  After he has incorporated
it, I will re-clone my krovak code from his.  Otherwise things will just
change again when my low level projections code is discarded in 
favor of his. 

If you can convince Gerald that is good enough for me!

Best regards,

On Wed, 30 Mar 2005 09:04:59 +0200, Radim Blazek <blazek@itc.it> wrote:

>> It seems that he used '+czech' switch:
>> 
>> Gerald I. Evenden wrote:
>>  > <S-JTSK> +proj=kocc +ellps=bessel +czech
>>  >      +lon_0=24d50 +lat_0=49d30
>>  >      +lat_t=59d42'42.69689 +lat_1=78d30 +k_0=.9999
>>  >
>>  > To generate Czech grid system values the flag option +czech
>>  > must be used to shift math-x to Y and math-y to X, changing
>>  > the sign to positive in Czech area.  That is plus y west of
>>  > Helsinki and plus x south of Helsinki.
>> 
>> As far as I understand, without '+czech' it gives 'math' axis (necessary
>> for any GIS), i.e. the same as current PROJ.4 + my patch. With 'czech+'
>> option it gives X = -y, Y = -x, i.e. current PROJ.4 + switched axis.

Radim,

So is that OK with you?  The default would be the same as with your
patch, right?  If that is good for you, I will adopt Gerald's current 
code and we won't have to worry about future proofing things
any further.   We will deal with any backward compatibility issues,
but as you note those are perhaps entirely hypothetical. 

Best regards,

PROJ mailing list / březen 2006

Jachym Cepicky jachym.cepicky at centrum.cz
Wed Mar 15 18:08:54 EST 2006

Hallo, 
after successful compilation of proj on ms windows (thank you, Mateusz),
I try to get "expected" results from UMN MapServer.

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

this corresponds more or less with the output from GRASS GIS. Parameters
for projection informations are taken from PROJ_INFO file from my GRASS
locations.

Now I would like to share some data via UMN MapServer as WMS. What is
the right way, how to define the projection "krovak" in Mapfile when I
would like view the data also in other projections?

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

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

So I have 
    
    EXTENT -603055.58 -1163534.63 -578334.19 -1138794.16
    WEB
        "wms_srs" "EPSG:4326"
    END
    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
    
In my mapfile

However GetCapabilities gives me strange result:
     <LatLonBoundingBox minx="37.9432" miny="69.5973" maxx="38.6756" maxy="69.8515" />

Using GRASSs  g.region w=-603055.58 s=-1163534.63 e=-578334.19 n=-1138794.16 -gb

returns

n=-1138794.16
s=-1163534.63
w=-603055.58
e=-578334.19
nsres=4.96696848
ewres=4.96712678
rows=4981
cols=4977
LL_W=16.50846607
LL_E=16.88255667
LL_N=49.41034697
LL_S=49.16549325

So what to do? How to setup my mapfile correctly? Why is actualy version
of Krovak in EPSG as it is (2065)? I know lot of people and companyes, who are having their
maps in "Krovakgis" projection (-x,-y), but none, who is using some
S-JTSK - Ferro X,Y version. This significantly complicates usage of e.g.
grass, qgis and other programs.

Results can be viewed on 

GetCapabilities: http://indica.mendelu.cz/cgi-bin/mapserv.exe?map=c:\mapservdata\data\data_sl
\krtiny.map&version=1.1.1&service=WMS&request=GetCapabilities
GetMap (default): http://indica.mendelu.cz/cgi-bin/mapserv.exe?map=c:\mapservdata\data\data_slp
\krtiny.map&version=1.1.1&service=WMS&request=GetMap&layers=landuse&format=image/png

Proj compilation and installation:

windows: > cp proj-4.9.src c:\proj
windows: > cd c:\proj\src
cygwin:  > patch PJ_krovak.c PJ_krovak.c.patch
windows: > nmake /f makefile.vc all
cygwin:  > cd /cygdrive/c/proj/
cygwin:  > ./configure --prefix=/cygdrive/c/proj
cygwin:  > make
cygwin:  > make install

Thank you for any reply 

Jachym
-- 
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

Oscar van Vlijmen ovv at hetnet.nl
Fri Mar 17 05:00:28 EST 2006

 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?

Clifford J Mugnier cjmce at lsu.edu
Fri Mar 17 09:10:03 EST 2006

The comment below regarding transformed ellipsoid height is getting quite
far from the general topic of cartography, GIS, and map projections.  When
concern is voiced over ellipsoid heights in combined datum shift and
projection transformations, the user is now deep into geodesy.

In fact, the topic is termed PHYSICAL Geodesy, and the appropriate
application of "GEOIDS" is in order.

The Bessel ellipsoid was cooked up in 1841.  The Krovak projection was
cooked up in 1922.  The Czech Republic is still referenced to the
Hermannskoegel Datum of 1871 when using the Krovak projection.  Note that
these are systems with origins from two centuries ago!  Trying to reconcile
these ancient systems with a GPS system is mostly a waste of time with
respect to ellipsoid heights.  That old datum does NOT have a geoid
associated with it.  It does not even have an old-timey astrogeodetic
deflection geoid computed muchless a geoid that could be comparable to the
EGM96 geoid.

The PROJ4 collection of cartographic and geometric utilities is not a
panacea to all problems geodetic!  If the Czech Republic is prepared to
sink several hundred million dollars into defining a geoid for its old
datum, then ellipsoid heights will make then make (more) sense to the
neophyte GIS user.  Until that happens, don't expect miracles from PROJ4.
You will need a Hermannskoegel 1871 Datum Geoid for that to happen, and
that's not going to happen.  It would be a waste of money.

I believe the Czech Republic government is working on a local national
geoid to be compatible with a new inertial datum that is associated with
the International Terrestrial Reference Frame (ITRF); a proper use of
public money.

Clifford J. Mugnier, C.P., C.M.S.
Chief of Geodesy,
CENTER FOR GEOINFORMATICS
Department of Civil Engineering
CEBA 3223A
LOUISIANA STATE UNIVERSITY
Baton Rouge, LA  70803
Voice and Facsimile:  (225) 578-8536 [Academic]
Voice and Facsimile:  (225) 578-4474 [Research]
======================================================
http://www.asprs.org/resources/GRIDS/
http://www.cee.lsu.edu/facultyStaff/mugnier/index.html
======================================================



> 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?

Jachym Cepicky jachym.cepicky at centrum.cz
Tue Mar 21 13:16:33 EST 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

Reakce Franka Warmerdama k vydání PROJ.4 4.5.0 / duben 2006

Maciek Sieczka wrote:
> On Tue, 25 Apr 2006 10:48:29 +0200
> "Martin Landa" <landa.martin@gmail.com> wrote:
>
>> Dear PROJ.4 developers,
>>
>> trying to build proj-4.5.0b2 I found out, that there is nothing new
>> with Krovak projection. As I know Radim's patch [1] has not been
>> implemented yet...
>>
>> Please, see bugreports
>>
>> http://bugzilla.remotesensing.org/show_bug.cgi?id=1133
>>
>> and
>>
>> http://bugzilla.remotesensing.org/show_bug.cgi?id=147
>>
>> I am sorry for disturbing you but it is really *frustrating* for Czech
>> and Slovak users
>
> + 1 Polish who happens to do some GIS for Krkonose :)

Martin / Maciek,

I hear your pain.  It was my intention not to do another PROJ release
till I dealt with this issue, but then pressure on the meridian issues
forced by hand.   My intent is to resolve this by adding arbitrary
axis orientation flags for the coordinate system definition rather than
have an assortment of variations of projection names for different
orientations.

Best regards,

Další zmínky o Křovákově zobrazení

PROJ.4 mailing list 2003/02 – "Horrors of --- Krovak projection"

Bugreport

Pokud tedy chcete přispět k vyřešení problému, přidávejte své poznámky k bugreportu 1133 nebo 147.

Jak aplikovat patch

Pracovat musíte se zdrojovými texty knihovny PROJ.4.

cd /cesta/ke/zdroji/proj.4/
cd src/

Pro verzi 4.4.9 a nižší,

wget http://mpa.itc.it/radim/jtsk/PJ_krovak.c.patch

verze 4.5.0

wget http://gama.fsv.cvut.cz/~landa/proj4/s-jtsk/PJ_krovak.c.patch

nebo s přepínačem +czech

wget http://gama.fsv.cvut.cz/~landa/proj4/s-jtsk/PJ_krovak_czech.c.patch

patch PJ_krovak.c PJ_krovak.c.patch

nebo

 patch PJ_krovak.c PJ_krovak_czech.c.patch

cd ..
./configure
make
make install

Výsledkem bude modifikované Křovákovo zobrazení použitelné v GIS.

PROJ.4 4.5.0 příběh s "dobrým" koncem

Ve verzi 4.5.0 je začleněn patch PJ_krovak_czech.c.patch!

echo "15 50" | proj +proj=krovak
-703105.69      -1058219.60

echo "15 50" | proj +proj=krovak +czech
703105.69       1058219.60

Praktická ukázka

Na tomto místě si prakticky demonstrujeme, jak pracovat s knihovnou PROJ.4 a souřadnicovým systémem S-JTSK.

Testovací bod:

WGS-84		S-JTSK
zem. délka	zem. šířka	Y	X
12d48'25.15992	49d27'8.14571	868208.52	1095793.96

Poznámka: v GIS pracujeme se zapornými souřadnicemi.

Manuální definice

WGS-84 → S-JTSK:

echo "12d48'25.15992  49d27'8.14571" | cs2cs +proj=latlong +datum=WGS84 \
+to +proj=krovak +ellps=bessel +no_defs \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

-868208.53      -1095793.57 -46.96

Diference: <math>\Delta Y = 1cm; \Delta X = 39cm</math>

EPSG

# S-JTSK (Ferro) / Krovak
<2065> +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999
+x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs  <>

LatLong (Bessel/Ferro) → S-JTSK:

echo "30d28'28.099  49d27'10.879" | cs2cs +proj=latlong +ellps=bessel +pm=ferro \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to +init=epsg:2065 \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

-868208.54      -1095793.57 0.00

Bohužel chybí

# S-JTSK / Krovak
<xxxx> +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999
+x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs  <>

echo "12d48'25.15992  49d27'8.14571" | cs2cs +proj=latlong +datum=WGS84 +to +proj=krovak +lat_0=49.5 \
+lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel \
+units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

-868208.53      -1095793.57 -46.96

ESRI

# S-JTSK (Ferro) / Krovak
<2065> +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999
+x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs  <>

Hodnota +lon_0=24.83333333333333 je chybná. Základní poledník vztažen k Ferru má zeměpisnou délku +lon_0=42.5.

# S-JTSK Krovak
<102065> +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813975277778 +k=0.9999
+x_0=0 +y_0=0 +ellps=bessel +units=m  +no_defs <>

echo "12d48'25.15992  49d27'8.14571" | cs2cs +proj=latlong +datum=WGS84 \
+to +init=esri:102065 +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

-868208.53      -1095793.57 -46.96

Gdalwarp a ogr2ogr

Knihovny GDAL a OGR slouží pro práci s množstvím rastrových a vektorových formátů. S jejich pomocí lze data importovat a velké množství formátů i exportovat. Ve spolupráci s knihovnou PROJ.4 je také možné provádět jejich transformaci.

Mapová projekce krovak není použita ve verzi knihovny GDAL 1.3.2 a nižších. Proto je potřeba vytvořit textový soubor ve známém formátu, který tuto projekci definuje a použít tento soubor jako jeden ze vstupních parametrů:

Soubor krovak.prf:

PROJCS["Krovak",
GEOGCS["Bessel 1841",
DATUM["D_unknown",
SPHEROID["bessel",6377397.155,299.1528128]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Krovak"],
PARAMETER["latitude_of_center",49.5],
PARAMETER["longitude_of_center",24.83333333333333], 
PARAMETER["scale_factor",0.9999],
PARAMETER["false_easting",0], 
PARAMETER["false_northing",0],
UNIT["Meter",1]]

Poznámka: V souvislosti s charakterem vstupních dat a modifikovanou knihovnu PROJ.4 zaměňte řádek

 UNIT["Meter",1]]

za

UNIT["Meter",-1]]

Vlastní transformaci (zde do systému WGS-84) provedeme příkazem

gdalwarp -co "TFW=YES" -s_srs krovak.prf -t_srs "+proj=latlong +datum=WGS84" raster.tif raster-wgs.tif

V případě vektorových dat vypadá příkaz podobně:

ogr2ogr -s_srs krovak.prf -t_srs '+proj=latlong +datum=WGS84' -f "ESRI Shapefile" vektor-out vektor.shp

GRASS lokace S-JTSK

Definice projekce GRASS lokace v S-JTSK by měla vypadat přibližně takto:

GRASS> g.proj -p

-PROJ_INFO-------------------------------------------------
name       : Krovak
proj       : krovak
ellps      : bessel
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1.0

nebo (EPSG 2065)

-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

V případě potřeby transformace dat do jiných souřadnicových systémů (např. WGS-84, UTM, S-42) je nutné doplnit geodetické datum (parametry transformace elipsoidu), viz [1]. Tyto parametry byly určeny pro celé území ČR, resp. SR. V případě menšího zájmového uzemí je vhodné vypočítat parametry znovu, případná chyba se tak zmenší, další informace naleznete zde.

ČR

towgs84    : 570.8,85.7,462.8,4.998,1.587,5.261,3.56

SR (zdroj: Geodetický a kartografický ústav Bratislava):

towgs84    : towgs84=485.021,169.465,483.839,7.786342,4.397554,4.102655,0

Založení nové GRASS lokace

wxGUI - Location wizard

Krok 1	Krok 2
Krok 3	Krok 4
Krok 5	Krok 6

Textové rozhraní

Would you like to create location <t1> ? (y/n) [y] <ENTER>

To create a new LOCATION, you will need the following information:
 
1. The coordinate system for the database
        x,y (for imagery and other unreferenced data)
        Latitude-Longitude
        UTM
        Other Projection
2. The zone for the UTM database
   and all the necessary parameters for projections other than
   Latitude-Longitude, x,y, and UTM
3. The coordinates of the area to become the default region
   and the grid resolution of this region
4. A short, one-line description or title for the location
 
 Do you have all this information? (y/n) [y] <ENTER>

Please specify the coordinate system for location <t1>
 
A   x,y
B   Latitude-Longitude
C   UTM
D   Other Projection
RETURN to cancel

D <ENTER>

Other Projection coordinate system? (y/n) [y] <ENTER>

Please enter a one line description for location <t1>
> popis grass lokace
 
=====================================================
popis grass lokace
=====================================================
ok? (y/n) [y] <ENTER>

Please specify projection name
Enter 'list' for the list of available projections
Hit RETURN to cancel request
> krovak <ENTER>

Do you wish to specify a geodetic datum for this location?(y/n) [y] <ENTER>

Please specify datum name
Enter 'list' for the list of available datums
or 'custom' if you wish to enter custom parameters
Hit RETURN to cancel request
>hermannskogel <ENTER>

Now select Datum Transformation Parameters
Please think carefully about the area covered by your data
and the accuracy you require before making your selection.
 
Enter 'list' to see the list of available Parameter sets
Enter the corresponding number, or <RETURN> to cancel request
>2 (3 pro SR)

Plural form of units: meters <ENTER>

<ESC><ENTER>

Do you accept this region? (y/n) [y] > <ENTER>

Hit RETURN <ENTER>

<ESC><ENTER>

GRASS> g.proj -p
-PROJ_INFO-------------------------------------------------
name       : Krovak
datum      : hermannskogel
towgs84    : 570.8,85.7,462.8,4.998,1.587,5.261,3.56
proj       : krovak
ellps      : bessel
a          : 6377397.1550000003
es         : 0.0066743722
f          : 299.1528128000
-PROJ_UNITS------------------------------------------------
unit       : meter
units      : meters
meters     : 1.0

QGIS a projekt v S-JTSK

Program QGIS umožňuje definici vlastního zobrazení v "Nastavení" → "Uživatelské zobrazení".

Přidejte novou projekci

sjtsk

s parametry

+proj=krovak +ellps=bessel +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

a uložte.

V "Nastavení" → "Vlastnosti projektu" můžete v záložce "Zobrazení" zvolit "User defined coordinate system" a v něm by se vámi definované "sjtsk" mělo objevit.

POZNÁMKA: Je potřeba mít modifikovanou knihovnu PROJ.4!

POZNÁMKA: Ve verzi QGISu 0.8 z 22.8.2006 nešlo tímto způsobem uložit projekci "sjtsk" :-( Uvedený postup popisuje, jak by to asi chodit mělo, resp. jak to chodí, když se jako parametry zadá řetězez "+proj=latlong +ellps=WGS84", tak se to korektně uloží a zobrazí. 23.8. zkouším kompilovat čerstvé SVN a případně píšu do qgis-listu -- Jachym 10:02, 23. 8. 2006 (CEST)

Související články

S-JTSK (zastaralý článek)

Obsah

Křovákovo zobrazení

Knihovna PROJ.4

Velkopáteční odpověď Franka Warmerdama

PROJ mailing list / březen 2006

Reakce Franka Warmerdama k vydání PROJ.4 4.5.0 / duben 2006

Další zmínky o Křovákově zobrazení

Bugreport

Jak aplikovat patch

PROJ.4 4.5.0 příběh s "dobrým" koncem

Praktická ukázka

Gdalwarp a ogr2ogr

GRASS lokace S-JTSK

Založení nové GRASS lokace

wxGUI - Location wizard

Textové rozhraní

QGIS a projekt v S-JTSK

Související články

Osobní nástroje

Jmenné prostory

Varianty

Zobrazení

Akce

Hledat

Navigace

Dokumentace

Geodata

stránky

Nástroje