raster - Create Zonal Grid in R

Suppose I have a n*n raster, and I want to create k square blocks (k can be divided exactly by n*n ) for zonal analysis:

for example, when n = 4 and k = 4 a 4*4 raster is create with value

1 1 2 2
1 1 2 2
3 3 4 4
3 3 4 4

How can I do this in R?

Answer

It would be nice to have the code execute correctly even when the block dimensions do not evenly divide the raster dimensions. Indeed, it's not any harder to create rectangular blocks (say, k1 rows by k2 columns) within a rectangular grid (of, say, n1 rows by n2 columns). To solve the problem as stated, set k1 = k2 = sqrt(k) and n1 = n2 = n.

For example, let's set up to create a 15 million cell raster with large rectangular blocks:

k1 <- 897; k2 <- 654; n1 <- 3 * 10^3; n2  <- 5 * 10^3

The solution is still a one-liner:

x <- outer(1:n1, 1:n2, function(i,j) (i-1) %/% k1 * ((n2+1) %/% k2) + (j-1) %/% k2 + 1)

It's reasonably fast: the timing for this operation was 1.34 seconds, more than ten million cells per second.

You can offset the blocks by changing the two "-1"s in the code to other (integral) values.