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

[nikos_geo]

I have a fine resolution raster (100m pixel size) and I want to downsample it (aggregate) to 460m pixel size. The downsampling should be done using a Gaussian filter (local function). The scale (zoom) factor is 4.6, so non-integer. The sigma (std) parameter of the Gaussian filter should be expressed in cells (pixels).

To put it differently I am trying to simulate coarse data as though they were measured with a coarse point spread function (PSF). The PSF is assumed to be a Gaussian filter and the units should be in pixels as well.

So, in theory, for each coarse pixel I need to go to its center and calculate the weights (from the PSF/Gaussian filter) needed for each fine pixel surround it. To do this I need to apply a transfer function (TF; e.g., Gaussian filter) to the fine data, but with a very large width. But that's in theory and I couldn't find any package in R that does this. The OpenImageR package has a function called down_sample_image but the units are not in pixels.

The function aggregate from the terra package does not allow the fact argument to be non-integer, so I thought to create a template raster (a raster with no pixel values) from a coarse resolution raster at 460m spatial scale that I have.

All in all, I need to create a custom function which downsamples a raster using a Gaussian filter (units in pixels). I am stuck with this problem for some time now so any help would be nice.

So far I created the template raster:

library(terra)

fr = rast("path/tirs.tif") # raster to be downsampled
cr = rast("path/ntl.tif") # raster to be used as template

# create an empty raster with the same dim, extent and crs as the coarse res raster
(er <- rast(ext(cr), resolution=res(cr)))  # er = empty raster
crs(er) <- crs(cr)

The raster layers I am using:

fr = rast(ncols=216, nrows=417, nlyrs=1, xmin=582700, xmax=604300, ymin=1005700, ymax=1047400, names=c('B10_median'), crs='EPSG:7767')

cr = rast(ncols=48, nrows=91, nlyrs=1, xmin=582360, xmax=604440, ymin=1005560, ymax=1047420, names=c('avg_rad'), crs='EPSG:7767')

Any recommendations on how I can proceed?

[nikos_geo]

Because the scale factor is 4.6, I needed to smooth the image and then resample the smoothed image using nearest neighbor.