Add line zonal stats notebook

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
+}

docs/index.md

@@ -30,6 +30,14 @@ rasterize/geometry_mask
 rasterize/exactextract
 ```
 
+```{toctree}
+---
+caption: Examples
+hidden:
+---
+examples/line-zonal-stats
+```
+
 ```{toctree}
 ---
 caption: Reference

pyproject.toml

@@ -62,6 +62,7 @@ docs = [
   'furo',
   'numpydoc',
   'sphinxext-opengraph[social_cards]',
+  'xvec',
 ]