Following Paul Hiemstra's code Function to smooth grids using n x n window
moving average (low-pass
filter)<https://stat.ethz.ch/pipermail/r-sig-geo/2008-June/003794.html>
,
I adapted it so I could apply any function (mean, max, min, quantile, etc.
). Sadly, using apply make it ten times slower, but maybe someone can
improve it.

I don't know if it could be of any use, and if it's faster than SAGA's
focal.function suggested by Alexander Brenning as I haven't make the
comparison, but if someone have an idea how to improve it, don't hesitate to
share it.

functions not prefixed with apply are from Paul Hiemstra's code...

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))
}
apply_matrix = function(m, kernel_size = 3, FUN, ...) {
    km <- floor(kernel_size/2)
    row_nums = rep(km: -km, each = kernel_size)
    col_nums = rep(km: -km, kernel_size)
    vec = NULL
    for(i in 1:length(row_nums)) {
            vec = rbind(vec, as.vector(shift_matrix(m, row_nums[i],
col_nums[i])))
    }
    out <- apply(vec, 2,FUN,...)
    # return
    matrix(out,ncol=dim(m)[2])
}
smooth_grid = function(grd, zcol, kernel_size = 3, add_to_name =
paste("_smooth_",kernel_size,sep="")) {
    if(!fullgrid(grd))
        fullgrid(grd) = TRUE
    newBandName = paste(zcol,add_to_name, sep = "")
    grd[[newBandName]] = as.vector(smooth_matrix(as.matrix(grd[zcol]),
kernel_size = kernel_size))
    grd
}
apply_grid = function(grd, band, kernel_size=3, add_to_name = NULL, FUN,
...){
    if(!fullgrid(grd))
        fullgrid(grd) = TRUE
    if (is.null(add_to_name)) {
        newBandName = paste(band, "_",
deparse(substitute(FUN)),"_",kernel_size,sep="")
    }else{
        newBandName = paste(band, "_", add_to_name, sep="")
    }

    mat = apply_matrix(as.matrix(grd[band]), kernel_size = kernel_size,
FUN=FUN, ...)
    grd[[newBandName]] = as.vector(mat)
    grd
}

#Usage example (and comparison) :

library(rgdal)
grd = readGDAL(system.file("pictures/Rlogo.jpg", package = "rgdal")[1],
band=1)
summary(y)
## Mon exemple
# grd <- readGDAL(choose.files(getwd(),caption="Select raster",
multi=FALSE))
cat("smooth 3:")
system.time(grd <- smooth_grid(grd,zcol="band1", kernel_size=3))
# les bandes sont ajoutées (si le kernel est différent)
cat("smooth 5:")
system.time(grd <- smooth_grid(grd, zcol="band1", kernel_size=5))
cat("mean 3:")
system.time(grd<-apply_grid(grd,"band1", kernel_size=5,FUN=mean))
cat("mean 5:")
system.time(grd<-apply_grid(grd,"band1", kernel_size=3,FUN=mean))
# you can also avoid border effect by using na.rm=TRUE, however you get a
smaller kernel on the borders
cat("median 5:")
system.time(grd<-apply_grid(grd,"band1",kernel_size=5, FUN=median,
na.rm=TRUE))
# calculate difference
grd@data$dif_3 <- grd@data$band1_mean_3-grd@data$band1_smooth_3
grd@data$dif_5 <- grd@data$band1_mean_5-grd@data$band1_smooth_5
spplot(grd,col.regions=grey(1:25/25))
summary(grd@data$dif_5)
summary(grd@data$dif_3)

--
Etienne

	[[alternative HTML version deleted]]

