# Calculating statistics per area for categorical raster using R?

I am working with the area of New York City (NYC) in ArcGIS 10.1. I have a landcover raster map with several different land cover categories. What I am trying to do is to build up a table in which, per borough, I have the percentage occupied per land cover category.

My first approach was to reclassify the raster, setting every category to `NoData`, keeping only one category as `1`. The next step would be to perform a zonal statistics as table that provides me with the `count` of pixels per borough. Knowing the resolution of my raster, I can then infer the area of the borough covered by that category and then divide it by the borough's total area.

However, this process seems quite repetitive and slow, generating tons of intermediate data (up to 7 reclassify runs).

Can you think of a way to do this faster using R?

• ArcGIS 10.1 is quite old (in fact, it's in Retired status). If you aren't using the terminal service pack and haven't installed the scores of patches, you should do that first. Please take the Tour, which emphasizes the importance of asking one question per Question (offering the option of R makes this two questions). Please Edit the question to focus on the software stack with which you want a solution. Commented Oct 3, 2018 at 20:14
• I reduced the scope of this particular question to R but if you also want to pursue an ArcGIS Desktop solution then I encourage you to ask another question with that as its scope.
– PolyGeo
Commented Oct 3, 2018 at 20:50
• If you're dealing with high-res land cover for NYC (e.g., the 6 inch resolution data available from the city) you might use QGIS, any version >3.0 as R might be slow with such a big dataset. The 'Zonal Histogram' tool gives counts per pixel value by polygon. (On a side-note I do work with NYC Landcover data too; feel free to ping me if it might be worth connecting) Commented Oct 4, 2018 at 20:09

In R, you can extract the raster data for each polygon and then summarize it. First, lets create some data (FYI, you can read in a shapefile using `raster::shapefile` or `rgdal::readOGR` and a raster using `raster::raster`).

``````library(raster)
library(rgeos)
r <- raster::raster(nrows=180, ncols=360, xmn=571823.6, xmx=616763.6, ymn=4423540,
ymx=4453690, resolution=270, crs = CRS("+proj=utm +zone=12 +datum=NAD83
+units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
r[] <- rpois(ncell(r), lambda=1)
r <- calc(r, fun=function(x) { x[x >= 1] <- 1; return(x) } )
p <- rgeos::gBuffer(sampleRandom(r, 5, sp=TRUE), byid=TRUE, width=1620)
plot(r)
``````

Now we can extract the raster values for each polygon. You will see that, since it is aerial data with more than one raster value per polygon, the function returns a list object.

``````( e <- extract(r, p) )
``````

Using `lapply` we can summarize the raster data for each polygon, even returning proportions. The list stays ordered so, any results that you distill from the list objects can be directly related back to the corresponding polygons.

``````( class.counts <- lapply(e, table) )
( class.prop <- lapply(e, FUN = function(x) { prop.table(table(x)) }) )
``````

Here is an example of an unequal number of classes. First, we add an extra class into the results.

``````class.prop[[2]] <- c(class.prop[[2]][1], class.prop[[2]][2]-0.10, 0.10)
names(class.prop[[2]]) <- c("0","1","4")
``````

Then we specify a function that twill handle the unequal vectors and use it to create a data.frame that can be jointed back to the polygon data. This is a less complicated to use version of `plyr::rbind.fill` that does not require an additional package.

``````rbind.fill <- function(x) {
nam <- sapply(x, names)
unam <- unique(unlist(nam))
len <- sapply(x, length)
out <- vector("list", length(len))
for (i in seq_along(len)) {
out[[i]] <- unname(x[[i]])[match(unam, nam[[i]])]
}
setNames(as.data.frame(do.call(rbind, out), stringsAsFactors=FALSE), unam)
}

( p.prop <- rbind.fill(class.prop) )
p <- SpatialPolygonsDataFrame(p, p.prop) #add to polygons
p@data #display data.frame in polygon object
``````
• Thank you very much for this answer! it totaly works out. I am curious about whether there is a similar way of proceeding using the package sf when dealing with shapefiles, since I am much more used to the management of spatial data as data frames. If I get the time to come up with this I will share it too. Commented Oct 22, 2018 at 20:22
• The sf packages does not have support for raster class objects. You really do not want to force the issue (eg., coerce rasters into point sf objects) as you will negate the memory safe aspects of the raster package. Given that a raster can have many million cells, this is a notable consideration. Commented Oct 22, 2018 at 20:39
• Indeed. My consideration goes through using sf to handle shapefiles, but rasterizing them at the time of the data extraction. That way I just "stamp" the values of the raster on a new one rather than relying on the polygon's topologies and geometry, which usually drains the machine's memory. Commented Oct 22, 2018 at 20:42

Update: After dealing with it for a while, I remembered the function "zonal statistics" in ArcMap that basically does this work by returning a table with a count of the pixels with each category per entity. So, after all, the solution was pretty simple.

• Be very weary of the ArcGIS zonal functions. I have documented several bugs with ESRI that seems to keep cropping back up. It seems to be related to the way that the index large rasters and the cursor getting lost. Commented Oct 29, 2018 at 17:57
• I see. That makes sense with other issues I have had before with calculating zonal statistics of entities that overlapped, which I ended up patching with an R code too. Thank you for the warning. Commented Oct 29, 2018 at 20:17