¿Agrupando cadenas de líneas conectadas en PostGIS?

12

Tengo una tabla de calles que he seleccionado en función de un conjunto de atributos (digamos que es speed_limit < 25). Hay grupos de calles que son localmente contiguas; Me gustaría agrupar estos conjuntos de cadenas de líneas conectadas en GeometryCollections. En la imagen a continuación, habría dos GeometryCollections: una con las líneas rojas y otra con las azules.

ingrese la descripción de la imagen aquí

Traté de ejecutar un par de consultas de "disolver, desagregar" en la línea de:

SELECT (ST_Dump(st_union)).geom
FROM 
    (SELECT ST_Union(geom) FROM roads) sq

Con todo lo que he probado, termino con una sola característica ( ST_Union) o mi geometría original ( ST_Dumpde ST_Union).

¿Tal vez es posible hacer esto con algún tipo de WITH RECURSIVEmagia?

dbaston
fuente
Algo no se ve bien con "(ST_Dump (st_union)). Geom"
Martin F
Como no alias ST_Union (geom), el nombre de la nueva geom heredó el nombre de la función para convertirse en st_union. Es por eso que se ve un poco divertido
LR1234567

Respuestas:

19

Entonces, por ejemplo. Aquí hay una tabla simple con dos grupos de bordes conectados:

drop table lines;
create table lines ( id integer primary key, geom geometry(linestring) );
insert into lines (id, geom) values ( 1, 'LINESTRING(0 0, 0 1)');
insert into lines (id, geom) values ( 2, 'LINESTRING(0 1, 1 1)');
insert into lines (id, geom) values ( 3, 'LINESTRING(1 1, 1 2)');
insert into lines (id, geom) values ( 4, 'LINESTRING(1 2, 2 2)');
insert into lines (id, geom) values ( 11, 'LINESTRING(10 10, 10 11)');
insert into lines (id, geom) values ( 12, 'LINESTRING(10 11, 11 11)');
insert into lines (id, geom) values ( 13, 'LINESTRING(11 11, 11 12)');
insert into lines (id, geom) values ( 14, 'LINESTRING(11 12, 12 12)');
create index lines_gix on lines using gist(geom);

Ahora, aquí hay una función recursiva que, dada la identificación de un borde, acumula todos los bordes que se tocan:

CREATE OR REPLACE FUNCTION find_connected(integer) returns integer[] AS
$$
WITH RECURSIVE lines_r AS (
  SELECT ARRAY[id] AS idlist, geom, id
  FROM lines 
  WHERE id = $1
  UNION ALL
  SELECT array_append(lines_r.idlist, lines.id) AS idlist, 
         lines.geom AS geom, 
         lines.id AS id
  FROM lines, lines_r
  WHERE ST_Touches(lines.geom, lines_r.geom)
  AND NOT lines_r.idlist @> ARRAY[lines.id]
)
SELECT 
  array_agg(id) AS idlist
  FROM lines_r
$$ 
LANGUAGE 'sql';

Eso solo nos deja con la necesidad de encontrar, después de que cada grupo se acumula, la identificación de un borde que aún no es parte de un grupo. Lo cual, trágicamente, requiere una segunda consulta recursiva.

WITH RECURSIVE groups_r AS (
  (SELECT find_connected(id) AS idlist, 
          find_connected(id) AS grouplist, 
          id FROM lines WHERE id = 1)
  UNION ALL
  (SELECT array_cat(groups_r.idlist,find_connected(lines.id)) AS idlist,
         find_connected(lines.id) AS grouplist,
         lines.id
  FROM lines, groups_r
  WHERE NOT idlist @> ARRAY[lines.id]
  LIMIT 1)
)
SELECT id, grouplist
FROM groups_r;   

Lo que en conjunto devuelve un buen conjunto con el ID de semilla y cada grupo que acumuló. Lo dejo como un ejercicio para que el lector vuelva a convertir las matrices de id en una consulta para crear geometría para el mapeo.

 id |   grouplist   
----+---------------
  1 | {1,2,3,4}
 11 | {11,12,13,14}
