# Raster aggregation in ArcGIS

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?

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

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")
``````

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

-