OSGeo / grass

GRASS GIS - free and open-source geospatial processing engine
https://grass.osgeo.org
Other
848 stars 308 forks source link

[Feat] Provide better explanation for radius parameter for r.resamp.filter #2569

Open nikosGeography opened 2 years ago

nikosGeography commented 2 years ago

I want to resample a raster from 15m to 460m using a Gaussian filter.

The goal

I am having a coarse image which I want to downscale. I also have a fine resolution band to assist the downscaling. The downscaling method I am using is called geographically weighted area-to-point regression Kriging (GWATPRK). The method consists of two steps:

  1. GWR and,
  2. ATPK on the GWR's residuals.

In order to perform GWR using raster data, those needs to have the same pixel size. This means that, my fine resolution image needs to be upscaled to match the spatial resolution of the coarse band. This upscaling of the fine band needs to be done using a Gaussian kernel (i.e., the PSF). I have found that GRASS GIS has a tool called r.resamp.filter. I am trying to run the function but I am getting the following error(s):

  1. ERROR: Differing number of values for filter= and [xy_]radius=

This error occurs when I use two filter kernels (e.g., gauss + box, or gauss + bartlett). I am using two kernels because according to the Manual:

Kernels with infinite extent (Gauss, normal, sinc, Hann, Hamming, Blackman) must be used in conjunction with a finite windowing function (box, Bartlett, Hermite, Lanczos).

Doesn't matter what numbers I put in the Filter radius or Filter radius (horizontal) and Filter radius (vertical) (see image below), I tested A LOT of numbers. 8JFqN

  1. ERROR: At least one filter must be finite

This error occurs when I use one filter kernel, the gaussian (I am interested in applying a Gaussian filter, because I want to model the point spread function during downscaling satellite imagery).

The steps I followed were:

r.external to import the raster g.region where I set the region using my original fine resolution image BUT in Resolution tab I changed the 2D resolution into 460 r.resamp.filter and the errors I mentioned Ultimately, I want to apply a Gaussian filter with sigma (std) = 0.5 to my image while I am upscaling it.

Here the image I am using

marisn commented 2 years ago

1st error indicates that you only provided one value for radius parameter. Specify two comma separated values – one for each of filters. Use same value in case if radius should be equal for both filters.

nikosGeography commented 2 years ago

Can you please provide some example or even better, could you use the image I provided? I tried what you've told me and the output raster has pixel values ranging from -Inf to +Inf. Some of my tries included:

  1. 0.2,0.2 all the way to 0.9,0.9 with step 0.05
  2. I change the filters gaussian + any of the box, barlett, Hermite, Lanczos.
petrasovaa commented 2 years ago
  1. 0.2,0.2 all the way to 0.9,0.9 with step 0.05

I didn't test but I would try different radius for the box and the gaussian filter. Unfortunately the documentation is unclear what radius means here.

nikosGeography commented 2 years ago

I am guessing it means the width of the psf (i.e., the std) but again, as you said it's unclear what radius means. I'll try different numbers and I'll post the results here.

marisn commented 2 years ago

Looking at the code shows that radius is used to calculate raster rows and columns. Thus any small value not resulting in an increment of row/col count doesn't make a lot of sense. Try to increase radius to 15 or even 500 to see real results. I do agree that the documentation could be improved to explain better units (same as for coordinate system) and magnitude of radius + give some examples on usage. Unfortunately this is not my topic so I can not propose changes.

marisn commented 2 years ago

@metzm could you help with this one?

metzm commented 2 years ago

With r.resamp.filter, the radius is in map units. At least one radius must be larger than the resolution of the source raster, otherwise a filter can not be applied. I use quite often r.resample.filter filter=gauss,box radius=<2 x resolution>,<3 x resolution> with good results.

r.resamp.filter like most other r.resamp.* modules is suitable for downscaling (increasing the resolution), but not for upscaling.

In the case here, a raster needs to be upscaled, changing the resolution from 15m to 460m, that means several input cells will fall into one output cell. The appropriate module here is r.resamp.stats to aggregate several input values into one output value.

In the case here I would use r.resamp.filter on the GWR's residuals to increase the resolution from 460m to 15m.

nikosGeography commented 2 years ago

@metzm I'll try to set larger radius. As for the aggregation method the reason for choosing the r.resamp.filter is because it gives the opportunity to change the pixel size using a Gaussian filter. My goal is to model the point spread function (PSF).

To put it differently, I am trying to simulate coarse data as though they were measured with a coarse PSF.

To do this I need to apply a transfer function (TF; e.g., Gaussian) to the fine data, but with a very large width. This produces the coarse data (460m pixel size).

metzm commented 2 years ago

With a Gaussian filter, the points in the center get more weight than the points at the edges. I think this is rather arbitrary when simulating a coarse PSF, you could just as well randomly pick one of the 15m cells falling into a target 460m cell.

When using r.resamp.filter to upscale from 15m to 460m, you could try filter=gauss,box radius=250,250 to cover approximately one target cell with the filter.

nikosGeography commented 2 years ago

@metzm the filter radius values were correct. I got a coarser resolution raster of ~460m pixel size. Just one more question before I close the issue. Why the value 250 produces the 460m pixel size? How the math works? Just a remainder that my initial raster has a 15m pixel size.

marisn commented 2 years ago

Just one more question before I close the issue.

Do not close this issue as in its new form it is still valid – essence of this discussion should be added to the documentation.

myxBeni commented 1 year ago

Do not close this issue as in its new form it is still valid – essence of this discussion should be added to the documentation.

Totally agree. Crucial to me is (credits to @metzm )

Helpful for me would be to know that a Box-filter in r.resamp.filter and radius 3 equals a Box of about 6 by 6 (radius of 3 units equals 6 units box length).