[GIS] Counting raster values using another raster as mask in R

aggregationrrasterzonal statistics

I would like to count the frequency of values in one raster, using another raster (different resolution) as the zonal condition.

The premise is that the coarser resolution raster represents the area that I would like to count the frequency of the finer resolution raster;

# make 2 rasters of the same extent, different resolutions
ext <- extent(0,1000,0,1000)
r1 <- raster(nrows=1000, ncols=1000,ext)
r1[] <- sample(seq(from = 1, to = 6, by = 1), size = 1000000, replace = TRUE)
r2 <- raster(nrows=10, ncols=10,ext)
r2[] <- sample(seq(from = 0, to = 1, by = 0.05), size = 100, replace = TRUE)
# create areas of interest in coarser raster
r2[r2 < 0.9] <- NA

Now i can't seem to find a function (might be blind) that will tell me what i want to know straight off – i thought zonal{raster} might have but it wont return a count function. So this is my work around:

# disaggregate the coarser raster to the same res as the finer raster
r2 <- disaggregate(r2,fact=c(100,100))
# mask the fine by the coarse 
r3 <- mask(r1,r2)
# and return a frequency table of the finer resolution
freq(r3,useNA="no")

But this seems a little round the houses.

ISSUE 1: Is there a function to zonal count in R with differing resolution rasters?

ISSUE 2: I up-scaled the above method to 2 very large rasters but got the error "Error: Failure during raster IO", why so? ## FIXED: issue with PC, not R ##

ISSUE 3: what if I change each cell in the coarser resolution to have an ID instead of a value and I want to count the frequency of values in the finer res raster per ID?

(have also tried changing coarse res raster to binary 1/NA and multiplying but also get the same error as Issue 2 – im using a powerful computer that has worked with bigger data in big stacks, so the issue is not that).

Best Answer

You could approach this as a raster/vector integration problem, where your course resolution data are essentially sampling blocks represented as polygons.

First, let's create your example data. 

library(raster)
library(sp)
ext <- extent(0,1000,0,1000)
  r1 <- raster(nrows=1000, ncols=1000,ext)
  r1[] <- sample(seq(from = 1, to = 6, by = 1), size = 1000000, replace = TRUE)
  r2 <- raster(nrows=10, ncols=10,ext)
  r2[] <- sample(seq(from = 0, to = 1, by = 0.05), size = 100, replace = TRUE)
  r2[r2 < 0.9] <- NA

Now we can coerce the coarse resolution raster to a SpatialPolygonsDataFrame. 

r2 <- rasterToPolygons(r2)

Using the raster::extract function we create a list object of raster values for each polygon. We can then use apply calculate frequencies with table. 

    r2.dat <- extract(r1, r2)  

# Just count values
( f <- do.call(rbind, lapply(r2.dat, FUN = function(x) { return( table(x)) })) ) 

# or, proportions
( f <- do.call(rbind, lapply(r2.dat, FUN = function(x) { return( prop.table(table(x))) })) )

The results will be ordered the same as the polygons and can be joined back to the polygon object. The only issue here is that the code is assuming that all values will be present in the frequencies. If this is not the case the vectors resulting from apply will not be equal. All you need to do if this error occurs is iterate through the list and add NA values where expected values are missing. You would need to write a function that is passed to lapply. 

r2@data <- data.frame(r2@data, as.data.frame(f))
  names(r2@data)[2:ncol(r2)] <- paste( "value", colnames(f), sep=".") 
  r2@data
Related Question