[prev in list] [next in list] [prev in thread] [next in thread]
List: r-sig-geo
Subject: Re: [R-sig-Geo] Function to smooth grids using n x n window moving
From: ajblist-geo () yahoo ! de
Date: 2008-06-25 13:22:10
Message-ID: 807823.35603.qm () web52208 ! mail ! re2 ! yahoo ! com
[Download RAW message or body]
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@uva.nl>
An: Paul Hiemstra <p.hiemstra@geo.uu.nl>; r-sig-geo@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@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_MONETA \
RY=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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_________________________________________________
> //de.overview.mail.yahoo.com
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic