### Setup for xarray-spatial Migration Examples

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/source/user_guide/local-migration.md

Imports necessary libraries and creates sample xarray DataArrays and a Dataset. It then stacks the variables into a single DataArray for easier manipulation.

```python
import numpy as np
import xarray as xr

arr1 = xr.DataArray([[np.nan, 4, 2, 0],
                     [2, 3, np.nan, 1],
                     [5, 1, 2, 0],
                     [1, 3, 2, np.nan]], name="arr1")

arr2 = xr.DataArray([[3, 1, 1, 2],
                     [4, 1, 2, 5],
                     [0, 0, 0, 0],
                     [np.nan, 1, 1, 1]], name="arr2")

arr3 = xr.DataArray([[3, 3, 2, 0],
                     [4, 1, 3, 1],
                     [6, 1, 2, 2],
                     [0, 0, 1, 1]], name="arr3")

ds = xr.merge([arr1, arr2, arr3])

# Stack all variables into a single DataArray with a "var" dimension.
# This is the key building block for all the replacements below.
stacked = ds.to_array(dim="var")
```

--------------------------------

### Install ASV dependencies

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/benchmarks/README.md

Install the necessary tools to run the benchmarking suite.

```bash
pip install asv virtualenv
```

--------------------------------

### Download xarray-spatial examples only

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/README.md

Downloads only the example notebooks. Note that many examples may not run without the sample data.

```bash
xrspatial copy-examples
```

--------------------------------

### Download xarray-spatial examples and data

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/README.md

Downloads both the example notebooks and the associated sample data.

```bash
xrspatial examples
```

--------------------------------

### Initializing Flow Path Start Points

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/11_Hydrology.ipynb

Creates an empty DataArray to serve as the starting grid for flow path tracing.

```python
start_points = xr.DataArray(
    np.full((H, W), np.nan, dtype=np.float64),
    dims=dem.dims, coords=dem.coords,
)
```

--------------------------------

### Importing Required Libraries

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/1_Surface.ipynb

Initial setup for data manipulation, visualization, and spatial analysis.

```python
import numpy as np
import pandas as pd
import xarray as xr

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, hsv_to_rgb
from matplotlib.patches import Patch

import xrspatial
```

--------------------------------

### Install xarray-spatial via package managers

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/source/getting_started/installation.md

Use these commands to install the library in your local environment.

```bash
# via pip
pip install xarray-spatial

# via conda
conda install -c conda-forge xarray-spatial
```

--------------------------------

### Verify xarray-spatial Installation

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/0_Getting_Setup.ipynb

Open a Python session and import the package to verify that the installation was successful.

```python
import xrspatial
```

--------------------------------

### Update README GeoTIFF Usage Examples

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-03-30-geotiff-perf-controls.md

Add these examples to the README.md file to demonstrate GeoTIFF I/O functionalities. They cover specifying data types, chunking with dask, and controlling compression levels for writing GeoTIFFs, as well as streaming dask arrays to a VRT file.

```python
open_geotiff('dem.tif', dtype='float32')              # half memory
```

```python
open_geotiff('dem.tif', dtype='float32', chunks=512)   # dask + half memory
```

```python
to_geotiff(data, 'out.tif', compression_level=1)       # fast scratch write
```

```python
to_geotiff(data, 'out.tif', compression_level=22)      # max compression
```

```python
to_geotiff(dask_da, 'mosaic.vrt')                      # stream dask to VRT
```

--------------------------------

### Commit User Guide Notebook

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-03-30-geotiff-perf-controls.md

Stage and commit the newly created user guide notebook for GeoTIFF performance controls. This command adds the notebook file to the Git staging area and commits it with a message indicating the addition and the relevant issue number.

```bash
cd .claude/worktrees/issue-1083
git add examples/user_guide/46_GeoTIFF_Performance.ipynb
git commit -m "Add user guide notebook for geotiff performance controls (#1083)"
```

--------------------------------

### Setup Complex Terrain and Observer

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/source/user_guide/surface.ipynb

Prepares a terrain dataset and an observer location for advanced viewshed analysis.

