GDAL rasters from irregular point or ASCII data

The following is an excerpt from the book: Geospatial Power Tools – Open Source GDAL/OGR Command Line Tools by Tyler Mitchell.  The book is a comprehensive manual as well as a guide to typical data processing workflows, such as the following short sample… 

A common question that seems to come up on the GDAL discussion lists seems to be:

How do I take a file of text data or OGR readable point data and turn it into an interpolated regular grid?

Depending on what your source data looks like, you can deal with this through a three step process, but is fairly straightforward if you use the right tools. Here we stick to just using GDAL/OGR command line tools.

1. Convert XYZ data into CSV text

If the source data is in a grid text format with rows and columns, then convert it to the XYZ output format with gdal_translate:
gdal_translate -of XYZ input.txt grid.xyz
If your source data is already in an XYZ format, or after above conversion, then change the space delimiters to commas (tabs and semicolon also supported). This is easily done on the command line using a tool like sed or manually in a text editor:


sed 's/ /,/g' grid.xyz > grid.csv

The results would look something like:

-120.5,101.5,100
-117.5,101.5,100
-116.5,100.5,150
-120.5,97.5,200
-119.5,100.5,100

2. Create a virtual format header file

Step 3 will require access to the CSV grid data in an OGR supported format, so we create a simple XML text file that can then be referred to as a vector datasource with any OGR tools.
Create the file grid.vrt with the following content, note the pointers to the source CSV file and the fields which are automatically discerned from the CSV file on-the-fly:

<OGRVRTDataSource>
<OGRVRTLayer name="grid">
<SrcDataSource>grid.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>WGS84</LayerSRS>
<GeometryField separator=" " encoding="PointFromColumns" x="field_1" y="field_2" z="field_3"/>
</OGRVRTLayer>
</OGRVRTDataSource>

If the first row of the grid.csv file includes custom field names, e.g. X, Y, Z, then substitute them for the field_1, field_2 and field_3 settings above.

3. Convert using gdal_grid

Now for the following step, convert this OGR vector VRT and CSV format into a GDAL raster format, while also interpolating any “missing” points in the grid. Various options are available for interpolating. See gdal_grid command syntax for more information, a minimal example is shown here using all the defaults:

gdal_grid -zfield field_3 -l grid grid.vrt newgrid.tif

The output product will fill in any holes in the source grid. Several algorithms are available, the default uses the inverse distance to a power option.

Results shown in the following image.


© Tyler Mitchell / Spatialguru.com