Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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:

[Raster from using gdal.rasterizelayer]

What I want for a final product:

[Wanted raster, interpolated]

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!

share|improve this question
up vote 13 down vote accepted

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.

share|improve this answer
Thanks for the help! I am giving the Scipy.interpolate.griddata approach a whirl. I will post back my results. – Doug Dec 2 '11 at 16:57
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 Mar 13 '12 at 16:42
@Doug No worries - happy to healp! – om_henners Mar 13 '12 at 17:06
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 Sep 9 '14 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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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