Korzystam z PostGIS2.0 do wykonywania skrzyżowań rastra / wielokąta. Mam problem ze zrozumieniem, której operacji powinienem użyć i jaki jest najszybszy sposób jej wykonania. Mój problem jest następujący:
- Mam wielokąt i raster
- Chcę znaleźć wszystkie piksele wchodzące w skład wielokąta i uzyskać sumę wartości pikseli
- I (zaktualizowany problem): Otrzymuję ogromne wartości dla niektórych pikseli, które nie istnieją w oryginalnym rastrze podczas wykonywania zapytania
Mam problem ze zrozumieniem, czy powinienem użyć, ST_Intersects()
czy ST_Intersection()
. Nie wiem też, jakie jest najlepsze podejście do sumowania moich pikseli. Oto pierwsze podejście, które próbowałem (# 1):
SELECT
r.rast
FROM
raster as r,
polygon as p
WHERE
ST_Intersects(r.rast, p.geom)
Zwraca listę rast
wartości, z którymi nie jestem pewien, co zrobić. Próbowałem obliczać statystyki podsumowujące przy użyciu, ST_SummaryStats()
ale nie jestem pewien, czy jest to suma ważona wszystkich pikseli, które znajdują się w obrębie wielokąta.
SELECT
(result).count,
(result).sum
FROM (
SELECT
ST_SummaryStats(r.rast) As result
FROM
raster As r,
polygon As p
WHERE
ST_Intersects(r.rast, p.geom)
) As tmp
Inne podejście, które wypróbowałem (# 2), wykorzystuje ST_Intersection()
:
SELECT
(gv).geom,
(gv).val
FROM
(
SELECT
ST_Intersection(r.rast, p.geom) AS gv
FROM
raster as r,
polygon as p
WHERE
ST_Intersects(r.rast, p.geom)
) as foo;
Zwraca to listę geometrii, które dalej analizuję, ale zakładam, że jest to mniej wydajne.
Nie jestem pewien, która kolejność jest najszybsza. Czy zawsze powinienem wybrać raster, polygon
lub polygon, raster
przekształcić wielokąt w raster tak, aby był raster, raster
?
EDYCJA: Zaktualizowałem podejście nr 2 z kilkoma szczegółami z R.K.
linku.
Stosując podejście nr 2, zauważyłem następujący błąd w wynikach, który jest jednym z powodów, dla których nie zrozumiałem wyniku. Oto obraz mojego oryginalnego rastra i zarys wielokąta, który jest używany do jego przecięcia, nałożony na wierzch:
A oto wynik skrzyżowania z PostGIS:
Problem z wynikiem polega na tym, że zwracane są wartości 21474836, których nie ma w oryginalnym rastrze. Nie mam pojęcia, dlaczego tak się dzieje. Podejrzewam, że jest to gdzieś związane z małymi liczbami (dzielenie przez prawie 0), ale zwraca zły wynik.
ST_SummaryStats()
dla # 1, ale nie jestem pewien, jak to zrobić dla # 2.Odpowiedzi:
Znalazłem samouczek na temat przecinania się wielokątów wektorowych z dużym pokryciem rastrowym za pomocą PostGIS WKT Raster . Może znajdziesz odpowiedzi, których szukasz.
W samouczku wykorzystano dwa zestawy danych, plik kształtu punktu, który został zbuforowany w celu wytworzenia wielokątów i serię 13 rastrów elewacji SRTM. Między nimi było wiele kroków, ale zapytanie użyte do przecięcia rastra i wektora wyglądało tak:
Wartości zostały następnie podsumowane przy użyciu:
Naprawdę nie znam wystarczającej liczby PostGIS, aby to wyjaśnić, ale z pewnością brzmi to tak, jak próbowaliście osiągnąć. Samouczek powinien rzucić światło na pośrednie kroki. Powodzenia :)
źródło
W odniesieniu do punktu 2 w pierwotnym pytaniu - kilka wydań deweloperskich Postgis 2.0 używało wersji biblioteki GDAL, która rzutowała zmiennoprzecinkowe na ints. Jeśli twój raster ma w sobie wartości zmiennoprzecinkowe, a używasz wersji GDAL niższej niż 1.9.0 lub wersji wstępnej wersji PostGIS 2.0, która nie wywoływała poprawnie GDALFPolygonize (), możesz napotkać ten błąd. Bilety w modułach śledzenia błędów PostGIS i GDAL zostały złożone i zamknięte. Ten błąd był aktywny w czasie pierwotnego pytania.
Pod względem wydajności przekonasz się, że używanie
ST_Intersects(raster, geom)
jest znacznie szybsze niż używanieST_Intersects(geom, raster)
. Pierwsza wersja rasteryzuje geometrię i przecina raster-przestrzeń. Druga wersja wektoryzuje geometrię i przecina przestrzeń wektorową, co może być znacznie droższym procesem.źródło
Miałem też dziwne problemy
ST_SummaryStats
z używaniemST_Clip
. Zapytanie o dane w inny sposób powiedziało mi, że min. Wartość mojego rastra to 32, a następnie maks. 300, aleST_SummaryStats
zwracałem -32700 dla wartości pikseli w moim docelowym wielokącie.Skończyło się na tym, że hackowałem ten problem:
źródło