```python
from xrspatial import viewshed

x_range = (-20e6, 20e6)
y_range = (-20e6, 20e6)
terrain = xr.DataArray(np.zeros((H, W)))
terrain = generate_terrain(terrain, x_range=x_range, y_range=y_range)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how="linear")

illuminated = hillshade(terrain)

OBSERVER_X = 0.0
OBSERVER_Y = 0.0

cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=x_range, y_range=y_range)
observer_df = pd.DataFrame({"x": [OBSERVER_X], "y": [OBSERVER_Y]})
observer_agg = cvs.points(observer_df, "x", "y")
observer_shaded = dynspread(shade(observer_agg, cmap=["orange"]), threshold=1, max_px=4)

stack(
    shade(illuminated, cmap=["black", "white"], alpha=128, how="linear"),
    terrain_shaded,
    observer_shaded,
)
```

--------------------------------

### Import xarray-spatial GeoTIFF modules

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/39_GeoTIFF_IO.ipynb

Initial setup for using xarray-spatial GeoTIFF functionality.

```python
%matplotlib inline
import tempfile
import os

import numpy as np
import xarray as xr

import matplotlib.pyplot as plt

import xrspatial
from xrspatial.geotiff import open_geotiff, to_geotiff, write_vrt
```

--------------------------------

### Plotting Pathfinding Comparisons

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/7_Pathfinding.ipynb

Use this code to visualize and compare the output of different pathfinding algorithms. Ensure that 'surface', 'path_8', 'path_4', 'path_friction', 'start', and 'goal' variables are defined and properly formatted before execution. Coordinates for 'start' and 'goal' should be in (y, x) order.

```python
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
paths = [path_8, path_4, path_friction]
colors = ['dodgerblue', 'darkorange', 'mediumpurple']
titles = ['8-connectivity', '4-connectivity', 'Friction-weighted']

for ax, path, color, title in zip(axes, paths, colors, titles):
    surface.plot.imshow(ax=ax, cmap=ListedColormap(['#e0e0e0']),
                        add_colorbar=False)
    path.plot.imshow(ax=ax, cmap=ListedColormap([color]),
                     add_colorbar=False)
    ax.scatter([start[1]], [start[0]], c='orange', s=100, zorder=10,
               edgecolors='white', linewidths=1)
    ax.scatter([goal[1]], [goal[0]], c='dodgerblue', s=100, zorder=10,
               edgecolors='white', linewidths=1)
    ax.set_title(title, fontsize=13)
    ax.set_axis_off()
plt.tight_layout()
```

--------------------------------

### Create a traversable grid for pathfinding

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/7_Pathfinding.ipynb

Initializes a 20x20 raster grid where all cells are traversable (value 1). This setup is useful for visualizing pathfinding results clearly. Start and goal coordinates are defined in (y, x) order.

```python
W, H = 20, 20

surface = xr.DataArray(
    np.ones((H, W), dtype=np.int32),
    dims=['y', 'x'],
    coords={'y': np.arange(H - 1, -1, -1), 'x': np.arange(W)}
)

start = (17, 2)
goal = (2, 17)
```

--------------------------------

### Trace Flow Paths and Visualize

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/11_Hydrology.ipynb

This snippet demonstrates tracing flow paths from identified start points and visualizing them with a neon color palette on a hillshade background. It includes dilating paths for better visibility and adding start-point markers.

