
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@stat.math.ethz.ch on behalf of Paul Hiemstra
Sent: Mon 6/23/2008 2:58 PM
To: r-sig-geo@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@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


	[[alternative HTML version deleted]]


