Estadísticas zonales para millones de polígonos superpuestos

Question

Yo soy de ArcGIS 10.1 con 64 bits GP en Windows 8.1 x64.

Quiero calcular el zonal estadísticas de búfer basado en polígonos parcelas y un kernel de densidad de la trama. Los 4,7 millones de polígonos se solapan y se almacenan en una única clase de entidad, haciendo que el problema particularmente difícil. He intentado utilizar el suplementario de Spatial Analyst Tools desarrollado por ESRI en diferentes formas, pero el guión no es lo suficientemente rápida. He tratado de calcular el zonal estadísticas pieza por pieza, el uso de pequeñas áreas censales (con varios cientos de búfer de polígonos de cada uno) y los grandes (miles de búfer de polígonos). Ninguna de estas implementaciones realizar adecuadamente. He incluido las siguientes características en una secuencia de comandos de Python:

• la carga de buffer de polígonos en el in_memory área de trabajo antes de cargar en la herramienta ZS

• la creación de las capas de entidades donde sea necesario

• escribir los resultados en un único archivo de texto

El código para la realización de ZS un polígono en un momento se publica a continuación, sin embargo, la herramienta sólo se completa acerca de una parcela a cada minuto. Me pregunto si alguien tiene alguna idea sobre cómo mejorar zonal estadísticas de rendimiento o la secuencia de comandos en general. Gracias!

import arcpy, os, string, csv
from arcpy import env

# arcpy.ImportToolbox("D:\ClusterAnalysis_Deployment\SpatialAnalystSupplementalTools\Spatial Analyst Supplemental Tools.pyt")
arcpy.CheckOutExtension("Spatial")

# LOCATIONS #
# Local machine paths
gdb = r"D:\path\to\geodatabase.gdb"
results = r"\path\to\results"

# Virtual machine paths

# INPUT DATA #
b = gdb + "\\" + "ParcelBuffers"
kd = gdb + "\\" + "kdensity750"

env.workspace = gdb
env.overwriteOutput = True

# TEMPORARY DATA #
buffers = arcpy.MakeFeatureLayer_management(b, "bdata")
kdensity = arcpy.Raster(kd)

parcelList = []

cursor = arcpy.da.SearchCursor(buffers, ("PIN"))
for row in cursor:
    parcelList.append(row[0])
del cursor

countDict = {}
countDict["Count:"] = 0

print "Script setup... done"

# GEOPROCESSING #
for PIN in parcelList:
    parcel_ram = "in_memory" + "\\" + "parcel"
    zs_table = results + "\\" + "zs_table_" + str(PIN) + ".dbf"
    solution = results + "\\" + "ZS_solutions.txt"

    arcpy.SelectLayerByAttribute_management(buffers, "NEW_SELECTION", "\"PIN\" = \'" + str(PIN) + "\'")
    count = int(arcpy.GetCount_management(buffers).getOutput(0))
    arcpy.CopyFeatures_management(buffers, parcel_ram)
    arcpy.SelectLayerByAttribute_management(buffers, "CLEAR_SELECTION")
    arcpy.gp.ZonalStatisticsAsTable_sa(parcel_ram, "PIN", kdensity, zs_table, "DATA", "ALL")

    table = arcpy.gp.MakeTableView_management(zs_table, "zs_table_lyr")

    fields = arcpy.ListFields(table)
    field_names = map(lambda n:n.name, fields)
    header = string.join(field_names, "\t")

    with open(solution, 'w') as x:
        x.write(header + "\n")
        cursor = arcpy.SearchCursor(table)
        for row in cursor:
            values = string.join(map(lambda n: str(row.getValue(n)), field_names), "\t")
            x.write(values + "\n")
        del row, cursor

    countDict["Count:"] = countDict["Count:"] + count
    print "Zonal statistics computed and saved for parcel: " + str(countDict["Count:"])

    arcpy.Delete_management(parcel_ram)
    arcpy.Delete_management(zs_table)
    arcpy.Delete_management(table)
    del count

Después de algunas pruebas adicionales de Felix del script creo que la forma en que mi conjuntos de datos se configuran está mal, que conduce a la unión espacial característica de consulta de tirar de los polígonos que se sobreponen. Aquí están los conjuntos de datos y sus atributos relevantes:

• fishnet característica (ID d por entero largo "FnID")
• minBound característica (ID d por entero largo "FnID")
• parentLR (Unirse FID es FnID y de Destino de la FID es parID, donde parID es entero largo campo creado para el búfer de polígonos)

La secuencia de comandos:

##  Zonal statistics on a set of overlapping
import arcpy, traceback, os, sys, time
from arcpy import env

arcpy.CheckOutExtension("Spatial")

env.workspace = "C:\Geoprocessing\ClusterAnalysis"
env.overwriteOutput = True

# Field name of fishnet ID
parID="FnID"
parID2="FnID_1"

# Location of input map document containing the minBound (minimum bounding geometry)
# layer and parentLR (spatial join of fishnet to polygons) layer

