rmatrixfilterrastergaussianblur

Add a gaussian blur with specified sigma to a satellite image using R


I am trying to add a gaussian blur with sigma = 0.5 to a satellite image. I have found the package spatialEco has a function called gaussian.kernel which does exactly what I want. So basically I convert my raster image into a matrix, then I create a gaussian matrix with sigma = 0.5 and finally, I multiply those two matrices. After the multiplication I am receiving an error: Error in gm * m : non-conformable arrays. I am guessing it's because in the function gaussian.kernel I don't set the n parameter to equals the dimensions of my matrix. Can some1 help me? Here is the code:

library(raster)
library(spatialEco)

pan = raster("path/pan.tif")

m = as.matrix(pan)
dim(m) #2775 1435 are the dimensions of the matrix.

gm = spatialEco::gaussian.kernel(sigma = 0.5, n = ????)

wm = gm * m

Solution

  • Just use raster.gaussian.smooth directly on the RasterLayer object as in:

    smoothed = spatialEco::raster.gaussian.smooth(pan, sigma = 0.5)
    plot(pan)
    plot(smoothed)
    

    Gaussian filters do not work like you described. The dimension of the kernel is much smaller than the image. This is because submatrices from the initial image are convoluted with the kernel in so-called focal or moving window operations. Just give the source code a look:

    # This is the code for spatialEco::gaussian.kernel
    gaussian.kernel <- function(sigma=2, n=5) 
    {
      m <- matrix(ncol=n, nrow=n)
      mcol <- rep(1:n, n)
      mrow <- rep(1:n, each=n)
      x <- mcol - ceiling(n/2)
      y <- mrow - ceiling(n/2)
      m[cbind(mrow, mcol)] <- 1/(2*pi*sigma^2) * exp(-(x^2+y^2)/(2*sigma^2))
      m / sum(m)
    }
    
    # This is the code for spatialEco::raster.gaussian.smooth
    raster.gaussian.smooth <- function(x, sigma = 2, n = 5, type = mean, ...) {  
      if (!inherits(x, "RasterLayer")) stop("MUST BE RasterLayer OBJECT")
      gm <- gaussian.kernel(sigma=sigma, n=n)
      return(raster::focal(x, w = gm, fun = type, na.rm=TRUE, pad=FALSE, ...) )
    }  
    

    The code for raster::focal can be found at Raster's Github repository. Also, for your information, matrix multiplication in R is done with the %*% operator rather than simply *.