Downsample (aggregate) raster by a non-integer factor, using a Gaussian filter kernel

Question

The task is to downsample (aggregate) a raster from 100m pixel size to 460m. The aggregation should be performed using a Gaussian filter. To better understand the task, I am following the paper ‘The effect of the point spread function on downscaling continua’. Most of the paper is irrelevant to my task, all I care from this paper is this one step:

(Note: by upscaling the authors mean aggregation)

One of the authors is my supervisor and I asked him if I can blur my fine resolution raster and then aggregate it using a common interpolation algorithm (nearest neighbor, bilinear etc). This is not the way to go. The aggregation should be done using a Gaussian kernel filter (the point spread function is assumed to be Gaussian). Also, If I blur and then resample is like I add extra PSF effect apart from what my image already has.

There is a post on Reddit, where a person suggests (without sharing how to do it) that this a common computer vision task. I share his suggestion:

My supervisor told me that the way I should create the aggregated raster is by applying a gaussian kernel filter to the fine data, but with a very large width. This large width I think it determines the output pixel size (which as I said I want it to be 460m).I say that based on this post.

According to my supervisor: For each new coarse pixel go to its center and calculate the weights (from the PSF) needed for each fine pixel surrounding it (PSF = point spread function = Gaussian filter).

You can download my data from here, or if you use R:

fr = rast(ncols=108, nrows=203, nlyrs=1, xmin=583400, xmax=594200, ymin=1005700, ymax=1026000, names=c('B10_median'), crs='EPSG:7767') # fine resolution raster

cr = rast(ncols=23, nrows=43, nlyrs=1, xmin=583280, xmax=593860, ymin=1006020, ymax=1025800, names=c('coarse_image'), crs='EPSG:7767') # template (coarse) resolution raster

I shared a template raster because I want my aggregated raster to match the resolution (ncols and nrows) of a coarse resolution raster that I will use later in my analysis.

Lastly, the units of σ (sigma) and the Gaussian should be in pixels.

Any recommendations on how to proceed? Preferably in R but it doesn’t really matter.

Answer 1

This would sound like a relatively benign rational resampling – sample up by a factor of 5, down by a factor of 23.

Since you would need an anti-imaging filter after the upsampling, and an anti-aliasing filter before downsampling, just keep the "stricter one", the 1/23-band low-pass for anti-aliasing.

Choosing a Gaussian filter for that gives you Gaussian "interpolation"¹.

The downsampled raster should have exactly the same dimensions as the coarse resolution raster because after the downsampling I want to perform regression between those two raster layers.

That is a paradox. The whole point of downsampling is to change the dimensions. So, either you resample, and change the "raster" of the image, or you don't.

---

¹ quotation marks courtesy of the fact that interpolation requires the original pixels, should they happen to end up on an integer pixel coordinate after interpolation, to be precisely preserved; which a Gaussian could never achieve as filter, since it has no periodic zeros in spatial domain that would allow for an exact solution to the interpolation problem in the infinite case.

Answer 2

So the solution is the one that @Cris Luengo proposed, with the only difference that I used nearest neighbor resampling. That is because the target resolution is 460m and my input is 100m, I cannot avoid using some form of interpolation technique (when aggregating isn't called interpolation but that's besides the point). So what I did is I firstly convolved the image using a Gaussian filter and then I used nearest neighbor resampling to downsample my raster to 460m. The code in R is this:

library(raster)
library(gridkernel)
library(gridprocess)

fr = raster("path/fine_resolution_image.tif") # image to be blurred
cr = raster("path/coarse_resolution_image.tif") # another image at the target resolution. to be used for the resampling of the blurred fr image

g = as.grid(fr)
smoothed = gaussiansmooth(g, sd = 0.3 * 460, max.r = 100) # units in pixels
r <- raster(smoothed)
resample(fr, cr, method="ngb", filename="path/blurred_resampled.tif")

The other solution, would be to downsample my coarse resolution raster from 460m to 500m, 'offline', and then use the aggregate function to downsample my 100m fine resolution raster to 500m, using a custom function for the fact argument. The custom function would be the Gaussian filter.