Tutorial: Working with Shapefiles and NetCDF in Python¶

This tutorial introduces key geospatial analysis workflows using Python:

• Reading shapefiles and converting to shapely geometries
• Testing if coordinates fall inside polygons
• Visualising spatial data and geographic boundaries
• Reading and analysing gridded NetCDF rainfall datasets
• Creating spatial masks and averaging values over regions

We’ll focus on the Murray-Darling Basin (MDB) as a case study.

Colab Setup¶

If you use Google Colab, run the next cell before you start the lab. It installs the geospatial libraries and downloads the Murray-Darling Basin boundary files into the notebook workspace.

In [1]:

%pip install --quiet netCDF4 pyshp shapely
from pathlib import Path
if not Path("MDB_boundaries/MDB_north_boundary.shp").exists():
    !curl -L --fail --silent --show-error -o MDB_boundaries.zip https://data.gadopt.org/water-course/MDB_boundaries.zip
    !unzip -o -q MDB_boundaries.zip

%pip install --quiet netCDF4 pyshp shapely from pathlib import Path if not Path("MDB_boundaries/MDB_north_boundary.shp").exists(): !curl -L --fail --silent --show-error -o MDB_boundaries.zip https://data.gadopt.org/water-course/MDB_boundaries.zip !unzip -o -q MDB_boundaries.zip

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

Import Required Libraries¶

We begin by importing essential Python libraries:

• numpy for numerical operations
• netCDF4 for handling NetCDF files
• matplotlib.pyplot for plotting
• pyshp (via shapefile.Reader) for shapefiles
• shapely.geometry to convert and operate on vector shapes

In [2]:

solution

import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from shapefile import Reader
from shapely.geometry import Point, shape

import numpy as np import netCDF4 as nc import matplotlib.pyplot as plt from shapefile import Reader from shapely.geometry import Point, shape

Load and Convert Shapefiles to Geometry¶

For this tutorial, we will use the MDB boundaries shapefile. Download the shapefiles from here and unzip it in the current directory. Use Reader() to load two shapefiles for the MDB: one for the north and one for the south. Convert each to a shapely geometry object so we can do spatial queries.

This enables operations like containment tests (Point.within(polygon), i.e., if a point is inside a polygon).

In [3]:

solution

# Load the north and south MDB shapefiles and convert each to a shapely shape
NMDB = Reader("MDB_boundaries/MDB_north_boundary.shp")
SMDB = Reader("MDB_boundaries/MDB_south_boundary.shp")
# Convert each to a shapely shape
NMDB_shape = shape(NMDB.shape())
SMDB_shape = shape(SMDB.shape())

# Load the north and south MDB shapefiles and convert each to a shapely shape NMDB = Reader("MDB_boundaries/MDB_north_boundary.shp") SMDB = Reader("MDB_boundaries/MDB_south_boundary.shp") # Convert each to a shapely shape NMDB_shape = shape(NMDB.shape()) SMDB_shape = shape(SMDB.shape())

Check if a Point Lies Within the Basin¶

Check if the point (151°E, 29°S) lies in either basin. Point().within(polygon) returns True if the point lies inside the shape.

In [4]:

solution

print(Point([151, -29]).within(NMDB_shape))
print(Point([151, -29]).within(SMDB_shape))

print(Point([151, -29]).within(NMDB_shape)) print(Point([151, -29]).within(SMDB_shape))

True
False

Plot Basin Boundaries and Test Coordinate¶

Extract the coordinate arrays from the shapefiles and plot them. Overlay the test point to confirm visually.

In [5]:

solution

# Extract the coordinate arrays from the shapefiles and plot them.
NMDB_coords = np.array(NMDB.shape().points)
SMDB_coords = np.array(SMDB.shape().points)
# Plot the north and south boundaries
plt.figure(figsize=(5, 5))
plt.plot(SMDB_coords[:, 0], SMDB_coords[:, 1])
plt.plot(NMDB_coords[:, 0], NMDB_coords[:, 1])
plt.scatter(151, -29, c='red')
plt.xlabel("Longitude (°E)")
plt.ylabel("Latitude (°N)")
plt.title("North and South MDB Boundaries with Test Point")
plt.show()

# Extract the coordinate arrays from the shapefiles and plot them. NMDB_coords = np.array(NMDB.shape().points) SMDB_coords = np.array(SMDB.shape().points) # Plot the north and south boundaries plt.figure(figsize=(5, 5)) plt.plot(SMDB_coords[:, 0], SMDB_coords[:, 1]) plt.plot(NMDB_coords[:, 0], NMDB_coords[:, 1]) plt.scatter(151, -29, c='red') plt.xlabel("Longitude (°E)") plt.ylabel("Latitude (°N)") plt.title("North and South MDB Boundaries with Test Point") plt.show()

