Gridding Geospatial Points to Raster with GDAL.Grid and Python

I once did a project where I analyzed extreme precipitation at different duration and return periods (i.e., Depth Duration Frequency-DDF) at several sites. These sites locate close to each other. Generally, we provide our client the final outputs with tabular format and charts. However, the client really want to see the spatial variability of the DDF at a specific duration and return year across sites.

Under such an extra requirement, a basic procedure is to grid a DDF from the scattered sites to a regular raster, and then visualize the raster with some contour lines. Sure, this could definitely be done with the spatial analysis toolbox of ArcMap. However, sometimes we may prefer finishing the task by a pythonic way. At the first glance, we may say 'en', it is easy! There are so many interpolation packages available. Sure, we can try scipy.spatial.griddata, or the algorithms of natural neighbor, bilinear interpolation, Cressman, Barnes or Radial basis function interpolation, etc.

However, these methods did not work as expected for the current task. The raster generated by some of them can not cover the whole domain, while others may produce weird contour lines after interpolation. A possible reason might be that there are just few points that are not large enough find appropriate relationships among points for these interpolation methods.

Finally, we perfectly solved the problem with GDAL Grid, which provides the desirous gridding function. Gridding is a process of creating a regular grid (or call it a raster image) from the scattered data. Typically you have a set of arbitrary scattered over the region of survey measurements and you would like to convert them into the regular grid for further processing and combining with other grids.

Let's have a shortest test in this tutorial.

import numpy as np
import geopandas as gpd

from osgeo import gdal
import rioxarray as rio
import matplotlib.pyplot as plt
%matplotlib inline        

Preprocess

We presume you already transformed a DDF data for a specific duration and return period to a shapefile. This can be done by the steps:

• Using?pandas?to read the DDF data from a .csv file, which contains the coordinate in Longitude and Longitude for each site
• using geopandas to create a GeoDataFrame from a Pandas.DataFrame with coordinates
• Save the GeoDataFrame to a shapefile with geopandas.to_file('data/gdal_grid/shp/Sites_DDF.shp')

So our scatter points are the sites, whose name is?data/gdal_grid/shp/Sites_DDF.shp. By the way, we give the interpolated raster a name of?data/gdal_grid/InterpolatedDDF.tif.

shp_file = "data/gdal_grid/shp/Sites_DDF.shp"
rst_file = 'data/gdal_grid/InterpolatedDDF.tif'        

Check the Sites

It is quite simple to load the site shapefile into a GeoDataFrame. The column of?DDF?is the target to be interpolated.

gdf_ddf  = gpd.read_file(shp_file)        

Grid scattered sites to a regular raster with gdal.Grid

gdal.Grid creates regular grid (raster) from the scattered data read from the OGR datasource. Input data will be interpolated to fill grid nodes with values. More info about gdal.Grid() is available in the?GDAL/OGR Python API. You can also choose various?interpolation methods. Here, we select the Inverse Distance Weighted (IDW) interpolation method. In addition, some input parameters could be put into?gdal.GridOptions().

Here we explicitly set the inputs parameters as follows:

• output file - rst_file
• input file - shp_file
• output raster format): GTiff
• interpolation algorithm: invdist
• target zfield: DDF

rasterDs = gdal.Grid(rst_file, 
                     shp_file, 
                     format='GTiff',
                     algorithm='invdist', 
                     zfield='DDF')
rasterDs.FlushCache()        

Visualize

Now we have the interpolated raster, let's have spatial visualization to check the interpolation effects.

We put the following three data layers into a map:

• interpolated raster
• contour lines
• points of sites

da_ddf = rio.open_rasterio(rst_file)

fig, ax = plt.subplots()
ax.axis('off')

# plot interpolated raster
da_ddf.plot(ax=ax, add_colorbar=False)

# plot contour lines
xx, yy = np.meshgrid(da_ddf.x, da_ddf.y)
cs = ax.contour(xx, yy, da_ddf[0], levels=10, cmap='Reds')
ax.clabel(cs, inline=1, fontsize=10)

# plot raw points
gdf_ddf.plot(ax=ax, markersize=50, marker="*", color='y', edgecolor='purple');

# set descriptions
ax.set_title('Interpolate Geospatial Points to Raster\n', 
          fontweight='bold',  size=22)

_ = ax.text(0.5, 1.025, r"with GDAL Grid", 
            fontsize=16, 
            horizontalalignment='center',             
            verticalalignment='center',
            fontstyle='italic',
            transform=ax.transAxes)
        

You can try your own data with the above procedure and you will definitely find it is quite beautiful. Moreover, the contour lines are all quite smooth.

Summary

GDAL?is the Swiss Army knife of geographic information systems (GIS), letting users read, write and refine rasters and vectors in hundreds of different formats via a common abstraction.

GDAL is really really powerful and it provides many many useful tools and functions. Sometimes, only a single function or command line may implement a complicated geospatial processing task easily.

It is worth spending time to learn it and until to master it. Have a try!!!