¿Proximidad entre la capa de puntos y la capa de líneas?

Question

He hecho esta pregunta varias veces en stackoverflow e irc entre #qgis y #postgis y también he intentado codificarlo o implementarlo yo mismo en postgis sin respuesta real.

Mediante programación (preferiblemente python), me gustaría dibujar una línea desde una capa de puntos, hasta su proyección en la línea más cercana de una capa de líneas o polígonos.

Por ahora, la mayoría de mis datos están en los formatos shape y postgis de ESRI; sin embargo, prefiero mantenerme alejado de una solución postgis, ya que soy un usuario predominantemente de shp + qgis.

Una solución ideal sería implementar GDAL/OGR con bibliotecas de python o similares

• Utilizando las bibliotecas GDAL/OGR, ¿por dónde debería empezar? ¿Sería posible dar un plan de solución?
• ¿Puedo utilizar NetworkX para realizar el análisis del vecino más cercano?
• ¿Es esto realmente posible?

Si es más fácil, los puntos podrían conectarse al punto final del segmento en lugar de un punto proyectado

Answer 1

Esta pregunta ha resultado ser un poco más complicada de lo que pensaba. Hay muchas implementaciones de la distancia más corta en sí misma, como la proporcionada por Shapely distancia (de GEOS). Sin embargo, pocas de las soluciones proporcionan el punto de intersección en sí, sino sólo la distancia.

Mi primer intento amortiguaba el punto por la distancia entre el punto y el polígono, y buscaba las intersecciones, pero los errores de redondeo impiden que esto dé una respuesta exacta.

Aquí hay una solución completa usando Shapely, basada en estas ecuaciones :

#!/usr/bin/env python
from shapely.geometry import Point, Polygon
from math import sqrt
from sys import maxint

# define our polygon of interest, and the point we'd like to test
# for the nearest location
polygon = Polygon(((0, 0), (0, 1), (1, 1), (1, 0), (0, 0)))
point = Point(0.5, 1.5)

# pairs iterator:
# http://stackoverflow.com/questions/1257413/1257446#1257446
def pairs(lst):
    i = iter(lst)
    first = prev = i.next()
    for item in i:
        yield prev, item
        prev = item
    yield item, first

# these methods rewritten from the C version of Paul Bourke's
# geometry computations:
# http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/
def magnitude(p1, p2):
    vect_x = p2.x - p1.x
    vect_y = p2.y - p1.y
    return sqrt(vect_x**2 + vect_y**2)

def intersect_point_to_line(point, line_start, line_end):
    line_magnitude =  magnitude(line_end, line_start)
    u = ((point.x - line_start.x) * (line_end.x - line_start.x) +
         (point.y - line_start.y) * (line_end.y - line_start.y)) \
         / (line_magnitude ** 2)

    # closest point does not fall within the line segment, 
    # take the shorter distance to an endpoint
    if u < 0.00001 or u > 1:
        ix = magnitude(point, line_start)
        iy = magnitude(point, line_end)
        if ix > iy:
            return line_end
        else:
            return line_start
    else:
        ix = line_start.x + u * (line_end.x - line_start.x)
        iy = line_start.y + u * (line_end.y - line_start.y)
        return Point([ix, iy])

nearest_point = None
min_dist = maxint

for seg_start, seg_end in pairs(list(polygon.exterior.coords)[:-1]):
    line_start = Point(seg_start)
    line_end = Point(seg_end)

    intersection_point = intersect_point_to_line(point, line_start, line_end)
    cur_dist =  magnitude(point, intersection_point)

    if cur_dist < min_dist:
        min_dist = cur_dist
        nearest_point = intersection_point

print "Closest point found at: %s, with a distance of %.2f units." % \
   (nearest_point, min_dist)

Para la posteridad, se ve así Extensión de ArcView maneja este problema bastante bien, lástima que esté en una plataforma muerta escrita en un lenguaje muerto...

Answer 2

Una respuesta de PostGIS (para la multilínea, si la cadena de líneas, eliminar la función st_geometryn)

select t2.gid as point_gid, t1.gid as line_gid, 
st_makeline(t2.geom,st_line_interpolate_point(st_geometryn(t1.geom,1),st_line_locate_point(st_geometryn(t1.geom,1),t2.geom))) as geom
from your_line_layer t1, your_point_layer t2, 
(
select gid as point_gid, 
(select gid 
from your_line_layer
order by st_distance(your_line_layer.geom, your_point_layer.geom)
limit 1 ) as line_gid
from your_point_layer
) as t3
where t1.gid = t3.line_gid
and t2.gid = t3.point_gid

Answer 3

Esto es un poco viejo, pero hoy estaba buscando soluciones a este problema (punto --> línea). La solución más simple que he encontrado para este problema relacionado es:

>>> from shapely.geometry import Point, LineString
>>> line = LineString([(0, 0), (1, 1), (2, 2)])
>>> point = Point(0.3, 0.7)
>>> point
POINT (0.3000000000000000 0.7000000000000000)
>>> line.interpolate(line.project(point))
POINT (0.5000000000000000 0.5000000000000000)

Answer 4

Si te entiendo bien la funcionalidad que pides está incorporada en PostGIS.

Para obtener un punto proyectado sobre una línea se puede utilizar ST_Closestpoint (en PostGIS 1.5)

Puedes leer algunos consejos sobre su uso aquí: http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/

También se puede utilizar para encontrar el punto más cercano de un polígono a otro polígono, por ejemplo.

Si quieres la línea entre los dos puntos más cercanos de ambas geometrías puedes utilizar ST_Shortestline. ST_Closestpoint es el primer punto de ST_Shortestline

La longitud de ST_Shortestline entre dos geometrías es la misma que ST_Distance entre las geometrías.

Answer 5

Dependiendo de sus intereses y de su caso de uso, podría ser útil buscar "algoritmos de coincidencia de mapas". Por ejemplo, hay un proyecto RoadMatcher en el wiki de OSM: http://wiki.openstreetmap.org/wiki/Roadmatcher .