Pisahkan poligon berdasarkan persimpangan menggunakan PostGIS

36

Saya memiliki tabel poligon PostGIS di mana beberapa berpotongan satu sama lain. Inilah yang saya coba lakukan:

  • Untuk poligon tertentu yang dipilih oleh id, beri saya semua poligon yang berpotongan. Pada dasarnya,select the_geom from the_table where ST_Intersects(the_geom, (select the_geom from the_table where source_id = '123'))
  • Dari poligon-poligon ini, saya perlu membuat poligon baru sehingga persimpangan menjadi poligon baru. Jadi jika poligon A berpotongan dengan poligon B, saya akan mendapatkan 3 poligon baru: A minus AB, AB, dan B minus AB.

Ada ide?

atogle
sumber
1
Hmmm, saya pikir melihat apa yang Anda maksud tetapi grafik sederhana mungkin melakukan keajaiban untuk membantu saya (dan orang lain) melihat persis apa yang Anda inginkan.
Jason
Menambahkan beberapa jawaban saya.
Adam Matan

Jawaban:

29

Karena Anda mengatakan Anda mendapatkan sekelompok poligon berpotongan untuk setiap poligon yang Anda minati, Anda mungkin ingin membuat apa yang disebut sebagai "hamparan poligon".

Ini bukan apa yang dilakukan solusi Adam. Untuk melihat perbedaannya, lihat gambar persimpangan ABC ini:

Persimpangan ABC

Saya percaya solusi Adam akan menciptakan poligon "AB" yang mencakup area "AB! C" dan "ABC", serta poligon "AC" yang mencakup "AC! B" dan "ABC", dan sebuah " BC "polygon yaitu" BC! A "dan" ABC ". Jadi poligon output "AB", "AC", dan "BC" akan tumpang tindih dengan area "ABC".

Hamparan poligon menghasilkan poligon yang tidak tumpang tindih, sehingga AB! C akan menjadi satu poligon dan ABC akan menjadi satu poligon.

Membuat overlay poligon di PostGIS sebenarnya cukup mudah.

Pada dasarnya ada tiga langkah.

Langkah 1 adalah mengekstrak garis [Perhatikan bahwa saya menggunakan cincin eksterior poligon, itu menjadi sedikit lebih rumit jika Anda ingin menangani lubang dengan benar]:

SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines

Langkah 2 adalah "simpul" garis (menghasilkan simpul di setiap persimpangan). Beberapa perpustakaan seperti JTS memiliki kelas "Noder" yang dapat Anda gunakan untuk melakukan ini, tetapi di PostGIS fungsi ST_Union melakukannya untuk Anda:

SELECT ST_Union(the_geom) AS the_geom FROM (...your lines...) AS noded_lines

Langkah 3 adalah membuat semua kemungkinan poligon yang tidak tumpang tindih yang dapat berasal dari semua baris tersebut, dilakukan oleh fungsi ST_Polygonize :

SELECT ST_Polygonize(the_geom) AS the_geom FROM (...your noded lines...)

Anda bisa menyimpan output dari masing-masing langkah ke tabel temp, atau Anda bisa menggabungkan semuanya menjadi satu pernyataan:

