Robbi
I'm not able to download your sample data on this laptop, but I understand the gist. This problem lends itself to an elegant solution in a few lines of script. I assume the points are snapped exactly on line endpoints - if not, you may need to run the Snap tool on the points first because the XY locations have to be identical for it to work.
This is a standalone script (not in a Toolbox) so just change the hard-coded paths (for "linefc" and "pointfc") and field names ("FromPointID", "ToPointID", "ID") and hit run in a Python IDE or command prompt.
__author__ = "John K. Tran"
__contact__ = "[email protected]"
__credits__ = "http://gis.stackexchange.com/questions/158490/extracting-the-closest-point-to-a-line-within-a-polygon-using-arcgis-for-desktop"
import arcpy
linefc = r"C:\Users\Me\Desktop\Sample.gdb\SampleLines" # Change this to path to your polyline feature class.
pointfc = r"C:\Users\Me\Desktop\Sample.gdb\SamplePoints" # Change this to path to your point feature class.
assert arcpy.Exists(linefc) and arcpy.Exists(pointfc), "Please ensure the paths to the Line FC and Point FC are correct."
with arcpy.da.UpdateCursor(linefc, ["SHAPE@", "FromPointID", "ToPointID"]) as linecursor: # Change the "FromPointID" and "ToPointID" to the respective field names in your polyline FC.
for linerow in linecursor:
linegeom = linerow[0]
firstpoint = linegeom.firstPoint
lastpoint = linegeom.lastPoint
firstpointXY = (firstpoint.X, firstpoint.Y)
lastpointXY = (lastpoint.X, lastpoint.Y)
with arcpy.da.SearchCursor(pointfc, ["SHAPE@XY", "ID"]) as pointcursor: # Change the "ID" to the respective field name in your point FC, like "ObjectID" etc.
for pointrow in pointcursor:
pointXY = pointrow[0]
pointID = pointrow[1]
if pointXY == firstpointXY:
linerow[1] = pointID
if pointXY == lastpointXY:
linerow[2] = pointID
del pointcursor
linecursor.updateRow(linerow)
del linecursor
print "Done!"
Here is a screenshot of the results I got, with the IDs for the points labeled (black) and the From and To attributes on the lines labeled (blue).
Looks like you're pretty fluent in NumPy, so you shouldn't have a problem with ArcPy. Let me know if you have any issues!
*A note on the XY locations. Sometimes precision errors occur in the GIS such that the very last hundredth, millionth, etc. decimals are not equivalent. I didn't do it, but you could build a tolerance in the == equalities so that it compares whether the XY are roughly equal rather than exactly equal. Or I'm sure someone out there knows some Environment Variable in ArcGIS to control resolution/tolerance in geoprocessing.
---Response to Robbi's Comments---
Hi Robbi,
So I think there are some issues keeping the script from working that relates to precision. Shapefiles are notoriously outdated data formats. Below, I've done a Select by Location for the Points that "Intersect" the Lines, and you can see not all the points are coincident with their lines (not including the point in the middle that I assume you doctored for testing purpose).
If I zoom into one of the unselected points, you can see at 1:1 scale, the point is just off the line. This has to do with the limited floating-point precision that shapefiles can store XY coordinate data.
(unsnapped)
One way to solve this is to use the "Snap" tool in the Editing Toolbox.
(snapped)
This will resolve most of our issues with the points. However, you can still see gaps of *0*s for some of the fromID and toID values. This is due to lack of precision in the shapefile storage format itself (I really suggest migrating to file geodatabase storage if interoperability is not an issue (has better precision, more schema options, no storage limit, etc.), shapefiles are somewhat legacy these days).
To populate the remaining values, we have to edit the == equality and do some math to calculate when the XY coords are roughly equivalent rather than absolutely equivalent.
I replaced the previous code which looks for straight identicalness:
if pointXY == firstpointXY:
linerow[1] = pointID
if pointXY == lastpointXY:
linerow[2] = pointID
With one that looks for values within a range. Since you're coordinate system was WGS84 (decimal degrees), I chose a tolerance of 0.0001 DD as being equivalent (which is roughly 1.1132 m at equator). In a coordinate system based on meters e.g. NAD83, you would use something like 1 m for your tolerance.
if (firstpointXY[0]-0.00001 < pointXY[0] < firstpointXY[0]+0.00001) and (firstpointXY[1]-0.00001 < pointXY[1] < firstpointXY[1]+0.00001):
linerow[1] = pointID
if (lastpointXY[0]-0.00001 < pointXY[0] < lastpointXY[0]+0.00001) and (lastpointXY[1]-0.00001 < pointXY[1] < lastpointXY[1]+0.00001):
linerow[2] = pointID
After making this adjustment... tadaaa, the Attribute Table is all populated.
Getting to the tail-end of your question, to export the table out as a .csv, just go to Table Options in the top-left >>> Export >>> Choose "Save as Type: Text File" and name it something like "My_Results.csv".
*Last note about efficiency: So since the script involves two nested cursor iterations (for every line, loop through every point - which for 1,000 lines and 100,000 points is 100 million passes) one solution would be to use a Spatial Join to have an attribute that describes which LineID each point intersects.
Then pass this LineID as the whereclause in your point cursor.
with arcpy.da.UpdateCursor(linefc, ["SHAPE@", "LineID", "FromPointID", "ToPointID"]) as linecursor:
for linerow in linecursor:
arcpy.SetProgressorLabel("Processing Line {0}".format(count))
linegeom = linerow[0]
lineID = linerow[1]
...
whereclause = "\"{0}\" = '{1}'".format("LineID", str(lineID))
with arcpy.da.SearchCursor(pointfc, ["SHAPE@XY", "PointID"], whereclause) as pointcursor:
We've basically said "only consider the points whose LineID matches the LineID of the current polyline." This limits our comparisons each iteration originally from say 100,000 points to essentially 2 each time. Much faster...
