I've been running into all sorts of issues using ArcGIS ZonalStats and thought R could be a great way. I'm fairly new to R but I have a coding background.

The situation is that I have several rasters and a polygon shapefile with many features of different sizes (though all features are bigger than a raster cell and the polygon features are aligned to the raster). I've figured out how to get the mean value for each polygon feature using the raster library with extract:

#load packages required
# ---Set the working directory-------
datdir <- "/test_data/"

#Read in grid of water depth
ras <- raster("test_data/raster/pl_sm_rp1000/w001001.adf")

#read in polygon shape file
proxNA <- shapefile("test_data/proxy/PL_proxy_WD_NA_test") 

#calc mean depth per polygon feature
#unweighted - only assigns grid to district if centroid is in that district
proxNA$RP1000 <- extract(ras, proxNA, fun = mean, na.rm = TRUE, weights = FALSE)

#plot depth values 

The issue I have is that I also need an area based ratio between the area of the polygon and all non NA cells in the same polygon. I know what the cell size of the raster is and I can get the area for each polygon, but the missing link is the count of all non-NA cells in each feature. I managed to get the cell number of all the cells in the polygon proxNA@data$Cnumb1000 <- cellFromPolygon(ras, proxNA)and I'm sure there is a way to get the actual value of the raster cell, which then requires a loop to get the number of all non-NA cells combined with a count, etc.

Do you have an idea or can you point me in the right direction?


4 Answers 4


The example data from Jeffrey, now using "terra"

r <- rast(ncols=10, nrows=10)
x <- runif(ncell(r))
x[round(runif(25,1,100),digits=0)] <- NA
values(r) <- x

cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- vect(list(cds1, cds2), type="polygons", crs="+proj=longlat")


extract(r, polys, fun=\(x) length(na.omit(x))/length(x))
#  ID     lyr.1
#1  1 0.8333333
#2  2 0.6666667

But it may be faster like this

extract(not.na(r), polys, fun=mean)
#  ID     lyr.1
#1  1 0.8333333
#2  2 0.6666667

If you have many rasters, first combine them into a single SpatRaster (if they have the same extent and resolution)

To get the polygon area you can use expanse

expanse(polys, unit="km") 
#[1] 60913738 31611037

I am not sure if you want the ratio based on the "real value" of the polygon(s) areas or the areas of the cells intersecting them. Here is some example code that uses all cells intersecting the polygons (basically, ratio of NA cells to non-NA cells). It is a dummy example and you will need to write your own function.

    # Create some example data

    r <- raster(ncols=10, nrows=10)
      x <- runif(ncell(r))
        x[round(runif(25,1,100),digits=0)] <- NA
          r[] <- x
      cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
        cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
          polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                                   Polygons(list(Polygon(cds2)), 2)))
            polys <- SpatialPolygonsDataFrame(polys, data.frame(ID=sapply(slot(polys, "polygons"), 
                                              function(x) slot(x, "ID"))))
  plot(polys, add=TRUE)

You can use this code snippet to add an area column to your polygon data by extracting from the area slot. This could be used if you want to ratio using the "real" polygon area.

# Add area of polygon(s)
polys@data <- data.frame(polys@data, Area=sapply(slot(polys, 'polygons'), 
                         function(i) slot(i, 'area')))  

The most efficient, and considerable faster, alternative to for loops are "apply" like functions. There are a number of these available in R that are utilized for different object classes or data structures. In this case, since extract returns a list, we will use lapply (list apply). This is a way to apply a base or custom function to a list object. The the object class stored in the list is a vector, the function is quite straight forward. If you use extract on a brick or stack raster object the resulting objects stored in the list would be data.frame objects.

# On a single raster object, extract returns list object with stored vectors.                           
( vList <- extract(r, polys, na.rm=FALSE) )

# Use lapply to apply function that calculates ratio of NA to non-NA values
#   wrapping lapply in unlist() collapses result into a vector  
aRatio <- function(x) { if(length(x[is.na(x)]) > 0) (length(x[is.na(x)]) / length(x[!is.na(x)])) else 0 }  
  ( vArea <- unlist( lapply(vList, FUN=aRatio ) ) )

# Assign ordered vector back to polygons
polys@data <- data.frame(polys@data, NAratio=vArea)

I do not have access to your files, but based on what you described, this should work:

masked_img=mask(nonNA_raster,mask_layer) #based on centroid location of cells
nonNA_count=cellStats(masked_img, sum)

Since this answer was originally posted, there is now a much faster method using the exactextractr package.

exactextractr::exact_extract(r, polys, 'count')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.