Interpolate raster at points

This example demonstrates the 2D interpolation of raster values to points using interp_points().

We open an example raster, a digital elevation model in South America.

import geoutils as gu

filename_rast = gu.examples.get_path("exploradores_aster_dem")
rast = gu.Raster(filename_rast)
rast = rast.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 2000, rast.bounds.bottom + 2000])

# Plot the raster
rast.plot(cmap="terrain")

We generate a random subsample of 100 coordinates to interpolate.

import numpy as np

rng = np.random.default_rng(42)
x_coords = rng.uniform(rast.bounds.left + 50, rast.bounds.right - 50, 50)
y_coords = rng.uniform(rast.bounds.bottom + 50, rast.bounds.top - 50, 50)

pc = rast.interp_points(points=(x_coords, y_coords))

We plot the resulting point cloud

pc.plot(ax="new", cmap="terrain", marker="x", cbar_title="Elevation (m)")

Important

The interpretation of where raster values are located can differ. The parameter shift_area_or_point (off by default) can be turned on to ensure that the pixel interpretation of your dataset is correct.

Let’s look and redefine our pixel interpretation into "Point". This will shift interpolation by half a pixel.

rast.area_or_point
rast.area_or_point = "Point"

We can interpolate again by shifting according to our interpretation, and changing the resampling algorithm (default to “linear”).

pc_shifted = rast.interp_points(points=(x_coords, y_coords), shift_area_or_point=True, method="quintic")
np.nanmean(pc - pc_shifted)

np.float32(-0.048779298)

The mean difference in interpolated values is quite significant, with a 2-meter bias!