2

I have a polyline shapefile of rivers/stream segments, and a point shapefile of intersection points (nodes) between the segments (edges). Each stream segment starts and finishes at an intersection point. I wish to generate an edge list file (.csv or similar) which contains the start and end node IDs for each stream segment ID, in a format similar to "Edge ID, From ID, To ID" (see example image below).

enter image description here

Preferably I would like a solution to this within ArcGIS 10. So far I've experimented with a iterative approach (iterating through each edge segment and selecting nodes within a +-5m radius), but feel that this might not be optimal for large stream networks (100,000+ intersections).

I have uploaded a basic sample dataset (one polyline and point shapefile matching the image above) here.

  • Hey Robbi, I just read your edit. So it looks like you've already got the iterative Python approach going, but want something more efficient. I'll give it some thought - some preliminary ideas are using Feature Vertices to Points, then the Near tool and/or a Spatial or Table Join to accomplish what you want. Are you on an ArcGIS Advanced license? – John Aug 17 '15 at 13:17
  • Looks like with the above approach, I can figure out which points relate to which lines e.g. Line 0 spatially joins to Points 0 and 1. However, I can't figure out the order of those points e.g. Point 0 is the start node, Point 1 is the end node. – John Aug 17 '15 at 14:17
1

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__ = "jtran20@masonlive.gmu.edu"
__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).

Solution Screenshot

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

IMG1

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

One way to solve this is to use the "Snap" tool in the Editing Toolbox.

Snap Snapped(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).

Table1

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.

Table2

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".

Export1 Export2

*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.

LineID

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

  • Hi John, this looks to be exactly what I was looking for (I am familiar with Python but haven't had much experience with ArcPy). The only problem is that when run with my example dataset, I get no output in the selected from/to fields. I was wondering if you could either test the script on the example I posted above (I've fixed the link), or send me a copy of the datasets you used to test if it's something wrong on my end? – Robbi Bishop-Taylor Aug 18 '15 at 1:48
  • Can you email me the example to jtran20@masonlive.gmu.edu? @RobbiBishop-Taylor Unfortunately I can't access cloud sharing (e.g. DropBox, Google Drive, Microsoft SkyDrive, etc.) from this laptop - work policy, but email is fine. – John Aug 18 '15 at 12:50
  • @RobbiBishop-Taylor The answer is updated to explain why your sample data didn't work. Basically the coordinate precision in shapefiles is horrible and causes the equivalence == codes to return False because the hundredth/thousandth/millionth decimal place is not equal. A mix of using the Snap tool and modifying the equivalences to accept a range of values solves it. More thorough explanation is provided above. Thanks! – John Aug 20 '15 at 2:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.