ArcPy or ArcObject code to find the closest point on a line.

George_GShi · ‎11-24-2015

Hi All,

Can you give me some hints (code snippet in Arcpy or ArcObject) on:

for a given point to find the closest point location on a line feature.

Thanks in advance!

George

DanPatterson_Retired · ‎03-27-2016

Richard, it is in your array forms... don't know how x,y got passed into an array object... here is some code using your numbers. ps, no need to cast your array to a float dtype, since numpy figures it out from the inputs and only needs guidance when there might be some confusion

>>> orig = np.array([ 6305489.12520918,2169986.76545799])
>>> orig
array([ 6305489.13,  2169986.77])
>>> dest = np.array([[6305087.957613736, 2169963.2243613005],[6305515.50958015, 2169962.629874304]])
>>> diff = dest-orig
>>> diff
array([[-401.17,  -23.54],
       [  26.38,  -24.14]])

where the array represents the dx, dy pairs for each destination from the points. einsum is faster... I posted something on my blog using it when I was comparing numpy to scipy, pdist (or cdist)... you should be able to calculate the distance of millions of points quicly.

Not the arrays are lists of lists, in the case of destination, I noticed that your tuples were open ended with the extra comma after the braces... you can list(a tuple) if your acquisition is not standard. NumPy Repository houses other documents etc that I have been working on for some time, for those interested in numpy and its dressed up brethern, scipy and pandas which use numpy arrays natively.

RichardFairhurst · ‎03-27-2016

Dan:

The X,Y got passed into an array object by the FeatureClassToNumpyArray, which is the method you recommended for people to use in your original post to convert a polyline to an array of points. That method is used in line 17 with the explode_to_points option set to True. The first output list used the ["SHAPE@X", "SHAPE@Y"] fields and the second list which had the extra comma used the ["SHAPE@XY"] field. The arrays are the unprocessed native output of that method with the explode_to_point=True option and those field settings.

I found the answer on how to make the output of FeatureClassToNumpyArray method with the explode_to_point=True option work with the rest of the code. I changed line 17 to :

dests = np.array(list(list(x) for x in arcpy.da.FeatureClassToNumPyArray(addressLines, ["SHAPE@X","SHAPE@Y"], "", sr, True)),dtype="float64")

Now the code works and the Origin and Dest arrays print as:

Origin Point = [ 6305489.12520918 2169986.76545799]

Destination Points = [[ 6305087.95761374 2169963.2243613 ]

[ 6305515.50958015 2169962.6298743 ]

...

[ 6305742.79783948 2169860.23900296]]

The numpy array rounds the orignal numbers from 9 decimal places to 8 decimal places (i.e., 6305087.957613736 became 6305087.95761374), but since that is way below my XY resolution and tolerance that should not be a problem.

So the full working code is:

class SplitLineToolClass(object):  
    """Implementation for SplitLine_addin.SplitLinetool (Tool)"""  
    def __init__(self):  
        self.enabled = False  
        self.cursor = 3 # Crosshairs
  
    def onMouseDownMap(self, x, y, button, shift):  
        # Pass if not a Left Mouse Button Click  
        if button != 1:  
            pass  
            return  
        addressLines = pythonaddins.GetSelectedTOCLayerOrDataFrame()  
        ## Code that verifies a polyline layer is selected in the TOC  
        desc=arcpy.Describe(addressLines)  
        sr = arcpy.SpatialReference(desc.spatialReference.factoryCode)  
        origin = np.array([x,y],dtype="float64")  
        dests = np.array(list(list(x) for x in arcpy.da.FeatureClassToNumPyArray(addressLines, ["SHAPE@X","SHAPE@Y"], "", sr, True)),dtype="float64")
        deltas = dests - origin   
        distances = np.hypot(deltas[:,0], deltas[:,1])   
        min_dist = np.min(distances)  
        wh = np.where(distances == min_dist)   
        closest = dests[wh[0]]  
        print("Closest Vertice = {0}".format(closest))