```python
accum_vals = flow_accum.values.copy()
accum_vals[np.isnan(accum_vals)] = 0
trib_mask = (accum_vals >= 50) & (accum_vals <= 200)
trib_rows, trib_cols = np.where(trib_mask)

regions = [
    (slice(20, H // 3), slice(20, W // 3)),
    (slice(20, H // 3), slice(2 * W // 3, W - 20)),
    (slice(H // 3, 2 * H // 3), slice(W // 3, 2 * W // 3)),
    (slice(2 * H // 3, H - 20), slice(20, W // 3)),
    (slice(2 * H // 3, H - 20), slice(2 * W // 3, W - 20)),
]
label = 1
for ys, xs in regions:
    in_region = (
        (trib_rows >= ys.start) & (trib_rows < ys.stop) &
        (trib_cols >= xs.start) & (trib_cols < xs.stop)
    )
    if np.any(in_region):
        idx = np.where(in_region)[0][len(np.where(in_region)[0]) // 2]
        start_points.values[trib_rows[idx], trib_cols[idx]] = float(label)
        label += 1

paths = xrspatial.flow_path(flow_dir, start_points)

n_path_cells = int(np.sum(~np.isnan(paths.values)))
n_starts = int(np.sum(~np.isnan(start_points.values)))
print(f"Traced {n_path_cells} path cells from {n_starts} start points")

from scipy.ndimage import maximum_filter
path_vals = paths.values.copy()
dilated = np.full_like(path_vals, np.nan)
for lbl in range(1, n_starts + 1):
    mask = (path_vals == lbl)
    thick = maximum_filter(mask.astype(float), size=3) > 0
    dilated[thick] = float(lbl)
paths_thick = xr.DataArray(dilated, dims=paths.dims, coords=paths.coords)

neon_colors = ['#ff00ff', '#00ffff', '#39ff14', '#ff6600', '#ffff00']
neon_cmap = ListedColormap(neon_colors[:n_starts])

sp_rows, sp_cols = np.where(~np.isnan(start_points.values))
sp_x = dem.x.values[sp_cols]
sp_y = dem.y.values[sp_rows]

fig, ax = plt.subplots(figsize=(10, 7.5))
fig.patch.set_facecolor('#1a1a2e')
ax.set_facecolor('#1a1a2e')

hillshade.plot.imshow(ax=ax, cmap='gray', alpha=0.3, add_colorbar=False)

paths_thick.plot.imshow(ax=ax, cmap=neon_cmap, alpha=0.95, add_colorbar=False)

ax.scatter(sp_x, sp_y, c='white', s=250, zorder=6, edgecolors='black',
           linewidths=2, marker='o')
for i, (xi, yi) in enumerate(zip(sp_x, sp_y)):
    ax.scatter([xi], [yi], c=neon_colors[i], s=120, zorder=7,
               edgecolors='none', marker='o')

ax.legend(handles=[Patch(facecolor=c, label=f'Path {i+1}')
                   for i, c in enumerate(neon_colors[:n_starts])],
          loc='lower right', fontsize=10, framealpha=0.85,
          facecolor='#1a1a2e', labelcolor='white', edgecolor='#444')
ax.set_axis_off()
```

--------------------------------

### Import xrspatial Morphology Submodule

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/17_Morphological_Operators.ipynb

Imports necessary libraries including numpy, xarray, matplotlib, and specific functions from xrspatial.morphology. Ensure these are installed before running.

```python
%matplotlib inline
import numpy as np
import pandas as pd
import xarray as xr

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

import xrspatial
from xrspatial.morphology import (
    _circle_kernel,
    morph_black_tophat,
    morph_closing,
    morph_dilate,
    morph_erode,
    morph_gradient,
    morph_opening,
    morph_white_tophat,
)
```

--------------------------------

### Import xarray-spatial dasymetric modules

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/14_Dasymetric_Mapping.ipynb

Standard setup for dasymetric analysis including the required disaggregate, pycnophylactic, and validation functions.

```python
%matplotlib inline
import numpy as np
import pandas as pd
import xarray as xr

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

from xrspatial import disaggregate
from xrspatial.dasymetric import pycnophylactic, validate_disaggregation
```

--------------------------------

### Dask Parallel Rasterization Setup

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/28_Rasterize.ipynb

Introduces the concept of using Dask for parallel rasterization by specifying `chunks`. The result is a lazy dask-backed DataArray that needs to be computed. Requires `dask`.

```python
import dask.array as da
```

--------------------------------

### Import xarray-spatial MCDA modules

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/41_MCDA.ipynb

Initial setup for the MCDA workflow, including necessary imports for terrain generation, proximity analysis, and MCDA functions.

```python
%matplotlib inline
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

from xrspatial import generate_terrain, slope, aspect
from xrspatial.proximity import proximity
from xrspatial.mcda import (
    standardize, ahp_weights, rank_weights,
    wlc, wpm, owa, fuzzy_overlay, boolean_overlay,
    constrain, sensitivity,
)
```

--------------------------------

### A* Search: Unweighted vs. Weighted Paths

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/source/user_guide/pathfinding.ipynb

Compares paths found using A* search with and without friction weighting. Choose start and goal coordinates to compare geometric distance path with least-cost path considering friction.

```python
# Choose start and goal coordinates
start_coord = (400.0, 50.0)
goal_coord = (100.0, 450.0)

# Unweighted path (geometric distance only)
path_unweighted = a_star_search(terrain, start_coord, goal_coord)

# Weighted path (friction from slope)
path_weighted = a_star_search(terrain, start_coord, goal_coord, friction=friction)

print(f"Unweighted path cost at goal: {np.nanmax(path_unweighted):.2f}")
print(f"Weighted path cost at goal:   {np.nanmax(path_weighted):.2f}")
```

--------------------------------

### Import necessary libraries for GeoTIFF performance notebook

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-03-30-geotiff-perf-controls.md

