S-JTSK

Z GRASSwikiCZ

Článek popisuje základní práci s daty v souřadnicovém systému S-JTSK Systém Jednotné trigonometrické sítě katastrální ve FOSS GIS. Na příkladech ukazuje transformaci běžných geodat mezi vybranými (WGS84, S-JTSK, UTM 33N) souřadnicovými referenčními systémy angl. coordinate reference system (CRS) obvykle používané na území České, případně Slovenské republiky s využitím nástroje z balíků Proj.4 (cs2cs) a GDAL (ogr2ogr, gdalwarp). Dále je součástí článku ukázka nastavení souřadnicových referenčních systémů v GIS aplikacích GRASS a QGIS.

Druhy transformací podle přesnosti (WGS84/ERTS89->S-JTSK) ve 2D:

  • 100-110 metrů na severozápad - základní transformace bez transformačního klíče; nepoužívá se
  • < 1 metr - sedmiprvková transformace s klíčem "+towgs"; viz níže globální transformační klíč pro ČR/SR
  • < 0,1 metr - pomocí gridu (ČR); zvláštní článek


Obsah

Teorie Křovákova zobrazení

Souřadnicový referenční 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í s poledníkem Ferro (dnes se užívá definice k poledníku Greenwich

Reference: [1]

Definice CRS v Proj.4/GDAL

Pro práci s CRS používat:

  • EPSG a ESRI kódy odkazující do vnitřní knihovny programu
  • kombinace kódů a doplňujících parametrů (jen Proj.4)
  • definici CRS ze souboru


název CRS oficiální EPSG/ESRI kódy odkazy v Proj.4 odkazy v GDAL správná definice dle EPSG/ESRI (formát Proj)
S-JTSK East North (Greenwitch) ESRI:102067
EPSG:5514
ESRI:102067 EPSG:102067 viz níže
S-JTSK South West (Ferro) EPSG:2065 chybně EPSG:2065 (bug report) chybně EPSG:2065 (report) viz níže
WGS84 EPSG:4326 EPSG:4326 EPSG:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
WGS84 / UTM zone 33N EPSG:32633 EPSG:32633 EPSG:32633 +proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

Pozn: Pro méně přesné submetrové transformace je možno WGS84 považovat za identický s ETRS 89 z důvodu jejich velmi malých odlišností.

EPSG:5514/ESRI:102067

Definice S-JTSK v typické verzi pro GIS se zápornými souřadnicemi (osy: East-North, poledník: Greenwich) vč. transformačního klíče pro ČR ve formátu OGC WKT:

PROJCS["Czech GIS S-JTSK (Greenwich) / Krovak",
   GEOGCS["Czech S-JTSK (Greenwich)",
       DATUM["Czech S-JTSK",
           SPHEROID["Bessel 1841",6377397.155,299.1528128,
               AUTHORITY["EPSG","7004"]],
           TOWGS84[570.8,85.7,462.8,4.998,1.587,5.261,3.56]],
       PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Krovak",
       AUTHORITY["EPSG","9819"]],
   PARAMETER["latitude_of_center",49.5],
   PARAMETER["longitude_of_center",24.83333333333333],
   PARAMETER["azimuth",0],
   PARAMETER["pseudo_standard_parallel_1",0],
   PARAMETER["scale_factor",0.9999],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   UNIT["Meter",1]]

Formát Proj (další formáty viz Spatial Reference/102067):

+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=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

EPSG:2065

Definice méně užívané verze S-JTSK s kladnými souřadnicemi (osy: South-West, poledník: Ferro) vč. transformačního klíče pro ČR ve formátu OGC WKT:

PROJCS["Czech S-JTSK (Ferro) / Krovak",

   GEOGCS["Czech S-JTSK (Ferro)",
       DATUM["Czech S-JTSK",
           SPHEROID["Bessel 1841",6377397.155,299.1528128,
               AUTHORITY["EPSG","7004"]],
           TOWGS84[570.8, 85.7, 462.8, 4.998, 1.587, 5.261, 3.56]],
       PRIMEM["Ferro",-17.666666666666668,
           AUTHORITY["EPSG","8909"]],
       UNIT["degree",0.017453292519943295],
       AXIS["Geodetic longitude",EAST],
       AXIS["Geodetic latitude",NORTH]],
   PROJECTION["Krovak",
       AUTHORITY["EPSG","9819"]],
   PARAMETER["latitude_of_center",49.5],
   PARAMETER["longitude_of_center",42.5],
   PARAMETER["azimuth",30.288139722222223],
   PARAMETER["pseudo_standard_parallel_1",78.5],
   PARAMETER["scale_factor",0.9999],
   PARAMETER["false_easting",0.0],
   PARAMETER["false_northing",0.0],
   UNIT["Meter",-1]]

Formát Proj (další formáty viz Spatial Reference/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 +to_meter=-1 +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56

Transformace mezi CRS

V případě potřeby transformace geodat do jiných SRS (např. WGS84 → S-JTSK) je nutné doplnit dvojici CRS o zpřesňující transformační klíč [1]. Tyto parametry byly určeny pro celé území ČR, resp. SR. V případě menšího zájmového území je možné vypočítat lokální transformační klíč, nebo rovnou použít S-JTSK-Grid případná chyba se tak až řádově zmenší. Nepoužití klíče vede ke zhoršení přesnosti o dva řády, přibližně < 100 metrů do S-JTSK, 3 vteřiny do WGS84 v horizontálních souřadnicích. Ve výšce Z versus HAE je chyba právě separace geoidu. Použití zpřesňujícího transformačního klíče nemá vliv na rychlost transformace.

+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
+towgs84=485.0,169.5,483.8,7.786,4.398,4.103,0

Jednotlivé prvky globálního klíče lze definovat přesněji, přínos ve vyšší přesnosti je však minimální: ČR [4], SR [5] Od roku 2011 je v databázi EPSG:2065 díky ČÚZK modifikovaný globální klíč, ani jeho přínos není významný[6].

Testovací bod č. 109 sítě DOPNUL [7]

WGS84 S-JTSK
zem. délka (lat) zem. šířka (lon) výška nad elipsoidem (HAE) Y X Z
12d48'25.15992
12,8069888666667
49d27'8.14571
49,4522626972222
559.417 -868208.52 -1095793.96 512.30


Testované verze GDAL, Proj.4

Níže uvedené příklady byly testovány:

  • Ubuntu 9.04, Proj.4 distribuční verze 4.6.0 (2007/12/21) a GDAL distribuční verze 1.5.2 (2008/05/29) dne 8.-10. května 2009
  • Ubuntu 9.04, Proj.4 Les-ejk [8] verze 4.6.1 (2008/09/21) dne 11. května 2009
  • Ubuntu 8.04, Proj.4 distribuční verze 4.6.0 (2007/12/21) a GDAL distribuční verze 1.4.4 (2007/11/23) dne 11. května 2009

Transformace textových geodat

Příklady transformace textových geodat za pomoci sady Proj.4 nástroje cs2cs. Výchozí souřadnici Z nebo HAE lze pro výpočty v rovině vynechat.

WGS84 → S-JTSK

echo "12d48'25.15992 49d27'8.14571 559.417" | cs2cs +init=epsg:4326 +to +init=esri:102067 \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 

-868208.54	-1095793.58 512.46

Variantní zápis:

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

-868208.54	-1095793.58 512.46

Diference: Δ Y = 1cm; Δ X = 38cm; Δ Z = 16cm

S-JTSK → WGS84

echo "-868208.53 -1095793.57 512.30" | cs2cs +init=esri:102067 \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to +init=epsg:4326

12d48'25.16"E	49d27'8.146"N 559.261

Variantní zápis:

echo "-868208.53 -1095793.57 512.30" | cs2cs +proj=krovak +ellps=bessel \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to +proj=latlong +datum=WGS84 

12d48'25.16"E	49d27'8.146"N 559.261

Diference: Δ lat = -1.0E-007 stupně; Δ long = -0.7E-007 stupně; Δ Z = 21cm

S-JTSK → UTM 33N

echo "-868208.53 -1095793.57 512.30" | cs2cs +init=esri:102067 \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to +init=epsg:32633

341060.93	5480045.49 559.26

Variantní zápis:

echo "-868208.53 -1095793.57 512.30" | cs2cs +proj=krovak +ellps=bessel \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 \
+to +proj=utm +zone=33N +datum=WGS84

341060.93	5480045.49 559.26

UTM 33N → S-JTSK

echo "341060.93 5480045.49 559.26" | cs2cs +init=epsg:32633 +to +init=esri:102067 \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 

-868208.53	-1095793.57 512.30

Variantní zápis je obdobný jako inverzní transformace S-JTSK → UTM 33N.

WGS84 → UTM 33N

echo "12d48'25.15992 49d27'8.14571 559.417" | cs2cs +init=epsg:4326 +to +init=esri:32633

341060.93	5480045.47 559.42

Variantní zápis:

echo "12d48'25.15992 49d27'8.14571 559.417" | cs2cs +proj=latlong +datum=WGS84 \
+to +proj=utm +zone=33N +datum=WGS84

341060.93	5480045.47 559.42

UTM 33N → WGS84

echo "341060.93 5480045.47 559.42" | cs2cs +init=esri:32633 +to +init=epsg:4326

12d48'25.16"E	49d27'8.146"N 559.420

Variantní zápis je obdobný jako inverzní transformace WGS84 → UTM 33N.

Transformace vektorových geodat

Příklady převodu vektorových geodat za pomoci sady GDAL nástroje ogr2ogr. Příklad prezentován na vrstvě typu ESRI Shapefile z datové sady FreeGeodataCZ archiv [2] nebo [3].

S-JTSK → WGS84

ogr2ogr -s_srs "+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=greenwich +units=m +no_defs \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" -t_srs "epsg:4326" -f "ESRI Shapefile" cfm_wgs84.shp cfm.shp

Variantní zápis využívající uživatelem vytvořený soubor krovak.prf[9]:

ogr2ogr -s_srs krovak.prf -t_srs "epsg:4326" -f "ESRI Shapefile" cfm_wgs84.shp cfm.shp

WGS84 → S-JTSK

ogr2ogr -s_srs "epsg:4326" -t_srs "+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=greenwich +units=m +no_defs \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" -f "ESRI Shapefile" cfm_jtsk.shp cfm.shp

Variantní zápis využívající uživatelem vytvořený soubor krovak.prf[9]:

ogr2ogr -s_srs "epsg:4326" -t_srs krovak.prf -f "ESRI Shapefile" cfm_jtsk.shp cfm.shp

Transformace rastrových geodat

Příklady transformace rastrových geodat za pomoci sady GDAL nástroje gdalwarp. Příklad prezentován na rastru typu tiff + world file z datové sady FreeGeodataCZ archiv [4] nebo [5].

S-JTSK → WGS84

gdalwarp -s_srs "+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=greenwich +units=m +no_defs \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" -t_srs "epsg:4326" dem_gtopo30.tif dem_gtopo30_wgs84.tif

Variantní zápis využívající uživatelem vytvořený soubor krovak.prf[9]:

gdalwarp -s_srs krovak.prf -t_srs "epsg:4326" dem_gtopo30.tif dem_gtopo30_wgs84.tif

WGS84 → S-JTSK

gdalwarp -s_srs "epsg:4326" -t_srs "+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=greenwich +units=m +no_defs \
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56" dem_gtopo30.tif dem_gtopo30_jtsk.tif

Variantní zápis využívající uživatelem vytvořený soubor krovak.prf[9]:

gdalwarp -s_srs "epsg:4326" -t_srs krovak.prf dem_gtopo30.tif dem_gtopo30_jtsk.tif

Lokace (location) s CRS S-JTSK v GRASS

Platné pro verze Grass 6.4.0 a vyšší, případně 6.2.3.

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 transformační klíč (towgs84).

Založení nové GRASS lokace

V GRASS 6.4.0 wxGUI pomocí Location wizard

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

Pomocí 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

Volba projekce S-JTSK v QGIS

QGIS od verze 1.0.1 do 1.7.3 obsahuje korektní CRS pro S-JTSK S-JTSK (Greenwich) / Krovak EPSG:102067, ale ve verzi 1.7.4 a 1.8.0 byl dočasně vypuštěn

+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=greenwich +units=m +no_defs


Pro potřeby použití více druhů CRS v jednom mapovém projektu je třeba nadefinovat S-JTSK s konkrétním globálním transformačním klíčem S-JTSK +towgs84:

+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=greenwich +units=m +no_defs 
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
Ukázka definice upraveného CRS v QGIS

Mapový projekt S-JTSK s vrstvou WGS84

Pokud chcete použít transformaci on-the-fly např. pro projekt v S-JTSK s vloženou vrstvou ve WGS84, musíte definovat mapovému projektu nový CRS identický k ESRI:102067 s transformačním klíčem +towgs84.

Původní CRS ve vlastnostech mapového projektu QGIS
Upravený CRS ve vlastnostech mapového projektu QGIS


Mapový projekt WGS84 s vrstvou S-JTSK

Pokud chcete použít transformaci on-the-fly pro projekt v WGS84 s vloženou vrstvou ve S-JTSK, musíte této vrstvě definovat CRS identický k ESRI:102067 s transformačním klíčem +towgs84.

Původní CRS ve vlastnostech datové vrstvy v QGIS
Upravený CRS ve vlastnostech datové QGIS

Odkazy

Související články

Reference

  1. Jan Ježek: Vývoj programového modulu pro převod souřadnic mezi kartografickými zobrazeními, diplomová práce ČVUT, Fakulta stavební, Praha 2003
  2. Bundesamt für Kartographie und Geodäsie: Description of Transformation - CZ_S-JTSK to ETRS89, 6. prosinec 2004
  3. Bundesamt für Kartographie und Geodäsie: Description of Transformation - SK_S-JTSK to ETRS89, 6. prosinec 2004
  4. Pavel Mátl: Využití systému WGS84 pro katastrální mapování, str. 54, diplomová práce ZČÚ, Plzeň 2006
  5. Geodetický a kartografický ústav Bratislava: Transformácia medzi S-JTSK a ETRS89, SKPOS - Slovenská Priestorová Observacná Služba yužitia signálov GNSS, 7. dubna 2009
  6. Coordinate Transformation: S-JTSK to WGS 84 (5) in GeoRepository
  7. Zdeněk Hrdina: Transformace souřadnic ze systému WGS-84 do systému S-JTSK, str. 18-21, ČVUT, Praha 1997
  8. Jáchym Čepický: Les-ejk UbuntuGIS repository. 2009
  9. 9,0 9,1 9,2 9,3 ESRI:102607 in Spatial Reference

Externí odkazy

Osobní nástroje
Jmenné prostory

Varianty
Akce
Navigace
Dokumentace
Geodata
stránky
Nástroje