Smoothing/interpolating raster in Python using GDAL?

I'd take a look at NumPy and Scipy - there's a good example of interpolating point data in the SciPy Cookbook using the scipy.interpolate.griddata function. Obviously this requires that you have the data in a numpy array;

  • Using the GDAL python bindings you can read your data into Python using gdal.Dataset.ReadAsArray() for a raster.
  • With OGR you would loop through the feature layer and extracting point data from the shapefile (or better yet, write the shapefile to a CSV using GEOMETRY=AS_XYZ [see the OGR CSV file format] and read the csv into Python).

Once you've got a gridded output you can then use GDAL to write the resulting numpy array to a raster.

Lastly, if you don't have any luck with the Scipy interpolate library, you could always try scipy.ndimage as well.


Have a look at the GDAL gridding API. I don't know if that is exposed in the Python bindings, but if not, you call call the gdal_grid utility via the subprocess module.

GDAL grid API only uses Inverse Distance Weighting, Moving Average and Nearest Neighbour, it doesn't implement splines. Another option is to use Scipy.


A bit old to this thread but I have written a simple module that uses the KNN algorithm from sklearn called skspatial.

https://github.com/rosskush/skspatial

You can import a shapefile using geopandas and select a column and it will interpolate a surface that can be exported to a raster. It's very basic and probably not the best way to do it, but it keeps everything pure python at least.