Open the Rainfall Data over OPeNDAP¶

We use the Australian Water Outlook (AWRA-L v7) monthly rainfall product, served by NCI's THREDDS server. Instead of downloading a file, we hand the OPeNDAP URL straight to netCDF4.Dataset and the server sends us only the bytes we ask for.

The variable is rain_day, in millimetres, and each time step is a monthly total. (Its long_name says "Daily Rainfall" even in the monthly file — that attribute is simply wrong upstream. Trust units and the time axis, not long_name.)

⚠️ Check the units before you trust a number. The Australian Water Outlook also publishes a decile version of exactly the same variable name, rain_day. In that product the values are percentile ranks between 0 and 1 with units: relative — not millimetres. Averaging it gives a number that looks plausible and means nothing. Always read units off the variable before you plot it.

Opening the dataset is instantaneous because it only reads the metadata. The actual transfer happens later, when we slice the array.

We use netCDF4 here because working at the index level is what makes the descending-latitude problem below visible. Lab 3 opens this same file with xarray, which decodes the time axis and attaches the coordinates for you, so you can select by latitude=-35.3 instead of by array index.

In [6]:

solution

url = (
    "https://thredds.nci.org.au/thredds/dodsC/iu04/australian-water-outlook/"
    "historical/v1/AWRALv7/processed/values/month/rain_day.nc"
)
data = nc.Dataset(url)
print("units:", data["rain_day"].units)
print("long_name:", data["rain_day"].long_name, " <- misleading, ignore it")
print("grid:", data["rain_day"].shape, "(time, latitude, longitude)")

url = ( "https://thredds.nci.org.au/thredds/dodsC/iu04/australian-water-outlook/" "historical/v1/AWRALv7/processed/values/month/rain_day.nc" ) data = nc.Dataset(url) print("units:", data["rain_day"].units) print("long_name:", data["rain_day"].long_name, " <- misleading, ignore it") print("grid:", data["rain_day"].shape, "(time, latitude, longitude)")

units: mm
long_name: Daily Rainfall  <- misleading, ignore it
grid: (1387, 681, 841) (time, latitude, longitude)

Decode the Time Axis Properly¶

The raw time values are "days since 1900-01-01". You may see code that converts these with time / 365.25 + 1900. Don't: that ignores leap years and the result drifts by roughly a month over a century. Use netCDF4.num2date, which reads the units and calendar attributes and gives you real datetime objects.

In [7]:

solution

time_var = data["time"]
dates = nc.num2date(time_var[:], time_var.units, only_use_cftime_datetimes=False)
print("units attribute:", time_var.units)
print("record runs", dates[0].date(), "to", dates[-1].date(), f"({len(dates)} months)")

time_var = data["time"] dates = nc.num2date(time_var[:], time_var.units, only_use_cftime_datetimes=False) print("units attribute:", time_var.units) print("record runs", dates[0].date(), "to", dates[-1].date(), f"({len(dates)} months)")

units attribute: days since 1900-01-01
record runs 1911-01-31 to 2026-07-31 (1387 months)

Subset to the Basin — Mind the Descending Latitude¶

The full grid is 681 × 841 cells over all of Australia and 1386 months. Pulling all of it would be several gigabytes, so we cut it down to the MDB bounding box and a five-year window (2019–2023) before asking the server for anything.

⚠️ Latitude runs from −10 to −44, i.e. descending. This catches almost everyone. A slice has to follow the direction the axis is stored in. With xarray you would write sel(latitude=slice(-24, -38)); writing sel(latitude=slice(-38, -24)) returns an empty array with no error at all, and the failure only surfaces later as a nan. We are indexing by position here, so we sidestep it — but note below that lat_idx starts at the northern edge.

We chose 2019–2023 deliberately: 2019 is the driest year in the whole 115-year record, and the La Niña years that follow are among the wettest.

In [8]:

solution

lats_full = data["latitude"][:]
lons_full = data["longitude"][:]
print("latitude goes from", lats_full[0], "to", lats_full[-1], "-> descending")

# Index ranges covering the basin. Because latitude descends, the *first* index is the north edge.
lat_idx = np.where((lats_full <= -24.0) & (lats_full >= -38.0))[0]
lon_idx = np.where((lons_full >= 138.0) & (lons_full <= 153.0))[0]
time_idx = np.where([2019 <= d.year <= 2023 for d in dates])[0]

