Issue with goemetry axes order (again)

sirex commented 1 year ago

It seems, that axes order was mixed again!

Data providers should provide data using axes order as specified in EPSG:

    AXIS["northing (X)",north],
    AXIS["easting (Y)",east],

WGS84:

    AXIS["geodetic latitude (Lat)",north],
    AXIS["geodetic longitude (Lon)",east],

But since we use always_xy=True flag:

https://github.com/atviriduomenys/spinta/blob/984ecf0b920ad9eb3be3b561081deaabfc7da1a8/spinta/backends/postgresql/types/geometry/helpers.py#L25-L29

Order of axes is always changed lon, lat:

always_xy (bool, default=False) – If true, the transform method will accept as input and return as output coordinates using the traditional GIS order, that is longitude, latitude for geographic CRS and easting, northing for most projected CRS. -- https://pyproj4.github.io/pyproj/dev/api/transformer.html#pyproj.transformer.Transformer.from_crs

But currently we do not respect always_xy=True order:

https://github.com/atviriduomenys/spinta/blob/984ecf0b920ad9eb3be3b561081deaabfc7da1a8/spinta/backends/postgresql/types/geometry/helpers.py#L33-L34

And instead of lon, lat we use lat, lon.

Resources

https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent

from pyproj import CRS, Transformer
from shapely.ops import transform
from shapely.geometry import Point

def trans(x, y, s, t, **kw):
    transformer = Transformer.from_crs(
        crs_from=CRS(f'EPSG:{s}'),
        crs_to=CRS(f'EPSG:{t}'),
        **kw,
    )
    shape = transform(transformer.transform, Point(x, y))
    return shape.x, shape.y

# WGS84 (4326) -> LKS94 (3346) Bell tower of Vilnius Cathedral 
trans(54.68570, 25.28669, 4326, 3346)                  # (north, east)
#| (6061789.991147185, 582964.0913465106)              # (north, east)
trans(25.28669, 54.68570, 4326, 3346, always_xy=True)  # (east, north)
#| (582964.0913465106, 6061789.991147185)              # (east, north)

# WGS84 (4326) -> LKS94 (3346) Rubikiai railway station
trans(55.53124, 25.28203, 4326, 3346)                  # (north, east)
#| (6155887.721039969, 580936.2572643314)              # (north, east)
trans(25.28203, 55.53124, 4326, 3346, always_xy=True)  # (east, north)
#| (580936.2572643314, 6155887.721039969)              # (east, north)

# WGS84 (4326) -> WGS84 / Pseudo-Mercator (3857) Bell tower of Vilnius Cathedral
trans(54.68570, 25.28669, 4326, 3857)                  # (north, east)
#| (2814901.454647363, 7301104.288392755)              # (east, north)
trans(25.28669, 54.68570, 4326, 3857, always_xy=True)  # (east, north)
#| (2814901.454647363, 7301104.288392755)              # (east, north)

# WGS84 (4326) -> WGS84 / Pseudo-Mercator (3857) Rubikiai railway station
trans(55.53124, 25.28203, 4326, 3857)                  # (north, east)
#| (2814382.705820266, 7465659.17820406)               # (east, north)
trans(25.28203, 55.53124, 4326, 3857, always_xy=True)  # (east, north)
#| (2814382.705820266, 7465659.17820406)               # (east, north)

CRS(f'EPSG:4326')
#| Name: WGS 84
#| Axis Info [ellipsoidal]:
#| - Lat[north]: Geodetic latitude (degree)
#| - Lon[east]: Geodetic longitude (degree)

CRS(f'EPSG:3346')
#| Name: LKS94 / Lithuania TM
#| Axis Info [cartesian]:
#| - X[north]: Northing (metre)
#| - Y[east]: Easting (metre)

CRS(f'EPSG:3857')
#| Name: WGS 84 / Pseudo-Mercator
#| Axis Info [cartesian]:
#| - X[east]: Easting (metre)
#| - Y[north]: Northing (metre)

So it looks, that this is the correct way of converting coordinates:

# LKS94 (3346) -> WGS84 (4326)  Bell tower of Vilnius Cathedral
trans(6061789, 582964, 3346, 4326)                  # (north, east)
#| (54.68569111173754, 25.286688302053335)          # (north, east)

# WGS84 / Pseudo-Mercator (3857) -> WGS84 (4326)  Bell tower of Vilnius Cathedral
trans(2814901, 7301104, 3857, 4326)                 # (east, north)
#| (54.68569850243032, 25.28668591583325)           # (north, east)

Here coordinates are passed in the same order as specified in EPSG database and in the result is always presented as north, east even if input is east, north.

We can put this in OSM link.

sirex commented 2 months ago

Dar vienas paklausimas dėl ašių eiliškumo:

