Exploser chevauchant vers de nouveaux polygones non chevauchants?

Étant donné plusieurs polygones qui se chevauchent de plusieurs manières, je voudrais exporter de ces entités tous les polygones qui ne se chevauchent pas de manière itérative.

Le produit serait un certain nombre de fonctionnalités sans chevauchement qui, réunies, constituent l'original.

Les produits pourraient ensuite être utilisés comme entrées dans les statistiques zonales, ce qui serait beaucoup plus rapide que l'itération des statistiques zonales sur chaque polygone.

J'ai essayé de coder cela dans ArcPy sans succès.

Le code pour cela existe-t-il déjà?

Il s'agit d'un problème de coloration de graphique .

Rappelez-vous qu'un coloriage de graphe est une affectation d'une couleur aux sommets d'un graphe de telle sorte qu'aucun deux sommets qui partagent un bord n'auront également la même couleur. Plus précisément, les sommets (abstraits) du graphique sont les polygones. Deux sommets sont connectés avec une arête (non orientée) chaque fois qu'ils se croisent (sous forme de polygones). Si nous prenons une solution au problème - qui est une séquence de (disons k ) collections disjointes des polygones - et attribuons une couleur unique à chaque collection de la séquence, alors nous aurons obtenu un k- coloriage du graphique . Il est souhaitable de trouver un petit k .

Ce problème est assez difficile et reste non résolu pour les graphiques arbitraires. Considérez une solution approximative simple à coder. Un algorithme séquentiel devrait faire l'affaire. L'algorithme Welsh-Powell est une solution gourmande basée sur un ordre décroissant des sommets par degré. Traduit dans la langue des polygones d'origine, triez d'abord les polygones dans l'ordre décroissant du nombre d'autres polygones qu'ils chevauchent. En travaillant dans l'ordre, donnez au premier polygone une couleur initiale. À chaque étape successive, essayez de colorer le polygone suivant avec une couleur existante: c'est-à-dire, choisissez une couleur qui n'est pasdéjà utilisé par l'un des voisins de ce polygone. (Il existe plusieurs façons de choisir parmi les couleurs disponibles; essayez soit celle qui a été la moins utilisée, soit choisissez-en une au hasard.) Si le polygone suivant ne peut pas être coloré avec une couleur existante, créez une nouvelle couleur et coloriez-la avec celle-ci.

Une fois que vous avez réalisé une coloration avec un petit nombre de couleurs, effectuez zonalstats couleur par couleur: par construction, vous êtes assuré qu'il n'y a pas deux polygones d'une couleur donnée qui se chevauchent.

Voici un exemple de code dans R. (Le code Python ne serait pas très différent.) Premièrement, nous décrivons les chevauchements entre les sept polygones montrés.

edges <- matrix(c(1,2, 2,3, 3,4, 4,5, 5,1, 2,6, 4,6, 4,7, 5,7, 1,7), ncol=2, byrow=TRUE)

Autrement dit, les polygones 1 et 2 se chevauchent, tout comme les polygones 2 et 3, 3 et 4, ..., 1 et 7.

Triez les sommets par degré décroissant:

vertices <- unique(as.vector(edges))
neighbors <- function(i) union(edges[edges[, 1]==i,2], edges[edges[, 2]==i,1])
nbrhoods <- sapply(vertices, neighbors)
degrees <- sapply(nbrhoods, length)
v <- vertices[rev(order(degrees))]

Un algorithme de coloration séquentielle (brut) utilise la première couleur disponible qui n'est pas déjà utilisée par un polygone se chevauchant:

color <- function(i) {
  n <- neighbors(i)
  candidate <- min(setdiff(1:color.next, colors[n]))
  if (candidate==color.next) color.next <<- color.next+1
  colors[i] <<- candidate
}

Initialisez les structures de données ( colorset color.next) et appliquez l'algorithme:

colors <- rep(0, length(vertices))
color.next <- 1
temp <- sapply(v, color)

Divisez les polygones en groupes selon la couleur:

split(vertices, colors)

La sortie de cet exemple utilise quatre couleurs:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Il a divisé les polygones en quatre groupes qui ne se chevauchent pas. Dans ce cas, la solution n'est pas optimale ({{3,6,5}, {2,4}, {1,7}} est un tricolore pour ce graphique). En général, la solution qu'il obtient ne devrait cependant pas être trop mauvaise.

La méthodologie recommandée par #whuber m'a inspiré à prendre une nouvelle direction, et voici ma solution arcpy, en deux fonctions. Le premier, appelé countOverlaps, crée deux champs, "overlaps" et "ovlpCount" pour enregistrer pour chaque poly les polys qui se chevauchent avec lui, et combien de chevauchements se sont produits. La deuxième fonction, explodeOverlaps, crée un troisième champ, "expl", qui donne un entier unique à chaque groupe de polys non chevauchants. L'utilisateur peut alors exporter de nouveaux fc en fonction de ce champ. Le processus est divisé en deux fonctions car je pense que l'outil countOverlaps peut s'avérer utile en soi. Veuillez excuser la négligence du code (et la convention de dénomination imprudente), car c'est assez préliminaire, mais cela fonctionne. Assurez-vous également que le "idName" champ représente un champ avec des ID uniques (uniquement testé avec des ID entiers). Merci whuber de m'avoir fourni le cadre nécessaire pour aborder ce problème!