Import NumPy, Xarray, OS, and temporary file utilities, along with the open_geotiff and to_geotiff functions from xrspatial.geotiff. These are required for creating, manipulating, and writing GeoTIFF files within the user guide notebook.

```python
import numpy as np, xarray as xr, os, tempfile
from xrspatial.geotiff import open_geotiff, to_geotiff
```

--------------------------------

### Visualize the traversable grid with start and goal markers

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/7_Pathfinding.ipynb

Renders the traversable grid using matplotlib, highlighting the start and goal points with distinct markers. This visualization helps in understanding the initial setup for pathfinding.

```python
fig, ax = plt.subplots(figsize=(8, 8))
surface.plot.imshow(ax=ax, cmap=ListedColormap(['#e0e0e0']),
                    add_colorbar=False)
ax.scatter([start[1]], [start[0]], c='orange', s=200, zorder=10,
           edgecolors='white', linewidths=2)
ax.scatter([goal[1]], [goal[0]], c='dodgerblue', s=200, zorder=10,
           edgecolors='white', linewidths=2)
ax.legend(handles=[
    Patch(facecolor='orange', alpha=0.78, label='Start (17, 2)'),
    Patch(facecolor='dodgerblue', alpha=0.78, label='Goal (2, 17)'),
], loc='lower left', fontsize=11, framealpha=0.9)
ax.set_axis_off()
```

--------------------------------

### Install xarray-spatial via pip

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/README.md

Install the xarray-spatial library using pip. This is the standard method for installing Python packages.

```bash
pip install xarray-spatial
```

--------------------------------

### Preview benchmark results in browser

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/benchmarks/README.md

Generate and serve the benchmark results locally.

```bash
asv publish
asv preview
```

--------------------------------

### Execute benchmarks

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/benchmarks/README.md

Run the full benchmark suite against the master branch.

```bash
cd benchmarks
asv run
```

--------------------------------

### Import required libraries

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/pharmacy-deserts.ipynb

Initializes the environment with necessary data processing, geospatial, and visualization libraries.

```python
import datashader as ds
from datashader.colors import inferno
from datashader.transfer_functions import set_background
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
import geopandas as gpd
import numpy as np
import pandas as pd
from spatialpandas import GeoDataFrame

import xrspatial
```

--------------------------------

### Compare connectivity settings

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/48_Sieve_Filter.ipynb

Demonstrates the difference between 4-connectivity (rook) and 8-connectivity (queen) in identifying clumps.

```python
sieved_4 = sieve(raster, threshold=6, neighborhood=4)
sieved_8 = sieve(raster, threshold=6, neighborhood=8)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))
for ax, data, title in zip(
    axes,
    [raster, sieved_4, sieved_8],
    ['Original', '4-connectivity (threshold=6)', '8-connectivity (threshold=6)'],
):
    im = ax.imshow(data.values, cmap=cmap, vmin=0.5, vmax=3.5, interpolation='nearest')
    ax.set_title(title)
fig.colorbar(im, ax=axes, ticks=[1, 2, 3], shrink=0.8)
plt.tight_layout()
plt.show()
```

--------------------------------

### Install xarray-spatial via conda

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/README.md

Install the xarray-spatial library using conda from the conda-forge channel. This is recommended for users who manage their environments with conda.

```bash
conda install -c conda-forge xarray-spatial
```

--------------------------------

### User Guide Notebook Imports

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-03-24-dask-graph-utilities.md

Imports necessary libraries for the Fused Overlap user guide notebook, including xarray, dask, and xrspatial utilities.

```python
import numpy as np
import dask.array as da
import xarray as xr
import xrspatial
from xrspatial.utils import fused_overlap, multi_overlap
```

--------------------------------

### Add top-level export to xrspatial/__init__.py

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-03-24-hypsometric-integral.md

Exposing the hypsometric_integral function at the package level.

```python
from xrspatial.zonal import hypsometric_integral  # noqa
```

--------------------------------

### Compare Suitability Combination Methods

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/41_MCDA.ipynb

Visualizes multiple suitability models side-by-side to compare different overlay operators.

