«Жадные» линии отсечения с многоугольником

9

Я хочу обрезать набор полилиний (черные линии на изображении ниже) по внешней границе многоугольника. Любые пустоты внутри многоугольника должны игнорироваться. Мой идеальный выход - пунктирные желтые линии. Начальные линии могут быть или не быть прямыми. Изображение является упрощенным примером, в действительности многоугольник намного сложнее и содержит сотни линий. Я не думаю, что выпуклый корпус будет работать (но я могу ошибаться). Я открыт для решений в arcgis, qgis, arcpy, shapely и т. Д. Кодирование желательно выполнять на python, поскольку я открыт для других вариантов, если это необходимо. Arcgis также предпочтительнее, чтобы моим коллегам было проще делиться инструментом, но это не является обязательным требованием.

Лучшее, о чем я могу думать сейчас, - это пересечь отдельную линию с многоугольником, создавая набор точек на всех пересечениях границ. Сортировка точек по расстоянию до начала линии. Самая дальняя и ближайшая (FAC) точки будут внешней границей многоугольника. Затем используйте точки FAC, чтобы выбрать правильные вершины из исходной линии и создать желтую пунктирную линию из соответствующих точек. Это должно работать, но кажется более сложным, чем необходимо.

Несколько дополнительных мыслей:

  • Линии линейные «достаточно», чтобы простой расчет расстояния между точками работал, линейные ссылки не должны быть необходимыми.
  • Это было бы легко в arcpy, если бы был инструмент для разделения линии в точке, но я не могу найти ее.

Мысли кто-нибудь?

пример

Майк Баннистер
источник
+1, интересная проблема! Мне интересно посмотреть, какие решения доступны =)
Джозеф
Трудно достичь только вашей средней линии - верх и низ появляются из клипа после заполнения любых пустот. Следовательно, я думаю, что вам следует сосредоточиться на этом вопросе и сузить область его действия до ArcPy, если это ваш предпочтительный инструмент. Вы всегда можете спросить о другом инструменте, если он не дает решения.
PolyGeo
линии пересекают многоугольники?
Эмиль Брундейдж
Эмиль, давайте предположим, что линии могут пересекать несколько полигонов. Однако, кроме геометрии, между полигонами нет никакой разницы, поэтому они могут быть растворены, объединены в составной элемент и т. Д., Если это облегчает алгоритм. Линия, пересекающая несколько многоугольников, вероятно, будет редкой, и это может быть отмеченным случаем, который должен рассматриваться вручную при необходимости.
Майк Баннистер
Какой у вас уровень лицензии?
Эмиль Брюндейдж

Ответы:

4

Я хочу добавить свое решение pyQGIS, больше ничего.

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)

Мой тестовый пример - до отсечения: перед клипом

После отсечения:

после

Чтобы получить полный набор атрибутов исходных линий, я думаю, что было бы лучше объединить их с результатом. В противном случае они должны быть созданы в разделе подготовки и установлены в самом внутреннем цикле. Но я не проверял, проходят ли они процесс растворения или теряются, потому что в принципе они могут иметь разные значения.

Детлеф
источник
Очень лаконичный ответ. Как снимки экрана QGIS всегда выглядят как QGIS?
Майк Баннистер
3

Это было бы легко в arcpy, если бы был инструмент для разделения линии в точке, но я не могу найти ее.

Если вы запустите Интегрировать с полигонами и линиями в качестве входных данных, он добавит вершину к каждой точке пересечения. (Осторожно, поскольку Integrate изменяет входы вместо создания новых выходов.)

Убедившись, что есть совпадающие вершины, вы можете перебрать вершины линии и проверить, не касается ли каждая из них другой функции. Из упорядоченного списка вершин, которые касаются, возьмите минимум и максимум из набора. Затем сделайте две строки для каждого объекта: A: (начало, ..., мин) и B: (максимум, ..., конец).

Другой вариант, хотя я не уверен, что ArcPy сохранит упорядочение частей объекта на основе упорядочения вершин во входном объекте, это запустить клип как есть. Для средней линии в вашем примере это должно привести к составной функции из трех частей. В зависимости от порядка вы можете выполнить итерацию для каждой составной строки, созданной Clip, и удалить все элементы, кроме первой и последней, из выходной многокомпонентной функции.

Джон Рейзер
источник
3

В этом случае нужно решить три вопроса:

  • Отверстия
  • Линии между полигонами
  • Конечные строки

введите описание изображения здесь

Отверстия

Поскольку любая линия в отверстии будет сохранена, удалите отверстия из полигонов. В приведенном ниже сценарии я делаю это с помощью курсоров и геометрий.

Линии между полигонами

Линии, которые касаются двух полигонов, должны быть удалены. В приведенном ниже сценарии я делаю это, выполняя пространственное соединение one to manyс моими линиями в качестве входного класса объектов и полигонами в качестве класса объектов соединения. Любая линия, которая генерируется дважды, касается двух многоугольников и удаляется.

Конечные строки

Чтобы удалить линии, которые касаются только многоугольника на одном конце, я конвертирую линии в конечные точки. Затем я использую слои объектов и выборки, чтобы определить, какие конечные точки являются плавающими. Я выбираю конечные точки, которые пересекают полигоны. Затем я переключаю свой выбор. Это выбирает конечные точки, которые не пересекаются с полигонами. Я выбираю любую линию, которая пересекает эти выбранные точки и удаляю их.

Результат

введите описание изображения здесь

Предположения

  • Входные данные - это классы объектов файловой базы геоданных
  • Доступна расширенная лицензия ArcGIS (из-за eraseи feature vertices to points)
  • Непрерывные, соединенные линии - это одна особенность
  • Полигоны не перекрываются
  • Нет многочастных многоугольников

скрипт

Сценарий ниже выводит класс пространственных объектов с именем вашего класса линейных пространственных объектов плюс _GreedyClip, в той же базе геоданных, что и ваш линейный класс пространственных объектов. Расположение рабочей области также необходимо.

#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)  
Эмиль Брундейдж
источник
Отличное решение, Эмиль. Это меньше кода, чем я закончил.
Майк Баннистер