Tabulate intersection of ArcGIS for Desktop returns empty output?

Question

I am using a US federal lands polygon layer (http://coastalmap.marine.usgs.gov/GISdata/basemaps/boundaries/fedlands/fedlanp020.htm)

I downloaded the layer, reprojected as Albers Equal Area, and then created a 25 km buffer around the polygons. I have been trying to use Tabulate Intersection (ArcGIS 10.2 toolbox), with the buffer as the zones and the original federal lands polygon layer as the classes in order to calculate the percentage area within each buffer that is federal. However, output is empty.

I have tried the following:

1. Check geometry
2. Repair geometry
3. Place both in the same geodatabase
4. Export as individual shapefiles
5. subset the buffer feature class (10000 and 1000 rows)
6. Select all and then run tabulate intersection
7. Opened new .mxd with only the two input files included

Tabulate Intersection will work when a few records are highlighted, but won't work on any subsets > 1000 records. Overlapping polygons don't seem to be an issue. I know R, so am willing to give that a shot, but have found that it bogs down with this many polygons.

I also attempted to use the intersect tool and received the error messages "ERROR 999999 Error executing function . .Table was not found . . . Invalid Topology". Other stack exchange messages suggest this relates to memory.

Is there a way using Arcpy/Python to iterate over the buffer feature class and perform the calculation on the individual polygons and concatenate the results?

Answer 1

I recommend installing 64bit ArcGIS Background Processing. This is not installed with the vanilla ArcGIS Desktop setup. The file is named something like, ArcGIS_BackgroundGP_for_Desktop_1031_145711.exe. After the install, enable it in ArcMap by checking Geoprocessing Options, enable Background Processing. This should fix your problem.

Answer 2

First make sure you have the data locally, preferably on a machine with lots of memory and an SSD... that can help. Check the coordinate systems and make sure they're the same. I've had success using QGIS for big data sets rather than ArcGIS (processing time down to hours from weeks)... Alternatively you can loop the operation in python, but that may be more pain than just trying it in QGIS.

Good luck :-)

Answer 3

I agree that using the 64 bit background processing may be the best way to go, but I can't wait for our IT people. I ended up using a python script that looped through the feature class. Code below could definitely be improved, but it got the job done for me. Any suggestions on improvement would be appreciated.

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

I had the same problem with ArcGIS's tabulate Intersection function and decided to learn to use R to tackle the problem. Here's my code using multi-core processing with a more detailed answer here, to do this using 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)

Here's the function:

    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. 
  }
}

So to use it, adjust the single@data$OWNER to fit with the column name of your identifying polygon (guess that could have been built into the function...) and put in:

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