```python
results = {
    'WLC': suit_wlc,
    'WPM': suit_wpm,
    'OWA (conservative)': suit_conservative,
    'OWA (optimistic)': suit_optimistic,
    'Fuzzy AND': fuzzy_overlay(criteria, operator='and'),
    'Fuzzy Gamma 0.9': fuzzy_overlay(criteria, operator='gamma', gamma=0.9),
}

fig, axes = plt.subplots(2, 3, figsize=(16, 9))
for ax, (title, data) in zip(axes.flat, results.items()):
    data.plot.imshow(ax=ax, cmap='RdYlGn', vmin=0, vmax=1,
                     add_colorbar=True)
    ax.set_title(title)
    ax.set_axis_off()
plt.tight_layout()
```

--------------------------------

### Initialize environment and imports

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/10_Mahalanobis_Distance.ipynb

Standard imports and configuration for plotting and Xarray-Spatial functionality.

```python
%matplotlib inline
import numpy as np
import xarray as xr

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Patch, Rectangle

import xrspatial
from xrspatial import mahalanobis

np.random.seed(42)
```

--------------------------------

### Notebook Markdown Header

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-spatial-autocorrelation.md

Markdown content for the spatial autocorrelation user guide notebook.

```markdown
# Spatial Autocorrelation: Moran's I and LISA

This notebook demonstrates how to measure spatial autocorrelation in raster data using `morans_i` (global) and `lisa` (local).

**Spatial autocorrelation** measures whether nearby pixels tend to have similar values (positive autocorrelation) or dissimilar values (negative autocorrelation). It answers the question: "Is the spatial pattern in this raster clustered, dispersed, or random?"

- **Global Moran's I** produces a single statistic for the entire raster.
- **LISA (Local Indicators of Spatial Association)** produces a per-pixel statistic, identifying where clusters and outliers are.
```

--------------------------------

### Add Visibility Imports to __init__.py

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-multi-observer-viewshed.md

Add imports for cumulative_viewshed, line_of_sight, and visibility_frequency to the main xrspatial package initialization file.

```python
from xrspatial.visibility import cumulative_viewshed  # noqa
from xrspatial.visibility import line_of_sight  # noqa
from xrspatials.visibility import visibility_frequency  # noqa
```

--------------------------------

### Aggregate Data and Visualize

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/Pathfinding_Austin_Road_Network.ipynb

Creates a raster representation of the streets and visualizes the start and goal points.

```python
H, W = 600, 800
cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=xrange, y_range=yrange)

street_agg = cvs.line(streets_spd, geometry="geometry")
street_shaded = dynspread(shade(street_agg, cmap=["salmon"]))

# Pick two locations
start = (30.08214069, -97.73662282)
goal = (30.17656606, -97.63753489)

start_df = pd.DataFrame({"x": [start[1]], "y": [start[0]]})
start_agg = cvs.points(start_df, "x", "y")
start_shaded = dynspread(shade(start_agg, cmap=["red"]), threshold=1, max_px=5)

goal_df = pd.DataFrame({"x": [goal[1]], "y": [goal[0]]})
goal_agg = cvs.points(goal_df, "x", "y")
goal_shaded = dynspread(shade(goal_agg, cmap=["lime"]), threshold=1, max_px=5)

set_background(stack(street_shaded, start_shaded, goal_shaded), "black")
```

--------------------------------

### Run benchmarks in development mode

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/benchmarks/README.md

Execute all benchmarks once to check for errors.

```bash
asv dev
```

--------------------------------

### Full Pipeline Comparison: Standard vs. Vegetation-Aware

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/34_Vegetation_Flood_Modeling.ipynb

Executes and compares two full hydrological modeling pipelines: a standard one with uniform parameters and a vegetation-aware one using land cover data. Prints key metrics like runoff, time of concentration, and flooded cell counts.

```python
# Standard pipeline
runoff_std = curve_number_runoff(rainfall, curve_number=80.0)
tt_std     = travel_time(fl_raster, slope_raster, mannings_n=0.03)
depth_std  = flood_depth(hand_raster, water_level=10)

# Vegetation-aware pipeline
n_veg      = vegetation_roughness(nlcd_raster, mode='nlcd')
cn_veg     = vegetation_curve_number(nlcd_raster, sg_raster)
runoff_va  = curve_number_runoff(rainfall, curve_number=cn_veg)
tt_va      = travel_time(fl_raster, slope_raster, mannings_n=n_veg)
depth_va   = flood_depth_vegetation(hand_raster, slope_raster, n_veg,
                                    unit_discharge=2.0)

print("Metric                      Standard    Vegetation-aware")
print("-" * 60)

ro_std_mean = np.nanmean(runoff_std.values)
ro_va_mean = np.nanmean(runoff_va.values)
print(f"Mean runoff (mm)            {ro_std_mean:>8.1f}    {ro_va_mean:>8.1f}")

tc_std = np.nanmax(tt_std.values[np.isfinite(tt_std.values)])
tc_va = np.nanmax(tt_va.values[np.isfinite(tt_va.values)])
print(f"Time of concentration       {tc_std:>8.1f}    {tc_va:>8.1f}")

d_std_vals = depth_std.values[~np.isnan(depth_std.values)]
d_va_vals = depth_va.values[~np.isnan(depth_va.values)]
print(f"Flooded cells               {len(d_std_vals):>8d}    {len(d_va_vals):>8d}")
if len(d_std_vals) > 0 and len(d_va_vals) > 0:
    print(f"Median flood depth (m)      {np.median(d_std_vals):>8.2f}    {np.median(d_va_vals):>8.2f}")
```

