Using gdal to interpolate elevation from a GeoTIFF

Question

I am new to GDAL and am looking for pointers about the Correct™ way to accomplish a few tasks. I would like to do this in Python, not on the command line.

I have a GeoTIFF that shows elevation at 10m intervals over a given area. That is, every pixel is a 10m x 10m area. The GeoTIFF contains a single raster image, where elevation in meters is expressed as a Float32. The GeoTIFF seems to have tags defining the geoid/datum (which UTM zone etc).

I wish to open that image, and query what the elevation is at a certain point, such as in the following pseudocode:

elevation = myGeoTiffData.getValueAtPosition(Position('WGS84', 39.584588, -122.461320))

Further, I would like to have this elevation interpolated (linearly, spline-wise, however else) if the requested position does not fall on a 10m x 10m lattice point. I see there are some functions to do this (GDAL grid?), but I do not understand how to use them in this case.

I have tried to read the GDAL documentation, but I can't seem to figure this out. Is there any advice people here can offer?

Answer 1

A high-level rasterstats tool is built on GDAL to fetch data from raster and vector data sources.

import rasterstats
from shapely.geometry import Point

raster = "elevation.tif"
pt = Point(245000, 1000000)

# default is interpolate="bilinear"
elev = rasterstats.point_query([pt], raster)[0]

# or use nearest
elev = rasterstats.point_query([pt], raster, interpolate="nearest")[0]

Transformations to/from lat/long coordinates can be done using pyproj.

Answer 2

rioxarray has an interpolation method that will work really well for this. See the details here https://corteva.github.io/rioxarray/stable/examples/interpolate_na.html

import rioxarray

raster_file = './elevation.tif'

ds = rioxarray.open_rasterio(raster_file)
point_elevation = ds.interp(x= -122.461320, y= 39.584588, method='linear')

Answer 3

gdallocationinfo is a utility for getting a value from a raster band or bands given a point location. However it does not do interpolation. I would recommend resampling to a higher resolution with an appropriate resampling algorithm. gdalwarp is the standard GDAL tool for this, and includes several resampling algorithms but it is up to you to choose the appropriate one.

I'm unaware of any GDAL tool that can do an on-the-fly continuous point interpolation.