„Chciwe” linie przycinające z wielokątem

Chcę przyciąć zestaw polilinii (czarne linie na obrazku poniżej) do zewnętrznej granicy wielokąta. Wszelkie puste przestrzenie w obrębie wielokąta należy zignorować. Moje idealne wyjście to przerywane żółte linie. Linie początkowe mogą, ale nie muszą być proste. Obraz jest uproszczonym przykładem, w rzeczywistości wielokąt jest o wiele bardziej złożony i zawiera setki linii. Nie sądzę, by wypukły kadłub działał (ale mogę się mylić). Jestem otwarty na rozwiązania w arcgis, qgis, arcpy, foremne itp. Kodowanie najlepiej byłoby w pythonie, ponieważ jestem otwarty na inne opcje, jeśli to konieczne. Arcgis byłby również lepszy, aby ułatwić moim współpracownikom udostępnianie tego narzędzia, ale nie jest to wymagane.

Najlepsze, co mogę teraz wymyślić, to przecięcie pojedynczej linii z wielokątem, tworząc zestaw punktów na wszystkich przecięciach granicznych. Posortuj punkty według odległości do początku linii. Najdalsze i najbliższe punkty (FAC) będą zewnętrzną granicą wielokąta. Następnie użyj punktów FAC, aby wybrać odpowiednie wierzchołki z oryginalnej linii i utwórz żółtą linię przerywaną z odpowiednich punktów. Powinien działać, ale wydaje się bardziej skomplikowany niż to konieczne.

Kilka dodatkowych myśli:

• Linie są „wystarczająco” liniowe, aby wystarczyło proste obliczenie odległości między punktami, odniesienie liniowe nie powinno być konieczne.
• Byłoby to łatwe, gdyby istniało narzędzie do dzielenia linii w punkcie, ale nie mogę jej znaleźć.

Myśli ktoś?

Chcę wrzucić moje rozwiązanie pyQGIS, nic więcej.

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)

Mój przypadek testowy - przed wycięciem:

Po wycięciu:

Aby uzyskać pełny zestaw atrybutów oryginalnych linii, myślę, że najlepiej byłoby połączyć je z wynikiem. W przeciwnym razie należy je utworzyć w sekcji przygotowania i ustawić w najbardziej wewnętrznej pętli. Ale nie sprawdziłem, czy przejdą proces rozpuszczania, czy też się zgubią, ponieważ w zasadzie mogą mieć różne wartości.

Byłoby to łatwe, gdyby istniało narzędzie do dzielenia linii w punkcie, ale nie mogę jej znaleźć.

Jeśli uruchomisz Zintegruj z wielokątami i liniami jako danymi wejściowymi, doda wierzchołek do każdego miejsca, w którym się przecinają. (Ostrożnie, ponieważ funkcja Integrate modyfikuje dane wejściowe zamiast tworzyć nowe dane wyjściowe).

Gdy upewnisz się, że są zbieżne wierzchołki, możesz iterować wierzchołki linii i sprawdzić, czy oba dotykają drugiej cechy. Z uporządkowanej listy wierzchołków, które się dotykają, weź minimum i maksimum z zestawu. Następnie utwórz dwie linie z każdej operacji: A: (początek, ..., min) i B: (maks., ..., koniec).

Inną opcją, chociaż nie jestem pewien, czy ArcPy zachowuje porządkowanie części elementu na podstawie kolejności wierzchołków w obiekcie wejściowym, byłoby uruchomienie klipu w niezmienionej postaci. W środkowej linii w twoim przykładzie powinna to być funkcja wieloczęściowa z trzema częściami. W zależności od kolejności można iterować każdą linię wieloczęściową utworzoną przez Clip i usuwać wszystkie oprócz pierwszej i ostatniej części funkcji wieloczęściowej.

W tym przypadku istnieją trzy problemy:

• Otwory
• Linie między wielokątami
• Linie końcowe

Otwory

Ponieważ dowolna linia w otworze zostanie zachowana, usuń otwory z wielokątów. W poniższym skrypcie robię to za pomocą kursorów i geometrii.

Linie między wielokątami

Linie dotykające dwóch wielokątów muszą zostać usunięte. W poniższym skrypcie robię to, wykonując sprzężenie przestrzenne one to manyz moimi liniami jako moją klasą elementów wejściowych, a wielokątami jako moją klasą elementów złączonych. Każda linia, która zostanie wygenerowana dwukrotnie, dotyka dwóch wielokątów i jest usuwana.

Linie końcowe

Aby usunąć linie dotykające tylko wielokąta na jednym końcu, przekształcam linie w punkty końcowe. Następnie używam warstw obiektów i selekcji, aby ustalić, które punkty końcowe są zmiennoprzecinkowe. Wybieram punkty końcowe, które przecinają wielokąty. Następnie zmieniam wybór. To wybiera punkty końcowe, które nie przecinają wielokątów. Zaznaczam dowolną linię, która przecina te wybrane punkty i usuwam je.

Wynik

Założenia

• Dane wejściowe to klasy obiektów geobazy danych
• Dostępna jest zaawansowana licencja ArcGIS (ze względu na erasei feature vertices to points)
• Ciągłe, połączone linie są pojedynczą funkcją
• Wieloboki nie nakładają się
• Nie ma wieloczęściowych wielokątów

Scenariusz

Poniższy skrypt generuje klasę obiektów o nazwie twojej klasy obiektów liniowych plus _GreedyClip, w tej samej geobazie, co twoja klasa obiektów liniowych. Potrzebna jest również lokalizacja obszaru roboczego.

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