¿Dibujar una línea entre puntos a una distancia específica en PostGIS?

9

Tengo datos de puntos a lo largo de las calles, me gustaría convertir esos puntos en líneas de colores simples. ¿Algún indicador de cómo se puede llamar este problema o algún algoritmo que pueda ayudarme a resolver esto? Puntos a lo largo de la calle que me gustaría convertir en líneas.

Esperaba usar PostGISfunciones para hacer esto, pero estoy abierto a sugerencias, estos son datos de un .shparchivo.

Edit1: se actualizó la imagen para demostrar la solución ideal de este problema.

Dibujar la línea se basaría puramente en la distancia entre esos puntos, no hay nada más que pueda usar para agruparlos. Idealmente, ¿serían puntos a la distancia máxima especificada a lo largo de la línea proyectada? Y por línea proyectada me refiero a encontrar el primer punto, luego el próximo más cercano a él, luego proyectar una línea y verificar si hay puntos en esta línea a una distancia máxima de cualquiera de los que ya están en la línea.

Mahakala
fuente
1
¿Qué software estás planeando usar?
ArMoraer
¿Estás tratando de convertir esto en aceras?
DPSEspacial
Esperaba usar las funciones de PostGIS para hacer esto, pero estoy abierto a sugerencias, estos son datos de un archivo .shp.
Mahakala
1
¿Podría mostrar exactamente qué puntos desea conectar en su dibujo o en otro dibujo? ¿Son solo dos puntos a la vez? O tres? ¿La distancia entre los puntos que deberían estar conectados es siempre la misma o es "justo" por debajo de cierto umbral?
Peter Horsbøll Møller
1
Muchas gracias a @dbaston y MarHoff, no tendré tiempo para probar sus ideas hasta finales de abril, desearía poder dividir la recompensa entre ellos, pero tendré que otorgarles esto a uno de ustedes y dbaston me dio algunas consultas para investigar entonces aceptaré su respuesta. Gracias a todos los que se tomaron el tiempo para responder! Gran comunidad para ser parte de :-)
Mahakala

Respuestas:

8

Puede usar una consulta recursiva para explorar el vecino más cercano de cada punto a partir de cada extremo de líneas detectado que desea construir.

Prerrequisitos : prepare una capa postgis con sus puntos y otra con un solo objeto Multi-linetring que contenga sus carreteras. Las dos capas deben estar en el mismo CRS. Aquí está el código para el conjunto de datos de prueba que creé, modifíquelo según sea necesario. (Probado en postgres 9.2 y postgis 2.1)

WITH RECURSIVE
points as (SELECT id, st_transform((st_dump(wkb_geometry)).geom,2154) as geom, my_comment as com FROM mypoints),
roads as (SELECT st_transform(ST_union(wkb_geometry),2154) as geom from highway),

ingrese la descripción de la imagen aquí

