[R-sig-Geo] Function to smooth grids using n x n window moving average (low-pass filter)

[Paul Hiemstra]

Dear list,

I recently implemented a smoothing algorithm that could be interesting 
for other people. It smooths a grid by calculating an average for a n x 
n window. The input is a SpatialGrid/PixelsDataframe. For a 3x3 window 
the algorithm creates 8 shifted matrices in addition to the original 
matrix (derived from the grid-object). Adding up these 9 matrices and 
dividing them by 9 gives the smoothed matrix.  The output is a 
SpatialGridDataFrame. Because the algorithm uses just a few simple 
matrix operations it is very fast and scalable to large grid-objects. 
Maybe the code could be added to the sp-package (Edzer or Roger?)?

This is the code, just three functions. I also provided a small but 
functioning example at the bottom:

shift_matrix = function(m, row, col) {
    nrow = dim(m)[1]
    ncol = dim(m)[2]
    if(row > 0) m = rbind(matrix(rep(NA,abs(row)*ncol), 
abs(row),ncol),m[1:(nrow-row),])
    if(row < 0) m = 
rbind(m[-(row-1):nrow,],matrix(rep(NA,abs(row)*ncol), abs(row),ncol))
    if(col > 0) m = cbind(matrix(rep(NA,abs(col)*nrow), nrow, 
abs(col)),m[,1:(ncol-col)])
    if(col < 0) m = 
cbind(m[,-(col-1):ncol],matrix(rep(NA,abs(col)*nrow), nrow, abs(col)))
    m
}

smooth_matrix = function(m, kernel_size = 3) {
    row_nums = rep(floor(kernel_size/2): -floor(kernel_size/2), each = 
kernel_size)
    col_nums = rep(floor(kernel_size/2): -floor(kernel_size/2), kernel_size)
    out = matrix(rep(0,dim(m)[1]*dim(m)[2]), dim(m))
    for(i in 1:length(row_nums)) {
        out = out + shift_matrix(m, row_nums[i], col_nums[i])
    }
    return(out / (kernel_size * kernel_size))
}

smooth_grid = function(grd, zcol, kernel_size = 3, add_to_name = 
"_smooth") {
    if(!fullgrid(grd)) fullgrid(grd) = TRUE
    grd[[paste(zcol,add_to_name, sep = "")]] = 
as.vector(smooth_matrix(as.matrix(grd[zcol]), kernel_size = kernel_size))
    grd
}

## Example script
library(gstat)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y
meuse.grid$log.zinc = idw(log(zinc)~1, meuse, meuse.grid)$var1.pred
meuse.grid = smooth_grid(meuse.grid, "log.zinc", kernel_size = 3)
spplot(meuse.grid, c("log.zinc","log.zinc_smooth"))

cheers,
Paul

ps The system i run:
R version 2.7.0 (2008-04-22)
i486-pc-linux-gnu

locale:
LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] gstat_0.9-45 sp_0.9-25

loaded via a namespace (and not attached):
[1] grid_2.7.0     lattice_0.17-6 tools_2.7.0

-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: 	+31302535773
Fax:	+31302531145
http://intamap.geo.uu.nl/~paul