⚙️ How-To Coastal & marine Gridded & climate data Coastal & Marine · step 15 of 20 Updated

How to subset an ocean model NetCDF by depth and time

Problem statement

An ocean model output is a four-dimensional array — time, depth, latitude, longitude — and subsetting it is where most of the work and most of the errors are. The depth axis in particular has conventions that no other dimension has: it can be positive down or positive up, it can be a coordinate in metres or a dimensionless sigma level, and its values can be layer centres or layer interfaces.

The practical consequence is that ds.sel(depth=50) can return the layer at 50 m, the layer at −50 m, the 50th sigma level, or nothing at all, and all four look like a successful selection.

This guide subsets by time, depth and area in the right order, and covers the checks that confirm the selection is what you meant.

Quick answer

Inspect before selecting, and use method="nearest" with a tolerance:

import xarray as xr

ds = xr.open_dataset("ocean.nc")
print(ds.dims, list(ds.coords))
for name in ("depth", "lev", "z", "deptht"):
    if name in ds.coords:
        d = ds[name]
        print(name, d.attrs.get("positive"), d.attrs.get("units"),
              float(d[0]), "→", float(d[-1]))

sub = (ds["thetao"]
       .sel(time=slice("2026-01-01", "2026-03-31"))
       .sel(depth=50, method="nearest", tolerance=10)
       .sel(lat=slice(51, 56), lon=slice(1, 9)))
print(sub.sizes, float(sub["depth"]))

A tolerance is the difference between a selection and a silent substitution. Without it, method="nearest" will happily return the 500 m layer when you asked for 50 m and the model has none between.

Stack of the four axes of an ocean model output with the convention question attached to each.
Four axes; only the depth axis has three separate conventions to establish.

Step-by-step solution

1. Find the depth coordinate and read its attributes

It may be called depth, lev, z, deptht, st_ocean or something model-specific. The CF positive attribute says whether it increases up or down, and units says whether it is metres or something else.

2. Establish whether it is a real depth or a sigma level

A sigma or hybrid coordinate is dimensionless and follows the bathymetry, so level 5 is a different physical depth at every grid point. Converting it to metres needs the formula terms named in the formula_terms attribute and the bathymetry variable.

3. Subset in the cheapest order

Time first, then depth, then area. Each step reduces what the next has to read, and on a chunked file it reduces which chunks are touched at all.

4. Use sel with a tolerance, not isel with a guessed index

isel(depth=3) is the fourth layer, whatever that is. sel(depth=50, method="nearest", tolerance=10) is a statement about metres that raises when it cannot be satisfied.

5. Check the latitude direction before slicing

Many ocean grids store latitude descending. slice(51, 56) on a descending axis returns nothing, silently, and the result is an empty array rather than an error.

lat = ds["lat"].values
s = slice(51, 56) if lat[0] < lat[-1] else slice(56, 51)

6. Check the longitude convention

Models frequently use 0–360. A slice from −10 to 10 returns nothing on such a grid, and subsetting across the prime meridian needs two slices and a concatenation — or a roll into the other convention.

7. Confirm the selection before computing anything

Print the sizes, the selected depth and the time range. A subset with the right shape and the wrong layer produces a plausible map.

Flow from opening lazily through time, depth and area selection to a loaded array.
Time, depth, area, then load — each step reduces what the next reads.

Code examples

Example 1 — describe the file before touching it 

import xarray as xr, numpy as np

DEPTH_NAMES = ("depth", "deptht", "lev", "z", "st_ocean", "olevel")

def describe_ocean(path):
    ds = xr.open_dataset(path)
    out = {"dims": dict(ds.sizes), "data_vars": list(ds.data_vars)}

    depth_name = next((n for n in DEPTH_NAMES if n in ds.coords), None)
    if depth_name:
        d = ds[depth_name]
        out["depth"] = {
            "name": depth_name, "n": int(d.size),
            "units": d.attrs.get("units"), "positive": d.attrs.get("positive"),
            "first": float(d[0]), "last": float(d[-1]),
            "sigma": "formula_terms" in d.attrs,
        }
    if "lat" in ds.coords:
        lat = ds["lat"].values
        out["lat"] = {"ascending": bool(lat[0] < lat[-1]),
                      "range": (float(lat.min()), float(lat.max()))}
    if "lon" in ds.coords:
        lon = ds["lon"].values
        out["lon"] = {"convention": "0-360" if lon.max() > 180 else "-180-180",
                      "range": (float(lon.min()), float(lon.max()))}
    if "time" in ds.coords:
        t = ds["time"]
        out["time"] = {"n": int(t.size), "first": str(t.values[0]),
                       "last": str(t.values[-1]),
                       "calendar": t.encoding.get("calendar")}
    return out

Running this first is thirty seconds that saves an afternoon, and the sigma flag in particular changes what every later step means.

Example 2 — a subsetter that handles the conventions 

import xarray as xr, numpy as np

