CZEMU POLSKA MA 400 TYS. KM KWADRATOWYCH?

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 🙂

Related Posts