# "Greedy" clipping lines with polygon

I am wishing to clip a set of polylines (black lines in image below) to the outer boundary of a polygon. Any voids within the polygon should be ignored. My ideal output are the dashed yellow lines. The initial lines may or may not be straight. The image is a simplified example, in reality the polygon is much more complex and there are hundreds of lines. I don't think a convex hull would work (but I might be wrong). I'm open to solutions in arcgis, qgis, arcpy, shapely, etc. Coding would preferably be in python by I'm open to other options if necessary. Arcgis would also be preferable to make it easier for my coworkers to share the tool but is not a requirement.

The best i can think of right now is to intersect an individual line with the polygon creating a set of points at all boundary intersections. Sort the points by distance to the start of the line. The furthest and closest (FAC) points will be the outer boundary of the polygon. Then use the FAC points to select the proper vertices from the original line and create the yellow dashed line from the appropriate points. It should work but seems more complicated than necessary.

• The lines are linear "enough" that a simple distance calculation between points should work, linear referencing shouldn't be necessary.
• This would be easy in arcpy if there was a tool to split a line at a point but I cannot find one.

Thoughts anyone?

• +1, interesting problem! I'm keen to see what solutions are available =) Commented Jan 12, 2016 at 16:19
• Only your middle line is hard to achieve - the top and bottom just come from a clip after filling any voids. Consequently, I think you should focus your question on that and narrow its scope to just ArcPy if that is your preferred tool. You can always ask about another tool, if that does not yield a solution.
– PolyGeo
Commented Jan 12, 2016 at 19:57
• do lines cross multiple polygons? Commented Jan 13, 2016 at 21:43
• Emil, lets assume that lines may cross over multiple polygons. However, other than geometry there is no difference between polygons so they can be dissolved, merged into a multipart feature, etc. if that makes the algorithm easier. A line crossing over multiple polygons would likely be rare and that can be a flagged case to be dealt with by hand if necessary. Commented Jan 13, 2016 at 22:20

I want to throw in my pyQGIS solution, nothing else.

``````from PyQt4.QtCore import QVariant
from qgis.analysis import QgsGeometryAnalyzer

# get layers
lines = QgsMapLayerRegistry.instance().mapLayersByName('lines')[0]
clipper = QgsMapLayerRegistry.instance().mapLayersByName('clipper')[0]

# prepare result layer
clipped = QgsVectorLayer('LineString?crs=epsg:4326', 'clipped', 'memory')
clipped.startEditing()
fni = clipped.fieldNameIndex('fid')
clipped.commitChanges()

prov = clipped.dataProvider()
fields = prov.fields()

for line in lines.getFeatures():
# to increase performance filter possible clippers
clippers = clipper.getFeatures(QgsFeatureRequest().setFilterRect(line.geometry().boundingBox()))
for clip in clippers:
# split the line
line1 = line.geometry().splitGeometry(clip.geometry().asPolygon()[0], True)
feats = []
# get the split points
vertices = [QgsPoint(vert[0], vert[1]) for vert in line1[2]]
for part in line1[1]:
# for each split part check, if first AND last vertex equal to split points
if part.vertexAt(0) in vertices and part.vertexAt(len(part.asPolyline())-1) in vertices:
# if so create feature and set fid to original line's id
feat = QgsFeature(fields)
feat.setAttributes([line.id()])
feat.setGeometry(part)
feats.append(feat)

# expose layer
clipped.updateExtents()

# now dissolve lines having the same value in field fni: here original line's id
diss = QgsGeometryAnalyzer()
diss.dissolve(clipped, 'E:\\clipped.shp', uniqueIdField=fni)
``````

My test case - before clipping:

After clipping:

To get the full set of attributes of the original lines I think it would be the best to join them with the result. Otherwise the have to be created in the prepare section, and set in the most inner loop. But I haven't tested if they pass the dissolve process or if they get lost, because in principle they could have different values.

• Very concise answer. How do QGIS screen shots always look like QGIS? Commented Jan 26, 2016 at 15:31

This would be easy in arcpy if there was a tool to split a line at a point but I cannot find one.

If you run Integrate with the polygons and lines as inputs, it will add a vertex to each where they intersect. (Careful, as Integrate modifies inputs instead of producing new outputs.)

Once you are sure there are coincident vertices, you can iterate over the line's vertices and test to see if each touches the other feature. From the ordered list of vertices that do touch, take the minimum and the maximum from the set. Then, make two lines from each feature, A: (start, ..., min) and B: (max, ..., end).

Another option, though I'm not sure if ArcPy preserves the feature part ordering based on the ordering of vertices in the input object, would be to run the clip as-is. For the middle line in your example, it should result in a multipart feature with three parts. Depending on the ordering, you could iterate over every multipart line produced by Clip and remove all but the first and last part of the out multipart feature.

There are three issues to contend with in this case:

• Holes
• Lines between polygons
• End lines

Holes

Since any line within a hole will be maintained, remove holes from polygons. In the script below I do so by use of cursors and geometries.

Lines between polygons

Lines that touch two polygons need to be removed. In the script below I do so by performing a spatial join of `one to many`, with my lines as my input feature class and my polygons as my join feature class. Any line that gets generated twice touches two polygons and is removed.

End lines

To remove lines that only touch a polygon on one end, I convert lines to end points. I then make use of feature layers and selections to determine which end points are floaters. I select the end points that intersect the polygons. I then switch my selection. This selects end points that don't intersect polygons. I select any line that intersects these selected points and delete them.

Result

Assumptions

• Inputs are file geodatabase feature classes
• ArcGIS advanced license is available (due to an `erase` and a `feature vertices to points`)
• Continuous, connected lines are a single feature
• Polygons do no overlap
• There are no multipart polygons