--------------------------------

### Render Elevation Data

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/xarray-spatial_classification-methods.ipynb

Renders the loaded elevation raster data using matplotlib. Ensure matplotlib is installed.

```python
raster.plot.imshow(figsize=(15, 10));
```

--------------------------------

### Commit Visibility Analysis Notebook

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-multi-observer-viewshed.md

Stages and commits the visibility analysis user guide notebook to the git repository.

```bash
git add examples/user_guide/37_Visibility_Analysis.ipynb
git commit -m "Add visibility analysis user guide notebook (#1145)"
```

--------------------------------

### Initialize Global Moran's I test file

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-spatial-autocorrelation.md

Header comment for the Global Moran's I test suite.

```python
# ---------------------------------------------------------------------------
# Global Moran's I -- NumPy
```

--------------------------------

### Import necessary libraries

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/34_Vegetation_Flood_Modeling.ipynb

Imports required libraries for xarray-spatial, numpy, matplotlib, and rasterio. Ensure these libraries are installed before running.

```python
%matplotlib inline
import numpy as np
import xarray as xr

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, LogNorm, BoundaryNorm
from matplotlib.patches import Patch

import xrspatial
from xrspatial.flood import (
    flood_depth,
    inundation,
    curve_number_runoff,
    travel_time,
    vegetation_roughness,
    vegetation_curve_number,
    flood_depth_vegetation,
    NLCD_MANNINGS_N,
    NLCD_CURVE_NUMBER,
)
```

--------------------------------

### Import necessary libraries

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/47_Streaming_GeoTIFF_Write.ipynb

Imports required libraries for data manipulation, geospatial operations, and plotting. Ensure these are installed before running.

```python
%matplotlib inline
import tempfile
import os

import numpy as np
import xarray as xr
import dask.array as da
import matplotlib.pyplot as plt

from xrspatial.geotiff import open_geotiff, to_geotiff
```

--------------------------------

### Import necessary libraries for Visibility Analysis Notebook

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-multi-observer-viewshed.md

Import NumPy, Xarray, Matplotlib, and specific functions from xrspatial for visibility analysis.

```python
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

from xrspatial.terrain import generate_terrain
from xrspatial.hillshade import hillshade
from xrspatial.visibility import (
    cumulative_viewshed,
    visibility_frequency,
    line_of_sight,
)
```

--------------------------------

### Import necessary libraries

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/33_Porto_Taxi_Trajectories.ipynb

Imports standard libraries for data manipulation, geospatial operations, plotting, and Xarray-Spatial. Ensure these are installed before running.

```python
%matplotlib inline
import json
import numpy as np
import pandas as pd
import xarray as xr

import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from shapely.geometry import LineString

import xrspatial  # registers .xrs accessor
from xrspatial.utils import ngjit
```

--------------------------------

### Initialize Random Number Generator

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/17_Morphological_Operators.ipynb

Sets up a random number generator for creating noisy binary masks.

```python
rng = np.random.default_rng(42)
```

--------------------------------

### Run Visibility Frequency Tests

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-multi-observer-viewshed.md

Command to verify the visibility frequency implementation.

```bash
pytest xrspatial/tests/test_visibility.py::TestVisibilityFrequency -v
```

--------------------------------

### Helper function to create a dask raster

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-03-24-dask-graph-utilities.md

Creates a dask-backed xarray DataArray for testing purposes. Requires dask and xarray to be installed.

