Spatial Interpolation with GDAL in Python #2: IDW and Linear Interpolation

In this second interpolation tutorial, I talk about the Inverse Distance to a Power and Linear Interpolation algorithms available for the program gdal_grid (or gdal.Grid() if you use the GDAL Python API). I explain a little about the theoretical background of these spatial interpolation algorithms and then apply them to create a regular grid from scattered point data in Python.
gdal_grid documentation: gdal.org/programs/gdal_grid.html
GDAL Grid Tutorial: gdal.org/tutorials/gdal_grid_...
GDAL/OGR Python API: gdal.org/python/
Great playlist on Delaunay triangulation: • Delaunay Triangulation...
Chapters:
0:00 Recap
1:06 Inverse Distance to a Power
10:54 Linear Interpolation
Updated 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
pts = layer = None
nearest neighbor interpolation
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=500:radius2=800:angle=20",
outputBounds = [ulx,uly,lrx,lry],
width = xsize, height = ysize)
ma = None
inverse distance to a power
idw = gdal.Grid("invdist.tif", "points.shp", zfield = "elevation",
algorithm = "invdist:power=3:radius1=2000:radius2=2000",
outputBounds = [ulx,uly,lrx,lry],
width = xsize, height = ysize)
idw = None
linear interpolation
lin = gdal.Grid("linear.tif", "points.shp", zfield = "elevation",
algorithm = "linear",
outputBounds = [ulx,uly,lrx,lry],
width = xsize, height = ysize)
lin = None

[@user-qy8yq5ff7y]

Extremely helpful! Please keep updating this series!💕

[@dazhouwu9425]

Thanks for the whole series of gdal in python tutorials! You really make everything clear!

[@No0neIsInnocent]

Many thanks for all your videos about GDAL, it's highly appreciated. 👍

[@Antonio-lp8hx]

highly appreciate your work, especially those awesome 3D animations! (y)

[@kmehour]

We need more python/gdal tutorials !!

[@duncangraham9600]

This was so good, thank you!

[@ernestbayah1640]

Fantastic. Many thanks for this,

[@user-iz5kk4sy6m]

Perfect explanation, thank you!

[@rameshbanagar3946]

Wonderful! Thank you for sharing :)

[@saberchaghamirza5678]

Its great tanks for sharing

[@AJRG1410]

That was a great vid! May I ask how did you do the animations? Loved them!

[@navneetsrivastava688]

Please complete the gdal python series. Also it would be very much helpful if you could create a video on how to convert lat long ( spatial coordinates) to pixel coordinates (x,y) and vice versa, with introduction to different projection and crs. Thank you for the video though

[@raghwendranarayanshandilya4525]

Can you please make a video on Kriging interpolation using Python, where X and Y coordinates and Z values are stored in a shapefile? Or if you already have please suggest the link.

[@gz6616]

Hi, I'm wondering is it possible to provide the 2nd argument to the Grid() function, not using a "points.shp" file name, but some other format? I don't have a points.shp file, but only the x- and y- and z- coordinates. How do I provide those to the Grid() function?

[@draufunddran]

how did you create these 3d animations? i would love to do something similar but actually not animated. Is this donw directly in QGIS?

[@makingsenseremotely6207]

Yes, I used the Qgis2threejs plugin for this

[@jeewankbastola2272]

Please provide the files as well so it wont be confusing to practice the same