PostGIS - cara efisien ST_Union semua poligon yang tumpang tindih dalam satu tabel

13

Tujuan saya adalah mengambil satu tabel dan st_union semua poligon yang menyentuh atau berdekatan satu sama lain menjadi satu poligon

Saya seorang pengembang C # yang mulai belajar tentang PostGIS. Dengan menggunakan kode di bawah ini, saya dapat melakukan ini, tetapi tampaknya tidak efisien, dan ada banyak hal untuk PostGIS yang baru bagi saya.

Dari upaya awal saya (masih dalam komentar), saya dapat mengurangi iterasi dengan menggunakan array_agg dengan ST_UNION alih-alih menyatukan hanya polys pada satu waktu.

Saya berakhir dengan 133 polis dari asal saya 173.

sql = "DROP TABLE IF Exists tmpTable; create table tmpTable ( ID varchar(50), Geom  geometry(Geometry,4326), Touchin varchar(50) ); create index idx_tmp on tmpTable using GIST(Geom); ";
                CommandText = sql;
                ExecuteNonQuery();

                sql = "";
                for (int i = 0; i < infos.Count(); i++)
                {
                    sql += "INSERT INTO tmpTable SELECT '" + infos[i].ID + "', ST_GeomFromText('" + infos[i].wkt + "', 4326), '0';";
                }
                CommandText = sql;
                ExecuteNonQuery();

                CommandText = "update tmpTable set touchin = (select id from tmpTable as t where st_intersects(st_buffer(geom, 0.0001), (select geom from tmpTable as t2 where t2.ID = tmpTable.ID ) ) and t.ID <> tmpTable.ID limit 1)";
                ExecuteNonQuery();

                CommandText = "select count(*) from tmpTable where touchin is not null";
                long touching = (long)ExecuteScalar();
                string thisId = "";
                // string otherId = "";
                while (touching > 0)
                {
                    CommandText = "select touchin, count(*)  from tmpTable where touchin is not null group by touchin order by 2 desc limit 1";
                    //CommandText = "select id, touchin from tmpTable where touchin is not null";
                    using (var prdr = ExecuteReader())
                    {
                        CommandText = "";
                        if (prdr.Read())
                        {
                            thisId = prdr.GetString(0);
                             // otherID = prdr.GetString(1);
                            CommandText = @"update tmpTable set geom = st_union(unioned) 
                                from (select array_agg(geom) as unioned from tmpTable where touchin = '" + thisId + "' or id = '" + thisId + @"') as data
                                where id = '" + thisId + "'";
                             // CommandText = "update tmpTable set geom = st_union(geom, (select geom from tmpTable where ID = '" + otherId + "')) where id = '" + thisId + "'";
                        }
                    }

                    if (!string.IsNullOrEmpty(CommandText))
                    {
                        ExecuteNonQuery();
                        //CommandText = "update tmpTable set geom = null, touchin = null where ID = '" + otherId + "'";
                        CommandText = "update tmpTable set geom = null, touchin = null where touchin = '" + thisId + "'";
                        ExecuteNonQuery();                            
                    }

                    CommandText = "update tmpTable set touchin = (select id from tmpTable as t where st_intersects(st_buffer(geom, 0.0001), (select geom from tmpTable as t2 where t2.ID = tmpTable.ID ) ) and t.ID <> tmpTable.ID limit 1)";
                    ExecuteNonQuery();

                    CommandText = "select count(*) from tmpTable where touchin is not null";
                    touching = (long)ExecuteScalar();
                }

Berikut contoh dataset yang saya gunakan:

