⚠️ Our 0.1 release refactored several early-development functions for long-term stability, to update your code see here. ⚠️
Future changes will come with deprecation warnings! 🙂

Interpolation from grid

Interpolation from grid#

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.crop([rast.bounds.left, rast.bounds.bottom, rast.bounds.left + 2000, rast.bounds.bottom + 2000])

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

We generate a random subsample of 100 points to interpolate, and extract the coordinates.

import geopandas as gpd
import numpy as np

# Replace by Raster function once done (valid coords)
np.random.seed(42)
x_coords = np.random.uniform(rast.bounds.left + 50, rast.bounds.right - 50, 50)
y_coords = np.random.uniform(rast.bounds.bottom + 50, rast.bounds.top - 50, 50)

vals = rast.interp_points(points=list(zip(x_coords, y_coords)))

Replace by Vector function once done

ds = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=x_coords, y=y_coords), crs=rast.crs)
ds["vals"] = vals
ds.plot(column="vals", cmap="terrain", legend=True, vmin=np.nanmin(rast), vmax=np.nanmax(rast), marker="x")
interpolation
<Axes: >

Important

The interpretation of where raster values are located differ, see XXX. 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.

print(rast.tags["AREA_OR_POINT"])
rast.tags["AREA_OR_POINT"] = "Point"
Area

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

vals_shifted = rast.interp_points(points=list(zip(x_coords, y_coords)), shift_area_or_point=True, method="quintic")
np.nanmean(vals - vals_shifted)
0.469777

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

Total running time of the script: (0 minutes 0.279 seconds)

Gallery generated by Sphinx-Gallery