DanPatterson_Retired · ‎03-27-2016

Yes I don't suspect 10^-8 or so meters is going to make much of a difference float128 is around the corner and then there is always the Decimal module if you really want to kill speed but are working at the atomic level

DanPatterson_Retired · ‎03-27-2016

PS.... if you need closest point on a line code... I have it somewhere, It replicates the Near tool functionality, if that is what you are looking for. I can't remember if it is a pure python implementation of a numpy/arcpy blending. In any event, if this is what you are looing for and I will try to dredge it up when I get on an Arc* machine.

RichardFairhurst · ‎03-27-2016

Dan:

Thanks for the offer, but I actually am using the code I posted to compare the nearest edge point returned by either the distanceTo or the queryPointAndDistance method of the polyline. I need at minimum two of the items the queryPointAndDistance method returns and may need the other two latter. I need the side offset (sideDist) to see if the line edge is inside my split tolerance, and I need the distance along the line where the point fell as a percent (perc) to actually split the line.

pt = arcpy.PointGeometry(arcpy.Point(x,y), sr)

ptOnLine, perc, sideDist, rightSide = pLine.queryPointAndDistance(pt, True)

The output I am after is to see if the vertice is within a maximum distance specified by the user of the real nearest point on the line edge from where the user clicked to determine if I will split at the closest edge or the closest vertice. After getting the queryPointAndDistance percent along the line of either the point on the line vertice returned by your code or point on the line edge returned by either the distanceTo or queryPointAndDistance method I then use that percent to split the line using:

fromLine = pLine.segmentAlongLine(0, perc, True)

toLine = pLine.segmentAlongLine(perc, 100, True)

I also use the percent to split address ranges.

So in my example, the full data showed that if a user chose a maximum vertice/edge difference of 10 feet (my units), then the edge should split and not the vertice, since the vertice is more than 11 feet from the closest point on the line edge based on where the user clicked.

Origin Point = [ 6305489.12520918 2169986.76545799]

Destination Points = [[ 6305087.95761374 2169963.2243613 ]

[ 6305515.50958015 2169962.6298743 ]

...

[ 6305742.79783948 2169860.23900296]]

Closest Vertice = [[ 6305515.50958015 2169962.6298743 ]]

Vertice Distance 35.7583756816 is within Max Distance 50.0

Difference between nearest vertice and nearest edge is 11.6595012766

JeremyMullins4 · ‎04-14-2016

Dan Patterson,

I would be interested in seeing what you've come up with to replace the Near tool. I only have a standard license. I have been messing with arcpy geometry and distanceTo but I haven't gotten very far.

Basically, I have an input of address points and street centerlines, and basically need to find the nearest two polylines from a particular address point and pull their individual IDs or (even better) their street names into the address points feature class as cross streets.

George_GShi · ‎11-27-2015

Thank you all for all your response. Special thanks to Dan for your detailed solutions.

Again, I need to create a function with its:

Input: a line and a point

Output: the point location and distance

Your comments remind me of two cases I need to handle:

1) to find the closest vertex of the line to the given point location.

2) to find a point location “on” (or intersecting) the line close to the given point

DanPatterson_Retired · ‎11-27-2015

The first case, I have given the necessary code if you explode the featureclass to points.

The second case, depends upon the license level that you have. I guess Darren's answer twigged your memory, so I presume that you can get all you need from the Near tool since it allows you to create from-to coordinates in the table, which you can use in the XY to Line tool to produce the necessary vectors. If you want to see the points on the line, just add the near result table as an event theme with the destination points as your XY coordinates. In that way, you can have everything.

I have variants which use pure arcpy and arctoolboxes, and two variants that use a combination of arcpy and numpy. One of these is reported in an earlier blog post of mine and I will post the other variant in one of my upcoming blogs.

In the interim, I presume that you have the necessary license to complete the tasks.

NeilAyres · ‎12-14-2015

Although its fun to see Dan come up with a numpy solution...

Would it fit your needs to use the queryPointAndDistance geometry method?