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
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