CREATE TABLE my_poly_overlay AS
SELECT geom FROM ST_Dump((
    SELECT ST_Polygonize(the_geom) AS the_geom FROM (
        SELECT ST_Union(the_geom) AS the_geom FROM (
            SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
        ) AS noded_lines
    )
)

Saya menggunakan ST_Dump karena output dari ST_Polygonize adalah kumpulan geometri, dan (biasanya) lebih nyaman untuk memiliki tabel di mana setiap baris adalah salah satu poligon yang membentuk lapisan poligon.

Jeff
sumber
Perhatikan bahwa ST_ExteriorRingsetiap lubang jatuh. ST_Boundaryakan mempertahankan cincin interior, tetapi juga akan membuat poligon di dalamnya.
jpmc26
Persatuan cincin eksterior mogok ketika ada terlalu banyak poligon. Sayangnya solusi "langsung" ini hanya berfungsi untuk pertanggungan kecil.
Pierre Racine
13

Jika saya mengerti dengan benar, Anda ingin mengambil (A adalah geometri kiri, B adalah kanan):

Gambar dari A∪B http://img838.imageshack.us/img838/3996/intersectab1.png

Dan ekstrak:

A ∖ AB

Gambar dari A ∖ AB http://img830.imageshack.us/img830/273/intersectab2.png

AB

Gambar AB http://img828.imageshack.us/img828/7413/intersectab3.png

dan B ∖ AB

Gambar dari B ∖ AB http://img839.imageshack.us/img839/5458/intersectab4.png

Yaitu - tiga geometri yang berbeda untuk setiap pasangan berpotongan.

Pertama, mari kita buat tampilan semua geometri berpotongan. Dengan asumsi nama tabel Anda polygons_table, kami akan menggunakan:

CREATE OR REPLACE VIEW p_intersections AS    -- Create a view with the 
SELECT t1.the_geom as t1_geom,               -- intersecting geoms. Each pair
       t2.the_geom as t2_geom                -- appears once (t2.id<t2.id)
    FROM polygons_table t1, polygons_table t2  
         WHERE t1.id<t2.id AND t1.the_geom && t2.the_geom 
                           AND intersects t1.the_geom, t2.the_geom;

Sekarang kita memiliki pandangan (praktis, tabel read-only) yang menyimpan pasangan geom berpotongan, di mana masing-masing pasangan hanya muncul satu kali karena t1.id<t2.idkondisi tersebut.

Sekarang mari kita kumpulkan persimpangan Anda - A∖AB, ABdan B∖AB, gunakan SQL UNIONdi ketiga kueri:

--AB:
SELECT ST_intersection(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION 

--AAB:
SELECT ST_Difference(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION

--BAB:
SELECT ST_Difference(t2.the_geom, t1.the_geom) 
    AS geom 
    FROM p_intersections;

Catatan:

  1. The &&operator yang digunakan sebagai filter sebelum intersectsoperator, untuk meningkatkan kinerja.
  2. Saya telah memilih untuk membuat VIEWalih - alih satu permintaan raksasa; Ini hanya untuk kenyamanan.
  3. Jika yang Anda maksud ABadalah gabungan, bukan persimpangan, dari Adan B- Gunakan ST_Union alih-alih st_intersection pada UNIONkueri di tempat yang sesuai.
  4. The tanda adalah tanda unicode untuk Set perbedaan; hapus dari kode jika membingungkan database Anda.
  5. Gambar-gambar dari kategori teori Set Wikimedia Commons yang bagus .
Adam Matan
sumber
Tiket bug saya di meta: meta.gis.stackexchange.com/questions/79/…
Adam Matan
Penjelasan yang bagus! Hasilnya sama seperti dalam solusi scw, tetapi caranya harus lebih cepat (tidak menghitung / atau menyimpan / persimpangan tambahan A dan B)
stachu
Terima kasih! Saya pikir saya tidak menyimpan informasi tambahan, karena saya hanya membuat tampilan SQL, bukan tabel.
Adam Matan
Ya, itu benar, tetapi Anda menghitung titik-temu tambahan A dan B, yang tidak diperlukan
stachu
5
Gambar dalam jawaban ini tidak berfungsi lagi.
Fezter
8

Apa yang Anda gambarkan adalah cara kerja operator Union di ArcGIS, tetapi sedikit berbeda dari Union atau titik-temu di dunia GEOS. Manual Shapely memiliki contoh cara kerja set di GEOS . Namun, wiki PostGIS memang memiliki contoh yang baik menggunakan linework yang seharusnya bisa membantu Anda.

Atau, Anda dapat menghitung tiga hal:

  1. ST_Intersection (geom_a, geom_b)
  2. ST_Difference (geom_a, geom_b)
  3. ST_Difference (geom_b, geom_a)

Itu harus menjadi tiga poligon yang Anda sebutkan di poin kedua Anda.

scw
sumber
2
Contoh wiki PostGIS bagus
fmark
Tidakkah ST_Intersect mengembalikan boolean jika mereka berpotongan atau tidak? Saya pikir ST_Intersection akan berfungsi.
Jason
Ya, salah ketik di pihak saya - diperbaiki dalam aslinya sekarang, terima kasih Jason!
scw
-2

Sesuatu seperti:

INSERT INTO new_table VALUES ((pilih id, the_geom dari old_table di mana st_intersects (the_geom, (pilih the_geom dari old_table di mana id = '123')) = true

EDIT: Anda membutuhkan persimpangan poligon yang sebenarnya.

Masukkan ke dalam nilai new_table ((pilih id, ST_Intersection (the_geom, (pilih the_geom dari yang lama di mana id = 123))

lihat apakah itu berhasil.

George Silva
sumber