mapDocument = env.workspace + "\\" + "ClusterAnalysis.mxd"
raster = env.workspace + "\\" + "ClusterAnalysis.gdb\intKD750"
rast = env.workspace + "\\" + "ClusterAnalysis.gdb\parcelRaster"
kdensity = arcpy.Raster(raster)

dbf="stat.dbf" # output zonal statistics table
joinLR="SD.shp" # intermediate data set for determining non-overlapping polygons
subset="subset"

##try:
def showPyMessage():
    arcpy.AddMessage(str(time.ctime()) + " - " + message)
def Get_V(aKey):
    try:
        return smallDict[aKey]
    except:
        return (-1)
def pgonsCount(aLayer):
        result=arcpy.GetCount_management(aLayer)
        return int(result.getOutput(0))

# Define variables based on layers in ArcMap document
mxd = arcpy.mapping.MapDocument(mapDocument)
layers = arcpy.mapping.ListLayers(mxd)
minBound,parentLR=layers[0],layers[1]

# Query parentLR polygons based on fishnet grid of minBound polygon
with arcpy.da.SearchCursor(minBound, ("SHAPE@","FnID")) as clipper:
    for rcrd in clipper:
        feat=rcrd[0]
        env.extent=feat.extent
        fp='"FnID"='+str(rcrd[1])
        parentLR.definitionQuery=fp
        nF=pgonsCount(parentLR)
        arcpy.AddMessage("Processing subset %i containing %i polygons" %(rcrd[1],nF))
        arcpy.AddMessage("Defining neighbours...")
        arcpy.SpatialJoin_analysis(parentLR, parentLR, joinLR, "JOIN_ONE_TO_MANY")
        arcpy.AddMessage("Creating empty dictionary")

        dictFeatures = {}
        with arcpy.da.SearchCursor(parentLR, parID) as cursor:
            for row in cursor:
                dictFeatures[row[0]]=()
            del row, cursor

        arcpy.AddMessage("Assigning neighbours...")
        nF=pgonsCount(joinLR)
        arcpy.SetProgressor("step", "", 0, nF,1)
        with arcpy.da.SearchCursor(joinLR, (parID,parID2)) as cursor:
            for row in cursor:
                aKey=row[0]
                aList=dictFeatures[aKey]
                aList+=(row[1],)
                dictFeatures[aKey]=aList
                arcpy.SetProgressorPosition()
            del row, cursor
        arcpy.AddMessage("Defining non-overlapping subsets...")
        runNo=0
        while (True):
            parentLR.definitionQuery=fp
            toShow,toHide=(),()
            nF=len(dictFeatures)
            arcpy.SetProgressor("step", "", 0, nF,1)
            for item in dictFeatures:
                if item not in toShow and item not in toHide:
                    toShow+=(item,)
                    toHide+=(dictFeatures[item])
                arcpy.SetProgressorPosition()
            m=len(toShow)
            quer='"parID" IN '+str(toShow)+ " AND "+fp
            if m==1:
                quer='"parID" = '+str(toShow[0])+ " AND "+fp
            parentLR.definitionQuery=quer
            runNo+=1
            arcpy.AddMessage("Run %i, %i polygon(s) found" % (runNo,m))
            arcpy.AddMessage("Running Statistics...")
            arcpy.gp.ZonalStatisticsAsTable_sa(parentLR, "PIN", kdensity, dbf,"DATA", "MEAN")
            arcpy.AddMessage("Data transfer...")
            smallDict={}
            with arcpy.da.SearchCursor(dbf, (parID,"MEAN")) as cursor:
                for row in cursor:
                    smallDict[row[0]]=row[1]
                del row, cursor
            with arcpy.da.UpdateCursor(parentLR, (parID,"MEAN")) as cursor:
                for row in cursor:
                    aKey=row[0]
                    row[1]=Get_V(aKey)
                    cursor.updateRow(row)
                del row, cursor
            for item in toShow:
                del dictFeatures[item]
            m=len(dictFeatures)
            if m==0:
                break
        parentLR.definitionQuery=fp
        del smallDict, dictFeatures
parentLR.definitionQuery=''

Produce el zonal estadísticas de error: http://resources.arcgis.com/en/help/main/10.1/index.html#//00vq0000000s010423

Answer 1

1. Convertir raster a puntos
2. Intersectar puntos con polígonos, tipo de salida POINT
3. Estadísticas sobre el valor de ráster, utilizando polígono como campo de caso. Usted podría considerar para dividir el área, si raster demasiado grande

Answer 2

Pre-procesamiento:

1. Crear red - polígonos con medida=grado de superposición de polígonos. Calcular campo ID = FID.
2. Espacial unirse a la superposición de polígonos con fishnet polígonos utilizando HAVE_THEIR_CENTER_IN.
3. Disolver polígonos, caso de campo de IDENTIFICACIÓN de red
4. Uso Mínimo de Delimitación de la geometría en el paso 3 de salida.

Negro polígonos en la foto de abajo son las salidas de la etapa 4 (minBound en la secuencia de comandos). Color polígonos (unique ID de red) son las salidas de la etapa 2 (parentLR en la secuencia de comandos)

La secuencia de comandos se supone:

• minBound es 1ª capa en el TOC
• parentLR es 2ª capa
• Identificador único de parentLR es parID - entero largo
• Identificador único de minBound es ID - número entero largo

Cambiar el nombre de la trama de la muestra.

minBound se utiliza para navegar por área de interés y el estrecho de procesamiento de medida y el número de polígonos en proceso en un momento, por modifuing parentLR definición de la consulta. Parte fundamental es la selección de la no superposición de grupos de polígonos, llevada a cabo por el análisis espacial de unirse a la mesa con ONE_TO_MANY opción y modificar la consulta de definición de parentLR. Un rendimiento de más de 2000 polígonos por minuto, por lo menos de 2 días para el proceso de 4.7 millin polígonos:). Windows 7, memoria RAM de 8 GB, 64 bits.

##  Zonal statistics on a set of overlapping 
import arcpy, traceback, os, sys
from arcpy import env
env.overwriteOutput = True
parID="parID"
parID2="parID_1"
dem=r'C:\From_MXD\ARC2MDEM_Clip'
env.workspace = "C:\\From_MXD"
dbf="stat.dbf"
joinLR="SD.shp"
subset="subset"

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)
    def Get_V(aKey):
        try:
            return smallDict[aKey]
        except:
            return (-1)
    def pgonsCount(aLayer):
            result=arcpy.GetCount_management(aLayer)
            return int(result.getOutput(0))    
    mxd = arcpy.mapping.MapDocument("CURRENT")
    layers = arcpy.mapping.ListLayers(mxd)
    minBound,parentLR=layers[0],layers[1]
    with arcpy.da.SearchCursor(minBound, ("SHAPE@","ID")) as clipper:
        for rcrd in clipper:
            feat=rcrd[0]
            env.extent=feat.extent
            fp='"ID"='+str(rcrd[1])
            parentLR.definitionQuery=fp
            nF=pgonsCount(parentLR)  
            arcpy.AddMessage("Processing subset %i containing %i polygons" %(rcrd[1],nF))
            arcpy.AddMessage("Defining neighbours...")
            arcpy.SpatialJoin_analysis(parentLR, parentLR, joinLR,"JOIN_ONE_TO_MANY")
            arcpy.AddMessage("Creating empty dictionary")
            dictFeatures = {}
            with arcpy.da.SearchCursor(parentLR, parID) as cursor:
                for row in cursor:
                    dictFeatures[row[0]]=()
                del row, cursor
            arcpy.AddMessage("Assigning neighbours...")
            nF=pgonsCount(joinLR)
            arcpy.SetProgressor("step", "", 0, nF,1)
            with arcpy.da.SearchCursor(joinLR, (parID,parID2)) as cursor:
                for row in cursor:
                    aKey=row[0]
                    aList=dictFeatures[aKey]
                    aList+=(row[1],)
                    dictFeatures[aKey]=aList
                    arcpy.SetProgressorPosition()
                del row, cursor    
            arcpy.AddMessage("Defining non-overlapping subsets...")
            runNo=0
            while (True):
                parentLR.definitionQuery=fp
                toShow,toHide=(),()
                nF=len(dictFeatures)    
                arcpy.SetProgressor("step", "", 0, nF,1)
                for item in dictFeatures:
                    if item not in toShow and item not in toHide:
                        toShow+=(item,)
                        toHide+=(dictFeatures[item])
                    arcpy.SetProgressorPosition()
                m=len(toShow)
                quer='"parID" IN '+str(toShow)+ " AND "+fp
                if m==1:
                    quer='"parID" = '+str(toShow[0])+ " AND "+fp
                parentLR.definitionQuery=quer
                runNo+=1
                arcpy.AddMessage("Run %i, %i polygon(s) found" % (runNo,m))
                arcpy.AddMessage("Running Statistics...")
                arcpy.gp.ZonalStatisticsAsTable_sa(parentLR, parID, dem, dbf, "DATA", "MEAN")
                arcpy.AddMessage("Data transfer...")
                smallDict={}
                with arcpy.da.SearchCursor(dbf, (parID,"MEAN")) as cursor:
                    for row in cursor:
                        smallDict[row[0]]=row[1]
                    del row, cursor    
                with arcpy.da.UpdateCursor(parentLR, (parID,"MEAN")) as cursor:
                    for row in cursor:
                        aKey=row[0]
                        row[1]=Get_V(aKey)
                        cursor.updateRow(row)
                    del row, cursor
                for item in toShow:
                    del dictFeatures[item]
                m=len(dictFeatures)
                if m==0:
                    break
            parentLR.definitionQuery=fp
            del smallDict, dictFeatures
    parentLR.definitionQuery=''
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()            

Usted puede aumentar la velocidad de eliminación de múltiples líneas de progreso de la secuencia de comandos. Sugiero probarlo en un subconjunto más pequeño de los polígonos, por ejemplo, 100 miles. Yo también creo que el tamaño de la red tiene que estar optimizado. Tenga en cuenta que yo soy derivados justo valor de la media de la trama. He tratado de documentar secuencia de comandos para mi mejor