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(dict(ds.sizes), 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["water_temp"]
       .sel(time=slice("2024-01-01", "2024-03-31"))
       .sel(depth=50, method="nearest", tolerance=10)
       .sel(lat=slice(51, 56), lon=slice(1, 9)))
print(sub.sizes, float(sub["depth"]))

{'time': 2, 'depth': 40, 'lat': 126, 'lon': 163} ['time', 'depth', 'lat', 'lon']
depth down m 0.0 → 5000.0
Frozen({'time': 2, 'lat': 126, 'lon': 100}) 50.0

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.

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.

Code examples

The data used on this page

Two real model outputs with different vertical coordinates. The z-level one is HYCOM's global 1/12° analysis — sea water temperature at its 40 fixed depths, on 1 January and 1 February 2024, for 51–56° N and 4° W–9° E — fetched from its NetCDF Subset Service and kept in HYCOM's own conventions, depth positive down and longitude 0–360. The sigma one is one hour of NOAA's San Francisco Bay forecast (SFBOFS), an FVCOM model on 20 terrain-following layers, from the NOAA open-data bucket on AWS.

import os
import subprocess

import xarray as xr

NCSS = ("https://ncss.hycom.org/thredds/ncss/GLBy0.08/expt_93.0?var=water_temp&north=56&south=51"
        "&time_start=2024-01-01T00:00:00Z&time_end=2024-03-01T00:00:00Z&timeStride=248&accept=netcdf4")
for name, west, east in (("hycom_east.nc", -4, 9), ("hycom_west.nc", 356, 360)):
    if not os.path.exists(name):
        subprocess.run(["curl", "-s", "--max-time", "1200", "-o", name,
                        f"{NCSS}&west={west}&east={east}"], check=True)

# HYCOM's tau variable has units "hours since analysis", which xarray cannot decode
asked = xr.open_dataset("hycom_east.nc", drop_variables=["tau"])
print(f"asked for 4° W–9° E, the service returned {float(asked.lon.min()):.2f}–{float(asked.lon.max()):.2f} on a 0–360 axis")
west = xr.open_dataset("hycom_west.nc", drop_variables=["tau"])
ocean = xr.concat([west, asked], dim="lon").sortby("lon")
ocean.to_netcdf("ocean.nc")
ds = xr.open_dataset("ocean.nc")

if not os.path.exists("sfbofs.nc"):
    subprocess.run(["curl", "-s", "-o", "sfbofs.nc", "https://noaa-nos-ofs-pds.s3.amazonaws.com/sfbofs/netcdf/"
                    "2025/08/01/sfbofs.t03z.20250801.fields.n000.nc"], check=True)
print(dict(ds.sizes), "longitudes", float(ds.lon.min()), "…", float(ds.lon.max()))

asked for 4° W–9° E, the service returned 0.00–8.96 on a 0–360 axis
{'time': 2, 'depth': 40, 'lat': 126, 'lon': 163} longitudes 0.0 … 359.9200439453125

The first line is already the lesson of step 6: the service silently dropped the part of the box west of the meridian, because on a 0–360 axis "−4 to 9" is two ranges. The second request fetches 356–360 and the two are joined.

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

import json

print(json.dumps(describe_ocean("ocean.nc"), indent=1, default=str))

{
 "dims": {
  "time": 2,
  "depth": 40,
  "lat": 126,
  "lon": 163
 },
 "data_vars": [
  "water_temp"
 ],
 "depth": {
  "name": "depth",
  "n": 40,
  "units": "m",
  "positive": "down",
  "first": 0.0,
  "last": 5000.0,
  "sigma": false
 },
 "lat": {
  "ascending": true,
  "range": [
   51.0,
   56.0
  ]
 },
 "lon": {
  "convention": "0-360",
  "range": [
   0.0,
   359.9200439453125
  ]
 },
 "time": {
  "n": 2,
  "first": "2024-01-01T00:00:00.000000000",
  "last": "2024-02-01T00:00:00.000000000",
  "calendar": "gregorian"
 }
}

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, "water_temp", time=("2024-01-01", "2024-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")

# a box that crosses the prime meridian
across = subset(ds, "water_temp", time=("2024-01-01", "2024-03-31"), depth_m=50, bbox=(-4, 51, 9, 56))
print("-4..9 E:", dict(across.sizes))
try:
    subset(ds, "water_temp", depth_m=1100, bbox=(1, 51, 9, 56))
