KRÓTKI KURS POSTGISA

PostGIS to rozszerzenie przestrzenne do PostgreSQL, które pozwala na przechowywanie i obsługę danych przestrzennych w środowisku bazodanowym. Dodaje ono do bazy przestrzenne typy danych, a także zestaw funkcji do ich obsługi. Rozszerzenie to ma wielką zaletę – pozwala na manipulacje danymi przestrzennymi bezpośrednio w bazie. Takie rozwiązanie jest zazwyczaj szybsze i wydajniejsze niż stosowanie  middleware’u – serwera danych przestrzennych.

Aby przystąpić do kursu musisz mieć zainstalowany lub dostęp do PostgreSQL z PostGIS. Najlepiej jak najnowsze czyli aktualnie PostgreSQL w wersji 9.3 a PostGIS w wersji 2.0 lub 2.1. Tutaj powinieneś znaleźć wszystko co potrzeba. Instalacja nie jest jakoś szczególnie skomplikowana więc nie będę temu poświęcał czasu. W necie też jest to dobrze udokumentowane.

Narzędzie do obsługi bazy danych – używam pgAdmin. To narzędzie jest wykorzystywane do zarządzania bazą, wykonywania zapytań itd.

Kolejną rzeczą, której potrzebujesz, są dane. Przygotowałem zestaw danych – fragmenty warstw z podziału administracyjnego IIRP. Należy pobrać paczkę, wypakować skrypt, podłączyć się w pgAdminie do odpowiedniej bazy i uruchomić go. Jeśli wykonałeś skrypt i wszystko poszło ok, masz już w bazie dane przestrzenne.

Ostatnią rzeczą jest QGIS. Potrzebujemy go do wizualizacji naszych danych. Warto je czasem wyświetlać, żeby naocznie się przekonać czy np. wykonaliśmy operację poprawnie. Kontrola podstawą zaufania 🙂

 Podłączanie warstw w QGISie wygląda tak:

postgis_tutorial_qgis

W tym miejscu warto dodać, że QGIS na bieżąco wyświetla zmiany w naszych danych, chyba że modyfikujemy strukturę tabeli (np. dodając kolumnę), wtedy musimy ją wczytać jeszcze raz.

Zacznijmy więc.

Jak wspomniałem we wstępie PostGIS dodaje do bazy nowe typy danych. Aby nasza tabela mogła przechowywać dane przestrzenne – „geometrię” – musi, pośród innych kolumn, mieć także tę o typie „geometry”. W danych, które przygotowałem, kolumny zawierające dane przestrzenne zawsze nazywają się „geometria”.

Jeśli już zidentyfikowaliśmy kolumny z geometrią w naszych danych, musimy rozpoznać co w nich się znajduje. Możemy poznać typy obiektów jakie znajdują się w tabeli – poligony, multipoligony, linie itd. oraz układ współrzędnych w jakim są nasze dane.

--sprawdzenie typów geometrii
select distinct st_geometrytype(geometria) from public.powiaty_2rp

W wyniku dostajemy „ST_MultiPolygon” co oznacza, że w tej warstwie znajdują się multipoligony. Użycie „distinct” w zapytaniu zapewnia nam, że dostaniemy w wyniku tylko unikalne wartości. Jeśliby go nie było zapytanie zwróciłoby tyle razy „ST_MultiPolygon” ile wierszy w tabeli. Jeśli w tabeli mamy różne typy geometrii (co jest możliwe) to w wyniku będzie tyle wierszy ile typów.

--sprawdzenie w jakim układzie są nasze dane
select distinct st_srid(geometria) from public.woj_2rp

Wynikiem zapytania będzie „4326” czyli kod EPSG naszego układu współrzędnych. Wykorzystywane w Polsce układy współrzędnych mają takie kody:

Układ92 - EPSG:2180
Układ 2000 strefa 5 - EPSG:2176
Układ 2000 strefa 6 - EPSG:2177
Układ 2000 strefa 7 - EPSG:2178
Układ 2000 strefa 8 - EPSG:2179
Układ 65 strefa 1 - EPSG:3120
Układ 65 strefa 2 - EPSG:2172
Układ 65 strefa 3 - EPSG:2173
Układ 65 strefa 4 - EPSG:2174
Układ 65 strefa 5 - EPSG:2175
WGS84 - EPSG:4326
WGS84/PseudoMercator - EPSG:3857 *-w zastosowaniach internetowych

Dane geometryczne trzymane są polu w postaci binarnej czyli mało przyjazne dla człowieka. PostGIS udostępnia nam funkcje, które pozwalają zwrócić geometrię w formie jaką pożądamy – jako tekst, gml, geojson itd.

