Smoothing/interpolating raster in Python using GDAL?

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

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?

Answer

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.