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

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

- 3 people found this helpful
Near—Help | ArcGIS for Desktop Is a good place to start if you have the advanced level (which is why Darren asked). There are several samples of using python code.

if you don't have advanced, then it will take a bit more work.

- 2 people found this helpful
For a python approach, simply

- get the coordinates of your origin ( a variety of ways)
- get the coordinates of your polyline
- any license FeatureClassToNumPyArray—Help | ArcGIS for Desktop using the explode_to_points option
- advanced license Feature Vertices To Points—Help | ArcGIS for Desktop

- convert both to numpy arrays ( a variety of ways, but FeatureClassToNumPyArray is best
- use the code example below

Assume a start location at the origin, and projected coordinates (values have been translated to the origin to simplify the visuals and the math)

>>> import numpy as np >>> origin = np.array([0,0],dtype="float64") >>> dests = np.array([[1,2],[-3,3],[-4,-4],[1,1],[5,-5]],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]] >>> closest array([[ 1., 1.]]) >>>

line (02) that is the origin

line (03 ) the destinations of the points constituting the polyline

line (04) get the coordinate differences)

line (05) calculate the distances, note the np.hypot uses dy,dx as inputs

line (06) find the minimum distance

line(07) find where the min distance is

line(08) use the position of the min distance to slice the original destinations. ( note wh returns a tuple and you want the first element)

the closest point...obviously is 1,1

see other purely arcpy options.

I'm not sure, but I believe this will only get you the closest vertex, not the closest location on a line (which doesn't necessarily mean a vertex). You need to do some trig in between each vertex. But, I could be missing some of the numpy magic going on here.

I left that part out ... but I will post about it soon. In the interim, George found your solution with the wrong license level worked for him.

Well, I don't believe I answered the question, but maybe thinking about license-levels did the trick...? Anyhow, I'd be interested to see what this looks like through numpy.

will do, the example I gave was for non-sequential points... caught it too late. But if you want to experiment for sequential points consider a polyline consisting of points a,b,c,d and a source point X

- you just find the closest point on the line from your source point... lets say 'c' wins and 'b' and 'd' are closer than 'a'
- you now have two potential segments to examine b-c and c-d (this could be multiple, but lets keep it simple.
- project X onto the candidate segments (the geometry stuff) to get the point on the segment and the distance to X from the segments (could be multiple)
- the segment from above that has shortest distance is the correct line segment and you have the intersection point and closest point onto the polyline

I have simplified some checks, but you can incorporate extent checks to see which segments of the polyline contain space that X can project on to. Remember, this is projection onto a segment and not a imaginary line extending from it. So soon, it is an interesting comparison of numpy solution and Near tool

Dan:

I was trying to use your code, since I want to find the minimum distance to the actual vertice points in a polyline with the onMouseDownMap function of an addin tool. Here is the code I tried:

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 = arcpy.da.FeatureClassToNumPyArray(addressLines, ["SHAPE@X","SHAPE@Y"], "", sr, explode_to_points=True) 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))

However, my adaptation of your code is throwing an error on the line that subtracts the two arrays (deltas = dests - origin on line 18 of the code above). Here is the error:

Traceback (most recent call last):

File ..., line ..., in onMouseDownMap

deltas = dests - origin

TypeError: unsupported operand type(s) for -: 'numpy.ndarray' and 'numpy.ndarray'

When I added print lines for the origin and dest arrays, it showed that the Origin array created by line 16 is a list with a space a number a double space and a number, while the Dest array was a list of tuples containing two numbers separated by a comma and a space. What is the best way to get the two arrays compatible for line 18?

Here is how the arrays printed:

Origin Point = [ 6305489.12520918 2169986.76545799]

Destination Points = [(6305087.957613736, 2169963.2243613005)

(6305515.50958015, 2169962.629874304)

...

(6305742.797839478, 2169860.239002958)]

When I tried the "SHAPE@XY" field instead the result seemed worse:

Destination Point = [([6305087.957613736, 2169963.2243613005],)

([6305515.50958015, 2169962.629874304],)

...

([6305742.797839478, 2169860.239002958],)]

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.

- 1 person found this helpful
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))

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.

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

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.

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

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.

What license level is available?