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?

  • 1
    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, 2012 at 15:52

3 Answers 3


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, 2012 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, 2012 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

# Load raster library

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

# Inspect zones raster

The result for k = 8 looks like this:

Zonal raster for analysis

  • 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, 2012 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, 2012 at 14:23
  • I think the kronecker product will work and I uploaded my answer. Thank you very much though!
    – Seen
    Oct 4, 2012 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, 2012 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){

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

  zone[] <- m3

  • 1
    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, 2012 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, 2012 at 12:30
  • Please see the help page for kronecker.
    – whuber
    Oct 5, 2012 at 12:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.