¿El vecino más cercano entre la capa de puntos y la capa de líneas? [cerrado]

34

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

Utilizando programación (lo más 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 línea o capa de polígono.

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

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

  • Usando las bibliotecas GDAL / OGR, ¿dónde debería comenzar? ¿Sería posible dar un plan de solución?
  • ¿Puedo usar NetworkX para hacer 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

    
pregunta dassouki 27.07.2010 - 19:54

7 respuestas

31

Esta pregunta resultó ser un poco más complicada de lo que pensé para acertar. Hay muchas implementaciones de la distancia más corta en sí misma, como la distancia (desde GEOS) . Sin embargo, pocas de las soluciones proporcionan el punto de intersección en sí, pero solo la distancia.

Mi primer intento amortiguó el punto por la distancia entre el punto y el polígono, y busqué intersecciones, pero los errores de redondeo evitan que esto dé una respuesta exacta.

Aquí hay una solución completa que usa 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, parece que esta extensión ArcView maneja este problema muy bien, muy mal, está en un muerto Plataforma escrita en un lenguaje muerto ...

    
respondido por el scw 29.07.2010 - 04:45
8

Una respuesta PostGIS (para multilínea, si es una cadena de líneas, elimina 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
    
respondido por el eprand 09.09.2010 - 02:31
4

Esto es un poco viejo, pero hoy estaba buscando soluciones para este problema (punto - > línea). La solución más sencilla 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)
    
respondido por el Richard Law 03.01.2014 - 03:27
4

Si entiendo bien, la funcionalidad que solicitas está integrada en PostGIS.

Para obtener un punto proyectado en una línea, puede usar ST_Closestpoint (en PostGIS 1.5)

Algunos consejos sobre cómo usarlo se pueden leer aquí: enlace

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

Si desea que la línea entre los dos puntos más cercanos en ambas geometrías puede usar ST_Shortestline. ST_Closestpoint es el primer punto en ST_Shortestline

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

    
respondido por el Nicklas Avén 06.10.2010 - 15:00
3

Vea el comentario a continuación sobre cómo mi respuesta no debe considerarse una solución confiable ... Dejaré esta publicación original aquí para que otros puedan examinar el problema.

Si entiendo la pregunta, este procedimiento general debería funcionar.

Para encontrar la ruta más corta entre un punto (como se define por x, yo x, y, z) y una polina (como se define por un conjunto de conexión de x, yo x, y, z) dentro del espacio euclidiano:

1) Desde un punto definido por el usuario dado (lo llamaré pt0), encuentre el vértice más cercano de la polilínea (pt1). OGRinfo puede sondear los vértices de una polilínea, y luego se pueden hacer cálculos de distancia a través de métodos estándar. Por ejemplo, iterar sobre una distancia calculada como: distance_in_radians = 2 * math.asin (math.sqrt (math.pow ((math.sin ((pt0_radians-ptx_radians) / 2)), 2) + math.cos (pt0_radians) * math.cos (ptx_radians) * math.pow ((math.sin ((pt0_radians-ptx_radians) / 2)), 2)))

2) Almacene el valor de la distancia mínima asociada (d1) y (pt1)

3) observe los dos segmentos que se alejan de pt1 (en la cadena de líneas de ogrinfo, estos serán los vértices anterior y posterior). Registre estos vértices (n2 y n3).

4) cree y = mx + b fórmula para cada segmento

5) Relaciona tu punto (pt0) con la perpendicular para cada una de esas dos fórmulas

6) Calcular la distancia y las intersecciones (d2 y d3; pt2, pt3)

7) Compara las tres distancias (d1, d2, d3) para las más cortas. Su pt0 al nodo asociado (pt1, pt2 o pt3) es el enlace más corto.

Esa es una respuesta de la corriente de conciencia: con suerte, mi imagen mental del problema y la solución se ajusta.

    
respondido por el glennon 27.07.2010 - 23:25
3

Aquí hay una secuencia de comandos de Python para QGIS > 2.0 hecha a partir de las sugerencias y soluciones indicadas anteriormente. Funciona bien para una cantidad razonable de puntos y líneas. Pero no lo probé con gran cantidad de objetos.

