Lignes de détourage «gourmandes» avec polygone

Je souhaite couper un ensemble de polylignes (lignes noires dans l'image ci-dessous) à la limite extérieure d'un polygone. Tout vide dans le polygone doit être ignoré. Ma sortie idéale sont les lignes jaunes en pointillés. Les lignes initiales peuvent être droites ou non. L'image est un exemple simplifié, en réalité le polygone est beaucoup plus complexe et il y a des centaines de lignes. Je ne pense pas qu'une coque convexe fonctionnerait (mais je peux me tromper). Je suis ouvert aux solutions en arcgis, qgis, arcpy, shapely, etc. Le codage serait de préférence en python si je suis ouvert à d'autres options si nécessaire. Arcgis serait également préférable pour faciliter le partage de l'outil par mes collègues, mais ce n'est pas obligatoire.

Le mieux que je puisse penser en ce moment est d'intersecter une ligne individuelle avec le polygone en créant un ensemble de points à toutes les intersections de limites. Triez les points par distance jusqu'au début de la ligne. Les points les plus éloignés et les plus proches (FAC) seront la limite extérieure du polygone. Utilisez ensuite les points FAC pour sélectionner les sommets appropriés à partir de la ligne d'origine et créez la ligne pointillée jaune à partir des points appropriés. Cela devrait fonctionner mais semble plus compliqué que nécessaire.

Quelques réflexions supplémentaires:

• Les lignes sont "suffisamment" linéaires pour qu'un simple calcul de distance entre les points fonctionne, un référencement linéaire ne devrait pas être nécessaire.
• Ce serait facile dans arcpy s'il y avait un outil pour diviser une ligne en un point mais je n'en trouve pas.

Vous pensez à quelqu'un?

Je veux ajouter ma solution pyQGIS, rien d'autre.

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)

Mon cas de test - avant coupure:

Après la coupure:

Pour obtenir l'ensemble complet des attributs des lignes originales, je pense qu'il serait préférable de les joindre au résultat. Sinon, ils doivent être créés dans la section de préparation et définis dans la boucle la plus interne. Mais je n'ai pas testé s'ils réussissent le processus de dissolution ou s'ils se perdent, car en principe ils peuvent avoir des valeurs différentes.

Ce serait facile dans arcpy s'il y avait un outil pour diviser une ligne en un point mais je n'en trouve pas.

Si vous exécutez l'intégration avec les polygones et les lignes en tant qu'entrées, cela ajoutera un sommet à chacun où ils se coupent. (Attention, car Integrate modifie les entrées au lieu de produire de nouvelles sorties.)

Une fois que vous êtes sûr qu'il y a des sommets coïncidents, vous pouvez parcourir les sommets de la ligne et tester pour voir si chacun touche l'autre entité. Dans la liste ordonnée des sommets qui se touchent, prenez le minimum et le maximum de l'ensemble. Ensuite, faites deux lignes à partir de chaque entité, A: (début, ..., min) et B: (max, ..., fin).

Une autre option, bien que je ne sois pas sûr qu'ArcPy préserve l'ordre des pièces d'entité en fonction de l'ordre des sommets dans l'objet d'entrée, serait d'exécuter le clip tel quel. Pour la ligne médiane dans votre exemple, cela devrait se traduire par une fonctionnalité en plusieurs parties en trois parties. Selon l'ordre, vous pouvez parcourir chaque ligne en plusieurs parties produite par Clip et supprimer tout sauf la première et la dernière partie de la fonction en plusieurs parties.

Il y a trois problèmes à résoudre dans ce cas:

• des trous
• Lignes entre polygones
• Lignes de fin

des trous

Étant donné que toute ligne dans un trou sera conservée, supprimez les trous des polygones. Dans le script ci-dessous, je le fais en utilisant des curseurs et des géométries.

Lignes entre polygones

Les lignes qui touchent deux polygones doivent être supprimées. Dans le script ci-dessous, je le fais en effectuant une jointure spatiale de one to many, avec mes lignes comme classe d'entités en entrée et mes polygones comme classe d'entités de jointure. Toute ligne générée deux fois touche deux polygones et est supprimée.

Lignes de fin

Pour supprimer des lignes qui ne touchent un polygone qu'à une extrémité, je convertis les lignes en points d'extrémité. J'utilise ensuite des couches d'entités et des sélections pour déterminer quels points d'extrémité sont des flottants. Je sélectionne les points d'extrémité qui coupent les polygones. Je change ensuite ma sélection. Cela sélectionne les points d'extrémité qui n'intersectent pas les polygones. Je sélectionne n'importe quelle ligne qui coupe ces points sélectionnés et les supprime.

Résultat

Hypothèses

• Les entrées sont des classes d'entités de géodatabase fichier
• La licence avancée ArcGIS est disponible (en raison de an eraseet a feature vertices to points)
• Les lignes continues et connectées sont une caractéristique unique
• Les polygones ne se chevauchent pas
• Il n'y a pas de polygones en plusieurs parties

Scénario

Le script ci-dessous génère une classe d' _GreedyClipentités avec le nom de votre classe d' entités linéaires plus , dans la même géodatabase que votre classe d'entités linéaires. Un emplacement d'espace de travail est également nécessaire.

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