How to smooth/interpolate a raster in python using GDAL?

Question

This is my first post so please feel free to critique where necessary.

Background: I am developing in Python and using GDAL from OSGEO to manipulate and interact with rasters and shapefiles.

Problem: I want to take a shapefile that has point features and interpolate it into a surface raster. Right now I am using the method 'RasterizeLayer' which burns a value from the point feature into the raster (which is set with all nodata values) but leaves all untouched pixels as a 'nodata' value. I am therefore left with a checkerboard type raster.

What I have after using RasterizeLayer:

What I want for a final product:

I believe the function I am looking for is known as 'Spline_sa()' from the arcgisscripting import. Does GDAL have a similar function, or is there a different method to get my desired output? Thanks!

Answer 1

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.

Answer 2

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.