CZEMU POLSKA MA 400 TYS. KM KWADRATOWYCH?
Robiąc analizy przestrzenne możecie narobić dużo różnych głupot. Na przykład dojść do tego, że Polska ma 400 tys. km kwadratowych. To taki przykład z czapy, ale dobrze ilustruje to o czym dziś napiszę. Na ilustracji powyżej widzicie dwie warstwy – czerwone prostokąty i niebieskie kwadraty. Prostokątów jest dwa – o wymiarach 1km na 2km i 1km na 4km. Kwadratów mamy trzy – każdy o boku 1km. Pytanie ogólne jest takie: jaką powierzchnię prostokątów zajmują kwadraty?
Jest to bardzo zgeneralizowany przypadek gdzie mamy przeciąć dwie warstwy przy czym jedna jest poprawna topologicznie (prostokąty nie nachodzą na siebie), a druga nie (kwadraty się nakładają). Może zdarzyć się tak, że nasze warstwy obejmują dużą powierzchnię, np. kraj i na pierwszy rzut oka są poprawne – nic się na siebie nie nakłada. Warstwę niepoprawną możemy uznać za poprawną bo np. nakładania są bardzo małe i nie widać, albo np. sprawdziliśmy tylko fragment warstwy i akurat tam ich nie znaleźliśmy. Używając takich warstw do zadania jak wyżej napisane – zliczenia jaką powierzchnię jednej warstwy zajmuje druga – w wyniku możemy dostać wartość nie do końca oczekiwaną.
Zadanie policzymy w jedynie słusznym PostGISie.
Najpierw tworzymy i zasilamy tabele. Poniżej macie gotowy kod, wystarczy uruchomić w pgAdminie. Geometria w układzie 92.
drop table if exists public.tmp_kwadraty;
drop table if exists public.tmp_prostokaty;
create table public.tmp_kwadraty (id integer, geometria geometry(Polygon, 2180));
create table public.tmp_prostokaty (id integer, geometria geometry(Polygon, 2180));
insert into public.tmp_kwadraty values (1, st_geomfromewkt('SRID=2180;POLYGON((664000 378000,
665000 378000,
665000 379000,
664000 379000,
664000 378000))'));
insert into public.tmp_kwadraty values (2, st_geomfromewkt('SRID=2180;POLYGON((664500 378000,
665500 378000,
665500 379000,
664500 379000,
664500 378000))'));
insert into public.tmp_kwadraty values (3, st_geomfromewkt('SRID=2180;POLYGON((665000 378000,
666000 378000,
666000 379000,
665000 379000,
665000 378000))'));
insert into public.tmp_prostokaty values (1, st_geomfromewkt('SRID=2180;POLYGON((663000 378500,
667000 378500,
667000 379500,
663000 379500,
663000 378500))'));
insert into public.tmp_prostokaty values (2, st_geomfromewkt('SRID=2180;POLYGON((665000 377500,
667000 377500,
667000 378500,
665000 378500,
665000 377500))'));
Po wyświetleniu w QGISie powinno wyglądać jak w nagłówku posta. Jak już mamy warstwy gotowe to piszemy zapytanie, które nam policzy co trzeba.
select t2.id, sum(st_area(st_intersection(t1.geometria, t2.geometria))) as sumpow from public.tmp_kwadraty t1, public.tmp_prostokaty t2 where st_intersects(t1.geometria, t2.geometria) group by t2.id
W zapytaniu sumujemy powierzchnie przecięcia i agregujemy po id prostokąta. Czyli w wyniku dostaniemy dwa rekordy – id prostokąta i powierzchnia kwadratów na tym prostokącie. I teraz wracamy do pytania z początku: ile powinno wyjść? Sprawdźmy.
id;sumpow; 1;1500000 2;750000
Generalnie wyszła nam głupota. Poprawnymi wartościami powinno być 1000000 i 500000. W ten właśnie sposób możemy otrzymać powierzchnię Polski równą 400tys. km kwadratowych 🙂 Otrzymaliśmy niepoprawny wynik przez nakładanie się obiektów na warstwie z kwadratami. Nie byłoby problemu gdyby kwadraty nie nakładały się. Dla dużych warstw i małych nakładań dostalibyśmy błędny wynik nawet nie zdając sobie z tego sprawy. Dlatego należy zawsze stosować metodę starych chińskich geodetów – kontrola kontroli kontroli 🙂 Ale naprawmy nasze zapytanie.
select t2.id, sum(st_area(st_intersection(t1.geometria, t2.geometria))) as sumpow from (select st_union(geometria) as geometria from public.tmp_kwadraty) t1, public.tmp_prostokaty t2 where st_intersects(t1.geometria, t2.geometria) group by t2.id
Aby było poprawnie przed przecięciem należy złączyć geometrię kwadratów w jeden obiekt – używamy do tego funkcji st_union(). W ten sposób likwidujemy nakładania. Teraz nasze nowe zapytanie zwraca takie wyniki:
id;sumpow; 1;1000000 2;500000
Zgadza się. Prawda? Po zsumowaniu policzonych powierzchni okazuje się, że w pierwszym przypadku dostaliśmy wynik o 750 tys.m kw za duży…
Pamiętajcie zatem – kontrola podstawą zaufania 🙂
