Python GDAL – Smoothing and Interpolating Raster

gdalinterpolationpythonraster

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:

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?

Best Answer

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.

Related Solutions

[GIS] Using signed Bytes with GDAL

If you want to use signed 8bit integers for GeoTIFFs use datatype Byte together with creation options PIXELTYPE=SIGNEDBYTE

https://gdal.org/drivers/raster/gtiff.html

gdal_translate -ot Byte -a_nodata -1 -co PIXELTYPE=SIGNEDBYTE input.tif output.tif

[GIS] Extract raster values from point using GDAL

Working from your provided example (which is the standard way of doing this), you can collect all the info in a python list.

Let's say your point layer has a unique ID field (if it doesn't, create one, as it really should). For this example, let's call it "id_points". You complement the code in your link with:

li_values = list()
for feat in lyr:
    geom = feat.GetGeometryRef()
    feat_id = feat.GetField('id_points')
    mx, my = geom.GetX(), geom.GetY()

    px = int((mx - gt[0]) / gt[1])
    py = int((my - gt[3]) / gt[5])

    intval = rb.ReadAsArray(px, py, 1, 1)
    li_values.append([feat_id, intval[0]])

Original code attribution

This gives you a list of feature IDs and their associated raster values. You can then save it in a CSV (for example):

import csv

with open(r'csv/file/path.csv', 'wb') as csvfile:
    wr = csv.writer(csvfile)
    wr.writerows(li_values)

This will give you an output in the form of a table, which you can then open anywhere you want.

Best Answer

Related Solutions

[GIS] Using signed Bytes with GDAL

[GIS] Extract raster values from point using GDAL

Related Question