0

I am trying to aggregate a raster using a custom function and the terra package. For the custom function I am using a Gaussian smoothing. When I try to run the aggregate function I am getting this error: Error in get(as.character(FUN), mode = "function", envir = envir) : object 'fun' of mode 'function' was not found. What am I doing wrong?

library(terra)

fr = rast("path/fr.tif")

matg = focalMat(fr, 
                0.5 * 500, 
                "Gauss")

ra = aggregate(fr, 
               fact = 5, 
               fun = matg)

writeRaster(ra, 
            "path/fr_aggr.tif", 
            overwrite = TRUE)

My raster:

fr = rast(ncols=222, nrows=205, nlyrs=1, xmin=4831500, xmax=4853700, ymin=3856600, ymax=3877100, names=c('B10_median'), crs='EPSG:7755')

After I tried the code from the answer below, I am getting this error: Error: [focal] test failed

Here is what I tried:

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

lin.conv <- function (y,K) {
  if(inherits(y, "matrix"))
    y <- as.vector(y)  
  K <- K / sum(K)
  X <- as.matrix(Matrix::bandSparse(length(y), 
                                    k = seq(-(length(K)-1),0,1), 
                                    diag = t(replicate(length(y), rev(K))), 
                                    symm=FALSE))
  out <- X %*% as.matrix(y, ncol=1)
  return(as.vector(out))
}

fconv <- function(x, sdv=0.2*500) {
  n <- floor((sqrt(1 + 8 * length(x)) - 1)/2)-1
  y <- gaussian.kernel(sigma=sdv, s = n)
  lcv <- lin.conv(x,y)
  return( lcv[ceiling(length(lcv)/2)] )
}

rgs <- focal(r, matrix(1,5,5), fun=fconv) #here is the error
2
  • The focalMat is creating a Gaussian Kernel (matrix), not a function (eg., mean) which is what aggregate is expecting. I think that you want to apply your kernel function to focal and then aggregate the results to something approximating your kernel size. You could also write a function that applies the Gaussian weights to the pixels being aggregated but, it would have to match the number of pixels being considered but the specified fact and you will need to figure out what is happening within the function, is it operating on a vector or matrix, and tailor you function accordingly. Commented Nov 11, 2022 at 17:35
  • That's exactly what I want to achieve, i.e., to write a function that applies the Gaussian weights to the pixels being aggregated. Could you recommend on how to proceed?
    – Nikos
    Commented Nov 11, 2022 at 18:28

1 Answer 1

1

I would do this at the source resolution via a focal function and then aggregate using the mean (it will already be normal). I would make the argument that, from a smoothing perspective, one could achieve good results by skipping Gaussian decomposition and just apply a smoothing spline or local polynomial regression. This would be the nonparametric analog. Anyway, here is an approach that uses a linear matrix convolution for the Gaussian decomposition.

First, here are the needed packages, a function to perform the convolution and one to apply it within a focal function.

library(terra)
library(spatialEco)

lin.conv <- function (y,K) {
  if(inherits(y, "matrix"))
    y <- as.vector(y)  
    K <- K / sum(K)
    X <- as.matrix(Matrix::bandSparse(length(y), 
           k = seq(-(length(K)-1),0,1), 
           diag = t(replicate(length(y), rev(K))), 
           symm=FALSE))
    out <- X %*% as.matrix(y, ncol=1)
  return(as.vector(out))
}

fconv <- function(x, sdv=2) {
  n <- floor((sqrt(1 + 8 * length(x)) - 1)/2)-1
  y <- gaussian.kernel(sigma=sdv, s = n)
  lcv <- lin.conv(x,y)
  return( lcv[ceiling(length(lcv)/2)] )
}

Before we apply it within focal, let's take a look at what it is doing first. Here we create a 5x5 matrix (emulating what would be a single focal operation) and a Gaussian kernel of the same dimensions. We then apply the linear convolution and plot the results. The black point represents what would be the estimated focal value.

( x <- matrix(runif(25),ncol=5,nrow=5) )
y <- gaussian.kernel(sigma=2.5, s = 5)
lcv <- lin.conv(x,y)
( v <- lcv[ceiling(length(lcv)/2)] )
  dev.new(width=11,height=8.5)
    plot(as.vector(x), type="l", lty=3)
      lines(1:length(x), lcv, col="red")
      points(x=which(lcv %in% v), y=v, 
             pch=20, cex=2)

smoothed distribution

Now that we have an inkling that this is, in fact, applying a smoothing function we can apply it to a raster and plot the results.

r <- rast(system.file("ex/elev.tif", package="terra"))
  rgs <- focal(r, matrix(1,5,5), fun=fconv)
    plot(c(r, rgs))

smoothed elevation

Then, you can aggregate the smoothed results. Since you have applied a Gaussian smoothing, the aggregated pixels will be normal so, the mean should be a stable metric for aggregating to a coarser resolution. Since it is sort of functioning as a focal operator you can pass our focal decomposition function directly to aggregate as well (which is what I believe you are after).

agg.rgs <- aggregate(rgs, fact = 3, fun ="mean") # OR
agg.rgs <- aggregate(r, fact = 3, fun =fconv)
  plot(agg.rgs)

smoothed aggregated results

5
  • The reason I chose Gaussian smoothing over other smoothers (splines etc) is because I am trying to simulate coarse data with a point spread function, which I assume is Gaussian (according to the literature).
    – Nikos
    Commented Nov 11, 2022 at 20:32
  • Generally speaking, yes Gaussian is assumed however, this is not a given pubmed.ncbi.nlm.nih.gov/21164793 You are correct, for your application, other smoothing functions such as splines would not be appropriate. Commented Nov 11, 2022 at 22:25
  • I edited my question because I'm getting this error when I try your code: error: Error: [focal] test failed. Any idea why?
    – Nikos
    Commented Nov 12, 2022 at 0:10
  • @Nikos I am not getting this error here but am elsewhere and am trying to track down what is going on. I think it may be a bug in the test function used to evaluate if a passed function can run under focal. I am able to recreate very inconstant results where I can run a something outside of a function but, exactly the same code will not run inside a function and I get this error. I reran the example and am not getting it here, are you or only on your data? Commented Nov 12, 2022 at 0:23
  • I tested it with the raster I shared plus one other image. I got the same error.
    – Nikos
    Commented Nov 12, 2022 at 11:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.