[gdal-dev] downsampling geotiff with a low-pass filter

[Kay F. Jahnke]

On 16.11.18 16:38, JDA wrote:
> 
> There’s a wide array of smoothing options available If you’re willing to 
> work in python. Based on https://gis.stackexchange.com/a/10467, the 
> basic idea is to load the raster into a numpy array and then convolve it 
> with either a kernel of arbitrary size.
> 
> I’ve written a program implementing that idea with both a Gaussian and a 
> mean kernel here:
> https://github.com/jacobdadams/general_scripts/blob/master/raster_chunk_processing.py. 
> It only smooths the data without downsampling, but there may be some 
> python functions that downsample as well.
> 
> There’s a lot of code in there for dealing with massive rasters by 
> processing them in chunks in parallel, but the blur_mean and blur_gauss 
> functions are where the smoothing is done. I’ve written an installation 
> and usage guide at 
> https://gisjake.blogspot.com/2018/10/rasterchunkprocessingpy-installation.html?m=1. 

I cloned your git repo and fetched the additional python modules my 
system did not have with the packet management system, namely astropy 
and skimage. At first I did this for python2, but no luck here, running 
the sript would throw an exception:

Processing chunks...

--- Error ---
__exit__


Traceback (most recent call last):
   File "/home/kfj/src/general_scripts/raster_chunk_processing.py", line 
1074, in <module>
     num_threads, verbose)
   File "/home/kfj/src/general_scripts/raster_chunk_processing.py", line 
961, in ParallelRCP
     maxtasksperchild=10
AttributeError: __exit__

So I picked up the packages in python 3, and got it to work. I used this 
command line:

python3 raster_chunk_processing.py -p 1 -m blur_gauss -o 25 \
         -s 1500 --verbose -r 15 -d 1.5 original.vrt jacob.tif

where original.vrt had a sample population of two DEM tiles on the 
border of my target area, having no-data values outside the region's 
territory. When inspecting the result in QGIS, it looked like it had 
done everything I wanted it to do:

- the VRT file was taken as valid input
- the no-data value coded in the VRT file was properly honoured
- there were no problems at tile borders
- there was a clean edge, valid data inside and no-data outside
- the applied blur looked like what I'd expect from the parameters
- the resulting contours were smooth and as intended

So, compared to the attempts with a kernel-filtered source in a VRT file 
which I outlined in my previous post, I'd say that using your script is 
the winner, because the margin treatment is what one might expect and 
the parametrization is straightforward. I'll just have to tweak the 
parameters to get precisely the amount of smoothing I want. I'll store 
the smoothed data as GeoTIFF in the original resolution and then proceed 
to downsample from that, and my problem should be solved.

Thank you for your reply and for sharing your scripts!

Kay