Script

The script below outputs a feature class with the name of your line feature class plus `_GreedyClip`, in the same geodatabase as your line feature class. A workspace location is also needed.

``````#input polygon feature class
polyFc = r"C:\Users\e1b8\Desktop\E1B8\Workspace\Workspace.gdb\testPolygon2"
#input line feature class
lineFc = r"C:\Users\e1b8\Desktop\E1B8\Workspace\Workspace.gdb\testLine"
#workspace
workspace = r"in_memory"

print "importing"
import arcpy
import os

#generate a unique ArcGIS file name
def UniqueFileName(location = "in_memory", name = "file", extension = ""):
if extension:
outName = os.path.join (location, name + "." + extension)
else:
outName = os.path.join (location, name)
i = 0
while arcpy.Exists (outName):
i += 1
if extension:
outName = os.path.join (location, "{0}_{1}.{2}".format (name, i, extension))
else:
outName = os.path.join (location, "{0}_{1}".format (name, i))
return outName

#remove holes from polygons
def RemoveHoles (inFc, workspace):
outFc = UniqueFileName (workspace)
array = arcpy.Array ()
sr = arcpy.Describe (inFc).spatialReference
outPath, outName = os.path.split (outFc)
arcpy.CreateFeatureclass_management (outPath, outName, "POLYGON", spatial_reference = sr)
with arcpy.da.InsertCursor (outFc, "SHAPE@") as iCurs:
with arcpy.da.SearchCursor (inFc, "SHAPE@") as sCurs:
for geom, in sCurs:
try:
part = geom.getPart (0)
except:
continue
for pnt in part:
if not pnt:
break
polygon = arcpy.Polygon (array)
array.removeAll ()
row = (polygon,)
iCurs.insertRow (row)
del iCurs
del sCurs
return outFc

#split line fc by polygon fc
def SplitLinesByPolygon (lineFc, polygonFc, workspace):
#clip
clipFc = UniqueFileName(workspace)
arcpy.Clip_analysis (lineFc, polygonFc, clipFc)
#erase
eraseFc = UniqueFileName(workspace)
arcpy.Erase_analysis (lineFc, polygonFc, eraseFc)
#merge
mergeFc = UniqueFileName(workspace)
arcpy.Merge_management ([clipFc, eraseFc], mergeFc)
#multipart to singlepart
outFc = UniqueFileName(workspace)
arcpy.MultipartToSinglepart_management (mergeFc, outFc)
#delete intermediate data
for trash in [clipFc, eraseFc, mergeFc]:
arcpy.Delete_management (trash)
return outFc

#remove lines between two polygons and end lines
def RemoveLines (inFc, polygonFc, workspace):
#check if "TARGET_FID" is in fields
flds = [f.name for f in arcpy.ListFields (inFc)]
if "TARGET_FID" in flds:
#delete "TARGET_FID" field
arcpy.DeleteField_management (inFc, "TARGET_FID")
#spatial join
sjFc = UniqueFileName(workspace)
arcpy.SpatialJoin_analysis (inFc, polygonFc, sjFc, "JOIN_ONE_TO_MANY")
#list of TARGET_FIDs
targetFids = [fid for fid, in arcpy.da.SearchCursor (sjFc, "TARGET_FID")]
#target FIDs with multiple occurances
deleteFids = [dFid for dFid in targetFids if targetFids.count (dFid) > 1]
if deleteFids:
#delete rows with update cursor
with arcpy.da.UpdateCursor (inFc, "OID@") as cursor:
for oid, in cursor:
if oid in deleteFids:
cursor.deleteRow ()
del cursor
#feature vertices to points
vertFc = UniqueFileName(workspace)
arcpy.FeatureVerticesToPoints_management (inFc, vertFc, "BOTH_ENDS")
#select points intersecting polygons
arcpy.MakeFeatureLayer_management (vertFc, "vertLyr")
arcpy.SelectLayerByLocation_management ("vertLyr", "", polygonFc, "1 FEET")
#switch selection
arcpy.SelectLayerByAttribute_management ("vertLyr", "SWITCH_SELECTION")
arcpy.MakeFeatureLayer_management (inFc, "lineLyr")
#check for selection
if arcpy.Describe ("vertLyr").FIDSet:
#select lines by selected points
arcpy.SelectLayerByLocation_management ("lineLyr", "", "vertLyr", "1 FEET")
#double check selection (should always have selection)
if arcpy.Describe ("lineLyr").FIDSet:
#delete selected rows
arcpy.DeleteFeatures_management ("lineLyr")

#delete intermediate data
for trash in [sjFc, "vertLyr", "lineLyr"]:
arcpy.Delete_management (trash)

#main script
def main (polyFc, lineFc, workspace):

#remove holes
print "removing holes"
holelessPolyFc = RemoveHoles (polyFc, workspace)

#split line at polygons
print "splitting lines at polygons"
splitFc = SplitLinesByPolygon (lineFc, holelessPolyFc, workspace)

#delete unwanted lines
print "removing unwanted lines"
RemoveLines (splitFc, polyFc, workspace)

#create output feature class
outFc = lineFc + "_GreedyClip"
outFcPath, outFcName = os.path.split (outFc)
outFc = UniqueFileName (outFcPath, outFcName)
arcpy.CopyFeatures_management (splitFc, outFc)
print "created:"
print outFc
print
print "cleaning up"
#delete intermediate data
for trash in [holelessPolyFc, splitFc]:
arcpy.Delete_management (trash)

print "done"

if __name__ == "__main__":
main (polyFc, lineFc, workspace)
``````
• Nice solution Emil. That is less code than I ended up with. Commented Jan 26, 2016 at 15:30