xarray-spatial 0.9.5

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

# 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

✅ = native backend    🔄 = accepted (CPU fallback)

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

GeoTIFF / COG I/O

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

open_geotiff and to_geotiff auto-dispatch to the correct backend:

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, 'ortho.tif', compression='jpeg')    # JPEG for orthophotos
write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT

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

Compression codecs: Deflate, LZW (Numba JIT), ZSTD, PackBits, JPEG (Pillow), JPEG 2000 (glymur), 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: mean, nearest, min, max, median, mode, cubic
• Planar config, big-endian byte swap, PixelIsArea/PixelIsPoint

Read performance (real-world files, A6000 GPU):

Read performance (synthetic tiled, GPU shines at scale):

Write performance (synthetic tiled):

Reproject / Merge

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.

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.

Reproject performance (reproject-only, 1024x1024, bilinear, vs rioxarray):

Full pipeline (read 3600x3600 Copernicus DEM + reproject to EPSG:3857 + write GeoTIFF):

Merge performance (4 overlapping same-CRS tiles, vs rioxarray):

Same-CRS tiles skip reprojection entirely and are placed by direct coordinate alignment.

Utilities

Surface

Hydrology

Flood

Multispectral

Classification

Focal

Proximity

Zonal

Interpolation

Morphological

Fire

Raster / Vector Conversion

Multivariate

Multi-Criteria Decision Analysis (MCDA)

Pathfinding

Diffusion

Dasymetric

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, matplotlib, zstandard

Optional:

• 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

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. For uncompressed files it reads 5-7x faster than rioxarray; for compressed COGs it is comparable or faster with GPU acceleration.

Citation

Cite this code:

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