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?

4 Answers 4


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.

  • Thanks for the help! I am giving the Scipy.interpolate.griddata approach a whirl. I will post back my results.
    – Doug
    Commented Dec 2, 2011 at 16:57
  • 1
    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!
    – Doug
    Commented Mar 13, 2012 at 16:42
  • @Doug No worries - happy to healp!
    – om_henners
    Commented Mar 13, 2012 at 17:06
  • 1
    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.
    – Miro
    Commented Sep 9, 2014 at 0:15

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.


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.


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

  • 1
    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
    – dimfalk
    Commented Apr 9 at 14:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.