21

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?

Improve this question
19

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.

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
  • 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 Mar 13 '12 at 16:42
  • @Doug No worries - happy to healp! – om_henners Mar 13 '12 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 Sep 9 '14 at 0:15
12

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.

Improve this answer
2

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.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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