Create Zonal Grid in R

Question

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 1

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.

Answer 2

If I understood correctly you want to create a raster of a certain size with constant values in zones. You can do that in two steps. First create a matrix with your blocks. Then convert it to a raster. Assuming you want k*k blocks, here is a script that creates those zones:

# Set matrix dimension
k <- 8
n <- k^2

# Create matrix
zones <- matrix(NA, n, n)

# Loop over blocks and fill with counter value
counter <- 1
for(i in 1:k){
    for(j in 1:k){
        zones[1:k + (i-1)*k,1:k + (j-1)*k] <- counter
        counter <- counter + 1
    }
}

# Load raster library
library(raster)

# Create zones raster from matrix, setting extents
zones <- raster(zones, xmn = 0, xmx = 10, ymn = 0, ymx = 10, crs = NA)

# Inspect zones raster
plot(zones)

The result for k = 8 looks like this:

Answer 3

I think the kronecker product will work. I wrote a fnction like this, Yraster is a Raster* object:

CreateZonalRasterK <- function(Yraster,k){

  zone<-Yraster
  #define the block size
  mul <- ceiling(nrow(Yraster)/k)
  #use kronecker to create zonal mask
  m1<- matrix(1:k^2,nrow=k,byrow=TRUE)
  m2<-matrix(1,mul,mul)

  m3<-kronecker(m1,m2)
  zone[] <- m3

  return(zone)
}