Intersekcja warstw w PostgreSQL\PostGIS

Wiemy już jak zainstalować PostgreSQL z rozszerzeniem PostGIS i stworzyć nową bazę z obsługą warstw GIS. W dzisiejszym wpisie pokażemy jak napisać zapytania SQL zwracające wynik intersekcji dwóch warstw wektorowych. Użyjemy do tego wykorzystywanych wcześniej warstw z liniami kolejowymi i drogami (pobierz).

Wczytajmy warstwy do naszej bazy danych korzystając z narzędzia shp2psql. Otwieramy pgAdmina i sprawdzamy, czy dane są zapisane w bazie:

Wczytane dane możemy wyświetlić w QGIS-ie (jak to zrobić pisaliśmy w jednym z poprzednich postów):

W pgAdminie rozwińmy kolumny warstwy drogi:

Warstwa drogi ma wiele atrybutów opisowych oraz atrybut geom zawierający geometrię obiektów. Wartości atrybutu zapisanego w kolumnie geom wykorzystywane są do zapytań przestrzennych w SQL.

Do otrzymania informacji o intersekcji obiektów z dwóch warstw w języku SQL musimy wykorzystać funkcję zdefiniowaną w PostGIS-ie. Funkcja ta to ST_INTERSECTION(). Jako argumenty tej funkcji podajemy atrybuty geom z dwóch tabel (warstw). Zapytanie dające w wyniku geometryczną intersekcję warstw zapisaną w tabeli przejazdy będzie miało postać:

CREATE TABLE przejazdy AS
SELECT ST_INTERSECTION(drogi.geom,kolej.geom) AS geom
FROM drogi, kolej

W wyniku zapytania otrzymujemy tabelę zawierająca 248784 rekordów. Tak duża liczba rekordów jest związana z tym, że funkcja sprawdza intersekcje każdego obiektu z każdym i zwraca wynik nawet w przypadku gdy obiekty się nie przecinają – wtedy otrzymujemy pusty obiekt geometryczny. Takiej warstwy nie będziemy mogli wyświetlić w QGIS-ie. Puste obiety moźnaby usunąć z utworzonej tabeli, ale spróbujmy zrobić to w jednym zapytaniu do którego dodamy warunek, że wynik intersekcji nie może być pusty. Do tego wykorzystamy funkcję ST_ISEMPTY():

CREATE TABLE przejazdy AS
SELECT ST_INTERSECTION(drogi.geom,kolej.geom) AS geom
FROM drogi, kolej
WHERE ST_ISEMPTY(ST_INTERSECTION(drogi.geom,kolej.geom)) = FALSE

W wyniku otrzymaliśmy tabelę z 57 rekordami. Dodajmy wynik do QGIS-a:

Przeanalizujmy poprzednie zapytanie krok po kroku:

Tworzymy nową tabelę o nazwie przejazdy zawierającą

CREATE TABLE przejazdy AS

Wynik intersekcji kolumn z geometrią pod nazwą geom

SELECT ST_INTERSECTION(drogi.geom,kolej.geom) AS geom

Z tabel drogi i kolej

FROM drogi, kolej

Gdzie wynik intersekcji nie jest pusty

WHERE ST_ISEMPTY(ST_INTERSECTION(drogi.geom,kolej.geom)) = FALSE

Funkcja ST_INTERSECTION() zwraca wynik intersekcji – w naszym przypadku punkty przecięć dróg z koleją. W przypadku gdy chcemy wybierać drogi przecinające kolej skorzystamy z innej funkcji ST_INTERSECTS() zwracającej wartości prawda/fałsz. Nasze zapytanie możemy napisać tak:

CREATE TABLE drogi_przez_kolej AS
SELECT DISTINCT(drogi.geom)
FROM drogi, kolej
WHERE ST_INTERSECTS(drogi.geom,kolej.geom)

Zapytanie zwróciło nam tabelę drogi_przez_kolej z 57 rekordami. Wyświetlamy ją w QGIS-ie:

Niektóre linie kolejowe przebiegają przez jeden i ten sam odcinek drogi, co spowodowało że niektóre odcinki drogi zostały wybrane wielokrotnie. Funkcja ST_INTERSECTS zwróciła dla nich wartość PRAWDA. Wyeliminować powtórzenia można dodając do zapytania funkcję DISTINCT o której pisaliśmy już wcześniej. Po modyfikacji nasze zapytanie wygląda następująco:

CREATE TABLE drogi_przez_kolej AS
SELECT DISTINCT(drogi.geom)
FROM drogi, kolej
WHERE ST_INTERSECTS(drogi.geom,kolej.geom)

Zapytanie zwraca 37 nie powtarzających się rekordów. Spróbujmy jeszcze na koniec zastosować dwie poznane funkcje do wyznaczenia przejazdów:

CREATE TABLE przejazdy_2 AS
SELECT DISTINCT(ST_INTERSECTION(drogi.geom,kolej.geom)) AS geom
FROM drogi, kolej
WHERE ST_INTERSECTS(drogi.geom,kolej.geom)

Czas wykonania zapytania wyniósł 330 msek. W porównaniu do poprzedniego tworzącego tabelę przejazdy, które zajęło 12 sek jest to znaczące przyspieszenie. Na podstawie tego małego porównania widzicie, że do określonego wyniku można dojść wieloma drogami (zapytaniami). Niektóre z nich mogą być kręte i długie:) Szukacie zatem zawsze optymalnego rozwiązania, testując kilka możliwych:)