Reading NetCDF Climate Data with xarray

A NetCDF climate file opens in one line and then refuses to behave like a raster: the time axis will not decode, a European bounding box returns zero cells, and rio.to_raster() complains there is no CRS. This guide is for anyone pulling CMIP6, ERA5 or gridded observational data into a Python raster workflow — it sits under xarray & rioxarray Raster Cubes in Mastering Core Geospatial Python Libraries, and covers the path from open_dataset to a georeferenced GeoTIFF subset without guessing at any of the metadata along the way.

Why This Approach / What Goes Wrong

NetCDF is a self-describing container, not a raster format. It stores n-dimensional arrays plus a bag of attributes, and the Climate and Forecast (CF) conventions are what turn that bag into something spatial. xarray reads the container faithfully and applies the CF decoders it can — but where a file is silent or non-standard, xarray stays silent too. It will not invent a projection, it will not normalise your longitudes, and it will not guess which of your four dimensions are the spatial pair.

Three specific gaps account for almost every failure. First, time. CF encodes time as an integer offset plus a units string like days since 1850-01-01 and a calendar attribute. NumPy's datetime64[ns] can only represent roughly 1678–2262 on a proleptic Gregorian calendar, so a 360_day model calendar or an 1850 epoch with a long run pushes decoding outside what NumPy can hold, and open_dataset raises rather than silently truncating. Second, longitude convention. Global model output is very often written on a 0–360 grid, so sel(lon=slice(-10, 5)) for western Europe matches nothing at all and returns an empty array with no warning — label-based selection has no notion of wrap-around. Third, CRS. A plain NetCDF has no grid_mapping variable; the grid is implicitly geographic, but GDAL needs that stated explicitly before any georeferenced write.

Reading the repr before writing any selection code is the whole discipline. Everything you need to know — which names are dimensions, which are coordinates, what units each variable claims, whether a CRS was declared — is printed there.

Prerequisites

• xarray>=2025.1 — the Dataset model, label-based selection and the xarray.coders decoder API
• rioxarray>=0.17 — registers the .rio accessor that adds CRS, transform and raster I/O
• netCDF4>=1.7 — the default backend engine for classic and NetCDF-4 files
• rasterio>=1.4 — the GDAL binding rioxarray writes through, and the tool you verify the output with
• cftime>=1.6 — non-standard calendar objects (360_day, noleap, julian)
• dask>=2024.5 — optional, but required for lazy chunked reads of multi-gigabyte cubes

conda install -c conda-forge "xarray=2025.1.*" "rioxarray=0.17.*" "netcdf4=1.7.*" \
  "rasterio=1.4.*" "cftime=1.6.*" "dask=2024.5.*"

Install the whole set from conda-forge. xarray's backend, rioxarray and rasterio all bind the same HDF5, netCDF-C and GDAL shared libraries, and mixing pip wheels with conda builds is the usual source of HDF5 error on a file that opens fine elsewhere.

Step-by-Step Implementation

1. Open the file and read the repr before writing any selection code.

open_dataset is lazy — it reads headers and coordinate arrays, not the data variables. Printing the object costs nothing and tells you the dimension names, coordinate ranges and attribute set you are about to depend on.

import xarray as xr

ds = xr.open_dataset("cmip6_tas_monthly.nc")   # engine auto-detected: netcdf4
print(ds)

# For a cube too large for RAM, add chunks to defer everything to dask:
# ds = xr.open_dataset("cmip6_tas_monthly.nc", chunks={"time": 120})

2. Interrogate dims, coords and attributes separately.

On a Dataset, .dims is a set-like view of dimension names in current xarray — use .sizes when you want the name-to-length mapping. Variable-level attributes are where CF hides units, standard_name and _FillValue, and they are per-variable, not global.

print(dict(ds.sizes))            # {'time': 1980, 'lat': 180, 'lon': 360}
print(list(ds.coords))           # ['time', 'lat', 'lon']
print(list(ds.data_vars))        # ['tas', 'pr']
print(ds.attrs.get("Conventions"))   # 'CF-1.10'
print(ds["tas"].attrs)           # {'units': 'K', 'standard_name': 'air_temperature', ...}
print(float(ds.lon.min()), float(ds.lon.max()))   # 0.5 359.5  -> a 0–360 grid

3. Decode times explicitly instead of hoping the default works.

Open once with decode_times=False to see the raw encoding. If the calendar is non-standard or the epoch is far outside the datetime64[ns] window, ask for cftime objects — they represent any CF calendar exactly, and string-based .sel() still works against them.

