Lab: Gridded Data with xarray¶

This lab introduces xarray, the library that makes gridded scientific data pleasant to work with in Python:

• Opening a remote NetCDF file over OPeNDAP, and why that costs almost nothing until you slice
• Selecting data by label (.sel) rather than by array position (.isel)
• The descending-latitude trap, which silently returns nothing instead of raising an error
• Reducing, plotting and re-binning an array by dimension name rather than by axis number

This is the same monthly rainfall file that Lab 2 read with netCDF4, opened a level higher up. Lab 2 worked at the index level on purpose, because that is what makes the latitude problem visible; here we let xarray carry the coordinates and see what that buys us.

Colab Setup¶

If you use Google Colab, run the next cell before you start the lab. netCDF4 gives xarray the OPeNDAP backend that reads the NCI data service.

In [1]:

%pip install --quiet xarray netCDF4

%pip install --quiet xarray netCDF4

Note: you may need to restart the kernel to use updated packages.

1. Why xarray¶

A numpy array is just numbers. If you have a rainfall array of shape (1386, 681, 841), then rain[0, 493, 741] is a valid expression and nothing in the object tells you that those three numbers mean January 1911, 34.65°S and 149.05°E — you have to keep separate coordinate arrays around and do the lookup yourself. That is exactly what Lab 2 did.

An xarray.DataArray is a numpy array plus named dimensions, coordinate values along each dimension, and a dictionary of attributes such as units. Once the labels travel with the data you can write rain.sel(latitude=-35.3, longitude=149.1) instead of rain[:, 493, 741], take a mean over ("latitude", "longitude") instead of over axis=(1, 2), and get a plot with the axes already labelled.

We start, as always, by importing what we need.

• xarray opens the gridded rainfall data and keeps its coordinate labels.
• numpy calculates with arrays of numbers.
• matplotlib.pyplot makes the figures.

xr, np and plt are short names for these libraries. The rest of the lab uses these names.

In [2]:

solution

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

import xarray as xr import numpy as np import matplotlib.pyplot as plt

2. Opening a Remote File¶

We use the same Australian Water Outlook (AWRA-L v7) monthly rainfall product as Lab 2, served over OPeNDAP by NCI's THREDDS server. xr.open_dataset(url) reads only the metadata: the dimension sizes, the coordinate values, the attributes. It returns in about a second and downloads essentially nothing.

The actual transfer happens later, and only for the part you ask for. It is triggered by .load(), by .values, by .plot(), or by any reduction such as .mean(). This laziness is the whole reason the workflow is affordable: the full array is over 3 GB and we will move about 20 MB of it.

A Dataset is a container of one or more variables; ds["rain_day"] pulls out a single DataArray, which is the thing we actually work with. That is the only distinction between the two you need for this course.

⚠️ Read units, not long_name. As in Lab 2: the long_name attribute on this variable says "Daily Rainfall" even in the monthly file. It is simply wrong upstream. The values are monthly totals in millimetres, which is what units and the time axis tell you. The Water Outlook also publishes a decile product under an identical filename whose values are percentile ranks in [0, 1] with units: relative. Averaging that gives a plausible-looking number that means nothing.

In [3]:

solution

url = (
    "https://thredds.nci.org.au/thredds/dodsC/iu04/australian-water-outlook/"
    "historical/v1/AWRALv7/processed/values/month/rain_day.nc"
)
ds = xr.open_dataset(url, engine="netcdf4")
print(ds)

rain = ds["rain_day"]
print("\nunits:", rain.attrs["units"])
print("long_name:", rain.attrs["long_name"], " <- wrong upstream, ignore it")
print("sizes:", dict(rain.sizes))
print("record runs", str(rain["time"].values[0])[:10], "to", str(rain["time"].values[-1])[:10])
print(f"the whole array would be {rain.nbytes / 1e9:.2f} GB — we have not downloaded it")

