Spatial Interpolation with GDAL in Python #1: Nearest Neighbor and Moving Average

In this tutorial, I will give an introduction to the spatial interpolation algorithms nearest neighbor and moving average. We will use gdal.Grid() from the geospatial data abstraction library GDAL to create a regular grid from scattered point data in Python.

gdal_grid documentation: https://gdal.org/programs/gdal_grid.html
GDAL Grid Tutorial: https://gdal.org/tutorials/gdal_grid_...
GDAL/OGR Python API: https://gdal.org/python/

Chapters:
0:00 Introduction
1:33 Nearest Neighbor
14:14 Moving Average

Code:

from osgeo import gdal
from osgeo import ogr

pts = ogr.Open("points.shp", 0)
layer = pts.GetLayer()

for field in layer.schema:
print(field.name)

dem = gdal.Open("dem.tif")
gt = dem.GetGeoTransform()

ulx = gt[0]
uly = gt[3]
res = gt[1]

xsize = dem.RasterXSize
ysize = dem.RasterYSize

lrx = ulx + xsize * res
lry = uly - ysize * res

dem = None

nearest neighbor interpolation
pts = layer = None
nn = gdal.Grid("nearest.tif", "points.shp", zfield="elevation",
algorithm = "nearest", outputBounds = [ulx,uly,lrx,lry],
width = xsize, height = ysize)
nn = None

moving average
ma = gdal.Grid("average.tif", "points.shp", zfield="elevation",
algorithm = "average:radius1=1000:radius2=800:angle=20",
outputBounds = [ulx,uly,lrx,lry],
width = xsize, height = ysize)
ma = None