raw = xr.open_dataset("cmip6_tas_monthly.nc", decode_times=False)
print(raw["time"].attrs)
# {'units': 'days since 1850-01-01', 'calendar': '360_day', 'axis': 'T'}

coder = xr.coders.CFDatetimeCoder(use_cftime=True)
ds = xr.open_dataset("cmip6_tas_monthly.nc", decode_times=coder)
print(ds["time"].values[:2])
# [cftime.Datetime360Day(1850, 1, 16, 0, 0, 0, 0, has_year_zero=True) ...]

On xarray older than 2025.1 the equivalent call is xr.open_dataset(path, use_cftime=True); that argument still works but is deprecated in favour of the coder object. Note that mask_and_scale=True is already the default, so _FillValue, scale_factor and add_offset are applied on read — packed 16-bit reanalysis data comes back as floats with NaN in the gaps without you asking.

4. Select by label and slice a bounding box.

.sel() matches coordinate labels, and slice bounds must run in the same direction as the coordinate. Latitude in reanalysis products is very often stored north-to-south, so slice(35, 60) returns nothing while slice(60, 35) returns the band you wanted. Read the direction from the data rather than assuming it.

lat_ascending = bool(ds.lat[0] < ds.lat[-1])
lat_bounds = slice(35, 60) if lat_ascending else slice(60, 35)

decade = ds["tas"].sel(time=slice("2000-01", "2009-12"))
print(decade.sizes["time"])      # 120

# Nearest-label lookup for a single grid cell (London, on a 0–360 grid)
cell = ds["tas"].sel(lat=51.5, lon=359.9, method="nearest")
print(float(cell.isel(time=0)))  # 279.4

5. Roll longitudes from 0–360 to −180…180.

The roll is pure metadata: reassign the coordinate values, then re-sort so the axis is monotonic increasing again. Nothing is interpolated and no data cell moves except in sort order. Coordinate arithmetic drops attributes by default, so put units and standard_name back or CF-aware readers downstream will complain.

def to_180(obj, lon_name="lon"):
    """Convert a 0–360 longitude axis to −180…180 and re-sort it."""
    if float(obj[lon_name].max()) <= 180:
        return obj
    attrs = dict(obj[lon_name].attrs)
    rolled = obj.assign_coords({lon_name: (((obj[lon_name] + 180) % 360) - 180)})
    rolled = rolled.sortby(lon_name)
    rolled[lon_name].attrs = attrs or {"units": "degrees_east", "standard_name": "longitude"}
    return rolled

ds180 = to_180(ds)
print(float(ds180.lon.min()), float(ds180.lon.max()))   # -179.5 179.5

europe = ds180["tas"].sel(lat=lat_bounds, lon=slice(-10, 5))
print(dict(europe.sizes))    # {'time': 1980, 'lat': 50, 'lon': 15}

6. Attach the CRS the file never stated.

Importing rioxarray registers the .rio accessor. Tell it which coordinates are spatial, then write the CRS — write_crs adds a spatial_ref coordinate and a grid_mapping attribute, which is exactly the CF construct the source file was missing. EPSG:4326 nominally defines latitude first, but rioxarray and GDAL always treat the x dimension as longitude, which is the same always_xy=True convention explained in Coordinate Systems with PyProj.

import rioxarray  # noqa: F401 — importing registers the .rio accessor

tas = europe.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=False)
tas = tas.rio.write_crs("EPSG:4326", inplace=False)
tas = tas.rio.write_coordinate_system(inplace=False)   # CF axis + units attrs

print(tas.rio.crs)          # EPSG:4326
print(tas.rio.transform())  # | 1.00, 0.00,-10.00| ...

If the file does carry a grid_mapping variable — many regional model outputs do, on a rotated pole or Lambert conformal grid — ds.rio.crs picks it up automatically and write_crs would overwrite correct metadata. Check ds.rio.crs is None before assuming.

7. Export the subset to GeoTIFF.

GeoTIFF is 2D plus bands, so reduce to at most three dimensions first and drop any singleton plev or ensemble axis. Sort latitude descending so the transform is north-up, declare the nodata value, then write. A 3D (time, y, x) array becomes one band per timestep.

import numpy as np

window = tas.sel(time=slice("2000-01", "2000-12")).sortby("lat", ascending=False)
window = window.rio.write_nodata(np.nan, inplace=False)

window.rio.to_raster(
    "tas_europe_2000.tif",
    driver="GTiff",
    dtype="float32",
    compress="DEFLATE",
    tiled=True,
)

