xarray-spatial

Latest Release		Downloads
License		People
Build Status		Coverage

---

:round_pushpin: Fast, Accurate Python library for Raster Operations

:zap: Extensible with Numba

:fast_forward: Scalable with Dask

:desktop_computer: GPU-accelerated with CuPy and Numba CUDA

:confetti_ball: Free of GDAL / GEOS Dependencies

:earth_africa: General-Purpose Spatial Processing, Geared Towards GIS Professionals

---

Xarray-Spatial is a Python library for raster analysis built on xarray. It has 150+ functions for surface analysis, hydrology (D8, D-infinity, MFD), fire behavior, flood modeling, multispectral indices, proximity, classification, pathfinding, and interpolation. Functions dispatch automatically across four backends (NumPy, Dask, CuPy, Dask+CuPy). A built-in GeoTIFF/COG reader and writer handles raster I/O without GDAL.

Installation

# via pip
pip install xarray-spatial

# with plotting helpers (matplotlib)
pip install xarray-spatial[plot]

# with vector rasterization (shapely): rasterize, polygonize
pip install xarray-spatial[vector]

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

Downloading our starter examples and data

Once you have xarray-spatial installed in your environment, you can use one of the following in your terminal (with the environment active) to download our examples and/or sample data into your local directory.

xrspatial examples : Download the examples notebooks and the data used.

xrspatial copy-examples : Download the examples notebooks but not the data. Note: you won’t be able to run many of the examples.

xrspatial fetch-data : Download just the data and not the notebooks.

In all the above, the command will download and store the files into your current directory inside a folder named ‘xrspatial-examples’.

xarray-spatial grew out of the Datashader project, which provides fast rasterization of vector data (points, lines, polygons, meshes, and rasters) for use with xarray-spatial.

xarray-spatial does not depend on GDAL or GEOS. Raster I/O, reprojection, compression codecs, and coordinate handling are all pure Python and Numba – no C/C++ bindings anywhere in the stack.

API reference docs and 33+ user guide notebooks cover every module.

Raster-huh?

Rasters are regularly gridded datasets like GeoTIFFs, JPGs, and PNGs.

In the GIS world, rasters are used for representing continuous phenomena (e.g. elevation, rainfall, distance), either directly as numerical values, or as RGB images created for humans to view. Rasters typically have two spatial dimensions, but may have any number of other dimensions (time, type of measurement, etc.)

Supported Spatial Functions with Supported Inputs

