0
$\begingroup$

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:

aggregation

(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:

enter image description here

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.

width of the psf

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.

$\endgroup$
6
  • $\begingroup$ We need more detail. Please edit your question with additional information. This being Stackexchange, it's best to just make it look like you meant it that way all along (i.e., don't just tack on a section at the end titled "edit" -- fold the edits into your question). "To do this I need to apply a transfer function ... but with a very large width." What is "very large" in this context? Give us the numbers. "The OpenImageR ... down_sample_image". Provide a link, please -- those of us who don't know R but do know image processing can interpret this for you -- with a link. $\endgroup$
    – TimWescott
    Nov 1, 2022 at 15:06
  • 1
    $\begingroup$ Just apply a Gaussian filter, then resample using interpolation (preferably cubic spline or Lanczos) at the desired output grid coordinates. This is identical to sampling the input with a Gaussian centered at each of the output grid coordinates. $\endgroup$ Nov 1, 2022 at 20:33
  • $\begingroup$ So this is what I did few months ago. I applied a Gaussian filter and then I resampled the image but it is wrong. Why? Because the Gaussian filter is the resampling. The change in the pixel size is determined by the width of the filter. Also, if you blur and then resample the image it's like you add aliasing to the image whereas if you apply a Gaussian kernel to downsample an image you remove the aliasing. $\endgroup$
    – Nikos
    Nov 1, 2022 at 20:44
  • $\begingroup$ I edited my question and I explained why blurring and then resampling is not the right way. $\endgroup$
    – Nikos
    Nov 1, 2022 at 20:48
  • $\begingroup$ I’m not sure what you plotted, but the method I suggested and what you describe in your answer should produce the same results. Maybe you’re not using the right boundary extension when applying the Gaussian filter, maybe there’s some other issue in your code. $\endgroup$ Nov 1, 2022 at 21:19

2 Answers 2

0
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ maybe I didn't make myself clear enough. So After the downsampling I should get a coarse resolution raster which I intend to use it to perform linear regression using another coarse resolution raster. In order to perform regression those two rasters should have exactly the same dimensions (rows and columns) and thus pixel size. That's why I shared two rasters. In case, the Gaussian downsampling doesn't produce a coarse raster with the same dimensions (and pixel size) as my template raster, if there was a way to use the template raster to assist the downsampling, if that makes sense. $\endgroup$
    – Nikos
    Nov 1, 2022 at 16:24
  • $\begingroup$ no, sorry, doesn't make sense. the scale factor fully defines the size. $\endgroup$ Nov 1, 2022 at 17:07
  • $\begingroup$ Yes I don't disagree with that, the scale factor defines the width of the Gaussian. I'm just sceptical about the downsampled image resolution, because the scale factor is not an integer number. $\endgroup$
    – Nikos
    Nov 1, 2022 at 19:00
  • $\begingroup$ You said you needed to downsample by a factor of 4.6, based on the pixel size before and after. That defines what your raster is. You cannot "zoom out" and "keep the same size". $\endgroup$ Nov 1, 2022 at 19:05
  • $\begingroup$ Maybe I explain it poorly. Yes the scale factor is 4.6 because my target pixel size is 460m and my input is 100m. I think I found a way to create such function (actually some one helped to find a way but that's irrelevant). Should I post an answer using the way I found and continue from there? Because I'm not sure it's going to work $\endgroup$
    – Nikos
    Nov 1, 2022 at 19:12
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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