datasets/gov/rc/ar/graapskritis.csv

Ten nurodyta: geometry(polygon,3346)

Tada žiūrim į epsg.io: https://epsg.io/3346

Ten nurodyta:

Coordinate system: Cartesian 2D CS. Axes: northing, easting (X,Y). Orientations: north, east. UoM: m.

Kas reiškia, kad X ašis yra iš į šiaurę, o Y ašis yra į rytus.

EPSG įrašas atitinka Lietuvos teisės aktus, žiūrim į LKS teisės aktą:

https://e-seimas.lrs.lt/portal/legalAct/lt/TAD/TAIS.24435

Ten nurodyta:

stačiakampės koordinatės su x šiaurinės abscisės pradžia pusiaujuje ir y rytinės ordinatės reikšme 24°C meridiane

Žiūrim į skaiciuokle.lt:

https://skaiciuokle.lt/skaiciuokles/lks-wgs-konvertavimas

Ten prašyta:

PASTABA. Dokumente Platuma(X) ir Ilguma(Y) - turi būti nurodyta BŪTENT TOKIA TVARKA, priešingu atveju sumaišius X/Y koord. bus paskaičiuota neteisingai.

Patikslinu, kad Platuma (Latitude) yra geodezinė ašis nuo pusiaujo į šiaurę. Kol kas viskas sutampa.

Jei turim koordinates LKS-94: 375215.075, 6251267.08496094

Konvertuojam su skaiciuokle.lt į WGS-84, gaunam: 56.3782999, 21.9796704

Logiškai mąstant, LKS-94 koordinatėse, padidinus pirmą ašį, t.y. 375215.075 pakeitus į 475215.075, turėtume gauti tašką šiaurėje?

Konvertavus gaunam: 475215.075, 6251267.08496094 -> 56.3941092, 23.5985844

Kas žemėlapyje yra į rytus:

Žemėlapis pakeitus pirmą ašį

Tai reiškia, kad ašys sukeistos vietomis?

Sukeitus ašis vietomis, Saugykloje sukurtas įrankis, ant žemėlapio rodo teisingą tašką:

https://get.data.gov.lt/_srid/3346/6251267.08496094/375215.075

Saugykloje naudojama PROJ biblioteka koordinačių konvertavimui:

https://pyproj4.github.io/pyproj/stable/api/transformer.html

Kur yra toks argumentas, kurio mes nenaudojame:

always_xy (bool, default=False) – If true, the transform method will accept as input and return as output coordinates using the traditional GIS order, that is longitude, latitude for geographic CRS and easting, northing for most projected CRS.

Šio argumento nenaudojame, kadangi Saugykla tikisi, kad šaltinis ašis pateiks tokia tvarka, kokią nurodo authority, o authority šiuo atveju yra LKS teisės aktas.

PROJ dokumentacijoje yra plačiau paaiškinta apie tai, kad daugelis GIS priemonių nesilaiko standarto ir ašis pateikia rytai, šiaurė eiliškumu, nepriklausomai nuo to, kaip nurodyta konkrečioje koordinačių sistemoje, kur LKS nurodo, kad ašys turi būti atvirkščiau - šiaurė, rytai.

Toks klausimas jau buvo iškeltas, kur buvo svarstoma, ar laikomės tokios tvarkos, kurią naudoja GIS priemonės, ar laikomės standartų. Buvo pasilikta, prie standartų laikymosi.

Bet klausimas yra atviras, jei praktiškiau, nesilaikyti standarto, tai matyti galima migruoti, tik migracija nebus tokia paprasta ir užims laiko, nes tokį pokytį reikia suderinti su visomis institucijomis, kurios pateikia koordinates pagal standartą.

sirex commented 2 months ago

LKS teisės akte rašoma:

Lietuvos geodezijos praktikoje ir žemėlapių sudarymui turi būti naudojamos LKS-94 stačiakampės plokštuminės x, y koordinatės, paskaičiuotos sukertant elipsoidą GRS-80 su gulsčiu cilindru ir konformiškai suprojektuojant Lietuvos teritoriją taip, kad projekcijos iškraipymo mastelis 24°C pagrindiniame meridiane būtų 0,9998. Išvyniojus cilindrą, gaunamos stačiakampės koordinatės su x šiaurinės abscisės pradžia pusiaujuje ir y rytinės ordinatės reikšme 24°C meridiane 500 000 metrų. Skersinė kirstinė cilindrinė Merkatoriaus projekcija dvigubai sumažina 6Ų zonos linijines deformacijas, lyginant su liestine projekcija (LKS-94 stačiakampių koordinačių aprašymas pateiktas 1 priede).

