How to smooth / "reinterpolate" a raster with GDAL?

by Pablo   Last Updated July 05, 2018 22:22 PM

Short:
I want to change the resolution of a raster and smooth the gray ramp like shown in the images bellow. The preference is to use GDAL, PIL or Numpy.

Description:
I'm kriging points into rasters with an output resolution of 20 meters with the High Performance Geostatistical Library. I don't want to change the output resolution because the interpolation time increase exponentially.
With this resolution the output image is ugly (pixelated and aliased). I don't know if it is conceptually correct but I want the image to be smoother like in the example bellow. It's something like 'reinterpolating' the image into a better resolution one. I'm using python so my preferences are GDAL, Python Imaging Library or Numpy. The answer could be theoretical, like pointing out the algorithm name or the concept of this kind of operation.

Source:
enter image description here

Destination:
enter image description here

EDIT Results with gdalwarp cubic spline:
enter image description here

Tags : python raster gdal


Answers 4


1) The hard way: With a bit of coding it's (relatively) easy to implement bilinear interpolation to accomplish decent resampling.

2) The easy way: use GDAL as explained in this previous GISSE post, but in reverse (decreasing the pixel size).

WolfOdrade
WolfOdrade
July 31, 2012 15:50 PM

To smooth out the variations, you need a low-pass filter. You could write your own using GDAL, or there's one using GRASS. I haven't tried it, but here's a guide http://wiki.awf.forst.uni-goettingen.de/wiki/index.php/Exercise_31

You may want to up-sample your raster first before applying the low-pass filter to get better resolution output.

spatialthoughts
spatialthoughts
July 31, 2012 15:53 PM

Use GDALReprojectImage, which is exposed in Python:

from osgeo import gdal
help(gdal.ReprojectImage)

For the smooth interpolation, use bilinear or cubic methods. This function is awkward, since it doesn't take keyword arguments, thus you need to find the position:

gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Bilinear)

Probably the tricky part is setting up dst_ds, which needs to have a geotransform similar to src_ds, but with modified cell sizes.

Mike T
Mike T
July 31, 2012 19:20 PM

you can use a rank/median filter with radius=5, i.e kernel size size=11, (for each rgb channels).

yildirim
yildirim
July 05, 2018 21:50 PM

Related Questions





Using signed Bytes with GDAL

Updated January 23, 2018 11:22 AM