Por supuesto, tenía que copiarse en modo inactivo o cualquier otra "solución pythonic" y guardarlo como "closest.point.py".

En la caja de herramientas de QGIS, vaya para el script, herramientas, agregue un script, y elíjalo.

##Vector=group
##CLosest_Point_V2=name
##Couche_de_Points=vector
##Couche_de_Lignes=vector

"""
This script intent to provide a count as for the SQL Funciton CLosestPoint
Ce script vise a recréer dans QGIS la Focntion SQL : CLosest Point
It rely on the solutions provided in "Nearest neighbor between a point layer and a line layer"
  http://gis.stackexchange.com/questions/396/nearest-pojected-point-from-a-point-                               layer-on-a-line-or-polygon-outer-ring-layer
V2 du  8 aout 2016
[email protected]
"""
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import os
import sys
import unicodedata 
from osgeo import ogr
from math import sqrt
from sys import maxint

from processing import *

def magnitude(p1, p2):
    if p1==p2: return 1
    else:
        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.0001 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 QgsPoint(ix, iy)

layerP = processing.getObject(Couche_de_Points)
providerP = layerP.dataProvider()
fieldsP = providerP.fields()
inFeatP = QgsFeature()

layerL = processing.getObject(Couche_de_Lignes)
providerL = layerL.dataProvider()
fieldsL = providerL.fields()
inFeatL = QgsFeature()

counterP = counterL= nElement=0

for featP in layerP.selectedFeatures():
    counterP+=1
if counterP==0:
    QMessageBox.information(None,"information:","Choose at least one point from point layer_"+ str(layerP.name())) 

indexLine=QgsSpatialIndex()
for featL in layerL.selectedFeatures():
    indexLine.insertFeature(featL)
    counterL+=1
if counterL==0:
    QMessageBox.information(None,"information:","Choose at least one line from point layer_"+ str(layerL.name()))
    #QMessageBox.information(None,"DEBUGindex:",str(indexBerge))     
ClosestP=QgsVectorLayer("Point", "Projected_ Points_From_"+ str(layerP.name()), "memory")
QgsMapLayerRegistry.instance().addMapLayer(ClosestP)
prClosestP = ClosestP.dataProvider()

for f in fieldsP:
    znameField= f.name()
    Type= str(f.typeName())
    if Type == 'Integer': prClosestP.addAttributes([ QgsField( znameField, QVariant.Int)])
    if Type == 'Real': prClosestP.addAttributes([ QgsField( znameField, QVariant.Double)])
    if Type == 'String': prClosestP.addAttributes([ QgsField( znameField, QVariant.String)])
    else : prClosestP.addAttributes([ QgsField( znameField, QVariant.String)])
prClosestP.addAttributes([QgsField("DistanceP", QVariant.Double),
                                        QgsField("XDep", QVariant.Double),
                                        QgsField("YDep", QVariant.Double),
                                        QgsField("XProj", QVariant.Double),
                                        QgsField("YProj", QVariant.Double),
                                        QgsField("Xmed", QVariant.Double),
                                        QgsField("Ymed", QVariant.Double)])
featsP = processing.features(layerP)
nFeat = len(featsP)
"""
for inFeatP in featsP:
    progress.setPercentage(int(100 * nElement / nFeatL))
    nElement += 1
    # pour avoir l'attribut d'un objet/feat .... 
    attributs = inFeatP.attributes()
"""

