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

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:

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.

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

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.

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