def subset(ds, var, time=None, depth_m=None, bbox=None,
           depth_tolerance=10.0, depth_name=None):
    a = ds[var]
    depth_name = depth_name or next((n for n in DEPTH_NAMES if n in a.coords), None)

    if time is not None:
        a = a.sel(time=slice(*time))

    if depth_m is not None and depth_name:
        d = a[depth_name]
        target = depth_m
        if d.attrs.get("positive") == "up":
            target = -abs(depth_m)
        a = a.sel({depth_name: target}, method="nearest", tolerance=depth_tolerance)

    if bbox is not None:
        w, s, e, n = bbox
        lat = a["lat"].values
        lat_slice = slice(s, n) if lat[0] < lat[-1] else slice(n, s)
        lon = a["lon"].values
        if lon.max() > 180:
            w, e = w % 360, e % 360
        a = a.sel(lat=lat_slice, lon=slice(w, e))

    return a

sub = subset(ds, "thetao", time=("2026-01-01", "2026-03-31"),
             depth_m=50, bbox=(1, 51, 9, 56))
print(sub.sizes, "at depth", float(sub[DEPTH_NAMES[0]]) if "depth" in sub.coords else "n/a")

The positive: up branch is the one that matters: on such a file, 50 m below the surface is a coordinate value of −50, and asking for +50 either raises on the tolerance or silently picks the surface.

Example 3 — convert sigma levels to metres 

import xarray as xr, numpy as np

def sigma_to_depth(ds, sigma_name="lev", bathy_name="depth_bnds", eta_name=None):
    """Ocean sigma: z = eta + sigma * (depth + eta). CF formula_terms names the parts."""
    sigma = ds[sigma_name]
    terms = dict(t.split(": ") for t in sigma.attrs["formula_terms"].split())
    s = ds[terms["sigma"]]
    depth = ds[terms["depth"]]
    eta = ds[terms["eta"]] if "eta" in terms else 0.0
    return eta + s * (depth + eta)

z = sigma_to_depth(ds)
print("physical depth range:", float(z.min()), "→", float(z.max()))

On a sigma grid there is no single depth for a level, so "the 50 m layer" is not a selection you can make — you interpolate onto a z grid first, and that is a resampling decision with its own error.

Explanation

Why the depth axis has so many conventions

Ocean models solve on grids that suit their numerics. A z-level model has fixed depths; a sigma model follows the bathymetry so that the same number of layers spans a shelf and a trench; a hybrid model mixes the two. Each is stored with the convention that matches the model, and CF provides positive and formula_terms so a reader can tell — which only helps if the reader looks.

Why tolerance matters more than method

method="nearest" never fails: there is always a nearest value. On a coarse model whose layers are 0, 10, 20, 50, 100, 200, 500, 1000 and 2000 m, asking for 300 m returns the 200 m layer, which is a 100 m error presented as a successful selection. A tolerance turns that into an exception.

Why subsetting order affects performance

On a chunked or remote file, each sel narrows the set of chunks that later operations touch. Selecting three months of a ten-year run first reduces the time axis by a factor of forty before the spatial selection is evaluated. Selecting the area first would read every timestep for that area, which is the same volume for a different reason and much worse on a file chunked by time.

Why a descending latitude axis fails silently

slice(51, 56) on a descending array asks for values from index of 51 to index of 56 in the order the array is stored, finds that 51 comes after 56, and returns an empty selection. Empty is a valid array, so nothing raises, and the error appears later as an empty plot or a NaN statistic.

Triage of four ways a depth selection returns the wrong layer: a positive-up axis, no tolerance, a sigma coordinate and a guessed index.
Every one produces a plausible map of the wrong layer.

Edge cases or notes

  • Check time.encoding["calendar"]. A 360-day calendar breaks pandas indexing.
  • Depth bounds versus centres. depth_bnds are interfaces; depth is the centre.
  • Land is NaN or a fill value. Mask before averaging.
  • Staggered grids. Velocity components sit on different grid points from temperature.
  • open_mfdataset for a run split by file, with combine="by_coords".
  • Chunking dictates performance. Match the subset to the chunk layout.
  • Units may be Kelvin or Celsius. Read the attribute.
  • Record the exact selection in the output's metadata, including the depth found.

FAQ

How do I select a depth layer in xarray?

sel(depth=50, method="nearest", tolerance=10). The tolerance is what turns a silent substitution into an error.

Why does my depth selection return the wrong layer?

Either the axis is positive: up, so 50 m below the surface is −50, or the nearest layer is further away than you assumed and there was no tolerance to catch it.

What is a sigma level?

A dimensionless vertical coordinate that follows the bathymetry, so the same level is a different physical depth at every point. Convert it with the CF formula_terms before selecting by metres.

Why is my spatial subset empty?

Usually a descending latitude axis, so slice(51, 56) selects nothing, or a 0–360 longitude grid sliced with negative values.

In what order should I subset?

Time, then depth, then area, then load. Each step reduces what the next reads.

How do I know which variable is the depth coordinate?

Look for depth, deptht, lev, z, st_ocean or olevel in the coordinates, then read its units and positive attributes.