--geometria jako tekst
select st_astext(geometria) from public.powiaty_2rp where pow = 'wołkowyski'

W wyniku powyższego zapytania otrzymujemy coś takiego:

„MULTIPOLYGON(((24.4314369867809 52.737032274128,24.3844114893152 52.7630961708411,24.3408267075299 52.7571010491867,24.2928575705516 

….
52.7554494460334,24.5784145261147 52.7499455966525,24.495755474848 52.7444580708797,24.4314369867809 52.737032274128)))”

Duuużo współrzędnych. Podobny wynik ale z oznaczeniem układu w jakim są dane otrzymamy wykonując takie zapytanie:

--geometria jako well known text
select st_asewkt(geometria) from public.powiaty_2rp where pow = 'wołkowyski'

Wynik:

„SRID=4326;MULTIPOLYGON(((24.4314369867809 52.737032274128,24.3844114893152 52.7630961708411,24.3408267075299 52.7571010491867,24.2928575705516 52.7478207371638,24.2668113183597 52.7496764595448,24.2463887829129 52.7460788526398,24.2312056050685 52.7432858164406,24.2020465758344
….
52.7605281930232,24.6395751628742 52.7546697848853,24.6186991945716 52.7524244269488,24.6082925523745 52.7513051190292,24.593973122475 52.7554494460334,24.5784145261147 52.7499455966525,24.495755474848 52.7444580708797,24.4314369867809 52.737032274128)))”

Możemy też pokazać geometrię jako GML:

--geometria jako GML
select st_asgml(geometria) from public.siedziby_gmin_2rp where gmi = 'Lipsk'

Wynikiem jest punkt zapisany w formacie GML:

„<gml:Point srsName=”EPSG:4326″><gml:coordinates>23.396131334731624,53.725834102448815</gml:coordinates></gml:Point>”

Dla webmasterów może być przydatny format GeoJSON. Jest to pewna pochodna JavaScriptu, która pozwala na łatwe umieszczanie danych przestrzennych w mapach internetowych.

--geometria jako geojson
select st_asgeojson(geometria) from public.siedziby_gmin_2rp where gmi = 'Lipsk'

I wynik: „{„type”:”Point”,”coordinates”:[23.3961313347316,53.7258341024488]}”

Oczywiście działa to też w drugą stronę – możemy zamienić na geometrię format czytelny dla człowieka, np.:

--geometria z well known text
select st_geomfromewkt('SRID=4326;POINT(21.0 52.1)')

Powyższe zapytanie zamieni współrzędne punktu leżącego w okolicach Warszawy na geometrię. Wynik będzie taki:

„0101000020E61000000000000000003540CDCCCCCCCC0C4A40”

W ten sposób możemy dodać powyższy punkt do warstwy z siedzibami powiatów:

--insert geometrii
insert into public.siedziby_pow_2rp(id, gmi, geometria) values(100, 'test', st_geomfromewkt('SRID=4326;POINT(21.0 52.1)'))

Wiemy już jak się dobrać do geometrii, sprawdzić w jakim jest układzie, jakiego jest typu i wyświetlić ją w łatwo czytelnej formie.

W jednej z warstw, które przygotowałem tkwi błąd. Warstwa posiada geometrię w nieokreślonym układzie współrzędnych. Sprawdźmy to.

--sprawdzenie układu w warstwie woj_osm
select distinct st_srid(geometria) from woj_osm --0

Powyższe zapytanie zwróciło zero. Czyli nie wiemy w jakim układzie są te dane. Aby się tego dowiedzieć możemy zastosować inżynierię wsteczną i po wyglądzie współrzędnych dopasować układ. Możemy także wczytywać taką warstwę w QGISie za każdym razem ustawiając jej inny układ współrzędnych, aż nasza warstwa wyląduje w miejscu, w którym powinna być. Możemy także zapytać wróżki i trzymać kciuki żeby dane nie były w układzie lokalnym. Chciałbym w tym miejscu szczególnie pozdrowić PODGiKi – moje ulubione instytucje w tym kraju.

Wracając do tematu. Jeśli już domyśliliśmy się jaki to układ, możemy ustawić go w naszej warstwie. Robi się to poleceniem st_setsrid.

--ustawianie układu współrzędnych
update woj_osm set geometria = st_setsrid(geometria, 3857)

Jak się domyślacie drugi parametr w tym poleceniu to kod EPSG naszego układu (tego od wróżki).

