[R-sig-Geo] Function to smooth grids using n x n window moving average (low-pass filter)
ajblist-geo at yahoo.de
ajblist-geo at yahoo.de
Wed Jun 25 15:22:10 CEST 2008
Hi,
FYI, the RSAGA package provides a function called 'focal.function', which can apply an R function operating on circular or quadratic moving windows to an entire grid. This function is currently very inefficient - but flexible.
Cheers
Alex
--
Alexander Brenning
brenning -at- uwaterloo.ca - T +1-519-888-4567 ext 35783
Department of Geography and Environmental Management
University of Waterloo
200 University Ave. W - Waterloo, ON - Canada N2L 3G1
http://www.fes.uwaterloo.ca/geography/faculty/brenning/
----- Ursprüngliche Mail ----
Von: "Hengl, T." <T.Hengl at uva.nl>
An: Paul Hiemstra <p.hiemstra at geo.uu.nl>; r-sig-geo at stat.math.ethz.ch
Gesendet: Montag, den 23. Juni 2008, 22:03:17 Uhr
Betreff: Re: [R-sig-Geo] Function to smooth grids using n x n window moving average (low-pass filter)
Hi Paul,
Thanks for your info and script.
A year ago, I also felt that much of image processing is really missing in R. At the moment, I am increasingly using now RSAGA, which is efficient, fast and easy (+SAGA is still developing; +it works both on Linux and Windows OS). Here are few examples to run e.g. edge filter (the only extra step is to export the sp grid layers to SAGA GIS format, but we are working on improving that):
> library(RSAGA)
> rsaga.env(path="C:/Program Files/saga_vc") # Windows OS
> # export the map to SAGA format:
> rsaga.esri.to.sgrd(in.grids="DEM.asc", out.sgrds="DEM.sgrd"), in.path="D:/MAPS")
> rsaga.filter.simple(in.grid="DEM.sgrd", out.grid="edge_filter.sgrd", mode="circle", method="edge"), radius=1000)
> # you can also run the filter using the generic module runner:
> rsaga.get.modules("grid_filter")
$grid_filter
code name interactive
1 0 Simple Filter FALSE
2 1 Gaussian Filter FALSE
3 2 Laplacian Filter FALSE
4 3 Multi Direction Lee Filter FALSE
5 4 User Defined Filter (3x3) FALSE
6 5 Filter Clumps FALSE
> rsaga.get.usage("grid_filter", 0)
> rsaga.geoprocessor(lib="grid_filter", module=0, param=list(INPUT="DEM.sgrd",RESULT="edge_filter.sgrd",METHOD=2,RADIUS=1000))
> # import the map to R:
> rsaga.sgrd.to.esri(in.sgrds="edge_filter.sgrd", out.grids="edge_filter.asc", out.path="D:/MAPS", prec=3)
> library(rgdal)
> edgef = readGDAL("edge_filter.asc")
see also: https://stat.ethz.ch/pipermail/r-sig-geo/2008-January/003078.html
Tom Hengl
http://spatial-analyst.net
-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch on behalf of Paul Hiemstra
Sent: Mon 6/23/2008 2:58 PM
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Function to smooth grids using n x n window moving average (low-pass filter)
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
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_________________________________________________
://de.overview.mail.yahoo.com
More information about the R-sig-Geo
mailing list