url = ( "https://thredds.nci.org.au/thredds/dodsC/iu04/australian-water-outlook/" "historical/v1/AWRALv7/processed/values/month/rain_day.nc" ) ds = xr.open_dataset(url, engine="netcdf4") print(ds) rain = ds["rain_day"] print("\nunits:", rain.attrs["units"]) print("long_name:", rain.attrs["long_name"], " <- wrong upstream, ignore it") print("sizes:", dict(rain.sizes)) print("record runs", str(rain["time"].values[0])[:10], "to", str(rain["time"].values[-1])[:10]) print(f"the whole array would be {rain.nbytes / 1e9:.2f} GB — we have not downloaded it")

<xarray.Dataset> Size: 3GB
Dimensions:      (time: 1387, nv: 2, latitude: 681, longitude: 841)
Coordinates:
  * time         (time) datetime64[ns] 11kB 1911-01-31 1911-02-28 ... 2026-07-31
  * latitude     (latitude) float64 5kB -10.0 -10.05 -10.1 ... -43.95 -44.0
  * longitude    (longitude) float64 7kB 112.0 112.0 112.1 ... 153.9 153.9 154.0
Dimensions without coordinates: nv
Data variables:
    time_bounds  (time, nv) datetime64[ns] 22kB ...
    rain_day     (time, latitude, longitude) float32 3GB ...
Attributes:
    var_name:                        rain_day
    title:                           Australian Landscape Water Balance AWRA-...
    Description:                     This data is provided by the Bureau of M...
    summary:                         Data produced by Bureau of Meteorology A...
    source:                          AWRA-L
    date_created:                    2021-09-19T19:54:37
    Conventions:                     CF-1.6, ACDD-1.3
    institution:                     Bureau of Meteorology
    email:                           awrams@bom.gov.au
    url:                             http://www.bom.gov.au/water/landscape
    date_modified:                   2021-09-19T19:54:37
    DODS_EXTRA.Unlimited_Dimension:  time

units: mm
long_name: Daily Rainfall  <- wrong upstream, ignore it
sizes: {'time': 1387, 'latitude': 681, 'longitude': 841}
record runs 1911-01-31 to 2026-07-31
the whole array would be 3.18 GB — we have not downloaded it

3. .sel and .isel¶

There are two ways to pick out part of an array, and mixing them up is the most common beginner mistake with xarray.

• .isel(...) selects by position, exactly like numpy indexing. isel(time=0) is the first time step in the file.
• .sel(...) selects by label, using the coordinate values. sel(time="1911-01") is January 1911 regardless of where it sits in the array.

Position and label are not interchangeable. isel(latitude=0) is not the southern edge of Australia — it is −10.0°, the northern edge, because this file stores latitude in descending order. Hold on to that; §4 is entirely about it.

Two conveniences worth knowing:

• method="nearest" snaps to the closest coordinate value. Without it, sel(latitude=-35.3) raises a KeyError if −35.3 is not exactly on the grid, and it is easy to conclude that xarray is broken when in fact you asked for a latitude that does not exist.
• Partial date strings work: sel(time="2022-10") selects October 2022. That only works because xarray decoded the raw "days since 1900-01-01" numbers into real datetimes when it opened the file — Lab 2 had to call nc.num2date by hand to get the same thing.

⚠️ Do not mix the two in one call. method="nearest" applies to every coordinate you name in that .sel, including time. Each month here is stamped at its end, so sel(time="2022-10", method="nearest") parses the string as 2022-10-01, finds the nearest stamp to be 2022-09-30, and quietly hands you September. Do the nearest-neighbour lookup on the coordinates that need it, then select the time separately.

In [4]:

solution

print("isel(time=0)     ->", str(rain.isel(time=0)["time"].values)[:10])
print("sel(time='1911-01') ->", str(rain.sel(time="1911-01")["time"].values[0])[:10])

print("isel(latitude=0) is latitude", float(rain["latitude"][0]), "— the NORTH edge")
print("isel(latitude=-1) is latitude", float(rain["latitude"][-1]), "— the south")

# Nearest grid point first, then the month — keeping method="nearest" off the time axis.
cell = rain.sel(latitude=-35.3, longitude=149.1, method="nearest")
canberra = cell.sel(time="2022-10").squeeze()
print("\nnearest grid point:", float(canberra["latitude"]), float(canberra["longitude"]))
print(f"Canberra rainfall, October 2022: {float(canberra):.1f} mm")

