Gridded Raster Fields with xarray and rioxarray #

An interpolated environmental surface is worthless the moment it leaves memory as a bare NumPy array. Strip away the coordinate axes, the coordinate reference system, and the nodata convention, and a 500×500 grid of PM2.5 values is just numbers — it cannot be reprojected, overlaid on a basemap, clipped to a watershed, or opened in QGIS without a human guessing its extent and projection. The gap between a scientifically correct interpolation and a usable geospatial product is entirely a labelling and serialization problem: attaching real-world x/y coordinates, pinning a CRS, encoding which cells are “no measurement,” and writing a format that downstream tools trust. xarray supplies the labelled n-dimensional container and rioxarray supplies the geospatial accessor that bridges it to GDAL. Together they turn a raw interpolation output into a CRS-aware raster field, and they are the raster half of the broader Geospatial Data Storage, Interpolation & GIS Export pipeline.

This page covers the full path from a bare 2-D array to a validated GeoTIFF or Zarr store: constructing the coordinate grid, wrapping it in a DataArray, assigning CRS and nodata through the rio accessor, and writing tiled output that reads efficiently at scale.

Prerequisites #

These are the pins and upstream conditions that make raster serialization deterministic. Version skew in the GDAL stack is the single most common source of “works on my machine” raster bugs, so pin exactly.

• Python 3.11 with the geospatial stack pinned: numpy==1.26.4, xarray==2024.3.0, rioxarray==0.15.0, rasterio==1.3.9, pyproj==3.6.1. Put these in a comment at the top of every code block. rioxarray binds to the GDAL that ships inside the rasterio wheel; mixing a system GDAL with a wheel GDAL produces silent CRS-parsing failures.
• A regular, monotonic grid. Both coordinate axes must be evenly spaced. rioxarray derives the affine geotransform from np.diff of the coordinates and refuses to write a valid GeoTIFF if the spacing is irregular. Irregular sensor spacing is fine as input to interpolation, but the output grid must be regular.
• Interpolated values, not raw points. The 2-D value array is the product of an upstream interpolation step. If you are still producing the surface itself, build it first with spatial interpolation using kriging and IDW for sensor fields — this page assumes you already hold a filled ny × nx array plus its extent and cell size.
• A known CRS for the grid. The interpolation should have been performed in a projected CRS (meters), not in raw WGS 84 degrees, so that distance-weighted methods behave isotropically. Record that EPSG code; you will stamp it onto the raster. If your points arrived in a different projection, resolve that during ingestion CRS mapping before interpolating.
• A nodata convention. Decide up front whether unsampled or masked cells are np.nan (float rasters, recommended) or an integer sentinel. Mixing conventions across a tile set breaks mosaicking.

---

Step-by-Step Workflow #

Step 1 — Construct the Coordinate Grid #

The coordinate vectors are the contract between array indices and real-world position. GeoTIFF geotransforms reference the top-left corner of the top-left pixel, but xarray coordinates label pixel centers. Build center coordinates and let rioxarray handle the half-cell shift — placing coordinates on cell edges is the usual cause of a half-pixel offset in QGIS.

Note the y axis runs north-to-south (descending) so the first raster row is the northernmost, which is what every GeoTIFF reader expects.

# python 3.11 · numpy==1.26.4 · xarray==2024.3.0 · rioxarray==0.15.0 · rasterio==1.3.9 · pyproj==3.6.1
from __future__ import annotations
import numpy as np


