Perpotongan raster dengan poligon menggunakan PostGIS - kesalahan artefak

Saya menggunakan PostGIS2.0 untuk melakukan beberapa persimpangan raster / poligon. Saya mengalami kesulitan memahami operasi mana yang harus saya gunakan, dan apa cara tercepat untuk melakukan ini. Masalah saya adalah sebagai berikut:

• Saya punya poligon dan raster
• Saya ingin menemukan semua piksel yang termasuk dalam poligon, dan mendapatkan jumlah nilai piksel
• Dan (masalah yang diperbarui): Saya mendapatkan nilai besar untuk beberapa piksel yang tidak ada di raster asli ketika saya melakukan kueri

Saya mengalami kesulitan memahami apakah saya harus menggunakan ST_Intersects()atau ST_Intersection(). Saya juga tidak tahu apa pendekatan terbaik untuk menjumlahkan piksel saya. Inilah pendekatan pertama yang saya coba (# 1):

SELECT 
    r.rast 
FROM
    raster as r, 
    polygon as p
WHERE
    ST_Intersects(r.rast, p.geom)

Ini mengembalikan daftar rastnilai, yang saya tidak yakin apa yang harus dilakukan. Saya mencoba menghitung statistik ringkasan menggunakan ST_SummaryStats()tetapi saya tidak yakin apakah ini adalah jumlah-tertimbang dari semua piksel yang ada di dalam poligon.

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

Pendekatan lain yang saya coba (# 2) menggunakan 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;

Ini mengembalikan daftar geometri yang saya analisis lebih lanjut, tetapi saya menganggap ini kurang efisien.

Saya tidak jelas yang mana urutan operasi tercepat juga. Haruskah saya selalu memilih raster, polygonatau polygon, raster, atau mengubah poligon menjadi raster sehingga itu raster, raster?

EDIT: Saya memperbarui pendekatan # 2 dengan beberapa detail dari R.K.tautan.

Menggunakan pendekatan # 2, saya telah memperhatikan kesalahan berikut dalam hasil yang merupakan bagian dari alasan mengapa saya tidak mengerti hasilnya. Ini adalah gambar raster asli saya, dan garis besar poligon yang digunakan untuk memotongnya, dilapis di atas:

Dan di sini adalah hasil persimpangan menggunakan PostGIS:

Masalah dengan hasilnya adalah bahwa ada nilai 21474836 yang dikembalikan, yang tidak ada di raster asli. Saya tidak tahu mengapa ini terjadi. Saya menduga itu terkait dengan angka kecil di suatu tempat (membaginya dengan hampir 0), tetapi mengembalikan hasil yang salah.

Saya menemukan tutorial tentang memotong vektor vektor dengan cakupan raster yang besar menggunakan PostGIS WKT Raster . Mungkin ada jawaban yang Anda cari.

Tutorial menggunakan dua dataset, file bentuk titik yang disangga untuk menghasilkan poligon dan serangkaian raster elevasi 13 SRTM. Ada banyak langkah di antaranya tetapi kueri yang digunakan untuk memotong raster dan vektor tampak seperti ini:

 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;

Nilai-nilai itu kemudian dirangkum menggunakan yang berikut:

 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;

Saya tidak benar-benar tahu cukup PostGIS untuk menjelaskan hal ini, tetapi kedengarannya seperti apa yang Anda coba capai. Tutorial harus menjelaskan langkah-langkah perantara. Semoga berhasil :)

Berkenaan dengan poin 2 dalam pertanyaan awal - beberapa rilis pengembangan postgis 2.0 menggunakan versi pustaka GDAL yang menampilkan float ke int. Jika raster Anda memiliki nilai float di dalamnya, dan Anda menggunakan versi GDAL lebih rendah dari 1.9.0, atau versi pralelease PostGIS 2.0 yang tidak benar memanggil GDALFPolygonize (), maka Anda mungkin akan menemukan bug ini. Tiket di PostGIS dan GDAL pelacak bug diajukan dan ditutup. Bug ini aktif sekitar saat pertanyaan awal.

Dalam hal kinerja, Anda akan menemukan bahwa menggunakan ST_Intersects(raster, geom)jauh lebih cepat daripada menggunakan ST_Intersects(geom, raster). Versi pertama merasterisasi geometri dan melakukan persimpangan ruang-raster. Versi kedua membuat vektor geometri dan melakukan persimpangan ruang vektor, yang bisa menjadi proses yang jauh lebih mahal.

Saya juga mengalami masalah yang aneh menggunakan ST_SummaryStatsdengan ST_Clip. Meminta data dengan berbeda mengatakan kepada saya nilai min raster saya adalah 32, dan kemudian maks 300, namun ST_SummaryStatskembali -32700 untuk nilai-nilai piksel dalam poligon target saya.

Saya akhirnya meretas masalah ini sebagai berikut:

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