# What the one-call version would have given you, for comparison:
oops = rain.sel(latitude=-35.3, longitude=149.1, time="2022-10", method="nearest")
print(f"all in one .sel with method='nearest': {float(oops):.1f} mm"
      f"  <- that is {str(oops['time'].values)[:7]}, the WRONG month")

print("isel(time=0) ->", str(rain.isel(time=0)["time"].values)[:10]) print("sel(time='1911-01') ->", str(rain.sel(time="1911-01")["time"].values[0])[:10]) print("isel(latitude=0) is latitude", float(rain["latitude"][0]), "— the NORTH edge") print("isel(latitude=-1) is latitude", float(rain["latitude"][-1]), "— the south") # Nearest grid point first, then the month — keeping method="nearest" off the time axis. cell = rain.sel(latitude=-35.3, longitude=149.1, method="nearest") canberra = cell.sel(time="2022-10").squeeze() print("\nnearest grid point:", float(canberra["latitude"]), float(canberra["longitude"])) print(f"Canberra rainfall, October 2022: {float(canberra):.1f} mm") # What the one-call version would have given you, for comparison: oops = rain.sel(latitude=-35.3, longitude=149.1, time="2022-10", method="nearest") print(f"all in one .sel with method='nearest': {float(oops):.1f} mm" f" <- that is {str(oops['time'].values)[:7]}, the WRONG month")

isel(time=0)     -> 1911-01-31
sel(time='1911-01') -> 1911-01-31
isel(latitude=0) is latitude -10.0 — the NORTH edge
isel(latitude=-1) is latitude -44.0 — the south

nearest grid point: -35.3 149.1
Canberra rainfall, October 2022: 181.3 mm

all in one .sel with method='nearest': 78.2 mm  <- that is 2022-09, the WRONG month

4. The Trap: Slices Must Follow the Stored Direction¶

⚠️ Latitude in this file runs from −10 to −44, i.e. descending. A slice must run in the direction the axis is stored. Writing sel(latitude=slice(-38, -24)) looks correct — smaller number first, as you would write any other range — and it returns an empty array with no error at all. Nothing warns you. The failure only surfaces much later, when every mean you compute comes back as nan and you spend an afternoon looking for the bug in the wrong place.

The correct order here is sel(latitude=slice(-24, -38)): north first, because north is stored first.

The habit that saves you: print .sizes after every slice. A dimension of length 0 is instantly obvious, and it costs one line.

Longitude in this file ascends (112 → 154), so slice(138, 153) is the right way round for it. There is no universal rule; it depends on how the file was written, which is why you check.

In [5]:

solution

wrong = rain.sel(latitude=slice(-38, -24))
right = rain.sel(latitude=slice(-24, -38))
print("slice(-38, -24) ->", wrong.sizes["latitude"], "  <- SILENTLY EMPTY, no error raised")
print("slice(-24, -38) ->", right.sizes["latitude"])

wrong = rain.sel(latitude=slice(-38, -24)) right = rain.sel(latitude=slice(-24, -38)) print("slice(-38, -24) ->", wrong.sizes["latitude"], " <- SILENTLY EMPTY, no error raised") print("slice(-24, -38) ->", right.sizes["latitude"])

slice(-38, -24) -> 0   <- SILENTLY EMPTY, no error raised
slice(-24, -38) -> 281

5. Plotting Straight from the DataArray¶

Because the coordinates and attributes travel with the data, .plot() on a two-dimensional DataArray gives you a map with both axes labelled, a labelled colourbar, and latitude the right way up — despite the axis being stored descending. You did none of that work.

sel(time="2022-10") keeps a length-1 time dimension, which would leave the array three-dimensional. .squeeze() drops it, which makes the intent explicit and keeps the plot title clean.

And a warning that lands hard here: metadata-driven plotting is only as good as the metadata. The colourbar on this figure will say "Daily Rainfall [mm]", because .plot() believes long_name, and long_name is wrong. These are monthly totals.

In [6]:

solution

oct22 = rain.sel(time="2022-10").squeeze()
plt.figure(figsize=(7, 5))
oct22.plot(cmap="Blues")
plt.title("Australian rainfall, October 2022")
plt.show()

