Smoothing/interpolating raster in Python using GDAL?

Question

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:

[Raster from using gdal.rasterizelayer]

What I want for a final product:

enter image description here

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 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.

• Doug: Thanks for the help! I am giving the Scipy.interpolate.griddata approach a whirl. I will post back my results.
• Doug: I apologize for taking so long to get back to this post. The answer above is basically what I did to solve my problem. I used the Scipy interpolate library to fill in those nodata spaces and then wrote it back to the rasterband. Thanks for the help guys!
• om_henners: @Doug No worries - happy to healp!
• Miro: How fast is this solution? Can it be used for 10k x 10k grid where only every 100x100 is known value? I tried gdal_fillnodata which is incredibly fast compared to any interpolation but it doesn't work well for too sparse points. At this moment I am using triangulation from Saga but it is very slow for medium arrays and fail with big ones.

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.

Answer 3

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.

Answer 4

The simplest way as I know so far is, use QGIS Raster>Projections>Warp(reproject), then you can specify the "resampling method to use" as "Cubic Spline", and output resolution in smaller units, which could interpolate the raster data very smoothly.

See the result below:

enter image description here

• dimfalk: Since the question does not ask for a QGIS solution, but the idea is good, maybe better refer directly to gdalwarp -r cubicspline? c.f. gdal.org/programs/gdalwarp.html