Geodéziai vonatkozási rendszerek és vetületek#
Bevezetés#
Figyelem
Ennek a bevezetőnek nem szándéka hogy teljes legyen, csupán megismétli a vetülettanból, illetve felsőgeodéziából tanultakat, valamint rendszerezi azokat oly módon, hogy megértsük a térinformatikában alkalmazott vetület definíciókat.
Egy pontot a térben természetes módon megadhatunk X, Y, Z térbeli (descartes-i) koordinátákkal. A geodéziában a térbeli koordináta rendszert a Föld középpontjához és annak forgásához kötik, ekkor geocentrikus koordináta rendszerről beszélünk, vagy külföldi szakirodalomban gyakran (ECEF) koordináta rendszerként hivatkoznak rá.
A térbeli koordináta rendszerben adott pont koordinátái nagy számok, mely számok alapján nehéz megmondani a pont körübelüli helyzetét a Földön. Ezért a három X, Y, Z koordináta helyett, a gyakorlatban gyakran alkalmazunk szög alapú (hosszúság és szélesség) koordinátákat. Ezek megadásához azonban gömb vagy ellipszoid használata szükséges. Ahhoz hogy ezt megtehessük, meg kell adnunk az ellipszoid méretére és alakjára (pl. a fél nagytengely hossza és a az ellipszoid lapultsága) vonatkozó adatokat, továbbá szükségünk van az ellipszoidnak a korábban említett térbeli koordináta rendszerben történő elhelyezéséhez és tájékozásához szükséges paraméterekre. Az ellipszoid elhelyezése lehetséges oly módon, hogy az a világ minden pontjában a legjobban közelítse a geoidot (pl. WGS84, ITRS), vagy úgy is, történhet, hogy az ellipszoid a Föld bizonyos részeit (lokálisan) jobban közelíti (pl. ETRS89, NAD83). A magasságot ezután a pont ellipszoid feletti magasságaként adjuk meg. Az így kapott koordinátákat földrajzi koordinátáknak nevezzük.
Látható, hogy az földrajzi koordináták megadásához szükségünk van egy térbeli koordináta rendszerre, így a két fajta pont reprezentáció a gyakrolatban szorosan összefügg. A térbeli koordináta rendszer leírását, valamint az ellipszoid alakját és eléhelyzését megadó parametereket együttesen geodéziai dátumnak hívjuk. A WGS84 geocentrikus koordináta rendszerében adott koordinátákat a WGS84 referencia ellipszoid segítségével átszámíthatjuk földrajzi koordinátákra.
Végül a vetületi koordinátákat az ellipszoidról, mint alapfelületről valamilyen sík vetületre történő vetítéssel kapjuk. Egy példa a vetületi rendszerre, az UTM, mely esetén a vetítés egy transzverzális elhelyezésű hengervetületre történik. Így tehát a vetítéshez szükségünk van egy alapfületre, vagy ellipszoidra, melynek elhelyzésére szükségünk van egy térbeli koordináta rendszerhez. Például a gyarkorlatban az UTM vetülethez használhatjuk a WGS84, vagy az ETRS89 dátumokat. (Az utóbbi Norvégia és egyes német tartományok hivatalos vetületei).
Ezek alapján egy pont megadása a következő módokon adható meg:
Geocentrikus koordináta: három X, Y, Z descartes-i koordánátával adható meg egy pont. Példák EPSG:4978 (WGS84), EPSG:4936 (ETRS89)
Földrajzi koordináta: egy térbeli pont megadásához a hosszúság, szélesség és esetlegesen ellipszoid vagy geoid feletti magasságot használjuk. Ahhoz hogy a hosszúság és szélességet értelmezni tudjuk, egy ellipszoidot illesztenek egy adott térbeli koordináta rendszerhez. Példák EPSG:4326 (WGS84), EPSG:4937 (ETRS89)
Vetületi koordináta: a vetületi koordináták esetén az eredmény x, y koordináta pár, melyet úgy kapunk, hogy az alapfelületről valamilyen sík vetületre vetítünk. Példák az UTM és az EOV.
Az EOV esetén az ellipszoidról egy ferde tengelyű sülyesztett hengervetületre történik a vetítés. A számítás során, a vetítés két lépésben történik: először az ellipszoidról egy simuló gömbre (az ún. új Gauss gömbre), majd onnan a ferde tengelyű sülyesztett hengervetületre. A nemzetközi szakirodalom, illetve térinformatikai rendszerek az ilyen típusú vetületet Hotine néven ismerik, lásd Snyder, Map projections: A working manual, p. 70. Megjegyezzük azonban, hogy a Hotine vetület egyenletei nem teljesen azonosak az EOV egyenleteivel, így az azzal adott megoldások csak közelítőek.
Eddig nem ejtettünk szót a magasság kezeléséről. Ismert, hogy a geodéziában a vizszintes és magassági koordinátákat gyakran külön kezeltek azok gyakorlati különbözőségei miatt. Általánosságban elmondhatjuk hogy két fajta magassági dátumot különböztethetünk meg: ellipszoid magasság, illetve tengerszint feletti (geoidi) magasság. Az ellipszoidi magasság könnyen megadható a térbeli koordinátákból, azonban a tengerszint feletti magasság megadásához ismerni kell az adott pontban a geoid unduláció értékét.
Irodalom:
Geodéziai rendszerek megadása#
import numpy as np
Ebben a leckében a pyroj könyvtárat fogjuk használni, mely lehetővé tesz koordináta-rendszerek definiálást és azok közötti átszámításokat. A könyvtárat Colabben először telepítenünk kell a következő módon:
!pip install pyproj
Requirement already satisfied: pyproj in /home/zoltan/env/uni/lib/python3.8/site-packages (3.5.0)
Requirement already satisfied: certifi in /home/zoltan/env/uni/lib/python3.8/site-packages (from pyproj) (2023.5.7)
Ezután a következő típusokat fogjuk használni a könyvtárból:
from pyproj import CRS, Transformer, Proj
Geodéziai koordináta-rendszereket különböző szabvány vagy kvázi-szabványok, konvenciók alapján tudjuk definálni. Az egyik legfontosabb ezek közül az ún. EPSG szabvány mely szinte a világ összes koordináta rendszeréhez egy számot rendel. Így például a WGS84 geocentrikus koordinátákkal, azaz X, Y, Z alakban a 4978:
crs = CRS.from_epsg(4978) # wgs84 xyz
A fenti kódban létrehozott crs változó tartalmazza a koordináta rendszert.
print(crs)
print(type(crs))
EPSG:4978
<class 'pyproj.crs.crs.CRS'>
Magáról a koordinát rendszerről további információkat tudunk lekérdezni:
print(crs.is_geocentric)
print(crs.is_geographic)
print(crs.is_projected)
True
False
False
Koordináta rendszer definíció WKT formátumban:
print(crs.to_wkt(pretty=True))
GEODCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[Cartesian,3],
AXIS["(X)",geocentricX,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(Y)",geocentricY,
ORDER[2],
LENGTHUNIT["metre",1]],
AXIS["(Z)",geocentricZ,
ORDER[3],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4978]]
Koordináta rendszer definíció proj4 formátumban:
print(crs.to_proj4())
+proj=geocent +datum=WGS84 +units=m +no_defs +type=crs
/home/zoltan/env/uni/lib/python3.8/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj = self._crs.to_proj4(version=version)
Hozzuk létre az ETRS89 “koordináta rendszert, és nézzük meg hogy, az geocentrikus.
crs = CRS.from_epsg(4936) # ETRS89 XYZ
print(crs.is_geocentric)
print(crs.is_geographic)
print(crs.is_projected)
True
False
False
A WGS84 rendszer földrajzi koordinátákkal adott változatát a 4326-os kóddal jelöli az EPSG:
crs = CRS.from_epsg(4326) # wgs84 latlon
print(crs.is_geocentric)
print(crs.is_geographic)
print(crs.is_projected)
False
True
False
Nézzük meg az ehhez tartozó WKT definíciót és hasonlítsuk össze a geocentrikus 4978-as változattal:
print(crs.to_wkt(pretty=True))
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Az EOV-t a 23700-as kóddal találjuk:
crs_eov = CRS.from_epsg(23700) # eov
print(crs_eov.is_geocentric)
print(crs_eov.is_geographic)
print(crs_eov.is_projected)
False
False
True
A vetületi koordináta rendszerek esetén, létezik egy alapfelület, melyet le tudunk kérdeni a következő módon:
crs_eov.geodetic_crs
<Geographic 2D CRS: EPSG:4237>
Name: HD72
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Hungary.
- bounds: (16.11, 45.74, 22.9, 48.58)
Datum: Hungarian Datum 1972
- Ellipsoid: GRS 1967
- Prime Meridian: Greenwich
print(crs_eov.to_wkt(pretty=True))
PROJCRS["HD72 / EOV",
BASEGEOGCRS["HD72",
DATUM["Hungarian Datum 1972",
ELLIPSOID["GRS 1967",6378160,298.247167427,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4237]],
CONVERSION["Egyseges Orszagos Vetuleti",
METHOD["Hotine Oblique Mercator (variant B)",
ID["EPSG",9815]],
PARAMETER["Latitude of projection centre",47.1443937222222,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8811]],
PARAMETER["Longitude of projection centre",19.0485717777778,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8812]],
PARAMETER["Azimuth of initial line",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8813]],
PARAMETER["Angle from Rectified to Skew Grid",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8814]],
PARAMETER["Scale factor on initial line",0.99993,
SCALEUNIT["unity",1],
ID["EPSG",8815]],
PARAMETER["Easting at projection centre",650000,
LENGTHUNIT["metre",1],
ID["EPSG",8816]],
PARAMETER["Northing at projection centre",200000,
LENGTHUNIT["metre",1],
ID["EPSG",8817]]],
CS[Cartesian,2],
AXIS["easting (Y)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (X)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Engineering survey, topographic mapping."],
AREA["Hungary."],
BBOX[45.74,16.11,48.58,22.9]],
ID["EPSG",23700]]
print(crs_eov.to_proj4())
+proj=somerc +lat_0=47.1443937222222 +lon_0=19.0485717777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67 +units=m +no_defs +type=crs
Transzformációk#
A BME permanens mérőállomás koordinátáit fogjuk használni a továbbiakban, amelynek koordinátái a következőek:
# %%
# BUTE GPS Point: https://www.gnssnet.hu/pdf/BUTE.pdf
y_eov = 650684.479 # easting
x_eov = 237444.175 # northing
z_eov = 137.177 # EOMA height
Végezzünk transzformációt az EOV és a WGS84 rendszer között. Az eredményt földrajzi koordinátákban szereténk megkapni. Ehhez először definiáljuk a két rendszert a korábbiaknak megfelelően:
crs_eov = CRS.from_epsg(23700)
crs_wgs84_llh = CRS.from_epsg(4326)
Majd a pyroj Transformer osztályát használjuk hogy egy Transformer típust hozzunk létre:
transformer = Transformer.from_crs(crs_eov, crs_wgs84_llh)
print(type(transformer))
<class 'pyproj.transformer.Transformer'>
Ezt a típust használhatjuk utána a transzformáció elvégzésére:
wgs_xyz = transformer.transform(y_eov, x_eov, z_eov)
print(wgs_xyz)
(47.480942681093, 19.05652711329527, 137.177)
Az eredmény egy tuple típus, mely a WGS84 földrajzi koordinátákat tartalmazza.
Hasonló módon végezzük el a transzformációt az EOV és a ETRS89 rendszerek között. Adjuk meg a földrajzi és geocentrikus koordinátákat is.
crs_eov = CRS.from_epsg(23700) # HD72 / EOV
crs_etrs89_llh = CRS.from_epsg(4258) # ETRS89/GRS80, lat/lon
crs_etrs89_xyz = CRS.from_epsg(4936) # ETRS89, ECEF-XYZ
# %%
eov2etrs89_xyz = Transformer.from_crs(crs_eov, crs_etrs89_xyz)
eov2etrs89_llh = Transformer.from_crs(crs_eov, crs_etrs89_llh)
# %%
# BUTE GPS Point: https://www.gnssnet.hu/pdf/BUTE.pdf
y_eov = 650684.479 # easting
x_eov = 237444.175 # northing
z_eov = 137.177 # EOMA height
p_etrs89_xyz = eov2etrs89_xyz.transform(y_eov, x_eov, z_eov)
p_etrs89_llh = eov2etrs89_llh.transform(y_eov, x_eov, z_eov)
print(p_etrs89_xyz)
print(p_etrs89_llh)
(4081854.584760393, 1410001.4285596195, 4678167.186847403)
(47.48094314528847, 19.05652833851546, 137.177)
Műveletek az ellipszoidon#
from geopy.distance import geodesic
# %% Geodetic distance between two points defined with lat/lon
coords_1 = (52.2296756, 21.0122287)
coords_2 = (52.406374, 16.9251681)
sol_1 = geodesic(coords_1, coords_2, ellipsoid='GRS-80').m * 100
sol_2 = geodesic(coords_1, coords_2, ellipsoid='WGS-84').m * 100
print('Difference GRS-80 vs. WGS-84 [cm]: ', sol_1 - sol_2)
Difference GRS-80 vs. WGS-84 [cm]: 0.00028607621788978577
Transzformációs gridek#
# %% Comparison to eht2: http://eht.gnssnet.hu/index.php
p_eht_etrs89_xyz = np.array([4081882.378, 1410011.142, 4678199.391])
p_eht_etrs89_llh = np.array([47.4809437284, 19.0565297521, 180.811])
print(p_etrs89_xyz - p_eht_etrs89_xyz)
[-27.79323961 -9.71344038 -32.2041526 ]
# %%
print('Difference XY [cm]: ', geodesic(p_eht_etrs89_llh[0:2], p_etrs89_llh[0:2]).m * 100)
# Note: the Z of p_etrs89_llh[2] is the same as input, no Z coordinate transformation
# print('Difference Z [m]: ', (p_eht_etrs89_llh[2] - p_etrs89_llh[2]))
Difference XY [cm]: 12.471731150364025
# %% Original definiation of HD72/EOV in proj4 format
print('Original definiation:')
print(crs_eov.to_proj4())
Original definiation:
+proj=somerc +lat_0=47.1443937222222 +lon_0=19.0485717777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67 +units=m +no_defs +type=crs
/home/zoltan/env/uni/lib/python3.8/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj = self._crs.to_proj4(version=version)
# %%
# Grid letölthető: https://github.com/OSGeoLabBp/eov2etrs
# < transzformáció 1 cm-es középhibán belül azonos az EHT-vel
# az EHT transzformációja 10 cm-es középhibájú
nadgrids = './data/etrs2eov_notowgs.gsb'# XY grid
geoidgrids = './data/geoid_eht2014.gtx' # geoid model for Z transformation
crs_eov_grid = Proj(init='EPSG:23700', nadgrids=nadgrids, geoidgrids=geoidgrids) # HD72 / EOV
print(crs_eov_grid.to_proj4())
+proj=somerc +lat_0=47.1443937222222 +lon_0=19.0485717777778 +k_0=0.99993 +x_0=650000 +y_0=200000 +ellps=GRS67 +nadgrids=./data/etrs2eov_notowgs.gsb +units=m +geoidgrids=./data/geoid_eht2014.gtx +geoid_crs=WGS84 +vunits=m +no_defs
/home/zoltan/env/uni/lib/python3.8/site-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
in_crs_string = _prepare_from_proj_string(in_crs_string)
# %% XYZ transformation
eov2etrs89_xyz = Transformer.from_proj(crs_eov_grid, crs_etrs89_xyz)
p_etrs89_xyz = eov2etrs89_xyz.transform(y_eov, x_eov, z_eov)
print(p_etrs89_xyz)
print(p_eht_etrs89_xyz)
print('Difference [cm]: ', np.linalg.norm(p_etrs89_xyz - p_eht_etrs89_xyz) * 100)
(4081882.3770687627, 1410011.1433535905, 4678199.3908524215)
[4081882.378 1410011.142 4678199.391]
Difference [cm]: 0.16496028796703496
# %% LLH transformation:
eov2etrs89_llh = Transformer.from_proj(crs_eov_grid, crs_etrs89_llh)
p_etrs89_llh = eov2etrs89_llh.transform(y_eov, x_eov, z_eov)
print('Difference [cm]: ', geodesic(p_etrs89_llh[0:2], p_eht_etrs89_llh[0:2]).m * 100)
print('Difference Z [cm]: ', (p_etrs89_llh[2] - p_eht_etrs89_llh[2])*100)
Difference [cm]: 0.16408611991976702
Difference Z [cm]: -4363.4000000000015
print(p_eht_etrs89_llh)
[ 47.48094373 19.05652975 180.811 ]