def build_grid_axes(
    xmin: float, ymin: float, xmax: float, ymax: float, cell_size: float,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Build pixel-center coordinate vectors for a regular raster grid.

    Returns (y, x) where y is DESCENDING (north first) so the written raster
    is not vertically flipped.

    Time complexity:  O(nx + ny).
    Space complexity: O(nx + ny) — axes only, not the full grid.
    """
    nx = int(round((xmax - xmin) / cell_size))
    ny = int(round((ymax - ymin) / cell_size))
    # pixel-CENTER coordinates: first center is half a cell in from the edge
    x = xmin + (np.arange(nx) + 0.5) * cell_size
    y = ymax - (np.arange(ny) + 0.5) * cell_size   # descending: north -> south
    return y, x

Parameter notes: cell_size is in the units of the projected CRS (meters for UTM). Choosing it is a resolution-versus-file-size trade-off covered in the tuning table below. Complexity: O(nx + ny) to build axes; the grid itself is O(nx · ny).

---

Step 2 — Wrap Values in a Labelled DataArray #

An xarray.DataArray is a NumPy array plus named dimensions and coordinates. The dimension order for a single-band raster must be (y, x) — rioxarray treats the last two dims as spatial and the y-then-x order matches GDAL’s row-major raster layout. Attach descriptive attributes now; they survive into the GeoTIFF metadata and Zarr .zattrs and save every downstream consumer from guessing units.

# python 3.11 · numpy==1.26.4 · xarray==2024.3.0
import xarray as xr


def wrap_as_dataarray(
    values: np.ndarray,        # shape (ny, nx), float32 recommended
    y: np.ndarray,             # descending coordinate vector
    x: np.ndarray,             # ascending coordinate vector
    var_name: str = "pm25",
    units: str = "ug/m3",
) -> xr.DataArray:
    """
    Assemble an interpolated 2-D array into a labelled DataArray.

    Time complexity:  O(1) — wraps the existing buffer, no copy of values.
    Space complexity: O(ny * nx) — references the values array.
    """
    if values.shape != (y.size, x.size):
        raise ValueError(
            f"values shape {values.shape} != (ny, nx) ({y.size}, {x.size})"
        )
    da = xr.DataArray(
        data=values.astype("float32"),
        dims=("y", "x"),
        coords={"y": y, "x": x},
        name=var_name,
        attrs={"units": units, "long_name": f"interpolated {var_name} surface"},
    )
    return da

Parameter notes: casting to float32 halves file size versus float64 at a precision cost that is negligible for environmental concentrations (well below sensor accuracy). Keep the raw values buffer contiguous — xr.DataArray wraps it without copying, so an upstream view can leak.

---

Step 3 — Assign CRS and Nodata Through the rio Accessor #

This is the step that makes the array geospatial. rio.write_crs stamps the projection so any reader can place the grid on Earth, and rio.write_nodata encodes which cells carry no measurement. Do both before writing — a GeoTIFF without a CRS opens as an “unknown” layer in QGIS, and one without a nodata tag renders masked cells as spurious zeros or black.

# python 3.11 · rioxarray==0.15.0 · pyproj==3.6.1
import numpy as np
import rioxarray  # noqa: F401 — registers the .rio accessor on xarray objects


def georeference(da: xr.DataArray, epsg: int) -> xr.DataArray:
    """
    Attach a CRS and a NaN nodata convention to a labelled DataArray.

    `write_crs` returns a NEW object (rioxarray is functional here), so always
    rebind. `inplace=False` is the default and the safe choice.

    Time complexity:  O(1) — metadata only.
    """
    da = da.rio.write_crs(f"EPSG:{epsg}")
    # NaN is the unambiguous fill for float environmental fields
    da = da.rio.write_nodata(np.nan, encoded=True)
    # sanity: rioxarray must be able to derive an affine transform
    transform = da.rio.transform()
    if transform.a <= 0 or transform.e >= 0:
        # a>0 (x increasing east), e<0 (y decreasing south) is the valid form
        raise ValueError(f"Unexpected geotransform orientation: {transform}")
    return da

Parameter notes: write_crs and write_nodata return new objects — a frequent bug is calling da.rio.write_crs(...) and discarding the result. The transform check catches a flipped y axis (e >= 0) before you write a mirrored file. Complexity: O(1); these touch metadata only.

---

Step 4 — Write GeoTIFF and Zarr #

For a single field bound for desktop GIS, write a tiled, compressed GeoTIFF. Tiling (512×512 blocks) lets readers fetch a window without decoding the whole image, and DEFLATE/LZW compression shrinks the smooth interpolated surfaces dramatically. For a stack — hourly surfaces over a month, or a multivariate Dataset — write Zarr, which chunks along time and reads one slice without touching the rest.

# python 3.11 · rioxarray==0.15.0 · xarray==2024.3.0
from pathlib import Path


def write_geotiff(da: xr.DataArray, path: str | Path, tile: int = 512) -> Path:
    """
    Serialize a georeferenced DataArray to a tiled, compressed GeoTIFF.

    Time complexity:  O(ny * nx) — one pass over cells.
    Space complexity: O(tile^2) streamed per block, not the whole array.
    """
    path = Path(path)
    da.rio.to_raster(
        path,
        driver="GTiff",
        tiled=True,
        blockxsize=tile,
        blockysize=tile,
        compress="DEFLATE",
        predictor=2,          # horizontal differencing — helps smooth fields
        num_threads="ALL_CPUS",
    )
    return path


def write_zarr(stack: xr.DataArray, path: str | Path,
               t_chunk: int = 1, xy_chunk: int = 512) -> Path:
    """
    Serialize a (time, y, x) stack to a chunked Zarr store.

    Chunk time=1 so a single-timestep read touches one chunk per tile.
    """
    path = Path(path)
    chunked = stack.chunk({"time": t_chunk, "y": xy_chunk, "x": xy_chunk})
    chunked.to_dataset(name=stack.name or "field").to_zarr(
        path, mode="w", consolidated=True,
    )
    return path

Parameter notes: predictor=2 (horizontal differencing) improves DEFLATE ratios on smooth surfaces by 20–40%; leave it off for classified integer rasters. consolidated=True writes a single metadata file so opening a Zarr store over object storage is one request, not thousands. Complexity: O(ny · nx) for the GeoTIFF, streamed block by block so peak memory is one tile.

For the specialized case of a single field that must stream efficiently from cloud storage — HTTP range requests, embedded overviews, strict layout — promote the GeoTIFF to a Cloud-Optimized GeoTIFF, covered in depth in writing Cloud-Optimized GeoTIFFs from sensor grids.

---

Step 5 — Reopen and Validate #

Never trust a write you have not read back. Reopen the file and assert the CRS, transform, nodata, and value range survived serialization intact.

# python 3.11 · rioxarray==0.15.0 · numpy==1.26.4
import numpy as np
import rioxarray


def validate_written_raster(path: str, expected_epsg: int) -> dict:
    """Reopen a written raster and confirm its georeferencing is intact."""
    da = rioxarray.open_rasterio(path, masked=True)
    report = {
        "crs_ok": da.rio.crs.to_epsg() == expected_epsg,
        "has_nodata": da.rio.nodata is not None,
        "shape": tuple(da.shape),
        "north_up": da.rio.transform().e < 0,
        "finite_fraction": float(np.isfinite(da.values).mean()),
    }
    da.close()
    return report

Complexity: O(ny · nx) to scan the finite fraction; drop that check for very large stacks and validate a windowed sample instead.

---

Configuration and Tuning #

Cell size and output format are deployment decisions driven by sensor spacing and consumer. A grid finer than your sensor density fabricates detail the interpolation cannot support; coarser than the network’s effective resolution throws away signal.

Compression choice: use DEFLATE with predictor=2 for continuous float fields (temperature, concentration) and LZW for classified integer rasters. Never use lossy JPEG compression on measurement data — it introduces block artifacts that read as fake spatial gradients.

Chunk sizing for Zarr: chunk time=1 for per-timestep dashboards, or time=24 if you routinely compute daily statistics across the whole grid. The chunk that matches your dominant read pattern minimizes I/O amplification.

---

Validation #

A correctly written raster field has a small set of invariants. Check them before the file leaves your pipeline.

• Transform orientation: transform.a > 0 (x increases eastward) and transform.e < 0 (y decreases southward). A positive e means a vertically flipped image — the most common silent raster bug.
• CRS round-trips: da.rio.crs.to_epsg() returns the exact code you wrote. A None result means the CRS was lost or is a custom WKT with no EPSG equivalent.
• Nodata is honored: reopen with masked=True and confirm masked cells are NaN, not a raw sentinel leaking as 0.
• Value range is physical: for float fields the finite (non-nodata) values fall inside sensor-plausible bounds — PM2.5 in roughly 0–500 µg/m³, air temperature −50 to 60 °C. Interpolation overshoot (kriging can produce negative concentrations near steep gradients) shows up here.
• Extent matches source: da.rio.bounds() equals the interpolation extent to within one cell size.

Expected output shape: a 300×300 UTM grid at 50 m cell size written as a tiled DEFLATE GeoTIFF lands around 150–400 KB for a smooth field — an order of magnitude smaller than the uncompressed 360 KB float32 payload only if the surface is genuinely smooth. A file that fails to compress signals noise the interpolation should have suppressed.

Quality-flag distribution: track finite_fraction per raster. A field where more than ~40% of cells are nodata means the interpolation extent overreaches the sensor coverage — mask the convex hull of the stations before writing rather than shipping a mostly-empty raster.

---

Failure Modes and Edge Cases #

Discarded return values from write_crs / write_nodata. rioxarray’s accessor methods are functional by default — they return a new object and leave the original untouched. da.rio.write_crs("EPSG:32633") on its own line does nothing observable; you must rebind da = da.rio.write_crs(...). This silently produces CRS-less rasters.

Ascending y axis writes an upside-down image. If the grid was built with np.linspace(ymin, ymax, ny), the y coordinate ascends and the GeoTIFF writes with the south row first, appearing vertically mirrored in every viewer. Build y descending, or fix it with da.sortby("y", ascending=False) immediately before writing.

Interpolating in geographic degrees. Building the grid in raw WGS 84 (EPSG:4326) makes cells physically non-square — a 0.01° cell is ~1.1 km north-south but only ~0.7 km east-west at 50° latitude — so distance-weighted interpolation becomes anisotropic. Interpolate and grid in a projected CRS, then reproject the final raster with rio.reproject("EPSG:4326") only if a consumer specifically needs degrees.

float64 grids double storage for no benefit. Environmental measurement precision is far below float64 resolution. Carrying float64 through to disk doubles every file and every network transfer. Cast to float32 at the DataArray step; the truncation error (~1e-7 relative) is invisible against sensor noise.

Nodata as 0 corrupts statistics. Leaving masked cells as 0 instead of NaN biases every spatial mean, gradient, and zonal statistic downstream. Any real reading of exactly zero (a valid PM2.5 measurement) also becomes indistinguishable from “no data.” Always encode a nodata value distinct from the physical domain.

Unchunked stacks trigger OOM. Calling .values or an eager reduction on a multi-gigabyte raster stack materializes the whole cube. Open large stacks with chunks={"x": 512, "y": 512} so operations stay lazy under Dask, and only .compute() the final reduced result.

---

Integration #

A written raster field is a hand-off point, and it feeds — or is fed by — the neighboring stages of the Geospatial Data Storage, Interpolation & GIS Export pipeline:

1. Upstream — interpolation produces the values. The 2-D array wrapped here is the direct output of spatial interpolation using kriging and IDW for sensor fields. Grid resolution and extent decisions are shared with that stage: interpolate onto the exact axes this page constructs so no resampling round-trip is needed.
2. Sideways — the vector alternative. When the deliverable is discrete features (station points, threshold-exceedance polygons) rather than a continuous field, skip rasterization and export vectors through GeoJSON and QGIS export workflows instead. Many projects ship both: a raster surface plus a GeoJSON of the source stations styled on top.
3. Downstream — cloud delivery. For rasters served from object storage to web maps or partner pipelines, convert the tiled GeoTIFF into a validated Cloud-Optimized GeoTIFF via writing Cloud-Optimized GeoTIFFs from sensor grids, which adds internal overviews and a strict layout for HTTP range reads.

Carry the interpolation provenance — method, variogram model or IDW power, station count, timestamp — into the raster attributes so an auditor opening the GeoTIFF a year later can reconstruct exactly how the surface was produced.

---

FAQ #

Why is my GeoTIFF flipped vertically compared to the source grid?

GeoTIFFs expect the first raster row to be the northernmost. rioxarray infers scan direction from the sign of the y-coordinate spacing, so the y coordinate must run from high (north) to low (south). If you built the grid with np.linspace(ymin, ymax, ny) the y axis is ascending and the image writes upside down. Reverse it with da.sortby("y", ascending=False) or build y with linspace(ymax, ymin, ny) before writing.

Should I store an interpolated surface as GeoTIFF or Zarr?

Use GeoTIFF for a single 2-D field consumed by desktop GIS such as QGIS or ArcGIS, and for anything shared with external partners. Use Zarr when you have a stack — hourly PM2.5 surfaces over a month, or a multi-variable Dataset — because Zarr chunks along time and reads a single slice without deserializing the whole cube. A common pattern is Zarr as the working archive and per-timestep GeoTIFF exports on demand.

What nodata value should I use for an interpolated environmental field?

For float32 concentration or temperature fields use np.nan as the in-memory fill and write it through with rio.write_nodata(np.nan, encoded=True). NaN is unambiguous because no valid measurement equals it, and GDAL, QGIS, and rasterio all honor a NaN nodata tag. Reserve integer sentinels like -9999 only for integer rasters where NaN cannot be represented; if you use one, make sure it lies outside every physically plausible reading.

How do I keep xarray from loading a huge raster stack into memory?

Chunk the DataArray with .chunk() so it is backed by Dask, or open existing files with rioxarray.open_rasterio(path, chunks={"x": 512, "y": 512}). Operations then stay lazy until you call .compute() or write. For Zarr, set the chunk shape at write time so each read touches only the tiles it needs. Never call .values on an unchunked multi-gigabyte stack — it materializes the entire array.

My cells look offset by half a pixel in QGIS — what went wrong?

xarray coordinates label pixel centers, but a GeoTIFF geotransform references the top-left corner of the top-left pixel. rioxarray handles the half-cell shift automatically when the coordinates are evenly spaced, so a visible offset almost always means your x or y spacing is irregular or the coordinates were placed on cell edges rather than centers. Rebuild the axes as xmin + (i + 0.5) * cell_size and confirm np.diff is constant before writing.

---

Related #

• Geospatial Data Storage, Interpolation & GIS Export — parent overview covering storage, interpolation, and export across the full geospatial pipeline
• Writing Cloud-Optimized GeoTIFFs from Sensor Grids — promote a tiled GeoTIFF to a validated COG with overviews for cloud-native reads
• Spatial Interpolation with Kriging and IDW for Sensor Fields — produces the interpolated value grid this page wraps and writes
• GeoJSON and QGIS Export Workflows — the vector alternative when the deliverable is discrete features rather than a continuous raster field