¿Tabular la intersección de ArcGIS for Desktop devuelve un resultado vacío?

Question

Estoy utilizando una capa de polígonos de tierras federales de EE.UU. ( http://coastalmap.marine.usgs.gov/GISdata/basemaps/boundaries/fedlands/fedlanp020.htm )

He descargado la capa, la he reproyectado como Albers Equal Area y he creado un búfer de 25 km alrededor de los polígonos. He intentado utilizar Tabulate Intersection (caja de herramientas de ArcGIS 10.2), con el tampón como zonas y la capa original de polígonos de tierras federales como clases, para calcular el porcentaje de superficie federal dentro de cada tampón. Sin embargo, la salida está vacía.

He probado lo siguiente:

1. Comprobar geometría
2. Geometría de reparación
3. Colocar ambos en la misma geodatabase
4. Exportar como shapefiles individuales
5. subconjunto de la clase de características de la memoria intermedia (10000 y 1000 filas)
6. Seleccione todo y ejecute la intersección tabulada
7. Nuevo .mxd abierto con sólo los dos archivos de entrada incluidos

Tabular intersección funcionará cuando se resalten unos pocos registros, pero no funcionará en ningún subconjunto > 1000 registros. La superposición de polígonos no parece ser un problema. Conozco R, así que estoy dispuesto a intentarlo, pero he descubierto que se atasca con tantos polígonos.

También intenté utilizar la herramienta de intersección y recibí los mensajes de error "ERROR 999999 Error al ejecutar la función No se encontró la tabla ". Topología no válida". Otros mensajes de intercambio de pila sugieren que esto está relacionado con la memoria.

¿Existe alguna forma de utilizar Arcpy/Python para iterar sobre la clase de características del buffer y realizar el cálculo en los polígonos individuales y concatenar los resultados?

Answer 1

Recomiendo instalar ArcGIS de 64 bits de procesamiento en segundo plano. Esto no se instala con la vainilla ArcGIS Desktop configuración. El archivo se llama algo así como, ArcGIS_BackgroundGP_for_Desktop_1031_145711.exe. Después de la instalación, habilitarlo en ArcMap marcando Opciones de geoprocesamiento, habilitar el procesamiento en segundo plano. Esto debería solucionar el problema.

Answer 2

Primero asegúrate de que tienes los datos localmente, preferiblemente en una máquina con mucha memoria y un SSD... eso puede ayudar. Comprueba los sistemas de coordenadas y asegúrate de que son iguales. He tenido éxito usando QGIS para grandes conjuntos de datos en lugar de ArcGIS (el tiempo de procesamiento se redujo a horas de semanas) ... Alternativamente, usted puede hacer un bucle de la operación en python, pero que puede ser más dolor que simplemente intentarlo en QGIS.

Buena suerte :-)

Answer 3

Estoy de acuerdo en que utilizar el procesamiento en segundo plano de 64 bits puede ser la mejor manera de hacerlo, pero no puedo esperar a nuestra gente de TI. Terminé usando un script de python que hace un bucle a través de la clase de características. Código de abajo sin duda podría mejorarse, pero hizo el trabajo para mí. Cualquier sugerencia de mejora sería apreciada.

import arcpy
from arcpy import env

env.workspace = "C:/GIS/US/FederalLands/Federal_Lands.gdb" ## path to feature class
featureClass = "C:/GIS/US/FederalLands/Federal_Lands.gdb/Fed_25KmBuff" ## path to feature class
fieldList = arcpy.ListFields(featureClass) ### simply pointer to location

row = 1 ## initiate

arcpy.Delete_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out_Temp",data_type="#")
arcpy.Delete_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp",data_type="#")
arcpy.Delete_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp_Tab",data_type="#")
while row < 49779:
    print row
    arcpy.FeatureClassToFeatureClass_conversion(in_features="Fed_25KmBuff",out_path="C:/GIS/US/FederalLands/Federal_Lands.gdb",out_name="Temp",where_clause=("OBJECTID >= " + str(row) + " AND OBJECTID < " + str(row + 100) ))
    arcpy.TabulateIntersection_analysis(in_zone_features="Temp",zone_fields="ORIG_FID",in_class_features="Federal_Lands_aea",out_table="C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp_Tab",class_fields="#",sum_fields="#",xy_tolerance="#",out_units="SQUARE_KILOMETERS")
    if row == 1:
        arcpy.Rename_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp_Tab",out_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out_Temp",data_type="Table")
    elif row > 1:
        arcpy.Merge_management(inputs="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out2;C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp_Tab",output="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out_Temp")

    arcpy.Delete_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out2",data_type="#")
    arcpy.Rename_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out_Temp",out_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Out2")
    arcpy.Delete_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp",data_type="#")
    arcpy.Delete_management(in_data="C:/GIS/US/FederalLands/Federal_Lands.gdb/Temp_Tab",data_type="#")
    row = row + 100

Answer 4

Tuve el mismo problema con la función tabular Intersección de ArcGIS y decidí aprender a utilizar R para abordar el problema. Aquí está mi código utilizando el procesamiento multi-núcleo con un respuesta más detallada aquí para hacerlo utilizando R.

#initiate multicore cluster and load packages
library(foreach)
library(doParallel)
library(tcltk)
library(sp)
library(raster)

cores<- 7
cl <- makeCluster(cores, output="") #output should make it spit errors
registerDoParallel(cl)

Esta es la función:

    multicore.tabulate.intersect<- function(cores, polygonlist, rasterlayer){ 
      foreach(i=1:cores, .packages= c("raster","tcltk","foreach"), .combine = rbind) %dopar% {

    mypb <- tkProgressBar(title = "R progress bar", label = "", min = 0, max = length(polygonlist[[i]]), initial = 0, width = 300) 

    foreach(j = 1:length(polygonlist[[i]]), .combine = rbind) %do% {
      final<-data.frame()
      tryCatch({ #not sure if this is necessary now that I'm using foreach, but it is useful for loops.

        single <- polygonlist[[i]][j,] #pull out individual polygon to be tabulated

        dir.create (file.path("c:/rtemp",i,j,single@data$OWNER), showWarnings = FALSE) #creates unique filepath for temp directory
        rasterOptions(tmpdir=file.path("c:/rtemp",i,j, single@data$OWNER))  #sets temp directory - this is important b/c it can fill up a hard drive if you're doing a lot of polygons

        clip1 <- crop(rasterlayer, extent(single)) #crop to extent of polygon
        clip2 <- rasterize(single, clip1, mask=TRUE) #crops to polygon edge & converts to raster
        ext <- getValues(clip2) #much faster than extract
        tab<-table(ext) #tabulates the values of the raster in the polygon

        mat<- as.data.frame(tab)
        final<-cbind(single@data$OWNER,mat) #combines it with the name of the polygon
        unlink(file.path("c:/rtemp",i,j,single@data$OWNER), recursive = TRUE,force = TRUE) #delete temporary files
        setTkProgressBar(mypb, j, title = "number complete", label = j)

      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) #trycatch error so it doesn't kill the loop

      return(final)
    }  
    #close(mypb) #not sure why but closing the pb while operating causes it to return an empty final dataset... dunno why. 
  }
}

Así que para usarlo, ajuste el single@data$OWNER para que se ajuste con el nombre de la columna de su polígono de identificación (supongo que podría haber sido incorporado en la función ...) y poner en:

myoutput <- multicore.tabulate.intersect(cores, polygonlist, rasterlayer)