23

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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Doug
  • 617
  • 1
  • 6
  • 10

3 Answers3

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.

Glorfindel
  • 1,096
  • 2
  • 9
  • 14
om_henners
  • 15,642
  • 2
  • 46
  • 86
  • Thanks for the help! I am giving the Scipy.interpolate.griddata approach a whirl. I will post back my results. – Doug Dec 02 '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 09 '14 at 00: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.

user2856
  • 65,736
  • 6
  • 115
  • 196
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.

rosskush
  • 176
  • 2
  • 6