Temukan poligon hulu

8

Ini adalah pertanyaan lanjutan untuk pertanyaan ini .

Saya memiliki jaringan sungai (multiline) dan beberapa poligon drainase (lihat gambar di bawah). Tujuan saya adalah untuk memilih hanya poligon headwater (hijau).

masukkan deskripsi gambar di sini

Dengan solusi John, saya dapat dengan mudah mengekstraksi titik awal sungai (bintang). Namun, saya dapat memiliki situasi (poligon merah) di mana saya memiliki titik awal dalam poligon, tetapi poligon itu bukan poligon hulu, karena diterbangkan melalui sungai. Saya hanya ingin poligon hulu.

Saya mencoba memilihnya dengan menghitung jumlah persimpangan antara poligon dan sungai (dasar pemikiran: poligon hulu hanya memiliki 1 persimpangan dengan sungai)

SELECT 
    polyg.*
FROM 
    polyg, start_points, stream
WHERE 
    st_contains(polyg.geom, start_points.geom)
    AND ST_Npoints(ST_Intersection(poly.geom, stream.geom)) = 1

, di mana poylg adalah poylgons, start_points dari johns answer and stream adalah jaringan sungai saya.

Namun, ini berlangsung selamanya dan saya tidak menjalankannya:

"Nested Loop  (cost=0.00..20547115.26 rows=641247 width=3075)"
"  Join Filter: _st_contains(ezg.geom, start_points.geom)"
"  ->  Nested Loop  (cost=0.00..20264906.12 rows=327276 width=3075)"
"        Join Filter: (st_npoints(st_intersection(ezg.geom, rivers.geom)) = 1)"
"        ->  Seq Scan on ezg_2500km2_31467 ezg  (cost=0.00..2161.52 rows=1648 width=3075)"
"              Filter: ((st_area(geom) / 1000000::double precision) < 100::double precision)"
"        ->  Materialize  (cost=0.00..6364.77 rows=39718 width=318)"
"              ->  Seq Scan on stream_typ rivers  (cost=0.00..4498.18 rows=39718 width=318)"
"  ->  Index Scan using idx_river_starts on river_starts start_points  (cost=0.00..0.60 rows=1 width=32)"
"        Index Cond: (ezg.geom && geom)"

Jadi pertanyaan saya adalah: Bagaimana saya bisa meminta kuota hulu secara efisien?

Pembaruan: Saya menambahkan beberapa data sampel ke dropbox saya . Data berasal dari Jerman barat daya. Ini dua file bentuk - satu dengan aliran dan satu dengan poligon.

EDI
sumber
Jadi, hanya untuk menjadi jelas, Anda ingin poligon yang hanya berisi titik awal, bukan titik awal itu sendiri. Dan poin awal didefinisikan seperti pada pertanyaan Anda sebelumnya (yang saya jawab, dan sejauh yang saya tahu), dengan benar?
John Powell
Jupp, hanya poligon yang berisi titik awal DAN tidak dilewati sungai / hanya mulai dari sungai. Poligon merah di atas mengandung titik awal, tetapi BUKAN poligon hulu karena sungai mengalir melewatinya / tidak mulai dalam poligon ...
EDi
Jadi, Anda ingin himpunan polygonsyang hanya berisi titik-titik yang merupakan sumber sungai (dari pertanyaan sebelumnya) dan untuk mengecualikan di mana dua sungai bertemu. Maaf, untuk semua pertanyaan, hanya ingin memastikan.
John Powell
Tidak, misalnya di poligon hijau bawah juga dua sungai bertemu. Saya ingin mengecualikan mereka polygonsyang memiliki sungai yang lewat (sungai masuk dan meninggalkan poligon) dan menjaga mereka yang mulai (dan sungai hanya menyisakan poligon ini).
EDi
1
Saya tidak tahu PostGIS, jadi saya tidak dapat membantu dengan kode langsung, namun, di ArcGIS, saya akan mengikuti baris ini: (1) berpotongan antara garis dan poligon menjadi file titik. (2) menghapus (identik) titik-titik identik. (3) tambahkan bidang numerik ke parameter titik dengan nilai 1 untuk setiap titik. (4) spasial bergabung dengan poligon ke titik dan menggunakan jumlah bidang numerik untuk menunjukkan jenis drainase. Jumlah 1 berarti itu adalah tanjung. Lebih tinggi dari 1 berarti ada lebih dari 1 pintu masuk atau keluar.
Mikkel Lydholm Rasmussen

Jawaban:

4

Saya percaya garis besar umum (sebagian diuji sejauh ini) adalah:

  1. Temukan titik yang mewakili sumber aliran, seperti dalam jawaban ini .

  2. Berpotongan dengan tabel poligon untuk mendapatkan hitungan simpul sumber dengan poligon.

  3. Gunakan ST_DumpPoints dalam hubungannya dengan grup dengan geometri untuk mendapatkan hitungan setiap titik. Idenya adalah untuk menghitung berapa banyak sungai yang bertemu pada titik tertentu.

Contoh permintaan seperti itu:

SELECT count(geom), ST_AsText(geom) as wkt
FROM 
   (SELECT (ST_DumpPoints(foo.geom)).geom 
   FROM 
     (SELECT 
        ST_Collect(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(10,10)),
                   ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(20,20))
        ) as geom
     ) foo 
 ) bar 
 GROUP BY geom; 

