I'm trying to use GDAL in a python script to do an inverse distance interpolation. Within the C documentation API I found GDALGridInverseDistanceToAPower function but I don't know how/if I can call that from Python.

Does somebody know how to perform this analysis or call this GDAL function in python?


Yes, you can call this from a Python script. But you don't directly call the low-level C API.

First, take a look at the GDAL Grid Tutorial for background info. From the Python library, the relevant function is gdal.Grid(destName, srcDS, **kwargs). You can see some examples of how it's used in test_gdal_grid_lib.py (from the test suite). Or a made-up example:

from osgeo import gdal

ds = gdal.Grid('result.tif', 'points.shp', format='GTiff',
               outputBounds=[0.0, 0.0, 100.0, 100.0],
               width=10, height=10, outputType=gdal.GDT_Float32,

The full set of spatial interpolation options for invdist (Inverse Distance to a Power) are documented here, and are passed as a colon-separated string to the algorithm keyword. And there are other algorithms than invdist available in a similar manner.

You can also construct many of the other gridding options with gdal.GridOptions, which has all of the available keyword arguments, i.e.:

>>> help(gdal.GridOptions)
Help on function GridOptions in module osgeo.gdal:

GridOptions(options=[], format=None, outputType=0, width=0, height=0, creationOptions=None, outputBounds=None, outputSRS=None, noData=None, algorithm=None, layers=None, SQLStatement=None, where=None, spatFilter=None, zfield=None, z_increase=None, z_multiply=None, callback=None, callback_data=None)
    Create a GridOptions() object that can be passed to gdal.Grid()
    Keyword arguments are :
      options --- can be be an array of strings, a string or let empty and filled from other keywords.
      format --- output format ("GTiff", etc...)
      outputType --- output type (gdal.GDT_Byte, etc...)
      width --- width of the output raster in pixel
      height --- height of the output raster in pixel
      creationOptions --- list of creation options
      outputBounds --- assigned output bounds: [ulx, uly, lrx, lry]
      outputSRS --- assigned output SRS
      noData --- nodata value
      algorithm --- e.g "invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0"
      layers --- list of layers to convert
      SQLStatement --- SQL statement to apply to the source dataset
      where --- WHERE clause to apply to source layer(s)
      spatFilter --- spatial filter as (minX, minY, maxX, maxY) bounding box
      zfield --- Identifies an attribute field on the features to be used to get a Z value from. This value overrides Z value read from feature geometry record.
      z_increase --- Addition to the attribute field on the features to be used to get a Z value from. The addition should be the same unit as Z value. The result value will be Z value + Z increase value. The default value is 0.
      z_multiply - Multiplication ratio for Z field. This can be used for shift from e.g. foot to meters or from  elevation to deep. The result value will be (Z value + Z increase value) * Z multiply value.  The default value is 1.
      callback --- callback method
      callback_data --- user data for callback

It would seem someone has created a pythonic method into running the gdal_grid tool. You can read about it and download a zipped file containing the scripts from the following link:

Creating a grid from scattered data using inverse of the distance with python (gdal_grid approach)

  • 1
    I'll try it now! – Tim C. Feb 27 at 13:28

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.