Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
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 '12 at 15:52
add comment

3 Answers

up vote 4 down vote accepted

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.

Image

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

share|improve this answer
    
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
add comment

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:

Zonal raster for analysis

share|improve this answer
    
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
add comment

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)
}
share|improve this answer
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 '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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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