def countOverlaps(fc,idName):
    intersect = arcpy.Intersect_analysis(fc,'intersect')
    findID = arcpy.FindIdentical_management(intersect,"explFindID","Shape")
    arcpy.MakeFeatureLayer_management(intersect,"intlyr")
    arcpy.AddJoin_management("intlyr",arcpy.Describe("intlyr").OIDfieldName,findID,"IN_FID","KEEP_ALL")
    segIDs = {}
    featseqName = "explFindID.FEAT_SEQ"
    idNewName = "intersect."+idName

    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal] = []
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal].append(idVal)

    segIDs2 = {}
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        segIDs2[idVal] = []

    for x,y in segIDs.iteritems():
        for segID in y:
            segIDs2[segID].extend([k for k in y if k != segID])

    for x,y in segIDs2.iteritems():
        segIDs2[x] = list(set(y))

    arcpy.RemoveJoin_management("intlyr",arcpy.Describe(findID).name)

    if 'overlaps' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'overlaps',"TEXT")
    if 'ovlpCount' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'ovlpCount',"SHORT")

    urows = arcpy.UpdateCursor(fc)
    for urow in urows:
        idVal = urow.getValue(idName)
        if segIDs2.get(idVal):
            urow.overlaps = str(segIDs2[idVal]).strip('[]')
            urow.ovlpCount = len(segIDs2[idVal])
        urows.updateRow(urow)

def explodeOverlaps(fc,idName):

    countOverlaps(fc,idName)

    arcpy.AddField_management(fc,'expl',"SHORT")

    urows = arcpy.UpdateCursor(fc,'"overlaps" IS NULL')
    for urow in urows:
        urow.expl = 1
        urows.updateRow(urow)

    i=1
    lyr = arcpy.MakeFeatureLayer_management(fc)
    while int(arcpy.GetCount_management(arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",'"expl" IS NULL')).getOutput(0)) > 0:
        ovList=[]
        urows = arcpy.UpdateCursor(fc,'"expl" IS NULL','','','ovlpCount D')
        for urow in urows:
            ovVal = urow.overlaps
            idVal = urow.getValue(idName)
            intList = ovVal.replace(' ','').split(',')
            for x in intList:
                intList[intList.index(x)] = int(x)
            if idVal not in ovList:
                urow.expl = i
            urows.updateRow(urow)
            ovList.extend(intList)
        i+=1

Cela fait un moment, mais j'ai utilisé ce code pour ma propre application et cela fonctionne très bien - merci. J'ai réécrit une partie de celui-ci pour le mettre à jour, l'appliquer aux lignes (avec une tolérance) et l'accélérer considérablement (ci-dessous - je l'exécute sur 50 millions de tampons qui se croisent et cela ne prend que quelques heures).

def ExplodeOverlappingLines(fc, tolerance, keep=True):
        print('Buffering lines...')
        idName = "ORIG_FID"
        fcbuf = arcpy.Buffer_analysis(fc, fc+'buf', tolerance, line_side='FULL', line_end_type='FLAT')
        print('Intersecting buffers...')
        intersect = arcpy.Intersect_analysis(fcbuf,'intersect')

        print('Creating dictionary of overlaps...')
        #Find identical shapes and put them together in a dictionary, unique shapes will only have one value
        segIDs = defaultdict(list)
        with arcpy.da.SearchCursor(intersect, ['Shape@WKT', idName]) as cursor:
            x=0
            for row in cursor:
                if x%100000 == 0:
                    print('Processed {} records for duplicate shapes...'.format(x))
                segIDs[row[0]].append(row[1])
                x+=1

        #Build dictionary of all buffers overlapping each buffer
        segIDs2 = defaultdict(list)
        for v in segIDs.values():
            for segID in v:
                segIDs2[segID].extend([k for k in v if k != segID and k not in segIDs2[segID]])

        print('Assigning lines to non-overlapping sets...')
        grpdict = {}
        # Mark all non-overlapping one to group 1
        for row in arcpy.da.SearchCursor(fcbuf, [idName]):
            if row[0] in segIDs2:
                grpdict[row[0]] = None
            else:
                grpdict[row[0]] = 1

        segIDs2sort = sorted(segIDs2.items(), key=lambda x: (len(x[1]), x[0])) #Sort dictionary by number of overlapping features then by keys
        i = 2
        while None in grpdict.values(): #As long as there remain features not assigned to a group
            print(i)
            ovset = set()  # list of all features overlapping features within current group
            s_update = ovset.update
            for rec in segIDs2sort:
                if grpdict[rec[0]] is None: #If feature has not been assigned a group
                    if rec[0] not in ovset: #If does not overlap with a feature in that group
                        grpdict[rec[0]] = i  # Assign current group to feature
                        s_update(rec[1])  # Add all overlapping feature to ovList
            i += 1 #Iterate to the next group

        print('Writing out results to "expl" field in...'.format(fc))
        arcpy.AddField_management(fc, 'expl', "SHORT")
        with arcpy.da.UpdateCursor(fc,
                                   [arcpy.Describe(fc).OIDfieldName, 'expl']) as cursor:
            for row in cursor:
                if row[0] in grpdict:
                    row[1] = grpdict[row[0]]
                    cursor.updateRow(row)

        if keep == False:
            print('Deleting intermediate outputs...')
            for fc in ['intersect', "explFindID"]:
                arcpy.Delete_management(fc)

Dans ce cas, j'utilise généralement la méthode suivante:

• Passez la classe Feature à travers UNION; (Il brise les polygones dans toutes ses intersections)
• Ajoutez les champs X, Y et Zone et calculez-les;
• Dissolvez le résultat par les champs X, Y, Zone.

Je crois que le résultat sera celui que vous vouliez, et vous pouvez même compter le nombre de chevauchements. Je ne sais pas si en termes de performances ce sera mieux pour vous ou non.