Teraz zadanie podobne – transformacja do innego układu współrzędnych. Załóżmy, że chcemy przetransformować warstwę powiaty_2rp do układu 92. Żeby nie psuć aktualnej warstwy, przy okazji transformacji zapiszemy jej kopię.

--zapisujemy kopię warstwy i przy okazji robimy transformację
create table powiaty_2rp_92 as
select st_transform(geometria, 2180) as geometria, gmi, pow, woj, opis from powiaty_2rp;

Jak widać wybraliśmy wszystkie kolumny oprócz geometrii – tę przetransformowaliśmy do układu 92 (EPSG:2180) i nadaliśmy jej nową nazwę – geometria

--i sprawdzenie
select distinct st_srid(geometria) from powiaty_2rp_92 --2180;

Przeliczanie pomiędzy układami przydaje się gdy chcemy liczyć powierzchnie czy długości. Postgisowe funkcje podają miary w jednostkach układu w jakim jest warstwa. Jeśli więc mamy warstwę w WGS84 i chcemy policzyć powierzchnie obiektów to zostanie ona nam zwrócona w stopniach kwadratowych, chyba nie tego oczekujemy. Weźmy dla przykładu nasze powiaty.

--powierzchnia
select pow as powiat, st_area(geometria) as powierzchnia from powiaty_2rp order by 1
--"augustowski";0.281515874900772
--"bialski";0.277413565673557
--...

Jak widać nie wygląda to ciekawie. Na szczęście nie musimy transformować do nowego układu całej warstwy tylko po to, by policzyć powierzchnię w metrach. Możemy zrobić transformację w locie czyli zagnieździć w funkcji st_area funkcję transformującą.

--powierzchnia w hektarach, powiaty posortowane od największego
select pow as powiat, st_area(st_transform(geometria,2180))/10000 as powierzchnia from powiaty_2rp order by 2 desc
--"brzeski";464867.113018932
--"bielski";500510.119097317
--...

Wyliczmy także obwód. Weźmy jeden powiat.

--wyliczamy obwód powiatu wołkowyskiego
select pow as powiat, st_area(st_transform(geometria,2180))/10000 as powierzchnia, st_perimeter(st_transform(geometria, 2180))/1000 as obwod_w_km 
from powiaty_2rp
where pow = 'wołkowyski' 
--"wołkowyski";392767.516204334;407.861021219954

Stwórzmy teraz bufor o wielkości 1km wokół granic powiatu.

--bufor - dla wygody do osobnej warstwy
create table bufor as
select pow, st_buffer(st_transform(geometria, 2180), 5000) as geometria from powiaty_2rp where pow = 'wołkowyski'

Po wczytaniu do QGISa warstwy bufor:

postgis_tutorial_buffer

Jedną z przydatnych i często wykorzystywanych funkcji jest funkcja generująca centroidy obiektów. W PostGISie możemy wygenerować centroidy na dwa sposoby – jako zwykły centroid funkcją st_centroid lub jako centroid ale leżący zawsze wewnątrz obiektu – funkcją st_pointonsurface. Jeśli obiekt ma nieregularny kształt może się zdarzyć, że centroid będzie leżał poza obiektem. Jest to sytuacja niepożądana, wtedy należy skorzystać z drugiej funkcji – będziemy mieli pewność, że punkt leży wewnątrz obiektu.

--centroid i pointonsurface
select st_asewkt(st_centroid(geometria)) as centroid, 
st_asewkt(st_pointonsurface(geometria)) as point_on_surface 
from powiaty_2rp where pow = 'brzeski'
--"SRID=4326;POINT(23.7570801157896 52.0929433963877)";"SRID=4326;POINT(23.8564677631663 52.0580031578082)"

Jak widzimy współrzędne punktów wygenerowanych przez obie funkcje różnią się.

Na koniec sedno GISu czyli zależności przestrzenne warstw – przecinanie, nachodzenie, stykanie itd. Niech pierwsze zadanie będzie takie: znaleźć siedziby gmin przedwojennego powiatu augustowskiego, które aktualnie leżą poza granicami województwa podlaskiego (kraju).

postgis_tutorial_intersect

Aby rozwiązać to zadanie musimy znaleźć siedziby z powiatu augustowskiego i sprawdzić czy się przestrzennie przecinają z województwem podlaskim (czerwona granica).

--przecięcie - siedziby gmin z powiatu augustowskiego i województwo podlaskie
select t1.gmi 
from siedziby_gmin_2rp t1, woj_osm t2
where t1.pow = 'augustowski' and not st_intersects(t1.geometria, t2.geometria)

