Sunday 30 April 2017

raster - Python - gdal.Grid() correct use


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;




newpoints.csv
wkbPoint




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



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)



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...