except KeyError as refusal:
    print("1,100 m, tolerance 10:", type(refusal).__name__, str(refusal)[:60])
print("the layers either side of 1,100 m:", ds["depth"].sel(depth=slice(900, 1500)).values)
print("without a tolerance:", float(ds["depth"].sel(depth=1100, method="nearest")))

Frozen({'time': 2, 'lat': 126, 'lon': 100}) at depth 50.0
-4..9 E: {'time': 2, 'lat': 126, 'lon': 0}
1,100 m, tolerance 10: KeyError "not all values found in index 'depth'"
the layers either side of 1,100 m: [ 900. 1000. 1250. 1500.]
without a tolerance: 1000.0

The crossing box came back empty. On a 0–360 grid, −4 becomes 356, and slice(356, 9) on an ascending axis selects nothing — without an error. The fix is the one step 6 describes, and it belongs in the function:

def lon_select(a, w, e):
    """Select w..e on either convention, crossing the meridian if it has to."""
    lon = a["lon"].values
    if lon.max() <= 180:
        return a.sel(lon=slice(w, e))
    w, e = w % 360, e % 360
    if w <= e:
        return a.sel(lon=slice(w, e))
    return xr.concat([a.sel(lon=slice(w, 360)), a.sel(lon=slice(0, e))], dim="lon")


box = lon_select(ds["water_temp"].sel(depth=50).sel(lat=slice(51, 56)), -4, 9)
print(dict(box.sizes), f"{float(box.lon[0]):.2f} … {float(box.lon[-1]):.2f}")

{'time': 2, 'lat': 126, 'lon': 163} 356.00 … 8.96

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

import re


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]
    # "sigma: siglay eta: zeta depth: h" -> {"sigma": "siglay", ...}; tolerant of a missing space
    terms = dict(re.findall(r"(\w+):\s*(\w+)", sigma.attrs["formula_terms"]))
    s = ds[terms["sigma"]]
    depth = ds[terms["depth"]]
    eta = ds[terms["eta"]] if "eta" in terms else 0.0
    return eta + s * (depth + eta)

fv = xr.open_dataset("sfbofs.nc")
print(fv["siglay"].attrs["formula_terms"], "|", fv["siglev"].attrs["formula_terms"])
z = sigma_to_depth(fv, sigma_name="siglay")
print("physical depth range:", round(float(z.min()), 2), "→", round(float(z.max()), 2))
deep = int(fv["h"].values.argmax())
print(f"layer 10 of 20: {float(z.isel(siglay=9, node=0).squeeze()):.2f} m at a node {float(fv['h'][0]):.1f} m deep, "
      f"{float(z.isel(siglay=9, node=deep).squeeze()):.2f} m at one {float(fv['h'][deep]):.1f} m deep")

sigma: siglay eta: zeta depth: h | sigma:siglay eta: zeta depth: h
physical depth range: -107.48 → 0.37
layer 10 of 20: -7.85 m at a node 17.0 m deep, -52.20 m at one 110.2 m deep

The first version of this function parsed formula_terms with split(": ") on space-separated tokens, which cannot work — no token contains ": " — and it would also have tripped on this very file, whose siglev attribute reads sigma:siglay with no space. A pattern match handles both.

Layer 10 of SFBOFS's 20 sits 7.85 m down where the bay is 17 m deep and 52.20 m down where it is 110 m deep. 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. HYCOM's layers step from 1,000 m to 1,250 m, so asking it for 1,100 m returns the 1,000 m layer — a 100 m error presented as a successful selection. With tolerance=10 the same request raises a KeyError instead (Example 2).

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.

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.

Internal links

• How to open NetCDF files with xarray in Python — the basics of the container
• How to select time and location in xarray — the general selection patterns
• NetCDF values look wrong — scaling, fill values and units
• CF conventions explained — where positive and formula_terms are defined
• Bathymetry explained: depths, datums and grids — the same sign question in two dimensions
• How to open many NetCDF files at once in Python — a run split across files
• Lazy loading in xarray explained — why order matters
• cftime and datetime errors in xarray — non-standard calendars

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.