# Keep the slices contiguous — OPeNDAP is efficient for ranges, terrible for scattered indices.
lats = lats_full[lat_idx[0]:lat_idx[-1] + 1]
lons = lons_full[lon_idx[0]:lon_idx[-1] + 1]
times = dates[time_idx[0]:time_idx[-1] + 1]
rain = data["rain_day"][
    time_idx[0]:time_idx[-1] + 1,
    lat_idx[0]:lat_idx[-1] + 1,
    lon_idx[0]:lon_idx[-1] + 1,
]
print("downloaded", rain.shape, "=", rain.nbytes / 1e6, "MB")
print("first month:", times[0].date(), " last month:", times[-1].date())

lats_full = data["latitude"][:] lons_full = data["longitude"][:] print("latitude goes from", lats_full[0], "to", lats_full[-1], "-> descending") # Index ranges covering the basin. Because latitude descends, the *first* index is the north edge. lat_idx = np.where((lats_full <= -24.0) & (lats_full >= -38.0))[0] lon_idx = np.where((lons_full >= 138.0) & (lons_full <= 153.0))[0] time_idx = np.where([2019 <= d.year <= 2023 for d in dates])[0] # Keep the slices contiguous — OPeNDAP is efficient for ranges, terrible for scattered indices. lats = lats_full[lat_idx[0]:lat_idx[-1] + 1] lons = lons_full[lon_idx[0]:lon_idx[-1] + 1] times = dates[time_idx[0]:time_idx[-1] + 1] rain = data["rain_day"][ time_idx[0]:time_idx[-1] + 1, lat_idx[0]:lat_idx[-1] + 1, lon_idx[0]:lon_idx[-1] + 1, ] print("downloaded", rain.shape, "=", rain.nbytes / 1e6, "MB") print("first month:", times[0].date(), " last month:", times[-1].date())

latitude goes from -10.0 to -44.0 -> descending

downloaded (60, 281, 301) = 20.29944 MB
first month: 2019-01-31  last month: 2023-12-31

Create Spatial Masks for Each Basin¶

We want to isolate grid cells that fall within the MDB. Loop through the lat/lon grid and use Point.within() to assign True to cells inside the north or south basin.

This is about 84,000 containment tests and takes a few seconds. It is only affordable because we cropped to the bounding box first — over the full Australian grid it would be 570,000 tests.

In [9]:

solution

# Build boolean masks for NMDB and SMDB based on which grid points fall inside each
NMDB_mask = np.zeros((len(lats), len(lons)), dtype=bool)
SMDB_mask = np.zeros((len(lats), len(lons)), dtype=bool)
# Loop through the lat/lon grid and use `Point.within()` to assign `True` to cells inside the north or south basin.
for ilat in range(len(lats)):
    for ilon in range(len(lons)):
        pt = Point([lons[ilon], lats[ilat]])
        if pt.within(NMDB_shape):
            NMDB_mask[ilat, ilon] = True
        elif pt.within(SMDB_shape):
            SMDB_mask[ilat, ilon] = True

# Build boolean masks for NMDB and SMDB based on which grid points fall inside each NMDB_mask = np.zeros((len(lats), len(lons)), dtype=bool) SMDB_mask = np.zeros((len(lats), len(lons)), dtype=bool) # Loop through the lat/lon grid and use `Point.within()` to assign `True` to cells inside the north or south basin. for ilat in range(len(lats)): for ilon in range(len(lons)): pt = Point([lons[ilon], lats[ilat]]) if pt.within(NMDB_shape): NMDB_mask[ilat, ilon] = True elif pt.within(SMDB_shape): SMDB_mask[ilat, ilon] = True

Plot the Mask to Verify Coverage¶

Combine the north and south masks, and plot to confirm the shape of the coverage.

Because latitude is stored descending, row 0 is the northernmost row — which is exactly what imshow expects, so the map comes out the right way up. We pass extent so the axes carry real coordinates rather than array indices.

In [10]:

solution

# Combine the north and south masks, and plot to confirm the shape of the coverage.
MDB_mask = NMDB_mask | SMDB_mask
extent = [lons[0], lons[-1], lats[-1], lats[0]]
plt.figure(figsize=(5, 5))
plt.imshow(MDB_mask, extent=extent)
plt.xlabel("Longitude (°E)")
plt.ylabel("Latitude (°N)")
plt.title("Spatial Mask for North + South MDB")
print("grid cells inside the basin:", MDB_mask.sum())
plt.show()