print(f"national maximum: {float(oct22.max()):.1f} mm")
print("note the colourbar says 'Daily Rainfall' — these are MONTHLY totals")

oct22 = rain.sel(time="2022-10").squeeze() plt.figure(figsize=(7, 5)) oct22.plot(cmap="Blues") plt.title("Australian rainfall, October 2022") plt.show() print(f"national maximum: {float(oct22.max()):.1f} mm") print("note the colourbar says 'Daily Rainfall' — these are MONTHLY totals")

national maximum: 693.6 mm
note the colourbar says 'Daily Rainfall' — these are MONTHLY totals

6. One Contiguous Request, Then Reduce by Name¶

The full record over the whole country is more than 3 GB. Nobody wants that down a network connection, so the rule for OPeNDAP is: subset first, with contiguous slice() objects, then .load() once. Fancy indexing and boolean masks turn into a great many small requests and are slow enough to look like a hang.

We take the Murray-Darling Basin bounding box for 2019–2023, the same window as Lab 2: 2019 is the driest year of the whole record, and the La Niña years that follow are among the wettest. That is about 20 MB, which is a few seconds.

Once it is in memory, .mean(dim=("latitude", "longitude")) averages over space and leaves the time axis alone. Reducing by dimension name rather than by axis number is the reason to use xarray at all: you cannot get it the wrong way round, and the result keeps its time coordinate, so plotting it gives real dates on the x-axis.

.idxmax("time") returns the coordinate value at which the maximum occurs — the date — rather than an integer index. That is almost always what you want, and unlike a bare .argmax() it is unambiguous about which dimension you mean.

In [7]:

solution

box = rain.sel(
    time=slice("2019", "2023"),
    latitude=slice(-24, -38),
    longitude=slice(138, 153),
).load()
print("downloaded", dict(box.sizes), f"= {box.nbytes / 1e6:.1f} MB")

series = box.mean(dim=("latitude", "longitude"))

plt.figure(figsize=(8, 3))
series.plot(marker="o", markersize=3)
plt.ylabel("box-mean rainfall (mm/month)")
plt.title("Mean rainfall over the MDB bounding box, 2019-2023")
plt.grid(alpha=0.3)
plt.show()

print(f"mean over 2019-2023: {float(series.mean()):.1f} mm/month")
print("wettest:", str(series.idxmax("time").values)[:7], f"at {float(series.max()):.1f} mm")
print("driest: ", str(series.idxmin("time").values)[:7], f"at {float(series.min()):.1f} mm")

box = rain.sel( time=slice("2019", "2023"), latitude=slice(-24, -38), longitude=slice(138, 153), ).load() print("downloaded", dict(box.sizes), f"= {box.nbytes / 1e6:.1f} MB") series = box.mean(dim=("latitude", "longitude")) plt.figure(figsize=(8, 3)) series.plot(marker="o", markersize=3) plt.ylabel("box-mean rainfall (mm/month)") plt.title("Mean rainfall over the MDB bounding box, 2019-2023") plt.grid(alpha=0.3) plt.show() print(f"mean over 2019-2023: {float(series.mean()):.1f} mm/month") print("wettest:", str(series.idxmax("time").values)[:7], f"at {float(series.max()):.1f} mm") print("driest: ", str(series.idxmin("time").values)[:7], f"at {float(series.min()):.1f} mm")

downloaded {'time': 60, 'latitude': 281, 'longitude': 301} = 20.3 MB

mean over 2019-2023: 40.5 mm/month
wettest: 2022-10 at 120.7 mm
driest:  2023-09 at 6.5 mm

7. .weighted — Checking Lab 4's Formula¶

The average above treats every grid cell as equal, and they are not: a 0.05° cell at 24°S covers noticeably more ground than one at 38°S. Lab 4 derives the exact area by hand,

$$ A = R^2 \, \big|\sin\phi_2 - \sin\phi_1\big| \, \big|\Delta\lambda\big| $$

and xarray has weighting built in as da.weighted(w).mean(dim=...). This cell is not an alternative to Lab 4 — it is a check that the two agree. We compute the box mean three ways: unweighted, weighted by $\cos(\text{latitude})$, and weighted by Lab 4's exact areas. The two weighted answers should be numerically indistinguishable, which is the point: at 0.05° the small-cell approximation is not merely "close enough", it is exact to within rounding.