Aquí están los pasos :

  1. Genere para cada punto la lista de cada vecino y su distancia que cumpla estos tres criterios.

    • La distancia no debe exceder un umbral definido por el usuario (esto evitará la vinculación a un punto aislado) ingrese la descripción de la imagen aquí
      graph_full as (
      SELECT a.id, b.id as link_id, a.com, st_makeline(a.geom,b.geom) as geom, st_distance(a.geom,b.geom) as distance
      FROM points a
      LEFT JOIN points b ON a.id<>b.id
      WHERE st_distance(a.geom,b.geom) <= 15
      ),
    • El camino directo no debe cruzar una carretera ingrese la descripción de la imagen aquí
      graph as (
      SELECt graph_full.*
      FROM graph_full RIGHT JOIN
      roads ON st_intersects(graph_full.geom,roads.geom) = false
      ),
    • La distancia no debe exceder una relación definida por el usuario de la distancia desde el vecino más cercano (esto debería acomodarse mejor a la digitalización irregular que la distancia fija) Esta parte fue realmente demasiado difícil de implementar, pegada al radio de búsqueda fijo

    Llamemos a esta tabla "el gráfico"

  2. Seleccione el punto de final de línea uniéndose al gráfico y manteniendo solo el punto que tenga exactamente una entrada en el gráfico. ingrese la descripción de la imagen aquí

    eol as (
    SELECT points.* FROM
    points  JOIN
    (SELECT id, count(*) FROM graph 
    GROUP BY id
    HAVING count(*)= 1) sel
    ON points.id = sel.id),

    Llamemos
    fácil a esta tabla "eol" (final de línea) ? que la recompensa por hacer un gran gráfico pero que las cosas se vuelvan locas en el próximo paso

  3. Configure una consulta recursiva que irá en ciclo de vecinos a vecinos a partir de cada eol ingrese la descripción de la imagen aquí

    • Inicialice la consulta recursiva usando la tabla eol y agregando un contador para la profundidad, un agregador para la ruta y un constructor de geometría para construir las líneas.
    • Pase a la siguiente iteración cambiando al vecino más cercano usando el gráfico y verificando que nunca retroceda usando la ruta
    • Una vez finalizada la iteración, conserve solo la ruta más larga para cada punto de partida (si su conjunto de datos incluye una posible intersección entre las líneas esperadas, esa parte necesitaría más condiciones)
    recurse_eol (id, link_id, depth, path, start_id, geom) AS (--initialisation
    SELECT id, link_id, depth, path, start_id, geom FROM (
        SELECT eol.id, graph.link_id,1 as depth,
        ARRAY[eol.id, graph.link_id] as path,
        eol.id as start_id,
        graph.geom as geom,
        (row_number() OVER (PARTITION BY eol.id ORDER BY distance asc))=1 as test
        FROM eol JOIn graph ON eol.id = graph.id 
        ) foo
    WHERE test = true
    
    UNION ALL ---here start the recursive part
    
    SELECT id, link_id, depth, path, start_id, geom  FROM (
        SELECT graph.id, graph.link_id, r.depth+1 as depth,
        path || graph.link_id as path,
        r.start_id,
        ST_union(r.geom,graph.geom) as geom,
        (row_number() OVER (PARTITION BY r.id ORDER BY distance asc))=1 as test
        FROM recurse_eol r JOIN graph ON r.link_id = graph.id AND NOT graph.link_id = ANY(path)) foo
    WHERE test = true AND depth < 1000), --this last line is a safe guard to stop recurring after 1000 run adapt it as needed

    Llamemos a esta tabla "recurse_eol"

  4. Mantenga solo la línea más larga para cada punto de inicio y elimine todas las rutas duplicadas exactas Ejemplo: las rutas 1,2,3,5 Y 5,3,2,1 son la misma línea descubierta por sus dos "finales de línea" diferentes

    result as (SELECT start_id, path, depth, geom FROM
    (SELECT *,
    row_number() OVER (PARTITION BY array(SELECT * FROM unnest(path) ORDER BY 1))=1 as test_duplicate,
    (max(depth) OVER (PARTITION BY start_id))=depth as test_depth
    FROM recurse_eol) foo
    WHERE  test_depth = true AND test_duplicate = true)
    
    SELECT * FROM result
  5. Comprueba manualmente los errores restantes (puntos aislados, líneas superpuestas, calles con formas extrañas)


Actualizado como se prometió, todavía no puedo entender por qué a veces la consulta recursiva no da exactamente el mismo resultado al comenzar desde el extremo opuesto de una misma línea, por lo que puede haber algún duplicado en la capa de resultados a partir de ahora.

Siéntase libre de preguntar, entiendo totalmente que este código necesita más comentarios. Aquí está la consulta completa:

WITH RECURSIVE
points as (SELECT id, st_transform((st_dump(wkb_geometry)).geom,2154) as geom, my_comment as com FROM mypoints),
roads as (SELECT st_transform(ST_union(wkb_geometry),2154) as geom from highway),

graph_full as (
    SELECT a.id, b.id as link_id, a.com, st_makeline(a.geom,b.geom) as geom, st_distance(a.geom,b.geom) as distance
    FROM points a
    LEFT JOIN points b ON a.id<>b.id
    WHERE st_distance(a.geom,b.geom) <= 15
    ),

graph as (
    SELECt graph_full.*
    FROM graph_full RIGHT JOIN
    roads ON st_intersects(graph_full.geom,roads.geom) = false
    ),

eol as (
    SELECT points.* FROM
    points  JOIN
        (SELECT id, count(*) FROM graph 
        GROUP BY id
        HAVING count(*)= 1) sel
    ON points.id = sel.id),


