S-JTSK / Grid
Článek popisuje použití programu cs2cs z knihovny Proj.4 (lze analogicky využít i GDAL nebo GRASS, PostGIS) pro transformace pomocí GRIDu s centimetrovou přesností mezi:
- plošnými souřadnicemi ETRS89 a S-JTSK,
- výškou nad elipsoidem GRS80 a výškovým systémem "Balt po vyrovnání" (Bpv)
- a dále se zabývá transformačními vztahy mezi WGS84/ITRF a ETRS89.
Teorie transformace GRIDem
Převod 2D souřadnic:
- metoda popisovaná v tomto článku skrze Proj.4:
- S-JTSK→ETRS89 - GRIDem,
- ETRS89→WGS84 - 7mi prvkovou transformací.
- metoda používaná jako oficiální postup ČÚZK[1][2]:
- S-JTSK→S-JTSK/05 - GRIDem,
- S-JTSK/05→ETRS89 - 7mi prvkovou transformací,
- ETRS89→WGS84 - tímto převodem se ČÚZK nezabývá.
Převod výšky je v Proj.4 i v oficiálním postupu ČÚZK analogický:
- výška nad elipsoidem GRS-80→výška Balt po vyrovnání(Bpv) - GRIDem.
Dostupné GRIDy
Grid | Typ | Vztah mezi CRS | Binární soubor gridu 5) | Přesnost vůči ČUZK DATAZ [m] 2) | Zdrojový soubor (autorův link) | ||||
---|---|---|---|---|---|---|---|---|---|
<0,05 | <0,10 | <1,00 | stdev.p | max 3) | |||||
Seidl2016 | XY | ETRS89/ETRF2000 × S-JTSK | seidl16_jtskcz.llb | 97,8% | 99,8% | 100,0% | 0,017 | 0,966 | seidl16_jtskcz.lla (Seidl) |
Seidl2014 | Z | H_GRS80 × H_Bpv | seidl16_bpvcz.gtx | 84,2% | 97,6% | 100,0% | 0,032 | 0,893 | - (Seidl) |
Ježek2008 4) | XY | ETRS89/ETRF89 × S-JTSK | jezek08_jtskcz.llb | 81,3% 1) | 99,4% | 100,0% | 0,023 | 0,957 | jezek08_jtskcz.lla (Ježek) |
1) Zhoršená přesnost je pravděpodobně dána testovanými body z novější realizace ETRS89 tj. realizace ETRF2000, která se po formálním provedením v roce 2.1.2011 projevila posunem v řádech 2 cm[3]. Autor v roce 2009 uvádí RMSE(střední kvadratická chyba) transformace tímto GRIDem 1,8cm, maximální 17 cm[4] 2) Trigonometrické body (TB) a zhušťovací body (ZhB) k datu 2016-12-19 s obojími souřadnicemi ETRF2000 a S-JTSK, celkem využito 41 647 bodů[5]. |
Transformace mezi ETRF2000 a S-JTSK pomocí GRIDu
- Na binární soubor GRIDu je praktické odkazovat jen jeho názvem bez cesty k souboru, což je možné docílit jejich nakopírováním s právy "sudo" do adresáře /usr/share/proj/ (pro Linux) případně tam, kde jsou jiné binární soubory s gridy (např. soubor "conus"). Pokud chcete umístit soubor do jiného než standardního ('zakompilovaného') adresáře, je nutné na něj ukázat pomocí systémové proměnné PROJ_LIB.
- Použití GRIDu automaticky ignoruje výchozí nebo explicitně zadanou 7mi prvkovou transformaci tj. parameter +towgs84.
- Stažení GRIDu Seidl16:
wget http://gis.templ.net/grid/seidl16_jtskcz.llb http://gis.templ.net/grid/seidl16_bpvcz.gtx sudo cp seidl16_jtskcz.llb seidl16_bpvcz.gtx /usr/share/proj/
- Samotnou transformaci pomocí GRIDu v PROJ.4 (cs2cs, proj, invproj) definice S-JTSK:
"+init=epsg:5514 +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx"
nebo rozšířený zápis:
"+proj=krovak +ellps=bessel +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx"
ETRF2000 → S-JTSK
B,L,HGRS80 → X,Y,hBpv na příkladu souřadnic GNSS stanice CZEPOS TUBO (Y = 599 131,597 X = 1 159 442,044 h = 279,584)[7]:
echo "16d35'34.20396 49d12'21.21026 324.272" | cs2cs -f "%.3f" +init=epsg:4258 +to +init=epsg:5514 +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx -599131.612 -1159442.061 279.592
Diference: Δ Y = 1,5 cm; Δ X = 1,7 cm; Δ h = 0,8 cm
Variantní a delší zápis:
echo "16d35'34.20396 49d12'21.21026 324.272" | cs2cs -f "%.3f" +proj=lonlat +ellps=GRS80 +towgs84=0,0,0 +to +proj=krovak +ellps=bessel +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx -599131.612 -1159442.061 279.592
S-JTSK → ETRF2000
X,Y,hBpv → B,L,HGRS80 na příkladu souřadnic GNSS stanice CZEPOS TUBO (B = 49°12'21.21026" L = 16°35'34.20396" H = 324.272)[8]:
echo "-599131.597 -1159442.044 279.584" | cs2cs -f "%.9f" +init=epsg:5514 +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx +to +init=epsg:4258 16.592834613 49.205891902 324.26351100
Diference: Δ L = 0,00065"; Δ B = 0,00059"; Δ H = 0,9 cm
Variantní a delší zápis:
echo "-599131.597 -1159442.044 279.584" | cs2cs -f "%.9f" +proj=krovak +ellps=bessel +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx +to +proj=lonlat +ellps=GRS80 +towgs84=0,0,0 16.592834613 49.205891902 324.26351100
Transformace mezi WGS84/ITRS a ETRS89
Teorie vztahů mezi souřadnicovými systémy
Na základě znalosti o pohybu kontinentálních desek (na území ČR cca 30 mm/rok na SV) a jejich vlastní vnitřní dynamiky (na území ČR v řádech desetin mm/rok u stabilních stanic, výjimečně až 3 mm/rok v případě GNSS stanice GEONAS Sněžka CSNE[9][10]) bylo zavedeno několik typů souřadnicových referenčních systémů (CRS):
- Mezinárodní - např. ITRS, pro něhož je podmínkou celková nulová rotace stanic skrze niž je systém realizován;
- Regionální (kontinentální) - např. ETRS, který je spřažen s euroasijskou kontinentální deskou. Nemá konstantní polohu souřadnicových os, protože ty se v čase natáčejí dle pohybu euroasijské kontinentální desky;[11]
- Národní (lokální) - např. S-JTSK
Protože tyto souřadnicové systémy jsou vlivem geologických procesů v čase proměnné, popisujeme při přesnějších měřeních jejich stav ještě rámcem a epochou udávanou v letech, např. ETRF2000 (2005.0) nebo ITRF2000 (2005.0)[12]:
- Referenční systém - definice základních konstant a algoritmů např. elipsoid, počátek, osy, jednotky
- Referenční rámec - autoritativní fyzická realizace referenčního systému na konkrétních měřicích bodech a v konkrétním čase
- Epocha např. (2005.0) - vztažný rok provedeného měření
Přehled realizací vybraných pro území ČR:
CRS | S-JTSK | ETRS89 | ITRF | WGS 84 |
---|---|---|---|---|
typ | národní | evropský | mezinárodní | mezinárodní |
elipsoid | Bessel1841 | GRS80 | GRS80 | WGS 84 |
rámec | S-JTSK | ETRF89 1) | ITRF88 | WGS 84 (Doppler) 2) |
S-JTSK/05 | ETRF90 | ITRF89 1) | WGS 84 (G730) | |
ETRF91 | ITRF90 2) | WGS 84 (G873) | ||
ETRF92 | ITRF91 | WGS 84 (G1150) | ||
ETRF93 | ITRF92 | WGS 84 (G1674) | ||
ETRF94 | ITRF93 | WGS 84 (G1762) | ||
ETRF96 | ITRF94 | |||
ETRF97 | ITRF96 | |||
ETRF2000 | ITRF97 | |||
ITRF2000 | ||||
ITRF2005 | ||||
ITRF2008 | ||||
ITRF2014 | ||||
1) rámce ETRF89 a ITRF89 jsou v epoše 1989.0 totožné 2) systémy ITRF90 a WGS84 (Doppler) jsou navázány převodním algoritmem |
Systémy ETRF89 a ITRF89 byly v epoše 1989.0 totožné a pro jejich následné realizace jsou v příslušných epochách známé vztahy.[14] Přestože elipsoidy GRS80 a WGS84 se liší pouze o 0,1 mm v délce poloměru vedlejší poloosy[15], je přesný převodní vztah definován pouze mezi mezi rámci WGS84/Doppler a ITRS90[16], v ostatních realizacích WGS84 není přesný vztah popsán.[17][18] Změna epochy v daném rámci CRS je možná pouze tehdy, pokud známe nejen souřadnice ale i rychlost (mm/rok) měřícího bodu, případně za použití modelu (např. NUVEL-A). Anebo když bezpečně víme, že v daném CRS při zvolené přesnosti je pohyb zanedbatelný. Proto na území ČR v rámci ETRS lze pro centimetrové přesnosti považovat pohyb měřicích bodů za zanedbatelný a změnu epochy v rámci desetiletí dle potřeby případně provést.
Transformační klíče
Vztah rámců v dané epoše | Klíč s prvky Tx[m], Ty[m], Tz[m], Rx["], Ry["], Rz["], sf-1[ppm] | |||
---|---|---|---|---|
WGS84/Doppler (1989.0) a ETRF2000 (1989.0) | +towgs84="0.0285,-0.5506,-0.2475,0.0183,-0.0003,0.00706,-0.00755" | |||
ITRF2008 (2005.0) a ETRF2000 (2005.0) $ | +towgs84="-0.0526,-0.0498,0.0675,-0.001296,-0.00784,0.012672,-0.00174" | |||
$ Při přesnosti 0,1 m lze považovat WGS84 (G1674) (2005.0) a ITRF2008 (2005.0) za totožné[19]. Podobně jsou i rámce IGS08 a IGb08 odvozeny z ITRF2008 Zdrojová data k výpočtu klíčů (např. pro změnu epochy): itrs-etrs_shift.xls na základě ITRF_Transformation_Parameters.xlsx (komentář). |
S-JTSK → WGS84 (1989.0)
X,Y,hBpv → B,L,HWGS84 (1989.0) na příkladu souřadnic fiktivního bodu:
echo "-600000.000 -1100500.000 300.000" | cs2cs -f "%.8f" +init=epsg:5514 +nadgrids=./seidl16_jtskcz.llb +geoidgrids=./seidl16_bpvcz.gtx +to +init=epsg:4258 | cs2cs -f "%.8f" +init=epsg:4258 +towgs84="0.0285,-0.5506,-0.2475,0.0183,-0.0003,0.00706,-0.00755" +to +init=epsg:4326 16.49279389 49.73193867 344.49183605
WGS84 (1989.0) → S-JTSK
B,L,HWGS84 (1989.0) → X,Y,hBpv na příkladu souřadnic fiktivního bodu:
echo "16.49279389 49.73193867 344.49183605" | cs2cs -f "%.8f" +init=epsg:4326 +to +init=epsg:4258 +towgs84="0.0285,-0.5506,-0.2475,0.0183,-0.0003,0.00706,-0.00755" | cs2cs -f "%.3f" +init=epsg:4258 +to +init=epsg:5514 +nadgrids=./seidl16_jtskcz.llb +geoidgrids=./seidl16_bpvcz.gtx -600000.000 -1100500.000 300.000
S-JTSK → ITRF2008 (2005.0)
X,Y,hBpv → B,L,HGRS80 (2005.0) na příkladu souřadnic GOPE (Y = 719512.302 m X = 1065662.523 m hBpv = 547.58587 m)(X = 3979316.130 Y = 1050312.476 Z = 4857067.109 m )[20][21]:
echo "-719512.302 -1065662.523 547.58587" | cs2cs -f "%.9f" +init=epsg:5514 +nadgrids=./seidl16_jtskcz.llb +geoidgrids=./seidl16_bpvcz.gtx +to +init=epsg:4258 | cs2cs -f "%.9f" +init=epsg:4258 +towgs84="-0.0526,-0.0498,0.0675,-0.001296,-0.00784,0.012672,-0.00174" +to +init=epsg:4326 14.785622455 49.913704608 592.587991651 echo "14.785622455 49.913704608 592.587991651" | cs2cs -f "%.3f" +ellps=GRS80 +proj=lonlat +to +ellps=GRS80 +proj=geocent 3979316.125 1050312.484 4857067.091
Diference: Δ X = 0,5 cm; Δ Y = 0,8 cm; Δ Z = 1,8 cm
Pozn: Protože Proj.4 uplatňuje parametr +towgs84 pouze proti elipsoidu WGS 84, je nutné využít EPSG:4326 namísto očekávaného EPSG:7911 (ITRF2008 geographic).
ITRF2008 (2005.0) → S-JTSK
B,L,HGRS80 (2005.0) → X,Y,hBpv na příkladu souřadnic GOPE (Y = 719512.302 m X = 1065662.523 m hBpv = 547.58587 m)(X = 3979316.130 Y = 1050312.476 Z = 4857067.109 m )[22][23]:
echo "3979316.130 1050312.476 4857067.109" | cs2cs -f "%.9f" +ellps=GRS80 +proj=geocent +to +ellps=GRS80 +proj=lonlat 14.785622328 49.913704698 592.603432324 echo "14.785622328 49.913704698 592.603432324" | cs2cs -f "%.9f" +init=epsg:4326 +to +init=epsg:4258 +towgs84="-0.0526,-0.0498,0.0675,-0.001296,-0.00784,0.012672,-0.00174" | cs2cs -f "%.3f" +init=epsg:4258 +to +init=epsg:5514 +nadgrids=./seidl16_jtskcz.llb +geoidgrids=./seidl16_bpvcz.gtx -719512.310 -1065662.512 547.601
Diference: Δ Y = 0,8 cm; Δ X = 1,1 cm; Δ h = 1,5 cm
Pozn: Protože Proj.4 uplatňuje parametr +towgs84 pouze proti elipsoidu WGS 84, je nutné využít EPSG:4326 namísto očekávaného EPSG:7911 (ITRF2008 geographic).
Další použití a poznámky
Použití GRIDu v PostgreSQL/PostGIS
Transformace lze obdobně nastavit i třeba v PostGIS pomocí editace tabulky spatial_ref_sys, sloupce proj4text, řádku srid=2065 - do položky stačí zadat "+proj=krovak +ellps=bessel +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx". Na databázovém serveru je nutné mít nahraný soubor s GRIDem.
bash:
wget http://gis.templ.net/grid/seidl16_jtskcz.llb seidl16_jtskcz.llb http://gis.templ.net/grid/seidl16_bpvcz.gtx sudo cp seidl16_jtskcz.llb seidl16_bpvcz.gtx /usr/share/proj
postgres (testováno na Linux Mint 18.1, Postgresql v9.5, Postgis v2.2.1):
INSERT INTO spatial_ref_sys(srid, auth_name, auth_srid, srtext, proj4text)VALUES (999, 'seidl', 999, 'NULL', '+proj=krovak +ellps=bessel +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx');
Dotaz probehl úspešne: bylo ovlivnený 1 rádek
SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(-599131.597 -1159442.044 279.584)',999),4258)) As wgs_geom; "POINT Z (16.5928346131882 49.20589190164 324.26351099927)" SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT(16.592834433 49.205891739 324.272)',4258),999)) As wgs_geom; "POINT Z (-599131.612000511 -1159442.06056716 279.592488671903)"
CZEPOS TUBO[24]:
- S-JTSK: Y = 599 131,597 X = 1 159 442,044 h = 279,584
- ETRF2000: B = 49°12'21.21026" L = 16°35'34.20396" H = 324.272
- ETRF2000: B = 49.205891739 L = 16.592834433 H = 324.272)
Použití GRIDu v GDAL
Pro převod pomocí knihovny ogr2ogr lze obdobně využít:
wget http://gis.templ.net/grid/seidl16_jtskcz.llb http://gis.templ.net/grid/seidl16_bpvcz.gtx sudo cp seidl16_jtskcz.llb seidl16_bpvcz.gtx /usr/share/proj/ ogr2ogr -s_srs "+proj=krovak +ellps=bessel +nadgrids=seidl16_jtskcz.llb +geoidgrids=seidl16_bpvcz.gtx +wktext" -t_srs "EPSG:4258" ./body_etrs.shp ./body_krovak.shp
Vytvoření binárního souboru gridu
- Stáhněte si soubor lla: např. seidl16_jtskcz.lla
- Pomocí příkazu "nad2bin" z Proj.4 převeďte .lla soubor na binární podobu a uložte do adresáře s binárními soubory. Ukázka zápisu příkazu nad2bin na linuxu:
wget http://gis.templ.net/grid/seidl16_jtskcz.llb http://gis.templ.net/grid/seidl16_bpvcz.gtx nad2bin < seidl16_jtskcz.lla ./seidl16_jtskcz.llb
Problematika změny ETRS89 na ETRF2000
Při změně 2.1.2011 ČÚZK v realizaci ETRS89 v ČR z rámce ETRF1989 na ETRF2000 došlo k posunu souřadnic ETRS89 (souřadnice S-JTSK zůstaly nedotčeny). Tento posun byl v rozměrech 1 až 2 cm, viz například souřadnice CZEPOS do 2011 a dále, nebo souřadnice stanice GOPE v transformační aplikaci EUREF, nebo Kostelecký.
Poznámky
- Výše popsaný text lze použít i pro GDAL/OGR, GRASS a ostatní systémy založené na PROJ.4.
- Problém se znaménky a případným přehozením os je věc samotného zobrazení (+proj=krovak) a tedy stále zůstává. Pro souřadnice vztažené k osám ve směru 'sever východ' (záporné hodnoty) lze přidat parametr +czech a používat Proj.4 ve verzi větší než 4.5.0.
- Debuging je možné zapnout pomocí nastavení systémové proměnné PROJ_DEBUG=true. Poté se vypisují informace o použitých gridech apod.
- V případě chybného výsledku prosím pošlete mail autorům GRIDu na michal . seidl @ gmail.com nebo h.jezek @ centrum.cz. Přiložte prosím také informace o verzi Proj, případně celý výpis při nastavení sys. proměnné PROJ_DEBUG=true.
Literatura
- Jan Ježek: Precise transformation between S-JTSK and ETRS89 (WGS-84) in GIS, Geoinformatics Ostrava 2009
- Jan Ježek: Implementace transformací souřadnicových systémů v GIS, disertace, FS ČVÚT, Praha 2009
- Jan Ježek: Studijní materiály MK2, 18. října 2012
- Jan Ježek: Vývoj programového modulu pro převod souřadnic mezi kartografickými zobrazeními, DP, ČVUT, prosinec 2003
- Jan Ježek: Výpočet transformačnı́ch koeficinetů vybraných 2D transformacı́, 2009
- Jan Ježek: prezentace disertační práce z Geoinformatics FCE CTU 2008.
- Klokan Petr Přidal: Konverze do formátu vhodného pro PROJ.4, 2008
- M. Seidl, P. Trasak: Proj.4 cartographic projection library as 3d transformation tool in Czech context, abstract, SGEM 2014
- M. Seidl: Improved proj.4 grid files for region of Czech republic, abstract, SGEM 2016
Reference
- ↑ http://www.cuzk.cz/Zememerictvi/Geodeticke-zaklady-na-uzemi-CR/GNSS/etrf00-jtsk-v1012-a-etrf00-jtsk-v1203.aspx
- ↑ http://geoportal.cuzk.cz/Default.aspx?head_tab=sekce-01-gp&mode=TextMeta&text=wcts&menu=19
- ↑ http://hvr.cz/doc-352/
- ↑ http://gis.vsb.cz/GIS_Ostrava/GIS_Ova_2009/sbornik/Lists/Papers/088.pdf
- ↑ http://dataz.cuzk.cz/hledej.php
- ↑ http://www.fce.vutbr.cz/ged/GEDASP/GNSS/Prezentace2016/08_Nagl-prevodni%20tabulky.pdf
- ↑ http://czepos.cuzk.cz/_souradnice.aspx
- ↑ http://czepos.cuzk.cz/_souradnice.aspx
- ↑ ftp://olggps.oeaw.ac.at/pub/EUREF_camp/ETRS89-CZ_FinalReport.pdf
- ↑ http://www.epncb.oma.be/_productsservices/coordinates/posvel_map.php
- ↑ http://transformace.webst.fd.cvut.cz/Iframe/ETRS_iframe.htm
- ↑ http://transformace.webst.fd.cvut.cz/Iframe/ETRS_iframe.htm
- ↑ https://confluence.qps.nl/pages/viewpage.action?pageId=29855173
- ↑ http://etrs89.ensg.ign.fr/memo-V8.pdf#page=7
- ↑ https://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS_84
- ↑ ftp://itrf.ensg.ign.fr/pub/itrf/WGS84.TXT
- ↑ http://www.unoosa.org/pdf/icg/2012/3.pdf#page=71
- ↑ https://confluence.qps.nl/pages/viewpage.action?pageId=29855173
- ↑ ftp://itrf.ensg.ign.fr/pub/itrf/WGS84.TXT
- ↑ http://www.epncb.oma.be/_productsservices/coordinates/crd4station.php?station=GOPE00CZE
- ↑ http://oko.pecny.cz/vesog/stanice/gope.html#coord
- ↑ http://www.epncb.oma.be/_productsservices/coordinates/crd4station.php?station=GOPE00CZE
- ↑ http://oko.pecny.cz/vesog/stanice/gope.html#coord
- ↑ http://czepos.cuzk.cz/_souradnice.aspx