Jeśli wykonacie to zapytanie, czeka was niespodzianka. Zapytanie zwraca zero wierszy, chociaż jak widzimy na powyższej ilustracji, powinno zwrócić dwa. Baza nie pokazała też żadnego błędu. O co chodzi? Warstwy są w różnych układach dlatego się nie przecinają. Musimy to naprawić – obie warstwy powinny być w tym samym układzie.

--przecięcie - siedziby gmin z powiatu augustowskiego i województwo podlaskie
select t1.gmi 
from siedziby_gmin_2rp t1, woj_osm t2
where t1.pow = 'augustowski' and not st_intersects(t1.geometria, st_transform(t2.geometria, 4326))
--"Hołynka"
--"Sopoćkinie"

Zapytanie zwraca dwa rekordy, co jest zgodne z tym co widzimy na mapie. Parę słów o tym zapytaniu. Funkcja st_intersects zwraca wartość boolean – Prawda lub Fałsz. Jako parametry tej funkcji dajemy geometrie z obu warstw. Żeby wybrać siedziby gmin z powiatu augustowskiego w zapytaniu dajemy warunek pow = ‚augustowski’, póżniej warunek logiczny and bo muszą być spełnione oba warunki. Na koniec przed funkcją st_intersects dajemy zaprzeczenie bo przecież nie szukamy tych, które się przecinają, tylko takich które leżą poza granicami kraju.

Kolejne zadanie, funkcja podobna, jednak inaczej działająca. Znaleźć jaki procent powierzchni powiatu augustowskiego został w granicach Polski.

--przecięcie - jaki procent powiatu augustowskiego został w granicach kraju
select pow_czesc_wspolna_ha/pow_powiatu_ha*100 as powierzchnia_proc
from (
select st_area(st_transform(t1.geometria,2180))/10000 as pow_powiatu_ha, 
st_area(st_transform(st_intersection(t1.geometria, st_transform(t2.geometria, 4326)), 2180))/10000 as pow_czesc_wspolna_ha
from powiaty_2rp t1, woj_osm t2
where t1.pow = 'augustowski' and st_intersects(t1.geometria, st_transform(t2.geometria, 4326))
) t100
--78.8318697808748

Mamy tu dwa selekty, jeden zagnieżdżony w drugim. W selekcie wewnętrznym wybieramy odpowiedni powiat t1.pow = ‚augustowski’ i dodajemy warunek na przecinanie się geometrii – jak w poprzednim przykładzie. Dalej jesteśmy w selekcie wewnętrznym. Wyliczamy dwie wartości – pierwsza to powierzchnia powiatu – st_area(st_transform(t1.geometria,2180))/10000 as pow_powiatu_ha; druga – st_area(st_transform(st_intersection(t1.geometria, st_transform(t2.geometria, 4326)), 2180))/10000 as pow_czesc_wspolna_ha to powierzchnia części wspólnej powiatu augustowskiego i województwa podlaskiego. Należy tu wspomnieć, że funkcja st_intersection zwraca geometrię części wspólnej dwu obiektów, a nie wartość boolean jak funkcja st_intersects. Możemy to sobie zwizualizować.W tym celu stworzymy warstwę z częścią powiatu jaki pozostał w granicach kraju.

--część wspólna
create table augustowski_cz_wspolna as
select t1.pow, st_intersection(t1.geometria, st_transform(t2.geometria, 4326)) as geometria
from powiaty_2rp t1, woj_osm t2
where t1.pow = 'augustowski' and st_intersects(t1.geometria, st_transform(t2.geometria, 4326))

postgis_tutorial_intersect2

Wracając do zapytania z dwoma selektami. W zewnętrznym selekcie jest tylko wykorzystanie wyliczonych wartości w wewnętrznym selekcie. Policzony jest procent powierzchni zgodnie z zadaniem. Wynik – 78.8 procent powierzchni powiatu augustowskiego zostało w granicach Polski.

Pokazałem tu kilka najprostszych funkcji przestrzennych PostgreSQL. Kompletną ich listę wraz z instrukcjami ich użycia znajdziecie tutaj. Co do zależności przestrzennych – w PostGISie jest zaimplementowana cała gama tych funkcji, pokazałem tylko st_intersects gdyż jest to funkcja najbardziej ogólna, reszta jest jej podzbiorem.

Tak jak wspomniałem, jest to ułamek możliwości PostGISa. Osobiście uważam, że w 90% przypadków, swobodnie może zastąpić klasycznego desktopowego GISa. Myślę też, że będzie to środowisko wydajniejsze i szybsze, przy tym pozwalające na rozbudowę, skalowanie. Same plusy, tylko brać i klepać zapytania 🙂

Related Posts