"Greedy" clipping lines with polygon

Question

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.

A few additional thoughts:

• 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?

Answer 1

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()
clipped.addAttribute(QgsField('fid', QVariant.Int))
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)

            prov.addFeatures(feats)

# expose layer
clipped.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([clipped])

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

Answer 2

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.

Answer 3

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
                    array.add (pnt)
                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)