```python
import numpy as np
import xarray as xr
da = pytest.importorskip("dask.array")

def _make_dask_raster(shape=(64, 64), chunks=16, dtype=np.float32):
    data = da.from_array(
        np.random.RandomState(42).rand(*shape).astype(dtype), chunks=chunks
    )
    return xr.DataArray(data, dims=['y', 'x'])
```

--------------------------------

### Import Necessary Libraries

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/source/user_guide/multispectral.ipynb

Imports essential libraries for data manipulation, visualization, and geospatial processing. Ensure these are installed before running the code.

```python
import datashader as ds
from datashader.colors import Elevation
import datashader.transfer_functions as tf
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
from datashader.transfer_functions import dynspread
from datashader.transfer_functions import set_background
from datashader.transfer_functions import Images, Image
from datashader.utils import orient_array
import numpy as np
import xarray as xr
from xrspatial.geotiff import open_geotiff
```

--------------------------------

### Compute Pairwise Least-Cost Corridors

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/22_Corridor_Analysis.ipynb

Demonstrates calculating corridors for multiple source locations simultaneously and visualizing the resulting cost-distance surfaces.

```python
# Three habitat patches
positions = [(5, 5), (5, 25), (25, 15)]
sources = []
for r, c in positions:
    s = np.zeros((n, n))
    s[r, c] = 1.0
    sources.append(make_raster(s))

corridors = least_cost_corridor(friction, sources=sources, pairwise=True)

print("Corridor pairs:", list(corridors.data_vars))

fig, axes = plt.subplots(1, 3, figsize=(16, 5))
for ax, name in zip(axes, corridors.data_vars):
    corridors[name].plot.imshow(ax=ax, cmap='inferno', add_colorbar=True,
                                cbar_kwargs={'label': 'Cost'})
    for r, c in positions:
        ax.plot(c, r, 'c*', markersize=14)
    ax.set_title(name, fontsize=11)
    ax.set_axis_off()

axes[0].legend(handles=[Patch(facecolor='cyan', label='Source locations')],
               loc='lower right', fontsize=10, framealpha=0.9)
plt.tight_layout()

# Save preview image
import pathlib
pathlib.Path('images').mkdir(exist_ok=True)
fig.savefig('images/corridor_analysis_preview.png', bbox_inches='tight', dpi=120)
```

--------------------------------

### Tagging a New Release with Git

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/RELEASE.md

Use this command to create a signed tag for a new release. Ensure you have a GPG key configured with GitHub. Replace 'v0.1.2' with your actual version number.

```bash
git tag -a v0.1.2 -s -m "Version 0.1.2"
git push --tags
```

--------------------------------

### Clean up temporary directory

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/47_Streaming_GeoTIFF_Write.ipynb

Removes the temporary directory and all its contents created during the GeoTIFF writing examples. This is a cleanup step to free up disk space.

```python
import shutil
shutil.rmtree(tmpdir, ignore_errors=True)
```

--------------------------------

### Generate synthetic terrain data

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/4_Focal.ipynb

Create a synthetic elevation raster using `generate_terrain` for demonstrating focal operations. This raster is reused across examples.

```python
W = 800
H = 600
x_range = (-20e6, 20e6)
y_range = (-20e6, 20e6)

terrain = xr.DataArray(np.zeros((H, W)))
terrain = terrain.xrs.generate_terrain(x_range=x_range, y_range=y_range)
illuminated = terrain.xrs.hillshade()
```

--------------------------------

### Verify Imports

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/docs/superpowers/plans/2026-04-01-spatial-autocorrelation.md

Confirm that the `morans_i` and `lisa` functions are correctly imported and accessible from the top-level `xrspatial` package.

```bash
cd .claude/worktrees/issue-1135
python -c "from xrspatial import morans_i, lisa; print('OK')"
```

--------------------------------

### Detect Sinks in DEM

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/11_Hydrology.ipynb

Identify sink cells in a raw DEM to locate depressions. This example calculates the number of sink cells and visualizes their distribution.

```python
flow_dir_raw = xrspatial.flow_direction(dem)
sinks = xrspatial.sink(flow_dir_raw)

n_sink_cells = int(np.sum(~np.isnan(sinks.values)))
n_sink_groups = len(np.unique(sinks.values[~np.isnan(sinks.values)]))
print(f"{n_sink_cells} sink cells in {n_sink_groups} depressions "
      f"({100 * n_sink_cells / (H * W):.1f}% of grid)")

sink_overlay = xr.where(np.isfinite(sinks), 1.0, np.nan)

fig, ax = plt.subplots(figsize=(10, 7.5))
plot_basemap(ax)
sink_overlay.plot.imshow(ax=ax, cmap=ListedColormap(['darkorange']),
                         alpha=0.9, add_colorbar=False)
ax.legend(handles=[Patch(facecolor='darkorange',
                         label=f'Sink cells ({n_sink_groups} depressions)')],
          loc='lower right', fontsize=11, framealpha=0.9)
ax.set_axis_off()
```