INSERT INTO tmpTable SELECT '872538', ST_GeomFromText('POLYGON((-101.455035985 26.8835084441,-101.455035985 26.8924915559,-101.444964015 26.8924915559,-101.444964015 26.8835084441,-101.455035985 26.8835084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872550', ST_GeomFromText('POLYGON((-93.9484752173 46.0755084441,-93.9484752173 46.0844915559,-93.9355247827 46.0844915559,-93.9355247827 46.0755084441,-93.9484752173 46.0755084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872552', ST_GeomFromText('POLYGON((-116.060688575 47.8105084441,-116.060688575 47.8194915559,-116.047311425 47.8194915559,-116.047311425 47.8105084441,-116.060688575 47.8105084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872553', ST_GeomFromText('POLYGON((-116.043688832 47.8125084441,-116.043688832 47.8214915559,-116.030311168 47.8214915559,-116.030311168 47.8125084441,-116.043688832 47.8125084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872557', ST_GeomFromText('POLYGON((-80.6380222359 26.5725084441,-80.6380222359 26.5814915559,-80.6279777641 26.5814915559,-80.6279777641 26.5725084441,-80.6380222359 26.5725084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872558', ST_GeomFromText('POLYGON((-80.6520223675 26.5755084441,-80.6520223675 26.5844915559,-80.6419776325 26.5844915559,-80.6419776325 26.5755084441,-80.6520223675 26.5755084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872559', ST_GeomFromText('POLYGON((-80.6400224991 26.5785084441,-80.6400224991 26.5874915559,-80.6299775009 26.5874915559,-80.6299775009 26.5785084441,-80.6400224991 26.5785084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872560', ST_GeomFromText('POLYGON((-80.6530226307 26.5815084441,-80.6530226307 26.5904915559,-80.6429773693 26.5904915559,-80.6429773693 26.5815084441,-80.6530226307 26.5815084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872568', ST_GeomFromText('POLYGON((-90.7892258584 30.7365084441,-90.7892258584 30.7454915559,-90.7787741416 30.7454915559,-90.7787741416 30.7365084441,-90.7892258584 30.7365084441))', 4326), '0';
INSERT INTO tmpTable SELECT '872569', ST_GeomFromText('POLYGON((-90.7832259127 30.7375084441,-90.7832259127 30.7464915559,-90.7727740873 30.7464915559,-90.7727740873 30.7375084441,-90.7832259127 30.7375084441))', 4326), '0';
Carol AndorMarten Liebster
sumber
Apakah data aktual itu sendiri diperlukan dalam pertanyaan Anda?
Paul
@ Paul - tidak yakin apakah itu akan membantu atau tidak.
Carol AndorMarten Liebster

Jawaban:

20

Cara paling sederhana adalah dengan ST_Unionseluruh tabel:

SELECT ST_Union(geom) FROM tmpTable;

Ini akan memberi Anda satu besar MultiPolygon, yang mungkin bukan yang Anda inginkan. Anda bisa mendapatkan komponen terlarut individu dengan ST_Dump. Begitu:

SELECT (ST_Dump(geom)).geom FROM (SELECT ST_Union(geom) AS geom FROM tmpTable) sq;

Ini memberi Anda poligon terpisah untuk setiap set input menyentuh , tetapi kelompok input yang dipisahkan oleh jarak pendek akan tetap sebagai geometri yang terpisah. Jika Anda memiliki akses ke PostGIS 2.2.0rc1 , Anda dapat menggabungkan geometri yang berdekatan menjadi satu GeometryCollectionmenggunakan ST_ClusterWithin :

SELECT unnest(ST_ClusterWithin(geom, 0.0001)) AS grp FROM tmpTable;

Jika Anda ingin Polygonsdalam GeometryCollectionharus dibubarkan, Anda dapat menjalankan ST_UnaryUnionpada setiap GeometryCollectiondalam hasil, seperti:

SELECT ST_UnaryUnion(grp) FROM
(SELECT unnest(ST_ClusterWithin(geom, 0.0001)) AS grp FROM tmpTable) sq;
dbaston
sumber
Itu jelas jauh lebih cepat, tetapi ada 2 hal: (1) Dapatkah saya menyimpan bidang ID dalam hasil? Tidak penting yang mana, tapi saya perlu mengambil hasilnya dan mendapatkan data lain dari itu. (2) Apakah ada cara untuk menambahkan ST_Buffer kembali?
Carol AndorMarten Liebster
1
(1) Tidak mudah, tetapi cara sederhana untuk mendapatkan kembali atribut adalah dengan menggabungkan spasi hasil poligon Anda terhadap titik interior poligon input Anda. (2) Menambahkan beberapa penjelasan untuk menangani geometri yang dekat tetapi tidak menyentuh.
dbaston
Terima kasih atas bantuannya - Saat ini saya tidak memiliki 2.2, jadi saya harus mengunjungi kembali ini ketika saya meningkatkannya. Untuk saat ini pengecualian buffer bukan merupakan pemecah kesepakatan.
Carol AndorMarten Liebster
Saya setuju ini yang paling sederhana. Saya bertanya-tanya apakah akan ada cara untuk melakukan permintaan rekursif yang menemukan menyentuh geoms, tapi saya menyerah pada itu- postgresql.org/docs/current/static/queries-with.html
chrismarx
1
@ chrismarx , periksa gis.stackexchange.com/a/94228/18189 untuk beberapa ide tentang solusi rekursif.
dbaston