4

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

4

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)

4
  • that's surprisingly simple, will try it thanks! I just interpreted it as; in case the point is NOT in a triangle, only then the radius is taken into account. Sep 5 '17 at 15:00
  • still no success.. Tried multiple radii, all resulting in the same raster. INVDIST is working fine, together with its parameters. Should I use a CRS by definition? Although it also works without on the invdist? checkout: gist.github.com/willemvanopstal/… Sep 5 '17 at 15:12
  • 1
    I've just tried 'linear and 'linear:radius=50' and they work in both cases. Two observations: first, the GeometryType should be wkbPoint25D if you pass also the z coordinate, then defining also the LayerSRS is always a good design choice. What GDAL version are you using? Sep 5 '17 at 17:12
  • Revisiting the pointType resolves the problem.. confuses me, but hé it works. Thanks for investigating! Sep 5 '17 at 18:06

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.