Applying a circular moving window filter to raster data in R

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

7 thoughts on “Applying a circular moving window filter to raster data in R

  1. Amy Whitehead

    Thanks! 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.

    Reply
    1. scrogster Post author

      I’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.

      Reply
  2. L. Velasco

    Thanks 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

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

  4. betterworldenergy

    You 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.

    Reply
  5. afana72

    Thanks 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

    Reply

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s