[GIS] Finding distance of point along polyline using ArcGIS Desktop

arcgis-desktopdistancelinear-referencingpoint

We need to find the distance of a series of point along a given polyline.

We have a polyline dataset repsenting rivers in a catchment, broken up into river reaches (FID 1 through 35 for example).

We also have a point dataset of abstraction licenses within this catchment, some of which are situated along rivers – we can tell from the attribute table which these are as they are defined as groundwater or surface water licenses.

The two are not related – i.e. one has not been derived from the other.
Therefore not all the surface water licenses lie exactly over the river polylines, despite in 'reality' being on these rivers.

We need to be able to see which points are "upstream/downstream" of each other, and so wanted to have a list of which river reach the point is on, followed by the distance downstream from the start of this reach. We can then compare the downstream distances of two point to see which is upstream and which is downstream – the desired outcome!

To solve the problem of finding points on/very close to the rivers we were told to use a (small) buffer around each river reach and then to query the buffer to find the points within it.

This is fine but we are then left with an area (buffer), not a line, so we can't measure a distance along/within this.

Can we use a buffer as described to find out the points on/near the river reach in question and then, somehow (?!) query these points for distance along the river reach (polyline segment)?

Best Answer

This is a perfect task for the linear referencing capabilities in ArcGIS. See the help for Locating Features Along Routes and probe from there. The tools include the ability to turn a layer of points near a route (the river reaches) into "point event" descriptors, which name the route (the reach) and the distance from the beginning of the route. That's exactly what you want. Moreover, these tools allow for a user-settable "tolerance": the points do not have to lie exactly on the routes.

The workflow is:

Make the river reaches "measurable." This sets up the ability to find distances along each reach.
Represent the points as a point layer.
Use the tool to create an "event table" for the points relative to the routes.
Check your work, because you may need to reverse the directions of some of the river reaches, depending on how consistently they were digitized.

Related Solutions

[GIS] Distance between points along polyline ArcGIS for Desktop

Linear reference should do the job, but can be bulky. This is why I am using this script:

# Import arcpy module
import arcpy, os, traceback, sys,time
from arcpy import env
env.overwriteOutput = True

infc = arcpy.GetParameterAsText(0)
routeid = arcpy.GetParameterAsText(1)
outfc=arcpy.GetParameterAsText(2)

fields = [f for f in arcpy.ListFields(infc)]
fieldNames = [f.name for f in arcpy.ListFields(infc)]
# define number of lines
result=arcpy.GetCount_management(infc)
nF=int(result.getOutput(0))
d=arcpy.Describe(outfc)
SR_p=d.spatialReference
d=arcpy.Describe(infc)
SR=d.spatialReference
if SR_p.name<>SR.name:
    arcpy.AddError("\nDifferent projections. Quit\n")
    raise SystemExit('Quit...')

## initial point layer clean-up
fieldsP = [f for f in arcpy.ListFields(outfc)]
fieldNamesP = [f.name for f in arcpy.ListFields(outfc)]
dToLine,dAlong="D_TO_LINE","D_ALONG"
for ent in [routeid,dToLine,dAlong]:
    if ent in fieldNamesP:
        arcpy.DeleteField_management(outfc, ent)
outFolder, theFile="in_memory","intLines"
env.workspace = outFolder
try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    def CalcChainage():
        pC,lineID,dToLine,Chainage=theRow
        shp=dictFeatures[lineID]
        theRow[2]=pC.distanceTo(shp)
        theRow[3]=shp.measureOnLine(pC)
        return

    with arcpy.da.SearchCursor(infc, ("SHAPE@",routeid)) as rows:
        dictFeatures = {}
        for feat,theID in rows:
            dictFeatures[theID]=(feat)
    del rows

    arcpy.SpatialJoin_analysis(outfc, infc, "theJoin", "JOIN_ONE_TO_ONE", "KEEP_ALL", "", "CLOSEST","",dToLine)
    arcpy.AddField_management("theJoin", dAlong, "DOUBLE")

##  calculate chainage
    result=arcpy.GetCount_management("theJoin")
    nF=int(result.getOutput(0))    
    arcpy.AddMessage("\nComputing chainage...")
    arcpy.SetProgressor("step", "", 0, nF,1)
    with arcpy.da.UpdateCursor("theJoin", ("SHAPE@",routeid,dToLine,dAlong)) as rows:
        for theRow in rows:
            CalcChainage()
            rows.updateRow(theRow)
            arcpy.SetProgressorPosition()
    del rows

##defining the type of route id field
    n=fieldNames.index(routeid)
    fType=fields[n].type
    fLength=fields[n].length
    try:
        if fType=="String":
            arcpy.AddField_management(outfc, routeid, "TEXT","","",fLength)
        if fType=="Integer":
            arcpy.AddField_management(outfc, routeid, "LONG")
        if fType=="SmallInteger":
            arcpy.AddField_management(outfc, routeid, "SHORT")
    except: pass
    try:
        arcpy.AddField_management(outfc, dToLine, "DOUBLE")
        arcpy.AddField_management(outfc, dAlong, "DOUBLE")
    except: pass

    arcpy.AddMessage("Transferring calcs to points...\n")            
    vFT=arcpy.da.TableToNumPyArray("theJoin",(routeid,dToLine,dAlong))
    arcpy.SetProgressor("step", "", 0, nF,1)
    with arcpy.da.UpdateCursor(outfc,(routeid,dToLine,dAlong )) as rows:
        n=0                               
        for theRow in rows:
            theRow=vFT[n]
            rows.updateRow(theRow)
            arcpy.SetProgressorPosition()
            n+=1
    del vFT, rows


except NameError, theMessage:
    arcpy.AddMessage (theMessage)
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()

Script expects these parameters:

To get this:

Perhaps you'd like to merge (dissolve) your polylines into single line. Make sure it won't create multipart shape, i.e. ends of segments must be snapped. Script tested on shapefiles.

R sf Package – How to Split Polyline by Point

If you snap the point to the line then you can split the line and then extract the parts from the resulting collection. Use a really small tolerance, I don't know how small it needs to be...

> site_snap = st_snap(site, reach, tol=1e-9)
> parts = st_collection_extract(st_split(reach$geometry, site_snap$geometry),"LINESTRING")

parts is now two features, either side of your point:

> parts
Geometry set for 2 features 
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 688358.8 ymin: 171876.4 xmax: 688465 ymax: 172171.2
epsg (SRID):    NA
proj4string:    +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs
LINESTRING (688400.1 171876.4, 688394.4 171890....
LINESTRING (688448.1 172086.7, 688445.9 172109....

Your original site doesn't intersect the line, so your split returned nothing. Test with:

> st_intersects(site, reach, sparse=FALSE)
      [,1]
[1,] FALSE
> st_intersects(site_snap, reach, sparse=FALSE)
     [,1]
[1,] TRUE

Your original site is 10^-11 units from the line feature:

> st_distance(site$geom, reach$geom)
Units: m
             [,1]
[1,] 9.550647e-11

Best Answer

Related Solutions

[GIS] Distance between points along polyline ArcGIS for Desktop

R sf Package – How to Split Polyline by Point

Related Question