# Combine the north and south masks, and plot to confirm the shape of the coverage. MDB_mask = NMDB_mask | SMDB_mask extent = [lons[0], lons[-1], lats[-1], lats[0]] plt.figure(figsize=(5, 5)) plt.imshow(MDB_mask, extent=extent) plt.xlabel("Longitude (°E)") plt.ylabel("Latitude (°N)") plt.title("Spatial Mask for North + South MDB") print("grid cells inside the basin:", MDB_mask.sum()) plt.show()

grid cells inside the basin: 40366

Visualise Rainfall for a Single Month¶

Pick one month and hide everything outside the basin, so only rainfall within the MDB shows. We use np.where to put nan outside the mask; imshow leaves nan blank.

In [11]:

solution

# October 2022 — the third year of the triple-dip La Niña and the wettest October on record.
imonth = next(i for i, d in enumerate(times) if d.year == 2022 and d.month == 10)
field = np.where(MDB_mask, rain[imonth], np.nan)

plt.figure(figsize=(5, 5))
plt.imshow(field, extent=extent, cmap="Blues")
plt.colorbar(label="Monthly rainfall (mm)")
plt.xlabel("Longitude (°E)")
plt.ylabel("Latitude (°N)")
plt.title(f"MDB rainfall, {times[imonth]:%B %Y}")
plt.show()

# October 2022 — the third year of the triple-dip La Niña and the wettest October on record. imonth = next(i for i, d in enumerate(times) if d.year == 2022 and d.month == 10) field = np.where(MDB_mask, rain[imonth], np.nan) plt.figure(figsize=(5, 5)) plt.imshow(field, extent=extent, cmap="Blues") plt.colorbar(label="Monthly rainfall (mm)") plt.xlabel("Longitude (°E)") plt.ylabel("Latitude (°N)") plt.title(f"MDB rainfall, {times[imonth]:%B %Y}") plt.show()

Average Rainfall Over the Basin¶

For each month, take the mean of all grid cells inside the mask. That gives a time series of basin-average monthly rainfall, in millimetres.

This is a plain average over cells. It is not quite the true basin average, because a 0.05° cell in the north of the basin covers about 15 % more ground than one in the south, so the northern cells are being under-counted. Lab 4 derives the latitude-dependent cell area and shows how to weight by it. For a basin this compact the difference is small, but the habit matters — and if you convert to a volume it matters a great deal.

In [12]:

solution

MDB_mean_rain = np.zeros(len(times))
for i in range(len(times)):
    MDB_mean_rain[i] = rain[i][MDB_mask].mean()
print(f"mean over 2019-2023: {MDB_mean_rain.mean():.1f} mm/month")
print(f"wettest: {times[MDB_mean_rain.argmax()]:%b %Y} at {MDB_mean_rain.max():.1f} mm")
print(f"driest:  {times[MDB_mean_rain.argmin()]:%b %Y} at {MDB_mean_rain.min():.1f} mm")

MDB_mean_rain = np.zeros(len(times)) for i in range(len(times)): MDB_mean_rain[i] = rain[i][MDB_mask].mean() print(f"mean over 2019-2023: {MDB_mean_rain.mean():.1f} mm/month") print(f"wettest: {times[MDB_mean_rain.argmax()]:%b %Y} at {MDB_mean_rain.max():.1f} mm") print(f"driest: {times[MDB_mean_rain.argmin()]:%b %Y} at {MDB_mean_rain.min():.1f} mm")

mean over 2019-2023: 41.0 mm/month
wettest: Oct 2022 at 143.7 mm
driest:  Sep 2023 at 5.4 mm

Plot the Basin-Average Time Series¶

Finally, plot the series against real dates. The seasonal cycle is visible, and so is the jump out of the 2019 drought into the wet La Niña years.

In [13]:

solution

plt.figure(figsize=(8, 3))
plt.plot(times, MDB_mean_rain, linewidth=1, marker='o', markersize=3)
plt.xlabel("Date")
plt.ylabel("Basin-average rainfall (mm/month)")
plt.title("Basin-average monthly rainfall over the Murray-Darling Basin")
plt.grid(alpha=0.3)
plt.show()

plt.figure(figsize=(8, 3)) plt.plot(times, MDB_mean_rain, linewidth=1, marker='o', markersize=3) plt.xlabel("Date") plt.ylabel("Basin-average rainfall (mm/month)") plt.title("Basin-average monthly rainfall over the Murray-Darling Basin") plt.grid(alpha=0.3) plt.show()

✅ Extensions¶

To go further, try:

• Averaging rainfall separately for NMDB and SMDB and comparing them
• Extending the time window back to 1911 and looking at the seasonal cycle or long-term variability
• Weighting by grid-cell area (see Lab 4) and converting the basin mean into a volume in GL
• Exporting the rainfall time series to CSV for reporting