Reduce raster around points

This example demonstrates the reduction of windowed raster values around a point using reduce_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 extract.

import geopandas as gpd
import numpy as np

# Replace by Raster function once done
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.reduce_points((x_coords, y_coords))

/home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/latest/examples/handling/raster_point/reduction.py:32: DeprecationWarning: reduce_points() is deprecated; use reduce_at_points() instead.
  pc = rast.reduce_points((x_coords, y_coords))

We plot the resulting point cloud

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

By default, reduce_points() extracts the closest pixel value. But it can also be passed a window size and reductor function to extract an average value or other statistic based on neighbouring pixels.

pc_reduced = rast.reduce_points((x_coords, y_coords), window=5, reducer_function=np.nanmedian)

np.nanmean(pc - pc_reduced)

/home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/latest/examples/handling/raster_point/reduction.py:42: DeprecationWarning: reduce_points() is deprecated; use reduce_at_points() instead.
  pc_reduced = rast.reduce_points((x_coords, y_coords), window=5, reducer_function=np.nanmedian)

np.float64(0.1735498046875)

The mean difference in extracted values is quite significant at 0.3 meters! We can visualize how the sampling took place in the windows.

# Replace by Vector function once done
coords = rast.coords(grid=True)
x_closest = rast.copy(new_array=coords[0]).reduce_points((x_coords, y_coords), as_array=True).squeeze()
y_closest = rast.copy(new_array=coords[1]).reduce_points((x_coords, y_coords), as_array=True).squeeze()
from shapely.geometry import box

geometry = [
    box(x - 2 * rast.res[0], y - 2 * rast.res[1], x + 2 * rast.res[0], y + 2 * rast.res[1])
    for x, y in zip(x_closest, y_closest)
]
ds = gpd.GeoDataFrame(geometry=geometry, crs=rast.crs)
ds["vals"] = pc_reduced.data
ds.plot(column="vals", cmap="terrain", legend=True, vmin=np.nanmin(rast), vmax=np.nanmax(rast))

/home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/latest/examples/handling/raster_point/reduction.py:52: DeprecationWarning: reduce_points() is deprecated; use reduce_at_points() instead.
  x_closest = rast.copy(new_array=coords[0]).reduce_points((x_coords, y_coords), as_array=True).squeeze()
/home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/latest/examples/handling/raster_point/reduction.py:53: DeprecationWarning: reduce_points() is deprecated; use reduce_at_points() instead.
  y_closest = rast.copy(new_array=coords[1]).reduce_points((x_coords, y_coords), as_array=True).squeeze()

<Axes: xlabel='Easting [metre]', ylabel='Northing [metre]'>