点层和线层之间最近的邻居？ [关闭]

#1 楼

结果这个问题比我认为正确的要棘手。最短距离本身有很多实现方式，例如Shapely提供的距离（距GEOS）。但是，很少有解决方案提供交点本身，但是仅提供了距离。从给出确切答案。

以下是使用Shapely的完整解决方案，它基于以下方程式：很好，太糟糕了，它在用死语编写的死语平台上...

#2 楼

PostGIS答案（对于多线串，如果为线串，请删除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

#3 楼

这有点老了，但是我今天正在寻找解决这个问题的方法（点->行）。针对此相关问题，我遇到的最简单的解决方案是：

>>> 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)

#4 楼

如果我理解正确，那么您所要求的功能是内置于PostGIS中的。有关如何使用它的提示，可以在这里阅读：
http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis -part1 /

例如，也可以找到一个多边形上与另一个多边形的最近点。

如果要在两个几何上的两个最接近点之间建立直线，可以使用ST_Shortestline。 ST_Closestpoint是ST_Shortestline中的第一个点

两个几何之间的ST_Shortestline的长度与几何之间的ST_Distance相同。

#5 楼

请参阅以下有关如何将我的答案视为不可靠解决方案的评论...我将把此原始帖子保留在这里，以便其他人可以检查问题。该过程应该起作用。

查找点（由x，y或x，y，z定义）和一个polyine（由x，y或x的连接集定义）之间的最短路径，y，z）在欧几里得空间内：

1）从给定的用户定义点（我称其为pt0），找到折线（pt1）的最近顶点。 OGRinfo可以轮询折线的顶点，然后可以通过标准方法进行距离计算。例如，对距离calc进行迭代，例如：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）存储关联的最小距离值（d1）和（pt1）

3）观察远离pt1的两个片段（在ogrinfo线串中，这将是先前和随后的顶点）。记录这些顶点（n2和n3）。

4）为每个线段创建y = mx + b公式

5）将每个点的点（pt0）与垂直线相关联这两个公式中的一个

6）计算距离和交点（d2和d3； pt2，pt3）

7）比较三个距离（d1，d2，d3）最短的。您到关联节点（pt1，pt2或pt3）的pt0是最短的链接。 >

#6 楼

这是QGIS> 2.0的python脚本，由上面给出的提示和解决方案组成。对于合理数量的点和线，它可以正常工作。
但是我并没有尝试使用大量的对象。 py”。

在QGIS工具箱中查找脚本，工具，添加脚本，然后选择脚本。

 ##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
jean-christophe.baudin@onema.fr
"""
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()
 

!!!警告!!!
请注意，由于此行命令可能会产生一些“怪异” /错误的投影点：

设置返回的最近邻数。
实际上，每个点都应以距每个线对象尽可能短的距离进行投影；并且找到的最小距离将给出正确的直线和投影点，作为我们寻求的最接近的邻居。选择一个减少到1的counterSelec值将返回遇到的“第一个”对象（更确切地说是边界框），它可能不是正确的对象。为了确定最短距离，可能不得不选择3或5个不同的线尺寸对象，或者甚至更接近的对象。该值越高，花费的时间越长。随着数百个点和线的出现，与3个或5个最近的邻居开始变得非常慢，而成千上万个点和线可能会出现这种值。

#7 楼

根据您的兴趣和用例，研究“地图匹配算法”可能会很有用。例如，OSM Wiki上有一个RoadMatcher项目：http://wiki.openstreetmap.org/wiki/Roadmatcher。