⚠️ This is a bounding-box number, not a basin number. The box is a rectangle; the Murray-Darling Basin is not. The percentage you print below is the weighting effect over this rectangle over these five years, and it is a different number from the one you get over the basin mask. If you are doing Assignment I, do not carry this figure across into Q2 — compute it over your own mask.

Note also that .weighted().mean() returns a depth in mm. Converting to a volume in GL needs actual areas in m², which is Lab 4's formula, not a normalised weight.

In [8]:

solution

w_cos = np.cos(np.deg2rad(box["latitude"]))

# Lab 4's exact formula, for a 0.05 deg cell centred on each latitude.
Rearth = 6370e3
phi1 = np.deg2rad(box["latitude"] - 0.025)
phi2 = np.deg2rad(box["latitude"] + 0.025)
w_area = Rearth**2 * np.abs(np.sin(phi2) - np.sin(phi1)) * np.deg2rad(0.05)

plain = float(box.mean())
cosw = float(box.weighted(w_cos).mean(dim=("latitude", "longitude")).mean())
areaw = float(box.weighted(w_area).mean(dim=("latitude", "longitude")).mean())

print(f"plain (unweighted) mean : {plain:.4f} mm/month")
print(f"cos(latitude)-weighted  : {cosw:.4f} mm/month")
print(f"Lab 4 exact-area weighted: {areaw:.4f} mm/month")
print(f"\nunweighted is {100 * (plain - areaw) / areaw:+.3f} % off")
print(f"the two weightings differ by {100 * abs(cosw - areaw) / areaw:.2e} % — the same answer")
print("this is over the bounding BOX, not the basin — do not reuse the number")

w_cos = np.cos(np.deg2rad(box["latitude"])) # Lab 4's exact formula, for a 0.05 deg cell centred on each latitude. Rearth = 6370e3 phi1 = np.deg2rad(box["latitude"] - 0.025) phi2 = np.deg2rad(box["latitude"] + 0.025) w_area = Rearth**2 * np.abs(np.sin(phi2) - np.sin(phi1)) * np.deg2rad(0.05) plain = float(box.mean()) cosw = float(box.weighted(w_cos).mean(dim=("latitude", "longitude")).mean()) areaw = float(box.weighted(w_area).mean(dim=("latitude", "longitude")).mean()) print(f"plain (unweighted) mean : {plain:.4f} mm/month") print(f"cos(latitude)-weighted : {cosw:.4f} mm/month") print(f"Lab 4 exact-area weighted: {areaw:.4f} mm/month") print(f"\nunweighted is {100 * (plain - areaw) / areaw:+.3f} % off") print(f"the two weightings differ by {100 * abs(cosw - areaw) / areaw:.2e} % — the same answer") print("this is over the bounding BOX, not the basin — do not reuse the number")

plain (unweighted) mean : 40.4923 mm/month
cos(latitude)-weighted  : 40.1760 mm/month
Lab 4 exact-area weighted: 40.1760 mm/month

unweighted is +0.787 % off
the two weightings differ by 1.59e-13 % — the same answer
this is over the bounding BOX, not the basin — do not reuse the number

8. resample: Monthly Totals to Annual Totals¶

Going from monthly values to annual values is a re-binning of the time axis, and resample is xarray's tool for it. resample(time="YE") makes year-end bins; .sum() then adds the months inside each bin.

There is a trap built into this, and it is the same one Assignment I warns about. The record ends part-way through the current year, so the final bin holds only a handful of months and its "annual total" is far too small. A plot of annual rainfall with a mysterious collapse in the last year is almost always this.

The fix is to pair .sum() with .count() and keep only the bins that actually contain twelve months:

annual = point.resample(time="YE").sum()
count  = point.resample(time="YE").count()
annual = annual.where(count == 12, drop=True)

To keep the download small we do this for a single grid cell over the full record. The same calls work unchanged on the 3-D array; only the amount of data moved changes.

