Wydajność w obliczaniu statystyk rastrowych w PostGIS

Próbuję obliczyć statystyki rastrowe (min., Maks., Średnia) dla każdego wielokąta w warstwie wektorowej przy użyciu PostgreSQL / PostGIS.

W odpowiedzi GIS.SE opisano, jak to zrobić, obliczając przecięcie wielokąta i rastra, a następnie obliczając średnią ważoną: https://gis.stackexchange.com/a/19858/12420

Korzystam z następującego zapytania (gdzie demjest moim rastrem, topo_area_su_regionjest moim wektorem i toidjest unikalnym identyfikatorem:

SELECT toid, Min((gv).val) As MinElevation, Max((gv).val) As MaxElevation, Sum(ST_Area((gv).geom) * (gv).val) / Sum(ST_Area((gv).geom)) as MeanElevation FROM (SELECT toid, ST_Intersection(rast, geom) AS gv FROM topo_area_su_region,dem WHERE ST_Intersects(rast, geom)) foo GROUP BY toid ORDER BY toid;

To działa, ale jest zbyt wolne. Moja warstwa wektorowa ma 2489 tys. Funkcji, z których każda zajmuje około 90 ms - przetworzenie całej warstwy zajęłoby dni . Szybkość obliczeń nie wydaje się znacznie poprawiona, jeśli obliczę tylko min i max (co pozwala uniknąć wywołań ST_Area).

Jeśli wykonam podobne obliczenia przy użyciu Pythona (GDAL, NumPy i PIL), mogę znacznie skrócić czas przetwarzania danych, jeśli zamiast wektoryzacji rastra (używając ST_Intersection) rasteryzuję wektor. Zobacz kod tutaj: https://gist.github.com/snorfalorpagus/7320167

Naprawdę nie potrzebuję średniej ważonej - podejście „jeśli dotyka, jest w” jest wystarczająco dobre - i jestem całkiem pewien, że właśnie to spowalnia.

Pytanie : Czy istnieje sposób, aby PostGIS zachowywał się w ten sposób? tzn. aby zwrócić wartości wszystkich komórek z rastra, którego dotyka wielokąt, zamiast dokładnego przecięcia.

Jestem bardzo nowy w PostgreSQL / PostGIS, więc może coś jest nie tak. Korzystam z PostgreSQL 9.3.1 i PostGIS 2.1 w systemie Windows 7 (2,9 GHz i7, 8 GB pamięci RAM) i poprawiłem konfigurację bazy danych zgodnie z sugestią tutaj: http://postgis.net/workshop/postgis-intro/tuning.html

Masz rację, ST_Intersectionzauważalne spowolnienie zapytania.

Zamiast ST_Intersectiongo używać lepiej przyciąć ( ST_Clip) swój raster z wielokątami (twoje pola) i zrzucić wynik jako wielokąty ( ST_DumpAsPolygons). Tak więc każda komórka rastrowa zostanie przekształcona w mały prostokąt wielokąta o różnych wartościach.

Do otrzymywania min, maks lub średnich ze zrzutów możesz użyć tych samych instrukcji.

To zapytanie powinno załatwić sprawę:

SELECT 
    toid,
    Min((gv).val) As MinElevation,
    Max((gv).val) As MaxElevation,
    Sum(ST_Area((gv).geom) * (gv).val) / Sum(ST_Area((gv).geom)) as MeanElevation
FROM (
    SELECT 
        toid,
        ST_DumpAsPolygons(ST_Clip(rast, 1, geom, true)) AS gv
    FROM topo_area_su_region,dem 
        WHERE ST_Intersects(rast, geom)) AS foo 
            GROUP BY toid 
            ORDER BY toid;

W instrukcji ST_Clipdefiniujesz raster, pasmo rastrowe (= 1), wielokąt i czy przycinanie powinno mieć wartość PRAWDA czy FAŁSZ.

Poza tym możesz użyć avg((gv).val)do obliczenia wartości średniej.

EDYTOWAĆ

Wynik twojego podejścia jest dokładniejszy, ale wolniejszy. Wyniki kombinacji ST_Clipi ST_DumpAsPolygonsignorują komórki rastrowe, które przecinają się z mniej niż 50% (lub 51%) ich wielkości.

Te dwa zrzuty ekranu z skrzyżowania CORINE Land Use pokazują różnicę. Pierwsze zdjęcie z ST_Intersection, drugie z ST_Clipi ST_DumpAsPolygons.