Reproject only after subsetting — warping the whole global cube to burn out a European window wastes most of the work. When the target grid must match an existing raster exactly, use the alignment path in Reprojecting Raster Cubes with reproject_match rather than reprojecting each side independently.

Verification

Check the three things that fail silently: the CRS is present, the longitude axis is in the range and order you think it is, and the written file carries a north-up transform with the band count you expect.

import rasterio

# 1. CRS was actually written, not just requested
assert tas.rio.crs is not None, "no CRS — rio.to_raster would emit an ungeoreferenced TIFF"
assert tas.rio.crs.to_epsg() == 4326

# 2. Longitude axis rolled and sorted
assert float(tas.lon.min()) >= -180.0 and float(tas.lon.max()) <= 180.0
assert tas.lon.to_index().is_monotonic_increasing, "sortby('lon') was skipped"

# 3. The GeoTIFF on disk agrees
with rasterio.open("tas_europe_2000.tif") as src:
    print(src.crs, src.count, src.width, src.height)
    # EPSG:4326 12 15 50
    assert src.crs.to_epsg() == 4326
    assert src.transform.e < 0, "south-up raster — sort lat descending before writing"
    band = src.read(1, masked=True)
    print(f"valid cells: {band.count()} / {band.size}")
    # valid cells: 750 / 750

A band count of 12 confirms the twelve monthly steps survived as bands; src.transform.e < 0 confirms the pixel height is negative, which is what every downstream GIS reader assumes.

Edge Cases & Debugging

• ValueError: did not find a match in any of xarray's currently installed IO backends. The netcdf4 package is missing, or the file is NetCDF-4/HDF5 and the installed engine only handles classic. Install netcdf4 and h5netcdf, then retry with an explicit engine="h5netcdf".
• ValueError: unable to decode time units ... with calendar '360_day'. Reopen with decode_times=xr.coders.CFDatetimeCoder(use_cftime=True), or with decode_times=False if you want the raw offsets and will call xr.decode_cf(ds) yourself.
• A .sel() bounding box returns zero cells. Either the longitude axis is 0–360 (check float(ds.lon.max()) > 180) or the latitude axis descends and your slice bounds are the wrong way round. Both are label-direction problems, and neither raises.
• MissingSpatialDimensionError: x dimension not found. rioxarray could not guess the spatial pair from names like rlon/rlat or nav_lon/nav_lat. Call .rio.set_spatial_dims(x_dim=..., y_dim=...), or .rename({"rlon": "x", "rlat": "y"}) first.
• to_raster fails on a 4D array. A plev or ensemble-member dimension is still attached. Reduce it with .isel(plev=0) or .squeeze(drop=True) so the array is at most (band, y, x).
• Values sit around 1e20 instead of NaN. The writer stored _FillValue as a string, so mask_and_scale skipped it. Mask explicitly with ds["tas"] = ds["tas"].where(ds["tas"] < 1e19) after checking the raw attribute.

Frequently Asked Questions

When should I use open_mfdataset instead of open_dataset? Whenever the archive splits one logical cube across files — one file per year or per variable is the standard CMIP6 layout. xr.open_mfdataset("tas_*.nc", combine="by_coords", chunks={"time": 120}) concatenates them lazily into a single dask-backed Dataset. Decode the time axis the same way you would for a single file, because a decode failure in one member fails the whole open, and pass parallel=True only when you already have a dask client running.

Do I need rioxarray at all if I only want to slice and plot? No. Plain xarray handles opening, decoding, label selection, resampling and statistics — .rio adds nothing to any of that. Bring in rioxarray at the point where the data has to become a raster in the GIS sense: writing a CRS, clipping by a vector geometry, reprojecting, or exporting to GeoTIFF or COG. Keeping the boundary clear also keeps the dependency out of pipelines that never leave NetCDF.

Why write EPSG:4326 explicitly when the coordinates are obviously latitude and longitude? Because "obviously" is a human inference and GDAL does not make it. Without a grid_mapping, rio.to_raster() writes a TIFF with a valid affine transform and no projection block, and every consumer downstream — QGIS, a tile server, PostGIS — either refuses it or silently assumes something. write_crs also pins the datum: degrees on WGS 84 and degrees on NAD 83 look identical in the coordinate arrays and differ by up to a metre or two on the ground.

Should I reproject before or after subsetting the cube? After, almost always. Reprojection is the most expensive operation in the chain and its cost scales with the number of cells warped, so slicing the time range and bounding box first can cut it by orders of magnitude. The exception is when the source grid is curvilinear or rotated — there, label-based slicing on the native axes does not correspond to a rectangle on the ground, so warp to a regular grid first and subset afterwards.