yang mengembalikan:

count  |  wkt      
-------+--------------
 2     | POINT(0 0)
 1     | POINT(10 10)
 1     | POINT(20 20)
  1. Jalankan perpotongan 3terhadap tabel poligon, untuk mendapatkan hitungan (jumlah simpul) persimpangan sungai per poligon.

  2. Bergabung dengan poligon dari 2pada 4, menolak titik-titik di mana jumlah (jumlah simpul) titik di persimpangan lebih besar dari jumlah sumber sungai, diperoleh dengan menjumlahkan sumber dengan poligon dari langkah 1 dan 2. Jika kondisi ini berlaku, itu berarti bahwa setidaknya satu dari sungai bertemu di persimpangan, berasal dari luar poligon yang dimaksud.

Ini semua dapat digabungkan bersama dalam urutan CTE yang besar, kecuali beberapa tabel dibuat dari langkah-langkah yang melibatkan poin (dan diindeks).

Saya tidak tahu apa runtime ini pada set data lengkap, hanya menguji bagian ini pada subset, tetapi dengan indeks spasial pada tabel poligon, akan ada beberapa bantuan - jelas tidak mungkin untuk menerapkan indeks ke titik-titik yang muncul dari ST_DumpPoints, sehingga pemindaian penuh akan diperlukan di sana, meskipun mereka harus ada dalam memori saat itu.

Ini tidak diposting sebagai jawaban penuh , tetapi sebagai pekerjaan yang sedang berjalan, dan kesempatan untuk menemukan kekurangan dalam logika. Permintaan kerja segera hadir.

EDIT 1

Ini adalah pertanyaan yang saya buat, yang tampaknya berfungsi pada sebagian kecil dari data Anda, tetapi berjalan selama berjam-jam pada dataset lengkap.

CREATE TABLE good_polys as  
   WITH 
     rivers as 
       (SELECT (ST_DUMP(ST_LineMerge(geom))).geom as geom FROM streams),
     start_points as
       (SELECT ST_StartPoint(geom) as geom FROM rivers),
     end_points as 
        (SELECT ST_EndPoint(geom) as geom FROM rivers),
     junctions as 
        (SELECT (ST_DumpPoints(geom)).geom 
        FROM (SELECT geom FROM streams) s),
     source_polygons as 
        (SELECT 
            count(rivers.geom) as source_count, 
            polygons.geom, 
            polygons.gid 
         FROM rivers, polygons
         WHERE st_intersects(polygons.geom, rivers.geom) 
         GROUP BY polygons.geom, polygons.gid),
     junction_polygons as 
        (SELECT 
            count(junctions.geom) as junction_count, 
            polygons.geom, 
            polygons.gid 
         FROM junctions, polygons
         WHERE st_intersects(polygons.geom, junctions.geom) 
         GROUP BY polygons.geom, polygons.gid)
    SELECT 
       jp.gid 
    FROM 
       junction_polygons jp, source_polygons sp 
    WHERE ST_Equals(jp.geom, sp.geom) 
    AND junction_count <= source_count;

EDIT 2 . Sementara ini tampaknya menghasilkan jawaban yang benar pada subset kecil, waktu berjalan pada dataset lengkap mengerikan , mungkin karena permintaan terakhir melakukan perbandingan n ^ 2 dan tidak menggunakan indeks spasial. Solusi yang mungkin adalah memecah kueri dan membuat tabel dari titik awal dan titik dalam kueri poligon, yang kemudian dapat diindeks secara spasial sebelum langkah terakhir.

John Powell
sumber
Permintaan sedang berjalan di desktop saya. Saya tidak tahu berapa lama waktu yang dibutuhkan atau apakah itu benar, meskipun itu terlihat masuk akal dari sampel kecil data Anda. Apakah Anda tahu berapa banyak poligon yang memenuhi kriteria Anda?
John Powell
Saya akan menjalankan kueri di server. Saya pikir hanya sebagian kecil dari poligon yang akan memenuhi kriteria seleksi ...
EDi
Itulah yang saya temukan di subset. Saya akan memposting pertanyaan saya setelah selesai
John Powell
Penyederhanaan besok.
John Powell
Maaf, sangat sibuk hari ini. Saya pikir jawabannya adalah menjalankan kueri sumber dan persimpangan sungai queries terlebih dahulu, berpotongan dengan tabel poligon untuk mendapatkan jumlah per poligon, simpan ini sebagai tabel, dan kemudian indeks mereka. Kemudian jalankan langkah terakhir, di mana geometri sama, dan membandingkan jumlah titik dari dua tabel. Saya berharap ini kemudian akan menggunakan indeks daripada melakukan perbandingan n² sebagai hadiah. Akan dikirim kembali nanti.
John Powell
0

Dalam kode semu, ini harus bekerja:

  • pilih semua dari poligon
  • (FULL OUTER?) Gabung dengan titik pada poligon memotong titik
  • (FULL OUTER?) Gabungkan garis di mana poligon memotong garis
  • adalah line.riverid tidak sama dengan point.riverid
  • dikelompokkan berdasarkan poligonnya
  • count (pointid)> 0

Saya tidak begitu yakin bagaimana cara membangun kueri, dan saya tidak dapat mengujinya tanpa database untuk diuji. Itu pertanyaan yang cukup gila, saya pikir. Tapi itu seharusnya berhasil!

Alex Leith
sumber