S-JTSK / Grid

Z FreeGIS portál
Přejít na: navigace, hledání

Č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].
3) Ve zbylých 0,2-2,4% nelze rozhodnout zda jde o chybu v DATAZ nebo lokální deformace základů S-JTSK[6].
4) Na základě 3070 bodů ČÚZK vybrané údržby s obojími souřadnicemi.
5) Binární gridy testovány na Linux Mint 18, Proj.4 v4.9.2.

Zdrojová data k testu: dataz20161219-test_grid.xls

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
Pozn: Další podobnosti systémů WGS84 a ITRS jsou v řádech 0,1 m popsány v odkaze[13]

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

Reference

  1. ↑ http://www.cuzk.cz/Zememerictvi/Geodeticke-zaklady-na-uzemi-CR/GNSS/etrf00-jtsk-v1012-a-etrf00-jtsk-v1203.aspx
  2. ↑ http://geoportal.cuzk.cz/Default.aspx?head_tab=sekce-01-gp&mode=TextMeta&text=wcts&menu=19
  3. ↑ http://hvr.cz/doc-352/
  4. ↑ http://gis.vsb.cz/GIS_Ostrava/GIS_Ova_2009/sbornik/Lists/Papers/088.pdf
  5. ↑ http://dataz.cuzk.cz/hledej.php
  6. ↑ http://www.fce.vutbr.cz/ged/GEDASP/GNSS/Prezentace2016/08_Nagl-prevodni%20tabulky.pdf
  7. ↑ http://czepos.cuzk.cz/_souradnice.aspx
  8. ↑ http://czepos.cuzk.cz/_souradnice.aspx
  9. ↑ ftp://olggps.oeaw.ac.at/pub/EUREF_camp/ETRS89-CZ_FinalReport.pdf
  10. ↑ http://www.epncb.oma.be/_productsservices/coordinates/posvel_map.php
  11. ↑ http://transformace.webst.fd.cvut.cz/Iframe/ETRS_iframe.htm
  12. ↑ http://transformace.webst.fd.cvut.cz/Iframe/ETRS_iframe.htm
  13. ↑ https://confluence.qps.nl/pages/viewpage.action?pageId=29855173
  14. ↑ http://etrs89.ensg.ign.fr/memo-V8.pdf#page=7
  15. ↑ https://en.wikipedia.org/wiki/World_Geodetic_System#A_new_World_Geodetic_System:_WGS_84
  16. ↑ ftp://itrf.ensg.ign.fr/pub/itrf/WGS84.TXT
  17. ↑ http://www.unoosa.org/pdf/icg/2012/3.pdf#page=71
  18. ↑ https://confluence.qps.nl/pages/viewpage.action?pageId=29855173
  19. ↑ ftp://itrf.ensg.ign.fr/pub/itrf/WGS84.TXT
  20. ↑ http://www.epncb.oma.be/_productsservices/coordinates/crd4station.php?station=GOPE00CZE
  21. ↑ http://oko.pecny.cz/vesog/stanice/gope.html#coord
  22. ↑ http://www.epncb.oma.be/_productsservices/coordinates/crd4station.php?station=GOPE00CZE
  23. ↑ http://oko.pecny.cz/vesog/stanice/gope.html#coord
  24. ↑ http://czepos.cuzk.cz/_souradnice.aspx
Citováno z „http://freegis.fsv.cvut.cz/wiki/index.php?title=S-JTSK_/_Grid&oldid=5655“