--------------------------------

### Sobel X vs. Prewitt X Gradient Comparison

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/37_Edge_Detection.ipynb

Compares the Sobel X and Prewitt X gradient operators. Use this to visualize differences in gradient estimation, especially where gradients change rapidly.

```python
grad_max = max(abs(float(sx.min())), abs(float(sx.max())),
               abs(float(px.min())), abs(float(px.max())))

for ax, data, title in zip(axes[0], [sx, px], ['Sobel X', 'Prewitt X']):
    data.plot.imshow(ax=ax, cmap='RdBu_r', vmin=-grad_max, vmax=grad_max, add_colorbar=False)
    ax.set_title(title, fontsize=13)
    ax.set_axis_off()

# Bottom left: difference map
diff = sx - px
diff_max = max(abs(float(diff.min())), abs(float(diff.max())))
im = diff.plot.imshow(ax=axes[1][0], cmap='RdBu_r', vmin=-diff_max, vmax=diff_max, add_colorbar=False)
axes[1][0].set_title('Sobel X − Prewitt X', fontsize=13)
axes[1][0].set_axis_off()
fig.colorbar(im, ax=axes[1][0], shrink=0.7, label='Difference')

# Bottom right: cross-section
row = H // 2
ax = axes[1][1]
cols = np.arange(W)
ax.plot(cols, sx.values[row], label='Sobel X', color='#2166ac', linewidth=1.2)
ax.plot(cols, px.values[row], label='Prewitt X', color='#b2182b', linewidth=1.2, linestyle='--')
ax.set_title(f'Cross-section at row {row}', fontsize=13)
ax.set_xlabel('Column')
ax.set_ylabel('Gradient response')
ax.legend(fontsize=10)

plt.tight_layout()
```

--------------------------------

### Write and read back GeoTIFF

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/39_GeoTIFF_IO.ipynb

Demonstrates round-trip I/O using to_geotiff and open_geotiff.

```python
tmpdir = tempfile.mkdtemp(prefix='xrs_geotiff_nb_')
path = os.path.join(tmpdir, 'elevation.tif')

# Write
to_geotiff(da, path, compression='deflate')
print(f'Wrote {os.path.getsize(path):,} bytes')

# Read back
loaded = open_geotiff(path)

# Verify round-trip
print(f'Shape: {loaded.shape}')
print(f'CRS:   {loaded.attrs.get("crs")}')
print(f'Name:  {loaded.name}')
print(f'Match: {np.allclose(loaded.values, da.values)}')
```

--------------------------------

### Display benchmark results

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/benchmarks/README.md

View the timing results for the master branch.

```bash
asv show master
```

--------------------------------

### Import necessary libraries for xarray-spatial

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/7_Pathfinding.ipynb

Imports common libraries used in xarray-spatial examples, including numpy, pandas, xarray, and matplotlib for data manipulation and visualization.

```python
import numpy as np
import pandas as pd
import xarray as xr

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch

import xrspatial
```

--------------------------------

### Generate synthetic terrain data

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/user_guide/5_Classification.ipynb

Creates synthetic elevation data using xarray's generate_terrain method for classification examples. The zfactor stretches the elevation range.

```python
W = 800
H = 600
x_range = (-20e6, 20e6)
y_range = (-20e6, 20e6)

terrain = xr.DataArray(np.zeros((H, W)))
terrain = terrain.xrs.generate_terrain(
    x_range=x_range, y_range=y_range, seed=1, zfactor=4000
)
illuminated = terrain.xrs.hillshade()
```

--------------------------------

### Import Libraries for Dask and xarray

Source: https://github.com/xarray-contrib/xarray-spatial/blob/master/examples/dask/distributed_reprojection.ipynb

Imports necessary libraries for Dask distributed computing, xarray data handling, and xarray-spatial functionalities. Ensure Dask and xarray are installed.

```python
import numpy as np
import xarray as xr
import dask
import dask.array as da

from dask.distributed import Client, LocalCluster
from pathlib import Path

import xrspatial
from xrspatial import reproject
```