"Gierige" Schnittlinien mit Polygon

Ich möchte eine Reihe von Polylinien (schwarze Linien im Bild unten) an der Außengrenze eines Polygons befestigen. Alle Hohlräume innerhalb des Polygons sollten ignoriert werden. Meine ideale Ausgabe sind die gestrichelten gelben Linien. Die Anfangslinien können gerade sein oder nicht. Das Bild ist ein vereinfachtes Beispiel, in Wirklichkeit ist das Polygon viel komplexer und es gibt Hunderte von Linien. Ich glaube nicht, dass eine konvexe Hülle funktionieren würde (aber ich könnte mich irren). Ich bin offen für Lösungen in Arcgis, QGIS, Arcpy, Shapely usw. Die Codierung erfolgt vorzugsweise in Python, da ich bei Bedarf für andere Optionen offen bin. Arcgis wäre auch vorzuziehen, um meinen Mitarbeitern das Teilen des Tools zu erleichtern, ist jedoch keine Voraussetzung.

Das Beste, was ich mir derzeit vorstellen kann, ist, eine einzelne Linie mit dem Polygon zu schneiden und an allen Grenzüberschneidungen eine Reihe von Punkten zu erstellen. Sortieren Sie die Punkte nach Abstand zum Zeilenanfang. Der am weitesten und am nächsten gelegene Punkt (FAC) ist die äußere Grenze des Polygons. Verwenden Sie dann die FAC-Punkte, um die richtigen Scheitelpunkte aus der ursprünglichen Linie auszuwählen und die gelbe gestrichelte Linie aus den entsprechenden Punkten zu erstellen. Es sollte funktionieren, scheint aber komplizierter als nötig.

Ein paar zusätzliche Gedanken:

• Die Linien sind "genug" linear, damit eine einfache Abstandsberechnung zwischen Punkten funktioniert. Eine lineare Referenzierung sollte nicht erforderlich sein.
• Dies wäre in arcpy einfach, wenn es ein Werkzeug zum Teilen einer Linie an einem Punkt gäbe, aber ich kann keine finden.

Gedanken jemand?

Ich möchte meine pyQGIS-Lösung einwerfen, sonst nichts.

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)

Mein Testfall - vor dem Clipping:

Nach dem Abschneiden:

Um die vollständigen Attribute der ursprünglichen Zeilen zu erhalten, ist es meiner Meinung nach am besten, sie mit dem Ergebnis zu verbinden. Andernfalls müssen die im Vorbereitungsabschnitt erstellt und in die innerste Schleife gesetzt werden. Aber ich habe nicht getestet, ob sie den Auflösungsprozess bestehen oder ob sie verloren gehen, weil sie im Prinzip unterschiedliche Werte haben könnten.

Dies wäre in arcpy einfach, wenn es ein Werkzeug zum Teilen einer Linie an einem Punkt gäbe, aber ich kann keine finden.

Wenn Sie Integrieren mit den Polygonen und Linien als Eingaben ausführen, wird jedem, an dem sie sich schneiden, ein Scheitelpunkt hinzugefügt. (Vorsicht, da Integrate Eingaben ändert, anstatt neue Ausgaben zu erzeugen.)

Sobald Sie sicher sind, dass es übereinstimmende Scheitelpunkte gibt, können Sie über die Scheitelpunkte der Linie iterieren und testen, ob jeder das andere Merkmal berührt. Nehmen Sie aus der geordneten Liste der Scheitelpunkte, die sich berühren, das Minimum und das Maximum aus der Menge. Machen Sie dann zwei Zeilen aus jedem Feature, A: (Start, ..., Min) und B: (Max, ..., Ende).

Eine andere Option, obwohl ich nicht sicher bin, ob ArcPy die Reihenfolge der Feature-Teile basierend auf der Reihenfolge der Scheitelpunkte im Eingabeobjekt beibehält, besteht darin, den Clip unverändert auszuführen. Für die mittlere Zeile in Ihrem Beispiel sollte sich eine mehrteilige Funktion mit drei Teilen ergeben. Abhängig von der Reihenfolge können Sie jede von Clip erzeugte mehrteilige Zeile durchlaufen und alle bis auf den ersten und letzten Teil der mehrteiligen Out-Funktion entfernen.

In diesem Fall gibt es drei Probleme:

• Löcher
• Linien zwischen Polygonen
• Endzeilen

Löcher

Entfernen Sie Löcher aus Polygonen, da jede Linie innerhalb eines Lochs beibehalten wird. Im folgenden Skript verwende ich Cursor und Geometrien.

Linien zwischen Polygonen

Linien, die zwei Polygone berühren, müssen entfernt werden. Im folgenden Skript führe ich dazu eine räumliche Verknüpfung von durch one to many, wobei meine Linien als Eingabe-Feature-Class und meine Polygone als Join-Feature-Class dienen. Jede Linie, die zweimal generiert wird, berührt zwei Polygone und wird entfernt.

Endzeilen

Um Linien zu entfernen, die nur ein Polygon an einem Ende berühren, konvertiere ich Linien in Endpunkte. Ich benutze dann Feature-Layer und Auswahlen, um zu bestimmen, welche Endpunkte Floater sind. Ich wähle die Endpunkte aus, die die Polygone schneiden. Ich wechsle dann meine Auswahl. Dadurch werden Endpunkte ausgewählt, die Polygone nicht schneiden. Ich wähle eine Linie aus, die diese ausgewählten Punkte schneidet, und lösche sie.

Ergebnis

Annahmen

• Eingaben sind Feature-Classes für Datei-Geodatabase
• Die erweiterte ArcGIS-Lizenz ist verfügbar (aufgrund von a eraseund a feature vertices to points).
• Durchgehende, verbundene Leitungen sind ein einziges Merkmal
• Polygone überlappen sich nicht
• Es gibt keine mehrteiligen Polygone

Skript

Das folgende Skript gibt eine Feature-Class mit dem Namen Ihrer Line-Feature-Class plus _GreedyClipin derselben Geodatabase wie Ihre Line-Feature-Class aus. Ein Arbeitsbereich wird ebenfalls benötigt.

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