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.

I'd like to aggregate a raster from 1m to 6m in ArcGIS with a custom rule. For example, each cell value of aggregated raster = mean value of 5 highest cell values from 36 aggregated values. How do I do that?

share|improve this question
2  
(1) In aggregating 6:1, there will be 6^2 = 36 values rather than 25. (2) Using ArcGIS, this task is possible but not easy. Would you consider a Python or R solution? –  whuber Dec 17 '12 at 16:04
    
Thanks, 36 indeed :) Python, R or any other solution would be much appreciated, though Python seems a little bit more friendly for ArcGIS user. –  Tomek Dec 17 '12 at 20:45

1 Answer 1

up vote 3 down vote accepted

Here is an R solution. Other questions on this site demonstrate how to read and write raster files in R, so let's get right to the solution. First, some random data to illustrate the procedure:

n.rows <- 1000
n.cols <- 1000
r <- outer(1:n.rows, 1:n.cols, 
  function(x,y) sin(x/100) * cos(((100+x)/(1+(y/100)^2)))^2) + rnorm(n.rows*n.cols, 0, .10)

Here's the code. It applies a user-supplied function f to each 6 by 6 block of the raster r; f finds and averages the five largest values. The code uses sapply to break the array into columns of width 6 and calls aggregate to spit out a column of statistics within each such vertical band. (I coded it in this order--split the columns first--because R stores arrays in column-major order. This will maximize the locality of reference, which is important for speed on very large arrays.)

f <- function(x) {y <- as.vector(x); mean(y[order(y, decreasing=TRUE)[1:5]])}
aggregate <- function(x, k=dim(x)[2], fun) sapply(seq(1, dim(x)[1]+1-k, k), 
  function(i) fun(x[1:k+i-1, ]))
agg <- sapply(seq(1, dim(r)[2]+1-6, 6), function(j) aggregate(r[, 1:6+j-1], k=6, fun=f))

This takes about two seconds for the sample grid of a million cells. Plots help us compare them:

library(raster)
r.r <- raster(r); agg.r <- raster(agg)
par(mfrow=c(1,2))
image(r.r, zlim=c(min(r), max(r)), col=terrain.colors(300), main="Original")
image(agg.r, zlim=c(min(r), max(r)), col=terrain.colors(300), main="Aggregate")

Maps

It is evident that the aggregate tends to pick out the larger values in each block.

share|improve this answer

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.