¿Existe una forma más rápida de unir espacialmente dos clases de entidades?

Question

Tengo dos grandes clases de entidades (3.5 millones de puntos y 3.8 millones de polígonos) que me gustaría unir espacialmente.

¿Hay una forma más rápida de hacer esto en lugar de simplemente usar la herramienta clásica de unión espacial?

Supongo que esta es una pregunta repetida, simplemente no la puedo encontrar.

Cuando uno clases de entidades de este tamaño, generalmente solo dejo que se ejecute durante unas horas para completarse.

En el futuro, me gustaría encontrar una forma de convertir este tipo de procesos de horas a minutos.

Solo traigo un campo de identificación único para crear una clave externa y luego relaciono las dos tablas.

Solo buscando consejos o ideas.

Unión de uno a uno Intersección

Amigable con Arcpy ArcGIS 10.1 Procesador: Intel(R) Xeon(R) CPU E5-1607 0 @ 3.00 GHz RAM: 16 GB Tipo de sistema: Sistema operativo de 64 bits

Answer 1

Posible solución utilizando rasters:

• Divida la extensión total utilizando fishnet
• Establezca el tamaño de celda de entorno
• a) comience con el primer mosaico, es decir, n = 1
• b) Convertir polígonos y puntos a raster
• Calcular estadísticas zonales (cuadrícula de puntos = zonas, cuadrícula de polígonos = valores)
• Contar el número de registros en la tabla de estadísticas. Si es menor que el recuento de puntos, disminuya el tamaño de la celda y vaya a b), de lo contrario transfiera los datos de la tabla de estadísticas a la tabla de puntos, continúe con el siguiente mosaico, es decir, n+=1 ir a a)

Probé el script a continuación utilizando 2.4 millones de puntos y la misma cantidad de polígonos. Utilicé una malla igual a 5*5 mosaicos. Esto significa que se procesaron aproximadamente 100 mil puntos a la vez. Esta es una captura de la ventana de resultados:

Probé la unión espacial en los mismos conjuntos de datos y cuando llegó al 5%, lo cancelé. Para este momento, la herramienta estaba en ejecución durante 8 minutos. Esto significa que con conjuntos de datos verdaderamente masivos (algunos millones), el enfoque de raster podría reducir considerablemente el tiempo de procesamiento.

Tenga en cuenta que utilicé raster para crear mis puntos y polígonos, es decir, los polígonos probados eran los cuadrados más básicos. Espero que la ventaja del enfoque de raster pueda ser aún mayor con formas más complejas. Podría reducir aún más el tiempo de procesamiento, utilizando una malla con menos mosaicos, pero no quería arriesgarme a la corrupción del espacio in_memory.

SCRIPT:

import arcpy, time, os
from arcpy import env
from arcpy.sa import ZonalStatisticsAsTable as ZS
env.overwriteoutput=True
t0=time.time()
## PARÁMETROS
workGDB=r'..\Scratch.gdb'
fish_net=r'D:\Scratch\fish_net.shp'
points=r'..\Scratch.gdb\points'
polygons=r'..\Scratch.gdb\polygons'
pointsGrid=r'in_memory\pntsGR'
polygonsGrid=r'in_memory\pgonsGR'
stats=r'D:\Scratch\stats.dbf'
##allStats=workGDB+os.sep+'allStats'

G=arcpy.Geometry()
tiles=arcpy.CopyFeatures_management(fish_net, G)
nTot=0
cSize=4
aDict={}
for n,tile in enumerate(tiles):
    arcpy.AddMessage("Procesando mosaico n. %i de %i" %(n+1,len(tiles)))
    env.extent=tile.extent
    result=arcpy.GetCount_management(points)
    nPoints=int(float(result.getOutput(0)))
    while True:
        env.cellSize = cSize
        arcpy.arcpy.PointToRaster_conversion(points, "POINTID",pointsGrid)
        arcpy.FeatureToRaster_conversion(polygons, "ID",polygonsGrid)
        outZSaT = ZS(pointsGrid, "Value",polygonsGrid,stats,"DATA", "MINIMUM")
        result=arcpy.GetCount_management(stats)
        nStats=int(float(result.getOutput(0)))
        if nStats>=nPoints:break
        cSize=float(cSize)/2
    arcpy.AddMessage("Actualizando búsqueda ....")
    arcpy.SetProgressor("step", "", 0, nPoints)
    with arcpy.da.SearchCursor(stats,("Value","MIN")) as cursor:
        for v,m in cursor:
            aDict[v]=m
            arcpy.SetProgressorPosition()
    nTot+=nPoints
    arcpy.AddMessage("Procesados %i puntos" %nTot)

arcpy.Delete_management("in_memory")
arcpy.AddMessage("Transfiriendo resultados ...")
arcpy.SetProgressor("step", "", 0, nTot)
with arcpy.da.UpdateCursor(points,("pointid","PGONID")) as cursor:
    for row in cursor:
        row[1]=aDict[row[0]]
        cursor.updateRow(row)
        arcpy.SetProgressorPosition()
arcpy.AddMessage("Tiempo total %i seg" %int(time.time()-t0))