A single grid cell is one 5 km square of the Australian landscape, not a region. The Canberra cell below sits in wet foothills country and comes out near 625 mm/yr. The basin-wide average used in Assignment I is around 463 mm/yr. Both are correct; they are answers to different questions, and this one is not a substitute for the other.

In [9]:

solution

point = rain.sel(latitude=-35.3, longitude=149.1, method="nearest").load()

annual = point.resample(time="YE").sum()
count = point.resample(time="YE").count()
print("year-end bins in the record:", annual.sizes["time"])
print("months in the final bin:", int(count.isel(time=-1)), "-> incomplete, drop it")

annual = annual.where(count == 12, drop=True)
years = annual["time"].dt.year
print(f"complete calendar years kept: {annual.sizes['time']} "
      f"({int(years[0])}-{int(years[-1])})")

plt.figure(figsize=(8, 3))
annual.plot.step()
plt.ylabel("annual rainfall (mm/yr)")
plt.title("Annual rainfall at one grid cell near Canberra")
plt.grid(alpha=0.3)
plt.show()

print(f"\nmean: {float(annual.mean()):.1f} mm/yr, std {float(annual.std()):.1f}")
print("wettest:", str(annual.idxmax("time").values)[:4], f"at {float(annual.max()):.1f} mm")
print("driest: ", str(annual.idxmin("time").values)[:4], f"at {float(annual.min()):.1f} mm")
print("\nthis is ONE 5 km cell, not a basin average")

point = rain.sel(latitude=-35.3, longitude=149.1, method="nearest").load() annual = point.resample(time="YE").sum() count = point.resample(time="YE").count() print("year-end bins in the record:", annual.sizes["time"]) print("months in the final bin:", int(count.isel(time=-1)), "-> incomplete, drop it") annual = annual.where(count == 12, drop=True) years = annual["time"].dt.year print(f"complete calendar years kept: {annual.sizes['time']} " f"({int(years[0])}-{int(years[-1])})") plt.figure(figsize=(8, 3)) annual.plot.step() plt.ylabel("annual rainfall (mm/yr)") plt.title("Annual rainfall at one grid cell near Canberra") plt.grid(alpha=0.3) plt.show() print(f"\nmean: {float(annual.mean()):.1f} mm/yr, std {float(annual.std()):.1f}") print("wettest:", str(annual.idxmax("time").values)[:4], f"at {float(annual.max()):.1f} mm") print("driest: ", str(annual.idxmin("time").values)[:4], f"at {float(annual.min()):.1f} mm") print("\nthis is ONE 5 km cell, not a basin average")

year-end bins in the record: 116
months in the final bin: 7 -> incomplete, drop it
complete calendar years kept: 115 (1911-2025)

mean: 625.6 mm/yr, std 165.9
wettest: 1950 at 1081.5 mm
driest:  1944 at 276.4 mm

this is ONE 5 km cell, not a basin average

Summary¶

This lab covered:

• xr.open_dataset(url, engine="netcdf4") on an OPeNDAP endpoint — metadata only, so it is nearly free
• laziness: the transfer happens at .load(), .values, .plot() or a reduction, not at open
• .attrs, .sizes and .nbytes for finding out what you have before you fetch it
• .sel (by label) versus .isel (by position), and method="nearest"
• the descending-latitude trap, and printing .sizes after every slice to catch it
• .plot() for a map and for a time series, straight off the array
• one contiguous .load(), then reducing by dimension name with .mean(dim=(...))
• .weighted(), and .resample(time="YE") with a .count() == 12 completeness check

Two things you will want next.

A calendar-month climatology — the average January, the average February, and so on — is the same idea as resample with different bins: da.groupby("time.month").mean(). Where resample cuts the time axis into consecutive chunks, groupby gathers every January in the record into one group. If you can read the resample cell above, you can read that one.

Getting back to numpy. .values hands you the plain numpy array at any point, and from there numpy, pandas and scipy work as they always do. Assignment I does exactly this: it uses xarray to open the file and pull one contiguous subset, then drops to .values and does the arithmetic in numpy and pandas. Using xarray for access and numpy for computation is a perfectly ordinary way to work.

Three habits to keep:

1. Read units before you trust a number, and do not trust long_name.
2. Print .sizes after every slice.
3. Subset before you .load().