(2 rows)
Paul Ramsey
fuente
Creo que este código podría ser más simple, si el tipo de geometría que admite el hashing en PostgreSQL (cuando escribe un RCTE más simple que no implica la acumulación de matrices de identificadores, aparece un error "Todos los tipos de datos de columna deben ser hashable"), por lo que hay un pequeña solicitud de mejora para mí.
Paul Ramsey
Este es un enfoque realmente asombroso. Noto algunos resultados extraños cuando lo aplico a un conjunto de prueba más grande; Veré si puedo reducir el problema a un simple ejemplo. 100 líneas: 85 grupos, grupo más grande = 3, 0.03 s //// 200 líneas: 144 grupos, grupo más grande = 9, 0.08 s //// 300 líneas: 180 grupos, grupo más grande = 51, 0.16 s /// / 400 líneas: 188 grupos, grupo más grande = 41, 0.27 s //// 500 líneas: 176 grupos, grupo más grande = 112, 0.56 s //// 600 líneas: 143 grupos, grupo más grande = 449, 1.0 s // // 650 líneas: 133 grupos, grupo más grande = 7601, 6,8 s
dbaston
La adición de este a los datos de prueba causará IDs duplicados en la grouplistmatriz: insert into lines (id, geom) values ( 15, 'LINESTRING(0 0, 10 10)');. El cambio array_agg(id)en la función return to array_agg(DISTINCT id)parece resolver el problema.
dbaston
Esta es una buena solución, entonces, ¿cómo podemos almacenar las geometrías en una tabla para que podamos ver las líneas conectadas?
zakaria mouqcit
6

Aquí hay un enfoque que utiliza una tabla temporal para agregar clústeres de forma incremental. Realmente no me importa el enfoque de la tabla temporal, pero parece funcionar bastante bien a medida que aumenta el número de líneas (tengo 1.2 M líneas en mi entrada).

DO
$$
DECLARE
this_id bigint;
this_geom geometry;
cluster_id_match integer;

id_a bigint;
id_b bigint;

BEGIN
DROP TABLE IF EXISTS clusters;
CREATE TABLE clusters (cluster_id serial, ids bigint[], geom geometry);
CREATE INDEX ON clusters USING GIST(geom);

-- Iterate through linestrings, assigning each to a cluster (if there is an intersection)
-- or creating a new cluster (if there is not)
FOR this_id, this_geom IN SELECT id, geom FROM lines LOOP
  -- Look for an intersecting cluster.  (There may be more than one.)
  SELECT cluster_id FROM clusters WHERE ST_Intersects(this_geom, clusters.geom)
     LIMIT 1 INTO cluster_id_match;

  IF cluster_id_match IS NULL THEN
     -- Create a new cluster
     INSERT INTO clusters (ids, geom) VALUES (ARRAY[this_id], this_geom);
  ELSE
     -- Append line to existing cluster
     UPDATE clusters SET geom = ST_Union(this_geom, geom),
                          ids = array_prepend(this_id, ids)
      WHERE clusters.cluster_id = cluster_id_match;
  END IF;
END LOOP;

-- Iterate through the clusters, combining clusters that intersect each other
LOOP
    SELECT a.cluster_id, b.cluster_id FROM clusters a, clusters b 
     WHERE ST_Intersects(a.geom, b.geom)
       AND a.cluster_id < b.cluster_id
      INTO id_a, id_b;

    EXIT WHEN id_a IS NULL;
    -- Merge cluster A into cluster B
    UPDATE clusters a SET geom = ST_Union(a.geom, b.geom), ids = array_cat(a.ids, b.ids)
      FROM clusters b
     WHERE a.cluster_id = id_a AND b.cluster_id = id_b;

    -- Remove cluster B
    DELETE FROM clusters WHERE cluster_id = id_b;
END LOOP;
END;
$$ language plpgsql;
dbaston
fuente
funciona perfectamente
zakaria mouqcit
@zakariamouqcit ¡Me alegra que esto haya funcionado para ti! Escribí esta respuesta antes de escribir la ST_ClusterIntersectingfunción en PostGIS. Si sus datos son lo suficientemente pequeños como para caber en la memoria, le sugiero que los revise para obtener una solución más eficiente.
dbaston
buscando esta pregunta me trajo aquí. Intenté la iterativa y la st_clusterintersecting pero encontré que st_clusterDBScan es la más adecuada. En caso de que alguien más sea traído aquí también. postgis.net/docs/manual-dev/ST_ClusterDBSCAN.html
D_C
De acuerdo, ST_ClusterDBSCAN es casi siempre la mejor
opción