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?

-
I can make sense of the phrase "k blocks" only by supposing each block has k elements. If by "k" you mean the number of blocks, then usually this problem has no solution (even when k divides n*n) and when it does it has multiple solutions (if we understand that a "block" might be rectangular). For instance, k=6 divides 6*6, but there is no way to decompose a 6*6 raster exactly into 6 square blocks. It can be decomposed into exactly 6 rectangular blocks of dimensions 1 x 6, 2 x 3, 3 x 2, or 6 x 1. – whuber Oct 4 '12 at 15:52

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.

-
I think your approach is very nice. The only concern for me is the smaller rectangles around the edge area. I think for ecological study, it requires almost all the grids to have similar area. – Seen Oct 5 '12 at 12:33
All you need to do is adjust the values of k1 and k2 to fill the region as closely as you can--the additional flexibility of this solution enables you to do that. – whuber Oct 5 '12 at 12:37

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
}
}

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:

-
Thanks for your answer @yellowcap. I am wondering if there is anyway to get rid of the double for loop. It could be very slow for a large matrix. – Seen Oct 4 '12 at 12:58
Sorry, no ideas on how to remove the double for loop. But I just tried it with `k=100` i.e. `n=10000` and the for-loop lasted about a second. – yellowcap Oct 4 '12 at 14:23
I think the kronecker product will work and I uploaded my answer. Thank you very much though! – Seen Oct 4 '12 at 14:34
`R` and loops don't work well together: `R` encourages you to use built-in functional programming constructs instead. In this case, the double looping is approximately 1,000 times slower than other approaches, so it's only good for small illustrations. – whuber Oct 4 '12 at 15:48

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)
}
``````
-
`kronecker` may be the basis of a good solution. But please run your code: it does not do what is described in the question. Note, too, that it may do some additional stuff you don't want--such as padding out the array--and it's built on `outer`, so it shouldn't be any more efficient than `outer` itself. – whuber Oct 4 '12 at 15:41
Thanks @whuber, I think my function just works for square grids which I need. Also, what do you man by "padding out the array" for this method? – Seen Oct 5 '12 at 12:30
Please see the help page for `kronecker`. – whuber Oct 5 '12 at 12:37