xarray-spatial

Xarray-Spatial is a library that implements common raster analysis functions using Numba. It is designed to be easy to install and extend, providing core GIS functions for developers and analysts.

Key characteristics:

• Performance: Uses Numba for high-performance raster analysis.
• Independence: Does not depend on GDAL or GEOS, making it fully extensible in Python, though this limits the breadth of operations compared to the full non-Python geo stack.
• Integration: Grew out of the Datashader project, which provides fast rasterization of vector data (points, lines, polygons, meshes, and rasters) for use with xarray-spatial.

The hydrology module provides high-level wrapper functions for various tasks (flow direction, accumulation, etc.). You select the routing algorithm using the routing keyword argument.

Supported routing algorithms:

• 'd8': Default algorithm.
• 'dinf': D-infinity routing.
• 'mfd': Multiple Flow Direction routing.

While you should use the top-level xrspatial functions for most tasks, the specific implementations for each algorithm are located in the xrspatial.hydro submodule (e.g., xrspatial.hydro.flow_direction_d8).

import xrspatial
# Example using 'dinf' routing
fdir = xrspatial.flow_direction(dem, routing='dinf')
acc = xrspatial.flow_accumulation(fdir, routing='dinf')

When using xarray-spatial with Dask-backed DataArray objects, the library attempts to maintain laziness to support out-of-core processing. However, different algorithms have different memory requirements. Understanding these levels helps you plan pipelines and avoid memory errors.

Laziness Levels

• Fully lazy: The function returns a Dask array without triggering computation. These are safe for arbitrarily large datasets.
• Partially lazy: The function computes small, bounded statistics (like scalars, quartiles, or a ~20K sample) during setup to configure the operation, then returns a Dask array for the main result. The heavy lifting remains lazy.
• Fully materialized: The algorithm requires random access to the full array and calls .compute() internally. These functions will load the necessary data into memory. Use caution with large inputs.

The reproject module splits work between pyproj for metadata and high-performance kernels (Numba/CUDA) for pixel-level operations.

• pyproj (Low Cost): Used once per call for CRS metadata parsing, EPSG code lookups, and estimating the output grid extent via Transformer.transform() on boundary points.
• Numba/CUDA (High Cost/Per-pixel): Handles the heavy lifting including coordinate transforms, resampling (bilinear, nearest-neighbor), datum grid interpolation (NTv2/NADCON), geoid undulation (EGM96/EGM2008), and Helmert datum shifts.
• Note on Cubic Resampling: Currently uses scipy.ndimage.map_coordinates and is restricted to the CPU (no Numba/CUDA kernel available).

All xarray-spatial functions assume the last two dimensions of a DataArray are ordered as (y, x). The library does not verify this via CRS attributes or CF conventions.

If your data is ordered as (x, y), you must transpose it before calling functions.

# Transpose from (x, y) to (y, x)
agg = agg.transpose('y', 'x')

The rasterizer benchmarks compare xarray-spatial (using both numpy and cupy backends) against other libraries including datashader, geocube, and rasterio.

Benchmarks are conducted across several dimensions:

• Geometry Types:
  • Polygons: circles, irregular, rectangles, stars, donuts, multipolygons.
  • Lines: lines, multilines.
  • Points: points, multipoints.
• Feature Counts (n): 50, 200, 1000, and 10000.
• Output Resolutions: 100px to 4000px wide.

Note that for Line and Point types, datashader uses a different API compared to the other libraries.

You are not required to use the provided example datasets. The library accepts any 2D xr.DataArray as input. You can provide data by:

• Creating a synthetic xr.DataArray manually.
• Reading local files or remote URLs (http(s), s3, gs, az) using the GeoTIFF I/O capabilities.
• Using data from other libraries converted to an xr.DataArray format.

All surface functions in xarray-spatial have the following characteristics:

• Output Dtype: All functions output float32 regardless of the input data type.
• Edge Handling: Cells at the edges within the 3x3 kernel radius are assigned NaN by default (controlled by the boundary='nan' parameter).

The merge() method is used to combine overlapping tiles into a single mosaic. Benchmarks comparing xarray-spatial against rioxarray using a 'first' strategy (where overlapping areas are resolved by taking the first available pixel) show:

• Small Tiles (256x256 to 512x512): xarray-spatial is significantly faster.
• Large Tiles (2048x2048): rioxarray performs better for very large merges.

The library uses NaN as the universal nodata sentinel. There is no option to use other values like -9999.

NaN Propagation Rules:

• Slope, aspect, hillshade: Any NaN neighbor produces a NaN output cell.
• Focal mean/std: NaN neighbors are excluded from the calculation; the cell itself is still computed.
• Zonal stats: NaN values are excluded from aggregation. An all-NaN zone returns NaN.
• Hydrology (flow direction) & Pathfinding (A):* NaN cells are treated as impassable barriers.

Edge Handling: When using boundary='nan', edge cells within the kernel radius of the raster boundary default to NaN. Use boundary='nearest' or boundary='reflect' to avoid this.