From a73c5aaa3c7ad1835594e2d59845deff2d82cdac Mon Sep 17 00:00:00 2001 From: Deepak Cherian Date: Mon, 9 Feb 2026 13:07:03 -0700 Subject: [PATCH] Add line zonal stats notebook --- docs/examples/line-zonal-stats.ipynb | 231 +++++++++++++++++++++++++++ docs/index.md | 8 + pyproject.toml | 1 + 3 files changed, 240 insertions(+) create mode 100644 docs/examples/line-zonal-stats.ipynb diff --git a/docs/examples/line-zonal-stats.ipynb b/docs/examples/line-zonal-stats.ipynb new file mode 100644 index 0000000..a9392e6 --- /dev/null +++ b/docs/examples/line-zonal-stats.ipynb @@ -0,0 +1,231 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Zonal stats along line segments\n", + "\n", + "This quick example shows off using {py:func}`~rasterix.rasterize.rasterize`, and Xarray's GroupBy to do zonal stats along a line." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": { + "tags": [ + "remove-cell" + ] + }, + "source": [ + "## Example dataset\n", + "\n", + "Here is an example dataset in the EPSG:3035 coordinate system with a coordinate system specified using a `\"GeoTransform'` attribute." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pyproj\n", + "import shapely\n", + "import xarray as xr\n", + "\n", + "ds = xr.Dataset(\n", + " coords={\n", + " \"spatial_ref\": (\n", + " (),\n", + " 0,\n", + " pyproj.CRS.from_epsg(3035).to_cf() | {\"GeoTransform\": \"100000 12000 0 60000 0 -12000\"},\n", + " )\n", + " },\n", + ")\n", + "ds[\"foo\"] = xr.DataArray(\n", + " np.random.default_rng(0).integers(0, 256, size=(100, 200), dtype=np.uint8),\n", + " dims=(\"y\", \"x\"),\n", + ")\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## Generate coordinates with `RasterIndex`\n", + "\n", + "We'll use {py:func}`~rasterix.assign_index` to assign lazy coordinate locations in `x` and `y`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "import rasterix\n", + "\n", + "ds = rasterix.assign_index(ds)\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## Make a geometry\n", + "\n", + "We will take a LineString and split it in to 80km segments.\n", + "\n", + "For simplicity the line will extend from the top-left corner to the bottom-right corner." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "import geopandas as gpd\n", + "\n", + "from rasterix.rasterize import rasterize\n", + "\n", + "# make a diagonal line across the raster\n", + "line = gpd.GeoDataFrame(\n", + " {\n", + " \"geometry\": [\n", + " shapely.LineString(\n", + " [\n", + " [ds.x[0].item(), ds.y[0].item()],\n", + " [ds.x[-1].item(), ds.y[-1].item()],\n", + " ]\n", + " )\n", + " ]\n", + " }\n", + ")\n", + "# 80km segments\n", + "line = line.segmentize(80_000)\n", + "coords = shapely.get_coordinates(line)\n", + "segment_coords = np.stack([coords[:-1], coords[1:]], axis=1) # shape (n-1, 2, 2)\n", + "segments = shapely.linestrings(segment_coords)\n", + "geoms = gpd.GeoDataFrame({\"geometry\": segments})\n", + "geoms" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## Rasterize the geometries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "res = (\n", + " rasterize(\n", + " ds,\n", + " geoms,\n", + " all_touched=True,\n", + " engine=\"rasterio\",\n", + " )\n", + " .compute()\n", + " .rename(\"segment\")\n", + ")\n", + "res.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "## Zonal stats with Xarray's GroupBy\n", + "\n", + "(faster when powered by [flox](https://flox.readthedocs.io/en/latest/)!)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "zonal = ds[\"foo\"].groupby(res).max()\n", + "# drop the last segment (background/NaN pixels)\n", + "zonal = zonal.isel(segment=slice(-1))\n", + "zonal" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "## Make a vector data cube with Xvec" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "import xvec # noqa: F401\n", + "\n", + "zonal_ds = geoms.assign(mean_foo=zonal.values).set_geometry(\"geometry\").set_index(\"geometry\").to_xarray()\n", + "zonal_ds = zonal_ds.xvec.set_geom_indexes(\"geometry\", crs=3035)\n", + "zonal_ds" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "## Plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "f, ax = zonal_ds[\"mean_foo\"].xvec.plot()\n", + "f.set_size_inches((8, 2))" + ] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/index.md b/docs/index.md index 5868eef..c66e7f5 100644 --- a/docs/index.md +++ b/docs/index.md @@ -27,6 +27,14 @@ rasterize/geometry_mask rasterize/exactextract ``` +```{toctree} +--- +caption: Examples +hidden: +--- +examples/line-zonal-stats +``` + ```{toctree} --- caption: Reference diff --git a/pyproject.toml b/pyproject.toml index 75aede2..3873442 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -62,6 +62,7 @@ docs = [ 'furo', 'numpydoc', 'sphinxext-opengraph[social_cards]', + 'xvec', ]