# Zonal statistics - computation time issue with assigning zones

I would like to upscale a 1/360 degree resolution raster (64800 x 129600) to a 1/12 resolution raster (2160 x 4320) by assigning to each cell in the low-res raster the number of cells in the high-res raster that are contained in that cell AND that have a certain value. In order to use the 'zonal' function, I first need to assign zones to the grid points of the high-res raster. Every 30x30 block of grid points should be one zone. (I.e. the zone raster should equal '1' in the top left 30x30 grid points, and 2160*4320=9331200 in the bottom right 30x30 grid points). Here is the code I am using to generate this zone raster ('raster_high' is the original high-res raster.)

``````zones <- raster_high;
for(i in 1:64800)
{
j <- floor(i/30);
zones[i,] <- rep(4320*j+1:4320*(j+1), each=30);
}
``````

However, even only the first step i=1 takes ages. Is there a quicker way to do this?

``````library(raster)

# MAke tests raster with integer values
r0 <- raster()
r0[] <- floor(rnorm(ncell(r), 30, 5))

# Define custom aggregation function
fun <- function(x, val = 28, ...) {
length(x[x == val])
}

# Aggregate raster using custom function for spatial aggregation
r1 <- aggregate(r0, 30, fun = fun)
``````

If you're concerned with `NAs`, you can use the following aggregating function instead.

``````fun2 <- function(x, val = 28, na.rm) {
sum(x[x == val], na.rm = na.rm)
}

fun(c(1,1,2,3,NA), 1)
# [1] 3

fun2(c(1,1,2,3,NA), 1, na.rm = TRUE)
# [1] 2
``````
• (+1) `aggregate` is the right tool. Consider using `sum(x==val)` for the custom function--it's clear and more efficient. In fact, applying `aggregate` (with the `sum` function) to the indicator grid `r0==val` can be quite a bit faster. I'm getting very bad scaling (with both methods) for huge grids, though, even when though they fit in RAM. Commented Jan 16, 2017 at 23:05