Koordinačių pradžia sutampa su ašinio meridiano ir ekvatoriaus susikirtimu projekcijos plokštumoje. Ašinio meridiano projekcija yra abscisių (x) ašis. Šios ašies teigiamoji kryptis nukreipta į šiaurę. Ordinačių (y) ašies teigiamoji kryptis nukreipta į rytus. Abscisių ašies ordinatė 500 000 m.

-- https://e-seimas.lrs.lt/portal/legalAct/lt/TAD/TAIS.24435

Lietuvos centro koordinatės:

x = 6,115,464.68 metrai į šiaurę, nuo pusiaujo
y =   431,194.04 metrai į rytus, nuo 24° - 500km

Apytiksliai x ašis nuo pusiaujo turi būti apie 6000 km:

EPSG duomenų bazėje nurodytas ašių eiliškumas:

    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],

sirex commented 2 months ago

Lietuvos centro taškai:

WGS84 (ašių eiliškumas: rytai, šiaurė)

24° į rytus, nuo meridiano
55° į šiaurę, nuo pusiaujo

LKS94 (ašių eiliškumas: šiaurė, rytai)

6,116,045 m į šiaurę, nuo pusiaujo
500,000 m į rytus, nuo 24° - 500km

Paveiksliuke žemiau pavaizduota, mėlyna linija, x ašis, nuo pusiaujo į šiaurę, raudona linija pažymėta y ašis, 24° - 500km.

download

Paveiksliukas, sugeneruotas, naudojant:

import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import LineString, Point
from geopy.distance import distance
from pyproj import Geod
from pyproj import CRS, Transformer

world = gpd.read_file('https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_countries.zip')

geod = Geod(ellps='WGS84')

b = (24, 0    )  # 24° į rytus, ties pusiauju (apačia)
t = (24, 56.45)  # 24° į rytus, ties šiauriausiu Lietuvos tašku (viršus)
c = [24, 55.18]  # 24° į rytus, ties Lietuvos centru (centras)
l = geod.fwd(24, 55.18, 90, -500_000)  # 24° - 500km, ties Lietuvos centru (keirė)
r = geod.fwd(24, 55.18, 90,  500_000)  # 24° + 500km, ties Lietuvos centru (dešinė)

x = LineString([b, t])
y = LineString([l, r])

fig, ax = plt.subplots(figsize=(10, 10))
world.plot(ax=ax, edgecolor='black', color='lightgray')

gpd.GeoSeries([y]).plot(ax=ax, color='red')
gpd.GeoSeries([x]).plot(ax=ax, color='blue')

WGS84 = 4326
LKS94 = 3346
transformer = Transformer.from_crs(WGS84, LKS94)
tt = transformer.transform(*t[::-1])

print(c)
print(transformer.transform(*c[::-1]))

plt.text(t[0] + 1, t[1], f'x ({tt[0]:,})', fontsize=12, color='blue')

ax.set_xlim(0, 40)
ax.set_ylim(-10, 70)

plt.show()

sirex commented 2 months ago

Aplinkos ministerijos išaiškinimas:

Dėl reikalavimų erdviniams duomenims, teikiamiems LKS94, turime vadovautis Lietuvos koordinačių sistemos apibrėžimu ir Valstybinės koordinačių sistemos tvarkos aprašu, kuriame koordinačių eiliškumas nurodomas: x (abscisė), y (ordinatė) arba B (platuma) L (ilguma). Geoportale erdviniai duomenys vaizduojami korektiškai pagal Lietuvos koordinačių sistemą (LKS’94) – t.y. x planimetrinė koordinatė kinta į Šiaurę, y planimetrinė koordinatė kinta į Rytus. Geoportalas sukurtas atsižvelgiant į INSPIRE direktyvos ir ją įgyvendinančių EK sprendimais priimamų techninių reglamentų reikalavimus. Lietuvos erdviniai duomenų rinkiniai sėkmingai perduodami į INSPIRE portalą ir kokybiškai vaizduojami skirtingose koordinačių sistemose, nes koordinačių sistemos apibrėžtis yra nurodyta metaduomenyse ir naudotojai (programinės priemonės) turi galimybę perkonvertuoti duomenis korektiškai tarp skirtingų koordinačių sistemų.

EK parengtas koordinačių sistemų reglamentavimo bendrinis techninis dokumentas - INSPIRE Data Specification on Coordinate Reference Systems – Technical Guidelines. Čia nurodytas Geodetic Coordinate Reference Systems (CRS) eiliškumas: „set of coordinates (X, Y, Z) and/or latitude (φ), longitude (λ) and either ellipsoidal (h) or gravity-related height (H)“. Dokumente išdėstyti kelių koordinačių sistemų taikymo principai, iš kurių Lietuvoje naudojame ETRS89 Transverse Mercator. Iš principo, eiliškumas tas pats - north (lat), east (lon).

Lietuvos atvirų duomenų portalas traktuotinas kaip Informacinė sistema su savo API reikalavimais, tad vertinant iš ISO taikymų, nėra griežto formato koordinačių eiliškumui. Pagal ISO 19168-2:2022 Geographic information – Geospatial API for features — Part 2: Coordinate Reference Systems by Reference nurodoma pačiame API koks koordinačių eiliškumas. Tai rašydamas noriu pasakyti, kad informacinių sistemų tvarkytojai turi turėti galimybę interpretuoti gaunamus duomenis pagal koordinačių sistemų apibrėžtis, o ne pagal pasirinktą programinę priemonę. Svarbu, kad metaduomenys yra parengti korektiškai ir yra aiškumas.

JustinasKen commented 1 month ago

Šiek tiek pasibandžiau lokaliai, issue gali būti, kad tiksliai negalėjom pasakyti po konvertacijos ar x yra lon ar lat,

uždėjus, always_xy, biblioteka po visų konversijų grąžina x - lon, y - lat eiliškumu. Ko mums ir reikia, kadangi openstreetmap yra privaloma nurodyti /lat/lon eiliškumu (dėl ko mums ir vis apsiversdavo, nes nežinodavom ar x yra lat ar lon).

Lokaliai pakeičiau always_xy = True ir nurodžiau lon, lat = centroid.x, centroid.y Paėmiau https://get-test.data.gov.lt/datasets/gov/rc/ar/gragatve/GraGatve/:format/json duomenis.

Darant prielaidą, kad duomenys yra sukelti tinkama tvarka (pagal Lietuvos koordinačių sistemą, kai x - North, y - East) padaviau užklausą: GET localhost:8000/_srid/3346/595690.698875/6044724.154125

ir gavau nuorodą: https://www.openstreetmap.org/?mlat=54.5301553077997&mlon=25.47843196015388#map=19/54.530155/25.478432

sirex commented 1 month ago

GET localhost:8000/_srid/3346/595690.698875/6044724.154125

Pirmasis skaičius 3346 (LKS) sistemai turi būti 6044724, nes tai yra kilometrai nuo pusiaujo į šiaurę.

sirex commented 1 month ago

https://get.data.gov.lt/_srid/4326/54.6981/25.2738 https://get.data.gov.lt/_srid/3346/6063156/582111 https://get.data.gov.lt/_srid/3857/2813472/7303494 https://get.data.gov.lt/_srid/4258/54.6981/25.2738

Visos šios nuorodos turėtu rodyti į #Tower pastato kampą:

JustinasKen commented 1 month ago

https://github.com/atviriduomenys/spinta/issues/737

užduotyje buvo paminėta, kad gerai generuojam žemėlapį (prie kodo pridėjau komentarą, kad 4326 yra taip pat north, east, eiliškumu todėl

lat, lon = centroid.x, centroid.y buvo teisingas pasirinkimas.

Problema, kadangi bounding box blogai tikrindavo įrašus, atsirado nemažai įrašų, kurie yra blogai sukelti, tai yra, turime įrašų kaip POINT (400000 6000000), kai turėtų būti POINT (6000000 400000)

sirex commented 1 month ago

Reikėtu identifikuoti visus duomenų rinkinius, kuriuose yra klaidingas eiliškumas ir perduoti @aurjas, susisiektu su įstaigų atstovais ir pataisytu.

Čia matyt reikalingas kažkoks skriptas, kuris pereitu per visas lentas, kur yra geometry laukas naudojamas, patikrintu kelis taškus, ar jie yra withing bounding box ir išvesti sąrašą modelių ir geometry lauko pavadinimą, jei koordinačių ašys yra sukeistos vietomis.

Bet pradžiai reikia išleisti bounding box tikrinimo pataisymą, o tada jau atskirai žiūrėti, kurie duomenis jo neatitinka.

JustinasKen commented 1 month ago

Reikėtu identifikuoti visus duomenų rinkinius, kuriuose yra klaidingas eiliškumas ir perduoti @aurjas, susisiektu su įstaigų atstovais ir pataisytu.

Čia matyt reikalingas kažkoks skriptas, kuris pereitu per visas lentas, kur yra geometry laukas naudojamas, patikrintu kelis taškus, ar jie yra withing bounding box ir išvesti sąrašą modelių ir geometry lauko pavadinimą, jei koordinačių ašys yra sukeistos vietomis.

Bet pradžiai reikia išleisti bounding box tikrinimo pataisymą, o tada jau atskirai žiūrėti, kurie duomenis jo neatitinka.

#139 užduotis tai atliks

Grigorjeva11 commented 1 week ago

@JustinasKen Kokia tolimesne eiga?

atviriduomenys / spinta

Issue with goemetry axes order (again) #410

Resources

Related