recurse_eol (id, link_id, depth, path, start_id, geom) AS (
    SELECT id, link_id, depth, path, start_id, geom FROM (
        SELECT eol.id, graph.link_id,1 as depth,
        ARRAY[eol.id, graph.link_id] as path,
        eol.id as start_id,
        graph.geom as geom,
        (row_number() OVER (PARTITION BY eol.id ORDER BY distance asc))=1 as test
        FROM eol JOIn graph ON eol.id = graph.id 
        ) foo
    WHERE test = true

UNION ALL
    SELECT id, link_id, depth, path, start_id, geom  FROM (
        SELECT graph.id, graph.link_id, r.depth+1 as depth,
        path || graph.link_id as path,
        r.start_id,
        ST_union(r.geom,graph.geom) as geom,
        (row_number() OVER (PARTITION BY r.id ORDER BY distance asc))=1 as test
        FROM recurse_eol r JOIN graph ON r.link_id = graph.id AND NOT graph.link_id = ANY(path)) foo
    WHERE test = true AND depth < 1000),

result as (SELECT start_id, path, depth, geom FROM
    (SELECT *,
    row_number() OVER (PARTITION BY array(SELECT * FROM unnest(path) ORDER BY 1))=1 as test_duplicate,
    (max(depth) OVER (PARTITION BY start_id))=depth as test_depth
    FROM recurse_eol) foo
WHERE  test_depth = true AND test_duplicate = true)

SELECT * FROM result
MarHoff
fuente
Hola @MarHoff, gracias por tu respuesta, tengo algo que buscar. No esperaba una solución completa, solo consejos sobre dónde buscar respuestas. Quiero entender esto más y seguiré cavando y probablemente tenga más preguntas más adelante. Necesito entender su algoritmo y esto me llevará algún tiempo de todos modos :)
Mahakala
Tengo un script de trabajo, previsualice aquí qgiscloud.com/MarHoff/test_qgiscloud_bis. Queda una pequeña advertencia para la desduplicación ... No más recompensas, no más presión, supongo, así que liberaré el lanzamiento cuando pueda. Sin embargo
MarHoff
gracias @MarHoff, si pudiera hubiera dividido esta recompensa, no puedo ver cómo puedo otorgarle ningún punto, pero muchas gracias por investigar esto y su prueba. Parece genuino :)
Mahakala
Hecho. Gracias por el rompecabezas, y perdón por despotricar. Si otra respuesta lo hizo por usted, entonces está totalmente bien, en algún momento simple es lo mejor ... Mi respuesta fue tal vez un poco pensar demasiado. Aunque es un buen ejemplo del uso de CTE + consulta recursiva + función de Windows + postgis en una sola consulta;)
MarHoff
8

Como señala @FelixIP, el primer paso es encontrar los puntos que formarán cada línea. Puede hacerlo llamando a ST_ClusterWithin con su distancia de separación máxima:

SELECT
  row_number() OVER () AS cid, 
  (ST_Dump(geom)).geom 
FROM (
  SELECT unnest(st_clusterwithin(geom, 0.05)) AS geom 
  FROM inputs) sq

Luego, necesitará usar algo de heurística para construir una línea a través de todos los puntos en cada grupo. Por ejemplo, si puede suponer que las líneas deseadas son monótonas en Y, puede ordenar los puntos en cada grupo y alimentarlos en ST_MakeLine . Combinar eso todos juntos se vería así:

SELECT 
  ST_MakeLine(geom ORDER BY ST_Y(geom)) AS geom
FROM (
  SELECT row_number() OVER () AS cid, 
  (ST_Dump(geom)).geom FROM (
    SELECT unnest(st_clusterwithin(geom, 0.05)) AS geom 
    FROM inputs) sq) ssq 
GROUP BY cid
dbaston
fuente
Camino a seguir, pero el enfoque Y-monótono (o incluso cambiar entre X / Y-monótono) no funcionará bien si el conjunto de datos contiene una carretera curva. Es el caso? Ordenar algoritmo es la parte más difícil de esta pregunta en mi humilde opinión.
MarHoff
@MarHoff: sí, las carreteras curvas serán un problema, pero estoy tratando de transformar la mayoría de los datos automáticamente y el descanso deberá hacerse manualmente. O seguiré profundizando más en el tema para encontrar una solución, pero puede llevar más tiempo que hacer que alguien arregle los datos restantes. Tendré que evaluar los resultados para poder decidir. Gracias por señalar esto!
Mahakala
Statut sintonizado Acabo de pensar en un truco que necesito revisar ...
MarHoff
¿Hay alguna forma sólida de hacer esto que no implique probar todos los ordenamientos posibles de puntos y encontrar cuál da la longitud total más corta?
dbaston
Si estos conjuntos de puntos siempre siguen las carreteras, proyecte la posición del punto en el segmento de la carretera (ST_Line_Locate_Point), luego ordene los puntos por el resultado.
travis