Each cell shows the feature tier for that function on that backend (see issue #2415). A blank cell means no implementation on that backend; a path that was previously documented as a CPU fallback is reported here as advanced (it works as a documented execution mode, but is not native-parity tested).

✅ stable · 🔼 advanced · 🧪 experimental · 🔧 internal · 🚫 unsupported

GeoTIFF / COG I/O · Classification · Diffusion · Focal · Morphological · Fire · Multispectral · Multivariate · MCDA · Pathfinding · Proximity · Reproject / Merge · Raster / Vector Conversion · Surface · Hydrology · Flood · Interpolation · Dasymetric · Zonal · Templates · Utilities

---

GeoTIFF / COG I/O

Native GeoTIFF and Cloud Optimized GeoTIFF reader/writer. No GDAL required.

VRT is supported as a conservative advanced feature for simple GeoTIFF mosaics, not as a full GDAL VRT replacement. See the VRT support matrix for the supported subset and what is out of scope.

Name	Description	NumPy	Dask	CuPy GPU	Dask+CuPy GPU	Cloud
open_geotiff	Read GeoTIFF / COG / VRT	✅	✅	🧪	🧪	🔼
to_geotiff	Write DataArray as GeoTIFF / COG / VRT	✅	✅	🧪	🧪	🔼

open_geotiff and to_geotiff select the backend from their parameters (gpu=, chunks=, .vrt path); GPU read/write is reached with gpu=True, not a separate function:

from xrspatial.geotiff import open_geotiff, to_geotiff

open_geotiff('dem.tif')                              # NumPy
open_geotiff('dem.tif', chunks=512)                  # Dask
open_geotiff('dem.tif', gpu=True)                    # CuPy (nvCOMP + GDS)
open_geotiff('dem.tif', gpu=True, chunks=512)        # Dask + CuPy
open_geotiff('https://example.com/cog.tif')          # HTTP COG
open_geotiff('s3://bucket/dem.tif')                  # Cloud (S3/GCS/Azure)
open_geotiff('mosaic.vrt')                           # VRT mosaic (auto-detected)

to_geotiff(cupy_array, 'out.tif')                    # auto-detects GPU
to_geotiff(data, 'out.tif', gpu=True)                # force GPU compress
to_geotiff(data, 'out.tif', compression='zstd')      # ZSTD for smaller files
to_geotiff(data, 'cog.tif', cog=True)                # COG with auto overviews
to_geotiff(data, 'cog.tif', cog=True,                # COG with explicit levels
           overview_levels=[2, 4, 8],
           overview_resampling='nearest')
to_geotiff(data, 'mosaic.vrt')                       # write a tiled VRT mosaic

open_geotiff('dem.tif', dtype='float32')             # half memory
open_geotiff('dem.tif', dtype='float32', chunks=512) # Dask + half memory
to_geotiff(data, 'out.tif', compression_level=1)     # fast scratch write
to_geotiff(data, 'out.tif', compression_level=22)    # max compression
to_geotiff(dask_da, 'out.tif')                       # stream Dask to single TIFF
to_geotiff(dask_da, 'mosaic.vrt')                    # stream Dask to VRT

# Accessor methods
da.xrs.to_geotiff('out.tif', compression='lzw')     # write from DataArray
ds.xrs.open_geotiff('large_dem.tif')                 # read windowed to Dataset extent

# xarray backend engine
import xarray as xr
xr.open_dataset('dem.tif', engine='xrspatial')   # open as a Dataset
xr.open_mfdataset('*.tif', engine='xrspatial',   # share one var name
                  backend_kwargs={'default_name': 'band_data'})

Compression codecs: Deflate, LZW (Numba JIT), ZSTD, PackBits, JPEG (Pillow, internal-only: requires allow_internal_only_jpeg=True and is not readable by GDAL), JPEG 2000 (glymur, experimental: requires allow_experimental_codecs=True), uncompressed

GPU codecs: Deflate and ZSTD via nvCOMP batch API; JPEG 2000 via nvJPEG2000; LZW via Numba CUDA kernels

Features:

• Tiled, stripped, BigTIFF, multi-band (RGB/RGBA), sub-byte (1/2/4/12-bit)
• Predictors: horizontal differencing (pred=2), floating-point (pred=3)
• GeoKeys: EPSG, WKT/PROJ (via pyproj), citations, units, ellipsoid, vertical CRS
• Metadata: nodata masking, palette colormaps, DPI/resolution, GDALMetadata XML, arbitrary tag preservation
• Cloud storage: S3 (s3://), GCS (gs://), Azure (az://) via fsspec
• GPUDirect Storage: SSD→GPU direct DMA via KvikIO (optional)
• Thread-safe mmap reads, atomic writes, HTTP connection reuse (urllib3)
• Overview generation (CPU and GPU): mean, nearest, min, max, median, mode, cubic
• Planar config, big-endian byte swap, PixelIsArea/PixelIsPoint

Consistency: 100% pixel-exact match vs rioxarray on all tested files (Landsat 8, Copernicus DEM, USGS 1-arc-second, USGS 1-meter).

---

Reproject / Merge / Resample

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Resample	Changes raster resolution (cell size) without reprojection. Nearest, bilinear, cubic, average, mode, min, max, median methods	Standard (interpolation / block aggregation)	✅	🔼	🔼	🔼
Reproject	Reprojects a raster to a new CRS with Numba JIT / CUDA coordinate transforms and resampling. Supports vertical datums (EGM96, EGM2008) and horizontal datum shifts (NAD27, OSGB36, etc.)	Standard (inverse mapping)	✅	🔼	🔼	🔼
Merge	Merges multiple rasters into a single mosaic with configurable overlap strategy. Same-CRS tiles skip reprojection and are placed by direct coordinate alignment	Standard (mosaic)	🔼	🔼	🔼	🔼

Built-in Numba JIT and CUDA projection kernels bypass pyproj for per-pixel coordinate transforms. pyproj is used only for CRS metadata parsing (~1ms, once per call) and output grid boundary estimation (~500 control points, once per call). Any CRS pair without a built-in kernel falls back to pyproj automatically.

Projection	EPSG examples	CPU Numba	CUDA GPU
Web Mercator	3857	✅️	✅️
UTM / Transverse Mercator	326xx, 327xx, State Plane	✅️	✅️
Ellipsoidal Mercator	3395	✅️	✅️
Lambert Conformal Conic	2154, 2229, State Plane	✅️	✅️
Albers Equal Area	5070	✅️	✅️
Cylindrical Equal Area	6933	✅️	✅️
Sinusoidal	MODIS grids	✅️	✅️
Lambert Azimuthal Equal Area	3035, 6931, 6932	✅️	✅️
Polar Stereographic	3031, 3413, 3996	✅️	✅️
Oblique Stereographic	custom WGS84	✅️	pyproj fallback
Oblique Mercator (Hotine)	3375 (RSO)	implemented, disabled	pyproj fallback

Vertical datum support: geoid_height, ellipsoidal_to_orthometric, orthometric_to_ellipsoidal convert between ellipsoidal (GPS) and orthometric (map/MSL) heights using EGM96 (vendored, 2.6MB) or EGM2008 (77MB, downloaded on first use). Reproject can apply vertical shifts during reprojection via the vertical_crs parameter.

Datum shift support: Reprojection from non-WGS84 datums (NAD27, OSGB36, DHDN, MGI, ED50, BD72, CH1903, D73, AGD66, Tokyo) applies grid-based shifts from PROJ CDN (sub-metre accuracy) with 7-parameter Helmert fallback (1-5m accuracy). 14 grids are registered covering North America, UK, Germany, Austria, Spain, Netherlands, Belgium, Switzerland, Portugal, and Australia.

ITRF frame support: itrf_transform converts between ITRF2000, ITRF2008, ITRF2014, and ITRF2020 using 14-parameter time-dependent Helmert transforms from PROJ data files. Shifts are mm-level.

---

Utilities

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Preview	Downsamples a raster to target pixel dimensions for visualization	Custom	✅	🔼	🔼	🔼
Rescale	Min-max normalization to a target range (default [0, 1])	Standard	✅	🔼	🔼	🔼
Standardize	Z-score normalization (subtract mean, divide by std)	Standard	✅	🔼	🔼	🔼
rechunk_no_shuffle	Rechunk dask arrays using whole-chunk multiples (no shuffle)	Custom	🔼	🔼	🔼	🔼
fused_overlap	Fuse sequential map_overlap calls into a single pass	Custom	🔼	🔼	🔼	🔼
multi_overlap	Run multi-output kernel in a single overlap pass	Custom	🔼	🔼	🔼	🔼
validate	Check a raster against the xarray-spatial input contract (.xrs.validate())	Custom	✅	🔼	🔼	🔼

---

Templates

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
from_template	Empty study-area grid for a named region (CONUS, NYC, …), a world city (London, Tokyo, … in its UTM zone), a country code, or a whole-world projection (web_mercator, wgs84/latlon, equal_earth); preserve='area'/'shape' picks an EPSG projection by property; list_templates() lists every accepted name	Custom	✅	🔼	🔼	🔼

---

Surface

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Aspect	Computes downslope direction of each cell in degrees	Horn 1981	✅	🔼	🔼	🔼
Northness	North-south component of aspect: cos(aspect) for linear models	Stage 1976	✅	🔼	🔼	🔼
Eastness	East-west component of aspect: sin(aspect) for linear models	Stage 1976	✅	🔼	🔼	🔼
Curvature	Measures rate of slope change (concavity/convexity) at each cell	Zevenbergen & Thorne 1987	✅	🔼	🔼	🔼
Hillshade	Simulates terrain illumination from a given sun angle and azimuth	GDAL gdaldem	✅	🔼	🔼	🔼
Roughness	Computes local relief as max minus min elevation in a 3×3 window	GDAL gdaldem	✅	🔼	🔼	🔼
Sky-View Factor	Measures the fraction of visible sky hemisphere at each cell	Zakek et al. 2011	✅	🔼	🔼	🔼
Slope	Computes terrain gradient steepness at each cell in degrees	Horn 1981	✅	🔼	🔼	🔼
Terrain Generation	Generates synthetic terrain from fBm or ridged fractal noise with optional domain warping, Worley blending, and hydraulic erosion	Custom (fBm)	🔼	🔼	🔼	🔼
TPI	Computes Topographic Position Index (center minus mean of neighbors)	Weiss 2001	✅	🔼	🔼	🔼
TRI	Computes Terrain Ruggedness Index (local elevation variation)	Riley et al. 1999	✅	🔼	🔼	🔼
Landforms	Classifies terrain into 10 landform types using the Weiss (2001) TPI scheme	Weiss 2001	✅	🔼	🔼	🔼
Viewshed	Determines visible cells from a given observer point on terrain	GRASS GIS r.viewshed	🔼	🔼	🔼	🔼
Cumulative Viewshed	Counts how many observers can see each cell	Custom	🔼	🔼	🔼	🔼
Visibility Frequency	Fraction of observers with line-of-sight to each cell	Custom	🔼	🔼	🔼	🔼
Line of Sight	Elevation profile and visibility along a point-to-point transect	Custom	🔼	🔼	🔼	🔼
Min Observable Height	Finds the minimum observer height needed to see each cell	Custom	🧪
Perlin Noise	Generates smooth continuous random noise for procedural textures	Perlin 1985	✅	🔼	🔼	🔼
Worley Noise	Generates cellular (Voronoi) noise returning distance to the nearest feature point	Worley 1996	✅	🔼	🔼	🔼
Hydraulic Erosion	Simulates particle-based water erosion to carve valleys and deposit sediment	Custom	🔼	🔼	🔼	🔼
Bump Mapping	Adds randomized bump features to simulate natural terrain variation	Custom	🔼	🔼	🔼	🔼

---

Hydrology

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Flow Direction	Direction of steepest descent out of each cell (D8 · Dinf · MFD via routing=)	O’Callaghan & Mark 1984; Tarboton 1997; Qin et al. 2007	✅	🔼	🔼	🔼
Flow Accumulation	Upstream cells or area draining through each cell (D8 · Dinf · MFD via routing=)	Jenson & Domingue 1988; Tarboton 1997; Qin et al. 2007	✅	🔼	🔼	🔼
Flow Length	Flow path length to the outlet or from the divide (D8 · Dinf · MFD via routing=)	Tarboton 1997; Qin et al. 2007	🔼	🔼	🔼	🔼
Flow Path	Traces downstream flow paths from start points (D8 · Dinf · MFD via routing=)	Tarboton 1997; Qin et al. 2007	🔼	🔼	🔼	🔼
Watershed	Labels each cell with the pour point it drains to (D8 · Dinf · MFD via routing=)	Tarboton 1997; Qin et al. 2007	✅	🔼	🔼	🔼
Stream Link	Assigns unique IDs to stream segments above a threshold (D8 · Dinf · MFD via routing=)	Tarboton 1997; Freeman 1991	✅	🔼	🔼	🔼
Stream Order	Strahler or Shreve stream ordering of the network (D8 · Dinf · MFD via routing=)	Strahler 1957, Shreve 1966	✅	🔼	🔼	🔼
HAND	Height Above Nearest Drainage (D8 · Dinf · MFD via routing=)	Nobre et al. 2011	✅	🔼	🔼	🔼
Fill	Fills depressions in a DEM using Planchon-Darboux iterative flooding (D8)	Planchon & Darboux 2002	✅	🔼	🔼	🔼
Sink	Identifies and labels depression cells (D8)	Standard (D8 tracing)	✅	🔼	🔼	🔼
Basin	Labels each cell with the outlet of the basin it drains to (D8)	Standard (D8 tracing)	✅	🔼	🔼	🔼
Snap Pour Point	Snaps pour points to the highest-accumulation cell within a search radius (D8)	Custom	✅	🔼	🔼	🔼
TWI	Topographic Wetness Index: ln(specific catchment area / tan(slope)) (D8)	Beven & Kirkby 1979	✅	🔼	🔼	🔼

---

Flood

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Flood Depth	Computes water depth above terrain from a HAND raster and water level	Standard (HAND-based)	✅	🔼	🔼	🔼
Inundation	Produces a binary flood/no-flood mask from a HAND raster and water level	Standard (HAND-based)	✅	🔼	🔼	🔼
Curve Number Runoff	Estimates runoff depth from rainfall using the SCS/NRCS curve number method	SCS/NRCS	✅	🔼	🔼	🔼
Travel Time	Estimates overland flow travel time via simplified Manning’s equation	Manning 1891	✅	🔼	🔼	🔼
Vegetation Roughness	Derives Manning’s roughness coefficients from NLCD land cover or NDVI	SCS/NRCS	✅	🔼	🔼	🔼
Vegetation Curve Number	Derives SCS curve numbers from land cover and hydrologic soil group	SCS/NRCS	✅	🔼	🔼	🔼
Flood Depth (Vegetation)	Manning-based steady-state flow depth incorporating vegetation roughness	Manning 1891	✅	🔼	🔼	🔼

---

Multispectral

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Atmospherically Resistant Vegetation Index (ARVI)	Vegetation index resistant to atmospheric effects using blue band correction	Kaufman & Tanre 1992	✅	🔼	🔼	🔼
Burn Area Index (BAI)	Spectral distance to charcoal reflectance point for burn scar detection	Chuvieco et al. 2002	✅	🔼	🔼	🔼
Enhanced Built-Up and Bareness Index (EBBI)	Highlights built-up areas and barren land from thermal and SWIR bands	As-syakur et al. 2012	✅	🔼	🔼	🔼
Enhanced Vegetation Index (EVI)	Enhanced vegetation index reducing soil and atmospheric noise	Huete et al. 2002	✅	🔼	🔼	🔼
Green Chlorophyll Index (GCI)	Estimates leaf chlorophyll content from green and NIR reflectance	Gitelson et al. 2003	✅	🔼	🔼	🔼
Modified Soil Adjusted Vegetation Index (MSAVI2)	Self-adjusting soil line vegetation index, no L parameter needed	Qi et al. 1994	✅	🔼	🔼	🔼
Normalized Burn Ratio (NBR)	Measures burn severity using NIR and SWIR band difference	USGS Landsat	✅	🔼	🔼	🔼
Normalized Burn Ratio 2 (NBR2)	Refines burn severity mapping using two SWIR bands	USGS Landsat	✅	🔼	🔼	🔼
Normalized Difference Built-up Index (NDBI)	Picks out built-up and urban areas from SWIR and NIR bands	Zha et al. 2003	✅	🔼	🔼	🔼
Normalized Difference Moisture Index (NDMI)	Detects vegetation moisture stress from NIR and SWIR reflectance	USGS Landsat	✅	🔼	🔼	🔼
Normalized Difference Snow Index (NDSI)	Separates snow and ice from clouds using green and SWIR bands	Hall et al. 1995	✅	🔼	🔼	🔼
Normalized Difference Water Index (NDWI)	Maps open water bodies using green and NIR band difference	McFeeters 1996	✅	🔼	🔼	🔼
Modified Normalized Difference Water Index (MNDWI)	Detects water in urban areas using green and SWIR bands	Xu 2006	✅	🔼	🔼	🔼
Normalized Difference Vegetation Index (NDVI)	Quantifies vegetation density from red and NIR band difference	Rouse et al. 1973	✅	🔼	🔼	🔼
Optimized Soil Adjusted Vegetation Index (OSAVI)	SAVI with fixed L=0.16, tuned for sparse vegetation	Rondeaux et al. 1996	✅	🔼	🔼	🔼
Soil Adjusted Vegetation Index (SAVI)	Vegetation index with soil brightness correction factor	Huete 1988	✅	🔼	🔼	🔼
Structure Insensitive Pigment Index (SIPI)	Estimates carotenoid-to-chlorophyll ratio for plant stress detection	Penuelas et al. 1995	✅	🔼	🔼	🔼
True Color	Composites red, green, and blue bands into a natural color image	Standard	✅	🔼	🔼	🔼

For a broader catalog of spectral indices and sensor-specific band combinations, see awesome-spectral-indices and its companion xarray library spyndex.

---

Classification

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Binary	Binarizes values by membership in a target set (1 if match, 0 otherwise)	Standard	✅	🔼	🔼	🔼
Box Plot	Classifies values into bins based on box plot quartile boundaries	PySAL mapclassify	✅	🔼	🔼	🔼
Equal Interval	Divides the value range into equal-width bins	PySAL mapclassify	✅	🔼	🔼	🔼
Head/Tail Breaks	Classifies heavy-tailed distributions using recursive mean splitting	PySAL mapclassify	✅	🔼	🔼	🔼
Maximum Breaks	Finds natural groupings by maximizing differences between sorted values	PySAL mapclassify	✅	🔼	🔼	🔼
Natural Breaks	Optimizes class boundaries to minimize within-class variance (Jenks)	Jenks 1967, PySAL	✅	🔼	🔼	🔼
Percentiles	Assigns classes based on user-defined percentile breakpoints	PySAL mapclassify	✅	🔼	🔼	🔼
Quantile	Distributes values into classes with equal observation counts	PySAL mapclassify	✅	🔼	🔼	🔼
Reclassify	Remaps pixel values to new classes using a user-defined lookup	PySAL mapclassify	✅	🔼	🔼	🔼
Std Mean	Classifies values by standard deviation intervals from the mean	PySAL mapclassify	✅	🔼	🔼	🔼

---

Focal

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Apply	Applies a custom function over a sliding neighborhood window	Standard	🔼	🔼	🔼	🔼
Hotspots	Identifies statistically significant spatial clusters using Getis-Ord Gi*	Getis & Ord 1992	✅	🔼	🔼	🔼
Emerging Hotspots	Classifies time-series hot/cold spot trends using Gi* and Mann-Kendall	Getis & Ord 1992, Mann 1945	🔼	🔼	🔼	🔼
Mean	Computes the mean value within a sliding neighborhood window	Standard	✅	🔼	🔼	🔼
Focal Statistics	Computes summary statistics over a sliding neighborhood window	Standard	✅	🔼	🔼	🔼
Bilateral	Feature-preserving smoothing via bilateral filtering	Tomasi & Manduchi 1998	🔼	🔼	🔼	🔼
GLCM Texture	Computes Haralick GLCM texture features over a sliding window	Haralick et al. 1973	🔼	🔼	🔼	🔼
Sobel X	Horizontal gradient via Sobel operator (detects vertical edges)	Sobel & Feldman 1968	🔼	🔼	🔼	🔼
Sobel Y	Vertical gradient via Sobel operator (detects horizontal edges)	Sobel & Feldman 1968	🔼	🔼	🔼	🔼
Laplacian	Omnidirectional second-derivative edge detector	Standard	🔼	🔼	🔼	🔼
Prewitt X	Horizontal gradient via Prewitt operator (detects vertical edges)	Prewitt 1970	🔼	🔼	🔼	🔼
Prewitt Y	Vertical gradient via Prewitt operator (detects horizontal edges)	Prewitt 1970	🔼	🔼	🔼	🔼

---

Proximity

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Allocation	Assigns each cell to the identity of the nearest source feature	Standard (Dijkstra)	✅	🔼	🔼	🔼
Balanced Allocation	Partitions a cost surface into territories of roughly equal cost-weighted area	Custom	🔼	🔼	🔼	🔼
Cost Distance	Computes minimum accumulated cost to the nearest source through a friction surface	Standard (Dijkstra)	✅	🔼	🔼	🔼
Least-Cost Corridor	Identifies zones of low cumulative cost between two source locations	Standard (Dijkstra)	🔼	🔼	🔼	🔼
Direction	Computes the direction from each cell to the nearest source feature	Standard	✅	🔼	🔼	🔼
Proximity	Computes the distance from each cell to the nearest source feature	Standard	✅	🔼	🔼	🔼
Surface Distance	Computes distance along the 3D terrain surface to the nearest source	Standard (Dijkstra)	🔼	🔼	🔼	🔼
Surface Allocation	Assigns each cell to the nearest source by terrain surface distance	Standard (Dijkstra)	🔼	🔼	🔼	🔼
Surface Direction	Computes compass direction to the nearest source by terrain surface distance	Standard (Dijkstra)	🔼	🔼	🔼	🔼

---

Zonal

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Apply	Applies a custom function to each zone in a classified raster	Standard	🔼	🔼	🔼	🔼
Clip Polygon	Clips a raster to an arbitrary polygon with masking	Standard	🔼	🔼	🔼	🔼
Crop	Extracts the bounding rectangle of a specific zone	Standard	✅	🔼	🔼	🔼
Regions	Identifies connected regions of non-zero cells	Standard (CCL)	✅	🔼	🔼	🔼
Trim	Removes nodata border rows and columns from a raster	Standard	✅	🔼	🔼	🔼
Zonal Statistics	Computes summary statistics for a value raster within each zone	Standard	✅	🔼	🔼	🔼
Zonal Cross Tabulate	Cross-tabulates agreement between two categorical rasters	Standard	🔼	🔼	🔼	🔼

---

Interpolation

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
IDW	Inverse Distance Weighting from scattered points (arrays or a GeoDataFrame) to a raster grid	Standard (IDW)	✅	🔼	🔼	🔼
Kriging	Ordinary Kriging with automatic variogram fitting (spherical, exponential, gaussian); accepts arrays or a GeoDataFrame	Standard (ordinary kriging)	🔼	🔼	🔼	🔼
Spline	Thin Plate Spline interpolation with optional smoothing; accepts arrays or a GeoDataFrame	Standard (TPS)	🔼	🔼	🔼	🔼

---

Morphological

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Erode	Morphological erosion (local minimum over structuring element)	Standard (morphology)	✅	🔼	🔼	🔼
Dilate	Morphological dilation (local maximum over structuring element)	Standard (morphology)	✅	🔼	🔼	🔼
Opening	Erosion then dilation (removes small bright features)	Standard (morphology)	✅	🔼	🔼	🔼
Closing	Dilation then erosion (fills small dark gaps)	Standard (morphology)	✅	🔼	🔼	🔼
Gradient	Dilation minus erosion (edge detection)	Standard (morphology)	✅	🔼	🔼	🔼
White Top-hat	Original minus opening (isolate bright features)	Standard (morphology)	✅	🔼	🔼	🔼
Black Top-hat	Closing minus original (isolate dark features)	Standard (morphology)	✅	🔼	🔼	🔼
Sieve	Remove small connected clumps from classified rasters	GDAL sieve	🔼	🔼	🔼	🔼

---

Fire

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
dNBR	Differenced Normalized Burn Ratio (pre minus post NBR)	USGS	✅	🔼	🔼	🔼
RdNBR	Relative dNBR normalized by pre-fire vegetation density	USGS	✅	🔼	🔼	🔼
Burn Severity Class	USGS 7-class burn severity from dNBR thresholds	USGS	✅	🔼	🔼	🔼
Fireline Intensity	Byram’s fireline intensity from fuel load and spread rate (kW/m)	Byram 1959	✅	🔼	🔼	🔼
Flame Length	Flame length derived from fireline intensity (m)	Byram 1959	✅	🔼	🔼	🔼
Rate of Spread	Simplified Rothermel spread rate with Anderson 13 fuel models (m/min)	Rothermel 1972, Anderson 1982	✅	🔼	🔼	🔼
KBDI	Keetch-Byram Drought Index single time-step update (0-800 mm)	Keetch & Byram 1968	✅	🔼	🔼	🔼

---

Raster / Vector Conversion

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Polygonize	Converts contiguous regions of equal value into vector polygons	Standard (CCL)	🔼	🔼	🔼	🔼
Contours	Extracts elevation contour lines (isolines) from a raster surface	Standard (marching squares)	🔼	🔼	🔼	🔼
Rasterize	Rasterizes vector geometries (polygons, lines, points) from a GeoDataFrame	Standard (scanline, Bresenham)	🔼		🔼

---

Kernel Density Estimation

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
KDE	Point-to-raster kernel density estimation (Gaussian, Epanechnikov, quartic); accepts arrays or a GeoDataFrame	Silverman 1986	🔼	🔼	🔼	🔼
Line Density	Line-segment-to-raster density estimation	Standard	🔼

---

Multivariate

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Mahalanobis Distance	Measures statistical distance from a multi-band reference distribution, accounting for band correlations	Mahalanobis 1936	✅	🔼	🔼	🔼

---

Multi-Criteria Decision Analysis (MCDA)

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Standardize	Converts criterion rasters to 0-1 suitability scale (linear, sigmoidal, gaussian, triangular, piecewise, categorical)	Standard	✅	🔼	🔼	🔼
AHP Weights	Derives criterion weights from pairwise comparisons using the Saaty eigenvector method with consistency ratio	Saaty 1980	✅	🚫	🚫	🚫
Rank Weights	Derives weights from a rank ordering (ROC, rank sum, reciprocal)	Standard	✅	🚫	🚫	🚫
WLC	Weighted Linear Combination (fully compensatory weighted sum)	Malczewski 2006	✅	🔼	🧪	🧪
WPM	Weighted Product Model (multiplicative, penalizes low scores)	Standard	✅	🔼	🧪	🧪
OWA	Ordered Weighted Averaging with tunable risk attitude	Yager 1988	✅	🔼	🧪	🧪
Fuzzy Overlay	Combines criteria using fuzzy set operators (AND, OR, sum, product, gamma)	Eastman 1999	✅	🔼	🧪	🧪
Boolean Overlay	Combines binary criterion masks using AND/OR logic	Standard	✅	🔼	🧪	🧪
Constrain	Masks exclusion zones from a suitability surface	Standard	✅	🔼	🧪	🧪
Sensitivity	Assesses weight stability via one-at-a-time or Monte Carlo perturbation	Standard	✅	🔼	🧪	🧪

---

Pathfinding

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
A* Pathfinding	Finds the least-cost path between two cells on a cost surface	Hart et al. 1968	✅	🔼	🔼	🔼
Multi-Stop Search	Routes through N waypoints in sequence, with optional TSP reordering	Custom	🔼	🔼	🔼	🔼

---

Diffusion

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Diffuse	Runs explicit forward-Euler diffusion on a 2D scalar field	Standard (heat equation)	✅	🔼	🔼	🔼

---

Dasymetric

Name	Description	Source	NumPy xr.DataArray	Dask xr.DataArray	CuPy GPU xr.DataArray	Dask GPU xr.DataArray
Disaggregate	Redistributes zonal totals to pixels using an ancillary weight surface	Mennis 2003	🔼	🔼	🔼	🔼
Pycnophylactic	Tobler’s pycnophylactic interpolation preserving zone totals via Laplacian smoothing	Tobler 1979	✅	🚫	🔼	🚫
Validate Disaggregation	Checks that disaggregated pixel sums match the original zone totals	Standard	✅	🔼	🔼	🔼

---

Usage

Quick Start

Importing xrspatial registers an .xrs accessor on DataArrays and Datasets, giving you tab-completable access to every spatial operation:

import xrspatial as xrs
from xrspatial.geotiff import open_geotiff, to_geotiff

# Read a GeoTIFF (no GDAL required)
elevation = open_geotiff('dem.tif')

# Surface analysis
slope = elevation.xrs.slope()
hillshaded = elevation.xrs.hillshade(azimuth=315, angle_altitude=45)
aspect = elevation.xrs.aspect()

# Reproject and write as a Cloud Optimized GeoTIFF
dem_wgs84 = elevation.xrs.reproject(target_crs='EPSG:4326')
to_geotiff(dem_wgs84, 'output.tif', cog=True)

# Classification
classes = elevation.xrs.equal_interval(k=5)
breaks = elevation.xrs.natural_breaks(k=10)

# Proximity
distance = elevation.xrs.proximity(target_values=[1])

# Multispectral
vegetation = nir.xrs.ndvi(red)
enhanced_vi = nir.xrs.evi(red, blue)

Dataset Support

The .xrs accessor works on Datasets too. Single-input functions apply the operation to each data variable. Multi-input functions (multispectral indices) accept string kwargs that map band aliases to variable names:

ds = xr.Dataset({'band_4': red, 'band_5': nir})

# Single-input: slope computed for each variable
slope_ds = ds.xrs.slope()

# Multi-input: map variable names to band parameters
ndvi_result = ds.xrs.ndvi(nir='band_5', red='band_4')

Function Import Style

All operations are also available as standalone functions:

import xrspatial as xrs

hillshaded = xrs.hillshade(elevation)
slope_result = xrs.slope(elevation)
vegetation = xrs.ndvi(nir, red)

Check out the user guide here.

---

Dependencies

Core: numpy, numba, scipy, xarray, zstandard

Optional:

• matplotlib — the .xrs.plot accessor helpers (pip install xarray-spatial[plot])
• shapely — the vector-to-raster paths, rasterize and polygonize (pip install xarray-spatial[vector])
• pyproj — WKT/PROJ CRS resolution
• cupy — GPU acceleration
• dask — out-of-core processing
• libnvcomp — GPU batch decompression (deflate, ZSTD)
• kvikio — GPUDirect Storage (SSD → GPU)
• fsspec + s3fs/gcsfs/adlfs — cloud storage

libnvcomp and kvikio are not pulled in by the gpu extra. They are runtime dependencies of the GeoTIFF GPU read path and must be installed separately (typically via conda from the rapidsai/nvidia channels), since libnvcomp ships as a system library and kvikio requires a matching CUDA toolkit.

Notes on GDAL

xarray-spatial does not depend on GDAL. The built-in GeoTIFF/COG reader and writer (xrspatial.geotiff) handles raster I/O natively using only numpy, numba, and the standard library. This means:

• Zero GDAL installation hassle. pip install xarray-spatial gets you everything needed to read and write GeoTIFFs, COGs, and VRT files.
• Pure Python, fully extensible. All codec, header parsing, and metadata code is readable Python/Numba, not wrapped C/C++.
• GPU-accelerated reads. With optional nvCOMP and nvJPEG2000, compressed tiles decompress directly on the GPU via CUDA – something GDAL cannot do.

The native reader is pixel-exact against rasterio/GDAL across Landsat 8, Copernicus DEM, USGS 1-arc-second, and USGS 1-meter DEMs.

Citation

Cite this code:

xarray-contrib/xarray-spatial, https://github.com/xarray-contrib/xarray-spatial, ©2020-2026.