Przecięcie rastra z wielokątem za pomocą PostGIS - błąd artefaktu

15

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ę rastwartoś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, polygonlub polygon, rasterprzekształ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:

wprowadź opis zdjęcia tutaj

A oto wynik skrzyżowania z PostGIS:

wprowadź opis zdjęcia tutaj

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.

djq
źródło
Jeśli chodzi o punkt drugi, chcesz uzyskać sumę wartości pikseli przecinających wielokąt?
RK
Tak, użyłem ST_SummaryStats()dla # 1, ale nie jestem pewien, jak to zrobić dla # 2.
djq
Wysłano link do referencji. Mam nadzieję że to pomogło.
RK
2
Maksymalna wartość skali na mapie to maksymalnie 32-bitowa liczba całkowita ze znakiem. Nie wiem, co to oznacza w twoim przypadku, ale może mieć to związek z wartościami NA. Zakres zapytania może mieć wartości zerowe, które nie są obsługiwane poprawnie. en.wikipedia.org/wiki/2147483647#2147483647_in_computing
yellowcap
6
FWIW, 21474836 nie jest maksymalną wartością 32-bitowego znaku int. Jednak 2 ^ 31-1 = 2147483647 to maksimum i zauważ, że 21474836 = 2147483647/100 (dzielenie liczb całkowitych). Wskazuje to, że wewnętrznie niektóre NA są generowane (być może wzdłuż komórek granicznych), są reprezentowane jako 2 ^ 31-1, a następnie kod „zapomina”, że są NA i (być może w procesie ponownego próbkowania?) Dzieli je przez 100.
whuber

Odpowiedzi:

6

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:

 CREATE TABLE caribou_srtm_inter AS
 SELECT id, 
        (gv).geom AS the_geom, 
        (gv).val
 FROM (SELECT id, 
              ST_Intersection(rast, the_geom) AS gv
       FROM srtm_tiled,
            cariboupoint_buffers_wgs
       WHERE ST_Intersects(rast, the_geom)
      ) foo;

Wartości zostały następnie podsumowane przy użyciu:

 CREATE TABLE result01 AS
 SELECT id, 
        sum(ST_Area(ST_Transform(the_geom, 32198)) * val) / 
        sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev
 FROM caribou_srtm_inter
 GROUP BY id
 ORDER BY id;

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 :)

RK
źródło
Dzięki @RK Przeczytałem ten samouczek. Myślę, że moje obliczenia są bardziej podstawowe, ale wciąż zastanawiam się nad tym podstawowym krokiem!
djq
2

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żywanie ST_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.

Pete Clark
źródło
0

Miałem też dziwne problemy ST_SummaryStatsz używaniem ST_Clip. Zapytanie o dane w inny sposób powiedziało mi, że min. Wartość mojego rastra to 32, a następnie maks. 300, ale ST_SummaryStatszwracałem -32700 dla wartości pikseli w moim docelowym wielokącie.

Skończyło się na tym, że hackowałem ten problem:

WITH first AS (
   SELECT id, (ST_Intersection(geom, rast)).val
   FROM raster_table
   INNER JOIN vector_table ON ST_Intersects(rast, geom)
)
SELECT id, COUNT(val), SUM(val), AVG(val), stddev(val), MIN(val), MAX(val)
FROM first
GROUP BY id
jczaplew
źródło