PostGIS - cómo eficientemente ST_Union todos los polígonos superpuestos en una sola tabla

13

Mi objetivo es tomar una sola tabla y unir todos los polígonos que se tocan o se unen entre sí en polígonos individuales

Soy un desarrollador de C # que está comenzando a aprender sobre PostGIS. Utilizando el código a continuación, pude lograr esto, pero parece ineficiente, y hay mucho en PostGIS que es nuevo para mí.

Desde mi intento inicial (todavía en los comentarios), pude reducir las iteraciones usando array_agg con ST_UNION en lugar de unir solo polys a la vez.

Terminé con 133 polys de mi origen 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();
                }

Aquí hay una muestra del conjunto de datos que estoy usando:

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
fuente
¿Se necesitan los datos reales en su pregunta?
Paul
@Paul: no estaba seguro de si sería útil o no.
Carol AndorMarten Liebster

Respuestas:

20

La forma más sencilla sería ST_Unionla tabla completa:

SELECT ST_Union(geom) FROM tmpTable;

Esto te dará uno enorme MultiPolygon, que probablemente no sea lo que quieres. Puede obtener los componentes individuales disueltos con ST_Dump. Entonces:

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

Esto le proporciona un polígono separado para cada conjunto de entradas táctiles , pero los grupos de entradas que estaban separados por una corta distancia permanecerán como geometrías separadas. Si tiene acceso a PostGIS 2.2.0rc1 , puede fusionar geometrías que están muy juntas en una sola GeometryCollectionutilizando ST_ClusterWithin :

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

Si desea que se disuelva Polygonsdentro del GeometryCollection, puede ejecutar ST_UnaryUnioncada uno GeometryCollectionen el resultado, como:

SELECT ST_UnaryUnion(grp) FROM
(SELECT unnest(ST_ClusterWithin(geom, 0.0001)) AS grp FROM tmpTable) sq;
dbaston
fuente
Eso es definitivamente mucho más rápido, pero 2 cosas: (1) ¿Puedo preservar el campo ID en los resultados? No es importante cuál, pero necesito tomar el resultado y obtener otros datos de él. (2) ¿Hay alguna forma de volver a agregar el ST_Buffer?
Carol AndorMarten Liebster
1
(1) No es fácil, pero una forma simple de recuperar los atributos es unir espacialmente sus polígonos de resultados contra un punto interior de sus polígonos de entrada. (2) Se agregó alguna explicación para manejar geometrías que están cerca pero no se tocan.
dbaston
Gracias por la ayuda: actualmente no tengo 2.2, así que tendré que volver a visitar esto cuando actualice a eso. Por ahora, la exclusión del búfer no es un factor decisivo.
Carol AndorMarten Liebster
Estoy de acuerdo en que este es el más simple. Me pregunto si habría una manera de hacer una consulta recursiva que encuentre gemas conmovedoras, pero me di por vencido- postgresql.org/docs/current/static/queries-with.html
chrismarx
1
@chrismarx, consulte gis.stackexchange.com/a/94228/18189 para obtener algunas ideas sobre la solución recursiva.
dbaston