Python – Correct Use of gdal.Grid() Function

gdalpythonrastervrt

Is there any good, complete API for the GDAL Python bindings?

I am interpolating and rasterizing points with elevation data using the gdal.Grid() method;

output = gdal.Grid('outcome.tif','newpoints.vrt')

In its most basic sense, it's working perfectly. But now, I'd like to try the linear approach;

output = gdal.Grid('outcome.tif','newpoints.vrt', algorithm = 'linear:radius=0')

That does generate a new .tif file (512kb), but is just one big black image. Also, not all options seem to be accepted.. for the invdist:power=4 for example it also results in a useless image. What is going wrong, am I missing something? Below is the .vrt I'm using;

<OGRVRTDataSource>
<OGRVRTLayer name="newpoints">
    <SrcDataSource>newpoints.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <GeometryField encoding="PointFromColumns" x="long" y="lat" z="depth"/>
</OGRVRTLayer>

EDIT
As a next example I used the gdal.GridOptions(), putting the parameters in that one. It is resulting in a correct file, but still ignores the options. How to actually pass the GridOptions onto the gdal.Grid() ?

edit
I found this question also; GDAL grid produces NoData
Is the linear interpolation just a bug in GDAL?

See my example workfiles;
https://gist.github.com/willemvanopstal/1eae4999fcfa62e89cca39a350b504d8

Best Answer

More info about gdal.Grid() and gdal.GridOptions() are available in the GDAL/OGR Python API. Instead, here's the available options of the interpolation algorithms. About the linear approach, this line (and more in detail the algorithm options 'linear:radius=0'):

output = gdal.Grid('outcome.tif','newpoints.vrt', algorithm = 'linear:radius=0')

returns always NoData because:

radius: In case the point to be interpolated does not fit into a triangle of the Delaunay triangulation, use that maximum distance to search a nearest neighbour, or use nodata otherwise. If set to -1, the search distance is infinite. If set to 0, nodata value will be always used. Default is -1.

(Source: http://www.gdal.org/gdal_grid.html#gdal_grid_algorithms)

So, I'd try to use a better value of radius.

Note that the "values to interpolate will be read from Z value of geometry record." otherwise you will need to specify the zfield option:

zfield --- Identifies an attribute field on the features to be used to get a Z value from. [...]

(Source: http://gdal.org/python/osgeo.gdal-module.html#GridOptions)

Related Question