for inFeatP in layerP.selectedFeatures():
    progress.setPercentage(int(100 * nElement / counterL))
    nElement += 1
    attributs=inFeatP.attributes()
    geomP=inFeatP.geometry()
    nearest_point = None
    minVal=0.0
    counterSelec=1
    first= True
    nearestsfids=indexLine.nearestNeighbor(geomP.asPoint(),counterSelec)
    #http://blog.vitu.ch/10212013-1331/advanced-feature-requests-qgis
    #layer.getFeatures( QgsFeatureRequest().setFilterFid( fid ) )
    request = QgsFeatureRequest().setFilterFids( nearestsfids )
    #list = [ feat for feat in CoucheL.getFeatures( request ) ]
    # QMessageBox.information(None,"DEBUGnearestIndex:",str(list))
    NBNodes=0
    Dist=DistT=minValT=Distance=0.0
    for featL in  layerL.getFeatures(request):
        geomL=featL.geometry()
        firstM=True
        geomL2=geomL.asPolyline()
        NBNodes=len(geomL2)
        for i in range(1,NBNodes):
            lineStart,lineEnd=geomL2[i-1],geomL2[i]
            ProjPoint=intersect_point_to_line(geomP.asPoint(),QgsPoint(lineStart),QgsPoint(lineEnd))
            Distance=magnitude(geomP.asPoint(),ProjPoint)
            toto=''
            toto=toto+ 'lineStart :'+ str(lineStart)+ '  lineEnd : '+ str(lineEnd)+ '\n'+ '\n'
            toto=toto+ 'ProjPoint '+ str(ProjPoint)+ '\n'+ '\n'
            toto=toto+ 'Distance '+ str(Distance)
            #QMessageBox.information(None,"DEBUG", toto)
            if firstM:
                minValT,nearest_pointT,firstM = Distance,ProjPoint,False
            else:
                if Distance < minValT:
                    minValT=Distance
                    nearest_pointT=ProjPoint
            #at the end of the loop save the nearest point for a line object
            #min_dist=magnitude(ObjetPoint,PProjMin)
            #QMessageBox.information(None,"DEBUG", " Dist min: "+ str(minValT))
        if first:
            minVal,nearest_point,first = minValT,nearest_pointT,False
        else:
            if minValT < minVal:
                minVal=minValT
                nearest_point=nearest_pointT
                #at loop end give the nearest Projected points on Line nearest Line
    PProjMin=nearest_point
    Geom= QgsGeometry().fromPoint(PProjMin)
    min_dist=minVal
    PX=geomP.asPoint().x()
    PY=geomP.asPoint().y()
    Xmed=(PX+PProjMin.x())/2
    Ymed=(PY+PProjMin.y())/2
    newfeat = QgsFeature()
    newfeat.setGeometry(Geom)
    Values=[]
    #Values.extend(attributs)
    fields=layerP.pendingFields()
    Values=[attributs[i] for i in range(len(fields))]
    Values.append(min_dist)
    Values.append(PX)
    Values.append(PY)
    Values.append(PProjMin.x())
    Values.append(PProjMin.y())
    Values.append(Xmed)
    Values.append(Ymed)
    newfeat.setAttributes(Values)
    ClosestP.startEditing()  
    prClosestP.addFeatures([ newfeat ])
    #prClosestP.updateExtents()
ClosestP.commitChanges()
iface.mapCanvas().refresh()

!!! ADVERTENCIA !!! Tenga en cuenta que algunos puntos proyectados "extraños" / incorrectos pueden producirse debido a este comando de línea:

nearestsfids=indexLine.nearestNeighbor(geomP.asPoint(),counterSelec)

El valor counterSelec en la misma establece la cantidad de devoluciones más cercanas. De hecho, cada punto debe proyectarse a la distancia más corta posible para cada objeto de línea; y la distancia mínima encontrada daría la línea correcta y el punto proyectado como los vecinos más cercanos que buscamos. Para reducir el tiempo de bucle, se utiliza el comando Vecino más cercano. Elegir un valor de counterSelec reducido a 1 devolverá el "primer" objeto encontrado (es más exacto el cuadro delimitador) y puede que no sea el correcto. Los objetos de diferente tamaño de línea pueden obligar a elegir pueden ser 3 o 5, o incluso objetos más cercanos para determinar la distancia más corta. Cuanto mayor sea el valor, más tiempo lleva. Con cientos de puntos y líneas, comienza a ser muy lento con 3 o 5 vecinos más cercanos, con miles de ellos puede funcionar con dichos valores.

    
respondido por el Jean-Christophe Baudin 08.08.2016 - 14:13
3

Dependiendo de sus intereses y caso de uso, podría ser útil buscar "algoritmos de correlación de mapas". Por ejemplo, hay un proyecto RoadMatcher en la wiki de OSM: enlace .

    
respondido por el underdark 27.07.2010 - 20:40

Lea otras preguntas en las etiquetas