Linhas de recorte “gananciosas” com polígono

9

Desejo cortar um conjunto de polilinhas (linhas pretas na imagem abaixo) no limite externo de um polígono. Quaisquer vazios dentro do polígono devem ser ignorados. Minha saída ideal são as linhas amarelas tracejadas. As linhas iniciais podem ou não ser retas. A imagem é um exemplo simplificado, na realidade o polígono é muito mais complexo e existem centenas de linhas. Não acho que um casco convexo funcionaria (mas posso estar errado). Estou aberto a soluções em arcgis, qgis, arcpy, bem torneado, etc. A codificação preferencialmente seria em python, pois estou aberto a outras opções, se necessário. A Arcgis também seria preferível para facilitar a partilha dos meus colegas de trabalho, mas isso não é um requisito.

O melhor que consigo pensar agora é cruzar uma linha individual com o polígono, criando um conjunto de pontos em todas as interseções de fronteira. Classifique os pontos por distância até o início da linha. Os pontos mais distantes e próximos (FAC) serão o limite externo do polígono. Em seguida, use os pontos FAC para selecionar os vértices adequados da linha original e crie a linha tracejada amarela a partir dos pontos apropriados. Deve funcionar, mas parece mais complicado do que o necessário.

Alguns pensamentos adicionais:

  • As linhas são lineares "o suficiente" para que um cálculo simples da distância entre pontos funcione; a referência linear não deve ser necessária.
  • Isso seria fácil no arcpy se houvesse uma ferramenta para dividir uma linha em um ponto, mas não consigo encontrar uma.

Pensamentos alguém?

Exemplo

Mike Bannister
fonte
+1, problema interessante! Estou ansioso para ver o que soluções são = disponíveis)
Joseph
Somente sua linha do meio é difícil de alcançar - as partes superior e inferior vêm de um clipe depois de preencher os vazios. Consequentemente, acho que você deve focar sua pergunta nisso e restringir seu escopo apenas ao ArcPy, se essa for sua ferramenta preferida. Você sempre pode perguntar sobre outra ferramenta, se isso não produzir uma solução.
PolyGeo
as linhas cruzam vários polígonos?
Emil Brundage
Emil, vamos supor que as linhas possam cruzar vários polígonos. No entanto, além da geometria, não há diferença entre os polígonos, para que possam ser dissolvidos, mesclados em um recurso de várias partes, etc., se isso facilitar o algoritmo. Uma linha cruzada sobre vários polígonos provavelmente seria rara e pode ser um caso sinalizado para ser tratado manualmente, se necessário.
Mike Bannister
Qual é o seu nível de licença?
Emil Brundage

Respostas:

4

Eu quero jogar na minha solução pyQGIS, nada mais.

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)

Meu caso de teste - antes do recorte: antes do clipe

Após o recorte:

depois de

Para obter o conjunto completo de atributos das linhas originais, acho que seria melhor juntá-las ao resultado. Caso contrário, eles devem ser criados na seção de preparação e configurados no loop mais interno. Mas não testei se eles passam no processo de dissolução ou se perdem, porque em princípio eles podem ter valores diferentes.

Detlev
fonte
Resposta muito concisa. Como as capturas de tela do QGIS sempre se parecem com o QGIS?
Mike Bannister
3

Isso seria fácil no arcpy se houvesse uma ferramenta para dividir uma linha em um ponto, mas não consigo encontrar uma.

Se você executar o Integrate com os polígonos e linhas como entradas, ele adicionará um vértice a cada onde eles cruzam. Cuidado, pois o Integrate modifica entradas em vez de produzir novas saídas.

Quando tiver certeza de que existem vértices coincidentes, você pode percorrer os vértices da linha e testar para ver se cada um deles toca o outro recurso. Na lista ordenada de vértices que tocam, tire o mínimo e o máximo do conjunto. Em seguida, faça duas linhas de cada recurso, A: (início, ..., min) e B: (máximo, ..., final).

Outra opção, embora não tenha certeza se o ArcPy preserva a ordem das peças do recurso com base na ordem dos vértices no objeto de entrada, seria executar o clipe como está. Para a linha do meio no seu exemplo, deve resultar em um recurso de várias partes com três partes. Dependendo da ordem, você pode iterar todas as linhas de várias partes produzidas pelo Clip e remover tudo, exceto a primeira e a última parte do recurso de saída múltipla.

John Reiser
fonte
3

Há três questões a serem enfrentadas neste caso:

  • Furos
  • Linhas entre polígonos
  • Linhas finais

insira a descrição da imagem aqui

Furos

Como qualquer linha dentro de um furo será mantida, remova os orifícios dos polígonos. No script abaixo, faço isso usando cursores e geometrias.

Linhas entre polígonos

As linhas que tocam em dois polígonos precisam ser removidas. No script abaixo, faço isso executando uma junção espacial de one to many, com minhas linhas como minha classe de recurso de entrada e meus polígonos como minha classe de recurso de junção. Qualquer linha que é gerada duas vezes toca em dois polígonos e é removida.

Linhas finais

Para remover linhas que tocam apenas um polígono em uma extremidade, eu converto linhas em pontos finais. Utilizo então as camadas e seleções de recursos para determinar quais pontos finais são flutuantes. Eu seleciono os pontos finais que cruzam os polígonos. Eu então mudo minha seleção. Isso seleciona pontos finais que não cruzam polígonos. Seleciono qualquer linha que cruze esses pontos selecionados e os apago.

Resultado

insira a descrição da imagem aqui

Premissas

  • As entradas são classes de recurso de geodatabase de arquivo
  • A licença avançada do ArcGIS está disponível (devido a um erasee a feature vertices to points)
  • Linhas conectadas contínuas são um recurso único
  • Os polígonos não se sobrepõem
  • Não há polígonos multipartes

Roteiro

O script abaixo gera uma classe de recurso com o nome da sua classe de recurso de linha plus _GreedyClip, no mesmo banco de dados geográfico da sua classe de recurso de linha. Um local de espaço de trabalho também é necessário.

#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)  
Emil Brundage
fonte
Solução agradável Emil. Isso é menos código do que eu acabei.
Mike Bannister