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
g...@data$dif_3 <- g...@data$band1_mean_3-grd@data$band1_smooth_3
g...@data$dif_5 <- g...@data$band1_mean_5-grd@data$band1_smooth_5
spplot(grd,col.regions=grey(1:25/25))
summary(g...@data$dif_5)
summary(g...@data$dif_3)
--
Etienne
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
