Counting raster values using another raster as mask in R?

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
# 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).

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
• thanks for your answer. I'm loathe to convert to vector when i'll be working with massive rasters, it's something i thought over in the original testing. – Sam Feb 10 '16 at 9:24
• however, if there is no direct function for zonal stats (count) for rasters of two resolutions, then perhaps i'll have to test my work around vs your workaround and see what will work better after upscaling. Yours has the advantage of working per ID as well. – Sam Feb 10 '16 at 9:25
• I looked at a way to keep the data raster but I believe that you will run into memory issues in short order. – Jeffrey Evans Feb 10 '16 at 17:10
• depends, i have a lot of memory at my disposal – Sam Feb 11 '16 at 12:20
• @JeffreyEvans: I know this is an old post, but I'm struggling with your last function to insert NAs where ther are missing values in the list. The line r2@data <- data.frame(r2@data, as.data.frame(f)) returns "Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, : arguments imply differing number of rows: 0, 285504" At first I thought "f" should be "r2.dat" since you mentioned iterating through the list, but that change also doesn't work. Any thoughts? – ecologist1234 Feb 23 '18 at 15:01