R sf Package – How to Split Polyline by Point

pointpolyline-creationrsfsplitting

I'm learning the sf package in R and would like to locate about 450 sites along the nearest reaches of a river network and to obtain the upstream and downstream distances along the reach where the site is located. I've include one simple feature POINT (site) and LINESTRING or polyline (reach) to provide a reproducible example at the bottom of my post.

Given a point located exactly along the polyline, how would I split the polyline at that point to obtain two polylines? I would like to obtain the lengths of the split polylines both upstream and downstream of the point in order to edit a graph representation of the whole river.

A similar question on SE is submitted here but involves splitting a polyline using polygons. My latest attempt with st_split() yields a GEOMETRYCOLLECTION, but it seems that the LINESTRING coordinates are not listed separately:

library(sf)
# Get the site/reach objects from dput()...
split.reach <- st_split(reach$geometry, site$geometry)
> split.reach
Geometry set for 1 feature 
geometry type:  GEOMETRYCOLLECTION
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
GEOMETRYCOLLECTION (LINESTRING (688400.1 171876...

To be clear, in this specific problem I'm interested in lengths along an individual river reach (i.e., an individual row or LINESTRING in an 'sf' object) that is split by a point and not along the whole river network. I would also like to avoid writing code that depends on a large number of spatial packages (i.e., stick with sf)

The dput for a site and reach as an example:

site <- structure(list(SiteId = "PF_00183", SampId = "PF_00183 _ 2014-10-15", 
    Longitude.In = 8.59713558862371, Latitude.In = 46.6953718769341, 
    X.In = 688524.841170469, Y.In = 172082.839248429, ReachId = 125010L, 
    Distance = 76.8220253856169, Longitude.Near = 8.59613333686445, 
    Latitude.Near = 46.6954168201561, geometry = structure(list(
        structure(c(688448.125, 172086.703125001), class = c("XY", 
        "POINT", "sfg"))), n_empty = 0L, precision = 0, crs = structure(list(
        epsg = NA_integer_, 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"), .Names = c("epsg", 
    "proj4string"), class = "crs"), bbox = structure(c(688448.125, 
    172086.703125001, 688448.125, 172086.703125001), .Names = c("xmin", 
    "ymin", "xmax", "ymax"), class = "bbox"), class = c("sfc_POINT", 
    "sfc"))), .Names = c("SiteId", "SampId", "Longitude.In", 
"Latitude.In", "X.In", "Y.In", "ReachId", "Distance", "Longitude.Near", 
"Latitude.Near", "geometry"), row.names = c(NA, -1L), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = c("SiteId", "SampId", "Longitude.In", 
"Latitude.In", "X.In", "Y.In", "ReachId", "Distance", "Longitude.Near", 
"Latitude.Near")), class = c("sf", "data.table", "data.frame"
))
reach <- structure(list(FNODE_ = 148558, TNODE_ = 148142, LENGTH = 343.00790968, 
    ReachId = 125010L, geometry = structure(list(structure(c(688400.125, 
    688394.375, 688385.875, 688375.6875, 688363.625, 688359.8125, 
    688358.8125, 688361.125, 688375.6875, 688389.375, 688411.1875, 
    688432.875, 688445.1875, 688448.125, 688445.875, 688446.625, 
    688448.125, 688455.375, 688465, 171876.40625, 171890.40625, 
    171906.203125, 171916.703125, 171929.09375, 171944.703125, 
    171952.59375, 171959.90625, 171980.296875, 171997.09375, 
    172017.796875, 172041.90625, 172063.5, 172086.703125, 172109.90625, 
    172129, 172138.5, 172153.5, 172171.203125), .Dim = c(19L, 
    2L), class = c("XY", "LINESTRING", "sfg"))), class = c("sfc_LINESTRING", 
    "sfc"), precision = 0, bbox = structure(c(688358.8125, 171876.40625, 
    688465, 172171.203125), .Names = c("xmin", "ymin", "xmax", 
    "ymax"), class = "bbox"), crs = structure(list(epsg = NA_integer_, 
        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"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("FNODE_", 
"TNODE_", "LENGTH", "ReachId", "geometry"), row.names = 125010L, class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_), .Names = c("FNODE_", 
"TNODE_", "LENGTH", "ReachId"), .Label = c("constant", "aggregate", 
"identity"), class = "factor"))

Best Answer

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

Related Solutions

[GIS] QGIS – split polyline

Everybody thanks a lot. Finally, I did it. Was to easy (when you know how). In QGIS 2.8.1 the tool is SPLIT LINES WITH LINES (on the V. Overlay, on Geoalgorithms, from Advanced Interface). I came there thanks to Chris W advice.

ArcPy – Standalone Python Script to Split a Polyline With a Point Layer

Here's a script that might do it for you. The native arcpy.Geometry class has a method called "cut" that will cut any feature using another polyline. Unfortunately, since you're using "points", we have to make "fake" lines out of these points. I essentially made a scratch polyline with the points [(Point.X+10, Point.Y+10), (Point.X-10, Point.Y-10)] - e.g. a diagonal line for every point to perform each cut. When you have more than one point intersecting a line, you have to re-cut the lines recursively smaller-and-smaller until they can't be cut up any further (without producing zero-length geometries).

A few caveats I should mention:

The script assumes there are no multi-part geometries in both feature classes. I always used Point.firstPoint for each point, so any additional multi-points thereafter will be left out. I haven't tested the effect of cutting a multi-part line, but I assume the descendants will also stay multi-part.
When you have a perfectly diagonal segment of a polyline, it won't get cut because our "fake cut line" is also diagonal. You may have to split these manually, if they exist.
The output won't have the attributes from the original polylines, but a Spatial Join can take care of that.
Lastly, of course, make sure both feature classes are in the same coordinate system.

You just have to make the skeleton tool in the toolbox that accepts 3 parameters:

Polyline Feature Class (Direction = Input, Type = Feature Class, Filter = Line)
Point Feature Class (Direction = Input, Type = Feature Class, Filter = Point)
Output Feature Class (Direction = Output, Type = Feature Class)

Here are some examples: enter image description here

And the code:

__author__ = "John K. Tran"
__contact__ = "jtran20@masonlive.gmu.edu"
__version__ = "4.0"
__created__ = "7/16/15"
__credits__ = "http://gis.stackexchange.com/questions/154708/standalone-python-script-to-split-a-polyline-with-a-point-layer"

"""
'Cuts', 'splits' or 'dices' a polyline feature class using a point feature class.
The output schema (i.e. list of fields) will be blank, but you can use a spatial join
to repopulate the attributes to the output feature class.
"""

import arcpy
import os
import sys

arcpy.env.overwriteOutput = True
arcpy.SetProgressor("default", "Firing up script...")

linefc = arcpy.GetParameterAsText(0) # Your polyline feature class
pointfc = arcpy.GetParameterAsText(1) # Your point feature class
output = arcpy.GetParameterAsText(2) # New output feature class.

spatialref1 = arcpy.Describe(pointfc).spatialReference
spatialref2 = arcpy.Describe(linefc).spatialReference
assert spatialref1.name == spatialref2.name, "Ensure both feature classes have the same projected coordinate system"

# If the Advanced license is available, just use the normal geoprocessing tools.
if arcpy.ProductInfo() in [u'ArcInfo', u'ArcServer']:
    arcpy.SetProgressorLabel("Splitting lines at points")
    arcpy.SplitLineAtPoint_management(linefc, pointfc, output, 1.0)
    arcpy.SetProgressorLabel("Deleting duplicate slices")
    outshapefieldname = arcpy.Describe(output).shapeFieldName
    arcpy.DeleteIdentical_management(output, [outshapefieldname])
    sys.exit(0)

# If not, proceed with script.
# Set up some preliminary variables to start.
arcpy.SetProgressorLabel("Gathering geometries")
points = [row[0] for row in arcpy.da.SearchCursor(pointfc, "SHAPE@")]
lines = [row[0] for row in arcpy.da.SearchCursor(linefc, "SHAPE@")]
cutlines = list()
iterations = 0

# Defining a function for the cut will let us reuse it again for subsequent cuts.
def CutLines():
    arcpy.SetProgressorLabel("Running cut: Iteration {0}".format(str(iterations)))
    for line in lines[:]:
        iscut = "Invalid"
        if line.length > 0.0: # Make sure it's not an empty geometry.
            iscut = "Not Cut"
            for point in points:
                if line.distanceTo(point) < 1.0: # Even "coincident" points can show up as spatially non-coincident in their floating-point XY values, so we set up a tolerance.
                    snappoint = line.snapToLine(point).firstPoint # To ensure coincidence, snap the point to the line before proceeding.
                    if not (snappoint.equals(line.lastPoint) and snappoint.equals(line.firstPoint)): # Make the sure the point isn't on a line endpoint, otherwise cutting will produce an empty geometry.
                        cutline1, cutline2 = line.cut(arcpy.Polyline(arcpy.Array([arcpy.Point(snappoint.X+10.0, snappoint.Y+10.0), arcpy.Point(snappoint.X-10.0, snappoint.Y-10.0)]), spatialref1)) # Cut the line.
                        if cutline1.length > 0.0 and cutline2.length > 0.0: # Make sure both descendents have non-zero geometry.
                            lines.insert(0, cutline1) # Feed the cut lines back into the "line" list to be recut.
                            lines.insert(0, cutline2) # The cut loop will only exit when all lines cannot be recut smaller and smaller (without producing zero-length geometries).
                            iscut = "Cut"
        if iscut == "Not Cut":
            cutlines.insert(0, line)
        lines.remove(line)

# Perform the pseudo-recurive Cut loop until the polyline FC can't be cut further.
while lines:
    CutLines()
    iterations += 1
    if iterations > 500: # Fail-safe to stop an infinite loop if something goes wrong (or if more than 500 points intersect a single line).
        break

# Create the output feature class.
arcpy.SetProgressorLabel("Creating output feature class")
arcpy.CreateFeatureclass_management(os.path.dirname(output), os.path.basename(output), "POLYLINE", None, None, None, spatialref2)

# Insert each cut line into the new feature class.
arcpy.SetProgressorLabel("Inserting cut lines")
with arcpy.da.InsertCursor(output, "SHAPE@") as insertcursor:
    for cutline in cutlines:
        insertcursor.insertRow((cutline,))

# Delete duplicates by comparing their geometry objects.
arcpy.SetProgressorLabel("Deleting duplicates")
crows = [row for row in arcpy.da.SearchCursor(output, ["OID@", "SHAPE@"])]
deleteOIDs = set() # Retain the OIDs for deletion, performed in the next cursor.
for drow in crows: # Compare each geometry to every other geometry in the output FC.
    dOID = drow[0]
    dgeom = drow[1]
    if dOID not in deleteOIDs:
        for srow in crows:
            sOID = srow[0]
            sgeom = srow[1]
            if sOID != dOID and sOID not in deleteOIDs: # Prevent redundant comparisons, only compare new combinations.
                dextent = dgeom.extent
                sextent = sgeom.extent
                if dextent.XMin == sextent.XMin and dextent.XMax == sextent.XMax and dextent.YMin == sextent.YMin and dextent.YMax == sextent.YMax: # See if their extents match first.
                    if dgeom.equals(sgeom): # Then see if they equal each other.
                        deleteOIDs.add(sOID) # If so, add the OID to the list of OIDs to delete.
with arcpy.da.UpdateCursor(output, "OID@") as deletecursor: # Delete the records here using your list of "DeleteOIDs".
    for delrow in deletecursor:
        delOID = delrow[0]
        if delOID in deleteOIDs:
            deletecursor.deleteRow()

# Clean up leftover variables.
namespace = dir()
for var in ["linecursor", "pointcursor", "insertcursor", "deletecursor"]:
    if var in namespace:
        exec "del " + var

arcpy.ResetProgressor()

If you get any errors, let me know.

Best Answer

Related Solutions

[GIS] QGIS – split polyline

ArcPy – Standalone Python Script to Split a Polyline With a Point Layer

Related Question