Comparing (gaussian) smoothing algorithms - spatialEco vs terra

[nikosGeography]

I am comparing smoothing algorithms in R. More specifically, the raster.gaussian.smooth from the spatialEcopackage and the focalfrom the terrapackage. Visually, the results look similar but the pixel values are way different. It seems that raster.gaussian.smooth rescales the data or something. In the picture below, the first layer is my original raster, the second is the result after using raster.gaussian.smooth and the third is the result from terra's focal:

How can I stop raster.gaussian.smooth from rescaling (or whatever it does) the data? is this possible?

Here is the code:

library(terra)
library(spatialEco)

r = rast("path/tirs.tif")

# spatialEco
smoothed = spatialEco::raster.gaussian.smooth(r, sigma = 0.5 * 350)

# terra focal
gf <- focalMat(r, 0.5*350, "Gauss")
r_gf <- focal(r, w = gf)

Here is the raster I am using:

tirs = rast(ncols=578, nrows=449, nlyrs=1, xmin=504400, xmax=562200, ymin=155700, ymax=200600, names=c('B10_median'), crs='EPSG:27700')

[jeffreyevans]

Nikolaos, You are using two different focal sizes for the Gaussian weights. The terra::focalMat function is resulting in an 11x11 matrix whereas, since you are not providing an “n” argument to spatialEco::raster.gaussian.smooth the default is using a 5x5 matrix. If you define the matrix size so that they match the results are quite a bit closer.

I do not know what the focalMat function is doing to derive the Gaussian matrix but, you do get different results with different rasters whereas, the spatialEco::gaussian.matrix function (under the hood of raster.gaussian.smooth) returns a consistent symmetric Gaussian matrix explicitly on based the defined sigma and size.

library(terra)
library(spatialEco)

 r <- rast(ncols=578, nrows=449, nlyrs=1, xmin=504400, xmax=562200, 
           ymin=155700, ymax=200600, names=c('B10_median'), 
           crs='EPSG:27700')
    r[] <- runif(ncell(r))

# Matrices based on how you are defining them
( gm <- round(focalMat(r, 350*0.5, "Gauss"),10) )
( gse <- round(gaussian.kernel(sigma = 0.5 * 350),10) )

# If raster.gaussian.smooth matrix is the same as focalMat with focal
spatialEco::raster.gaussian.smooth(r, sigma = 0.5 * 350, n=nrow(gm))
focal(r, w = gm)

Best, Jeff