The raster package for R provides a variety of functions for the analysis of raster GIS data. The *focal()* function is very useful for applying moving window filters to such data. I wanted to calculate a moving window mean for cells within a specified radius, but *focal()* did not provide a built-in option for this. The following code generates an appropriate weights matrix for implementing such a filter, by using the matrix as the *w* argument of *focal()*.

require(raster) #function to make a circular weights matrix of given radius and resolution #NB radius must me an even multiple of res! make_circ_filter<-function(radius, res){ circ_filter<-matrix(NA, nrow=1+(2*radius/res), ncol=1+(2*radius/res)) dimnames(circ_filter)[[1]]<-seq(-radius, radius, by=res) dimnames(circ_filter)[[2]]<-seq(-radius, radius, by=res) sweeper<-function(mat){ for(row in 1:nrow(mat)){ for(col in 1:ncol(mat)){ dist<-sqrt((as.numeric(dimnames(mat)[[1]])[row])^2 + (as.numeric(dimnames(mat)[[1]])[col])^2) if(dist<=radius) {mat[row, col]<-1} } } return(mat) } out<-sweeper(circ_filter) return(out) }

This example uses a weighs matrix generated by *make_circ_filter()* to compute a circular moving average on the Meuse river grid data. For a small raster like this, the function is more than adequate. For large raster datasets, it’s quite slow though.

#make a circular filter with 120m radius, and 40m resolution cf<-make_circ_filter(120, 40) #test it on the meuse grid data f <- system.file("external/test.grd", package="raster") r <- raster(f) r_filt<-focal(r, w=cf, fun=mean, na.rm=T) plot(r, main="Raw data") #original data plot(r_filt, main="Circular moving window filter, 120m radius") #filtered data

Advertisements

Amy WhiteheadThanks! I’ve been looking for something that I can use to calculate the percentage of a particular habitat within a specified radius of each cell and this looks like it will do the job nicely.

scrogsterPost authorI’m glad the code was useful to you Amy. That’s exactly the task that I use it for myself, though it is a little bit slow for large datasets.

Dan RosauerThanks Michael, this was really timely help! Precisely what I needed.

L. VelascoThanks for the code Michael. I do not quite understand the res parameter. Let’s suppose I have a raster with resolution 0.005 and epsg:4326 projection, and I would like to obtain a moving window filter with a radius of 20km, what values should take the parameters?

I would like to ask you as well whether you know how to calculate the length in meters of one cell of my raster taking into account the resolution and projection.

Many thanks in advance

Pingback: Clump and eliminate classes in a raster in ArcMap or R | Question and Answer

betterworldenergyYou could also use Factorial Kriging. This could be used to remove or ‘isolate’ various spatial components (removing nugget, long range only, short range only). Just convert raster to points at pixel centroids, then save as a raster to get the final image.

afana72Thanks for the code Michael. Really it’s very useful for raster applications, but I have a question for you. What if I need to filter a vector spatial data formed by x,y,z? do you have anything in mind to handle such problem? I tried to have a look at the available libraries and I couldn’t fund anything.

Ashraf, cheers