Reprojekcja warstw w PostGIS

W dzisiejszym poście zaprezentujemy Wam jak używać funkcji reprojekcji danych wektorowych w bazie PostgreSQL/PostGIS. Wykorzystamy do tego warstwę kolej z poprzednich postów, którą możecie pobrać tutaj. W przypadku gdy ten post jest Waszą pierwszą przygodą z PostgreSQL odsyłamy do naszych wcześniejszych wpisów, gdzie pokazujemy jak zainstalować bazę i wgrać do niej dane.

Po wgraniu warstwy do bazy w pierwszej kolejności sprawdzamy jaki układ został dla niej zdefiniowany. Służy do tego funkcja ST_SRID(), gdzie między nawiasami wstawiamy kolumnę z geometrią:

SELECT ST_SRID(geom) FROM kolej LIMIT 1

W wyniku zapytania otrzymujemy informację o układzie współrzędnych (kod EPSG) dopisanym do warstwy. W naszym przypadku jest to 2180, czyli układ PUWG 1992.

Dane dla układów współrzędnych w bazie PostgreSQL/PostGIS są przechowywane w tabeli spatial_ref_sys. Korzystając z zapytań SQL możemy przeglądać rekordy tabeli, np.:

SELECT * FROM spatial_ref_sys WHERE srid = 2180

w wyniku którego otrzymujemy rekord dotyczący PUWG 1992:

2180;"EPSG";2180;"PROJCS["ETRS89 / Poland CS92",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","890 (...)";"+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "

Zapis układu w standardzie proj4 jest przechowywany w kolumnie proj4text:

SELECT proj4text FROM spatial_ref_sys WHERE srid = 2180

"+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "

Wszystko jest w porządku do momentu, gdy układ jest przypisany poprawnie. Co zrobić, gdy nie mamy przypisanego układu (wartość 0) lub jest on błędny? Wiedząc jaki jest docelowy układ warstwy sprawa jest prosta, ale co zrobić jeśli nie wiemy nic o układzie naszych danych? Musimy skorzystać z kolumny geom gdzie zapisana jest geometria obiektów i na jej podstawie spróbować określić układ. Zapytanie SQL:

SELECT geom FROM kolej LIMIT 1

Zwraca nam jednak wartość obiektów geometryczny w postaci przechowywanej w bazie danych:

01050000208408000001000000010200000002000000F3A7AE9545F32341E09704439C971E418659758004F32341A0CCAA2150971E4

Czytelną wersję zapisu obiektu otrzymamy korzystając z komendy ST_AsText()

SELECT ST_AsText(geom) FROM kolej LIMIT 1

zwracającej obiekt geometryczny w formacie WKT (Well-known Text):

MULTILINESTRING((653730.792348145 501223.065447209,653698.250895307 501204.03287811))

Wiedząc już w jakim układzie współrzędnych przedstawione są nasze dane możemy przypisać im kod EPSG korzystając z funkcji ST_SetSRID(). Funkcja ta jedynie zmienia przypisanie układu do warstwy, a nie przeprowadza reprojekcji geometrii. Widać to najlepiej gdy dodatkowo wywołamy funkcję ST_SRID():

SELECT ST_SRID(ST_SetSRID(geom,4326)) FROM kolej LIMIT 1

SRID przypisany dla warstwy to 4326. Wywołując dodatkowo ST_AsText():

SELECT ST_AsText(ST_SetSRID(geom,4326)) FROM kolej LIMIT 1

widzimy że w wyniku działania funkcji ST_SetSRID() mamy współrzędne w układzie PUWG 1992:

"MULTILINESTRING((653730.792348145 501223.065447209,653698.250895307 501204.03287811))"

Do warstwy przypisany został tylko nowy SRID. Funkcja ta nadaje się tylko do poprawy błędnego przypisania układu. Chcąc przeprowadzić reprojekcję musimy skorzystać z funkcji ST_Transform():

SELECT ST_SRID(ST_Transform(geom,4326)) FROM kolej LIMIT 1

Sprawdźmy jak będą wyglądały współrzędne obiektów warstwy po zastosowaniu dodatkowo funkcji ST_AsText():

SELECT ST_AsText(ST_Transform(geom,4326)) FROM kolej LIMIT 1

Geometria została przetransformowana do nowego układu:

"MULTILINESTRING((21.2580670000033 52.3553704999872,21.2575807000033 52.3552085999872))"

Funkcje ST_SetSRID() i ST_Tranform() stosujemy w jednym zapytaniu w przypadku, gdy chcielibyśmy przeprowadzić reprojekcję warstwy dla której nie mamy zdefiniowanego SRID:

SELECT ST_Transform(ST_SetSRID(geom,2180),4326) FROM kolej LIMIT 1

Dla słabo zaznajomionych ze składnią języka SQL niezrozumiałym na pewno jest występujące często na końcu zapytania wyrażenie LIMIT 1. Używamy go do zwrotu przez zapytanie tylko jednego wyniku. Liczba 1 określa liczbę zwracanych wyników.

W poście zapoznaliście się z czterema nowymi funkcjami z których możecie skorzystać w PostgreSQL/PostGIS:

ST_SRID(geometry geom) – zwraca układ w postaci kodu EPSG przypisany do geometrii warstwy

ST_SetSRID(geometry geom, integer srid) – przypisuje nowy układ (srid) do geometrii warstwy

ST_Transform(geometry geom, integer srid) – zwraca nową geometrię z współrzędnymi transformowany do nowego układu (srid)

ST_AsText(geometry geom) – zwraca geometrię obiektów warstwy w formacie WKT (Well-known text)

Testujcie, a jednocześnie uczcie się nowo poznanych funkcji na Waszych danych:)