1. Python basics: what it is, how it runs, variables, lists, loops, functions, modules
2. Station data with pandas: reading a CSV, summarising, grouping
3. Time series plots with matplotlib
4. Gridded data with xarray: NetCDF, selecting, unit conversion, reductions
5. Spatial plots: maps, colour maps, regions, multi-panel figures
6. Saving results, the Climate Data Store, and where to go next

Two data files sit next to the notebook:

• Sampledata.csv: daily rainfall at three sites, 2005 to 2024
• fnl_africa_202302-202303.nc: NCEP FNL analysis cut to Africa, Feb to Mar 2023

1. Install Python and JupyterLab

Install Miniconda, open a terminal (Anaconda Prompt on Windows), then:

conda create -n climate -c conda-forge python=3.12 \
    xarray netcdf4 pandas matplotlib cartopy jupyterlab
conda activate climate

2. Get the files

Put the notebook and both data files in one folder. In the terminal, cd into that folder and run jupyter lab. It opens in your browser.

Jupyter in three keys

Stuck? Kernel → Restart, then run from the top.

No install possible? Google Colab runs notebooks in the browser. Upload the files and run !pip install xarray netcdf4 cartopy first.

A high-level programming language known for its simplicity and readability.

Created by Guido van Rossum, first released in 1991. Now on version 3.

Free and open source: anyone can read, modify and distribute it. That is why a huge scientific community has grown around it.

Why climate scientists use it

• Readable: code looks like the maths it describes
• One language for everything: reading files, statistics, maps, papers' figures
• The ecosystem: NumPy, pandas, xarray, matplotlib, cartopy, SciPy, all free
• Shared conventions: NetCDF and CF metadata work out of the box
• The community: CMIP, ERA5, Copernicus tutorials are written in Python

Python is interpreted: a program called the interpreter reads your code line by line and executes it. No compiling step.

You can run it three ways:

• Interactively, typing one line at a time
• As a script, a .py file run top to bottom
• In a notebook, cells of code and text mixed, run in any order

The notebook shows the value of the last line of each cell. Use print() to show anything else.

A variable is a name for a value. Create one with =.

n_years    = 30          # int
temp_mean  = 16.4        # float
station    = "Cape Town" # str
is_raining = False       # bool

temp_day = [18, 20, 24, 16]   # list
coords   = (-33.9, 18.4)     # tuple

Names: letters, numbers and underscores. Cannot start with a number. Temp and temp are different.

# Arithmetic:  +  -  *  /  **  %  //
temp_c = temp_k - 273.15
mm_day = prate * 86400     # kg m⁻² s⁻¹ → mm/day

# Strings and f-strings
print(f"{station}: {rain:.1f} mm")
Cape Town: 12.5 mm

# Lists count from 0
temp_day[0]      # 18, the first
temp_day[-1]     # 16, the last
temp_day[1:3]    # [20, 24], end not included
temp_day.append(22)

A precipitation rate in kg m⁻² s⁻¹ is the same as mm/s, because 1 kg of water over 1 m² is 1 mm deep.

Slicing with [start:stop] excludes stop. Remember this: xarray's slice() behaves differently.

We are given today's temperature and the past mean for that day. If today is 4 degrees or more above the mean, issue a warning.

temp_mean = 16
temp_day  = 20

if temp_day >= temp_mean + 4:
    print("HOT - Issue a warning")
else:
    print("No warning necessary")

More than one temperature? Loop over the list.

temp_day = [18, 20, 24, 16]

for t in temp_day:
    if t >= temp_mean + 4:
        print(t, "HOT - Issue a warning")
    else:
        print(t, "No warning necessary")

A function is a named piece of code that does one task. It may take inputs and may return a value.

def warn_temp(temp_day, temp_mean):
    if temp_day >= temp_mean + 4:
        return "HOT - Issue a warning"
    return "No warning necessary"

for t in temp_day:
    print(t, warn_temp(t, temp_mean))

A module is a file of functions someone else wrote. Load it with import.

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

temps = np.array([18, 20, 24, 16])
temps.mean()        # 19.5
temps + 273.15      # whole array at once

These four import lines start almost every climate notebook. The short names are a convention everyone uses.

import pandas as pd

data = pd.read_csv("Sampledata.csv", sep=";")
data.head()

Rows are days. Columns are variables. The numbers on the left are the index.

sep=";" because this file uses semicolons. The default is a comma.

A single column, data["Site_A_Mediterranean_mm"], is a Series.

data.head()       # first 5 rows
data.tail(3)      # last 3 rows
data.shape        # (7304, 6)
data.columns      # column names
data.describe()   # count, mean, std, min, max…

data["Site_A_Mediterranean_mm"].mean()
data["Site_B_TropicalMonsoon_mm"].max()

Methods need brackets (). Attributes such as shape and columns do not.

Build a real date from the three columns, so pandas and matplotlib understand time:

data["Date"] = pd.to_datetime(
    data[["year", "month", "day"]]
)

Keep only the rows that match a condition:

wet = data[data["Site_A_Mediterranean_mm"] > 20]
len(wet)

groupby splits the table into groups, applies a statistic to each, and joins the results back into a table.

sites = ["Site_A_Mediterranean_mm",
         "Site_B_TropicalMonsoon_mm",
         "Site_C_TemperateMaritime_mm"]

# mean annual cycle: one row per month
monthly_mean = data.groupby("month")[sites].mean()

# annual totals: one row per year
annual = data.groupby("year")[sites].sum()

import matplotlib.pyplot as plt

plt.figure(figsize=(12, 4))            # 1 figure

plt.plot(data["Date"],                   # 2 draw
         data["Site_A_Mediterranean_mm"])

plt.xlabel("Date")                       # 3 label
plt.ylabel("Rainfall (mm)")
plt.title("Daily rainfall - Site A")

plt.show()                               # 4 show

The figure is the whole image. Inside it, the axes is the drawing area. Every plt. call adds to the current axes.

• Several lines: call plt.plot once per line with a label=, then plt.legend()
• marker="o" puts a dot on each point
• plt.grid(True) adds a grid
• plt.axhline(0, linestyle="--") draws a zero line, useful for anomalies
• plt.bar(names, values, yerr=std) for bars with error bars

plt.figure(figsize=(10, 5))

for site in sites:
    plt.plot(annual.index, annual[site],
             marker="o", label=site)

plt.xlabel("Year")
plt.ylabel("Annual rainfall (mm)")
plt.legend()
plt.grid(True)

plt.savefig("annual_rainfall.png",
            dpi=150, bbox_inches="tight")
plt.show()

savefig must come before show. After show, the figure is gone and you save a blank image.

• dpi=150 for slides, dpi=300 for a paper
• bbox_inches="tight" trims the white margins
• .pdf or .svg give vector output for publication

The loop over sites replaces three nearly identical plt.plot lines. That is what loops are for.

A .nc file stores three things together:

• Dimensions: time, plev, lat, lon
• Coordinates: the actual values along each dimension
• Variables with attributes: units, long_name, standard_name

The CF conventions say how to name these, so any tool can read any CF file. xarray follows them, which is why your plots label themselves.

xarray's two objects:

NCEP FNL (Final) operational global analysis, February to March 2023, cut to Africa for today. The original global file is 1.3 GB with the same structure.

import xarray as xr

ds = xr.open_dataset("fnl_africa_202302-202303.nc")
ds                     # interactive summary

t = ds["t"]            # one DataArray
t.dims                 # ('time','plev','lat','lon')
t.shape                # (236, 3, 86, 91)
t.attrs["units"]       # 'K'

That means the order of work is:

1. Select the variable, level, region, time you need
2. Convert units
3. Reduce: mean, max, sum along a dimension
4. Plot

# one level, first time step
t.sel(plev=100000).isel(time=0)

# nearest grid point to Cape Town, all times
t.sel(plev=100000, lat=-33.9, lon=18.4,
      method="nearest")

# a region: southern Africa
t.sel(lat=slice(0, -40), lon=slice(10, 52))

# a range of dates
t.sel(time=slice("2023-02-01", "2023-02-28"))

# one moment
t.sel(time="2023-02-15T12")

• .sel uses coordinate values: a real latitude, a real date
• .isel uses positions: 0, 1, 2 … like a list
• method="nearest" when your value is off the grid
• slice(start, stop) includes both ends, unlike a Python list

t_c = t - 273.15
t_c.attrs["units"] = "°C"

pr = ds["prate"] * 86400
pr.attrs["units"] = "mm/day"

# collapse time → a map
t_map = t_c.sel(plev=100000).mean(dim="time")

# collapse space → a time series
pr_sa = pr.sel(lat=slice(0, -40), lon=slice(10, 52))
pr_series = pr_sa.mean(dim=("lat", "lon"))

Arithmetic applies to every value at once. Attributes are dropped by arithmetic, so set units again and the plots stay labelled.

dim= names what you collapse. What is left decides what you get:

# 6-hourly → daily means
pr_daily = pr_series.resample(time="1D").mean()

# monthly mean maps: one per calendar month
pr_monthly = pr_sa.groupby("time.month").mean(dim="time")
pr_monthly.sel(month=3) - pr_monthly.sel(month=2)

• resample changes the time step: "1D" daily, "MS" month start, "YS" year start
• groupby("time.month") pools every February together, every March together: the way to make a climatology
• Also time.year, time.season, time.dayofyear

Same ideas as pandas, because xarray borrowed them.

plt.figure(figsize=(12, 5))
t_map.plot()
plt.title("Mean 1000 hPa temperature, Feb–Mar 2023")
plt.show()

# control the colours
pr.mean(dim="time").plot(cmap="YlGnBu", vmin=0, vmax=15)

# a difference: diverging, centred on zero
pr_diff.plot(cmap="BrBG", center=0)

Two dimensions left → a filled map with a colour bar. One dimension left → a line. Axis labels and the colour bar label come from the long_name and units attributes.

• cmap= the colour map
• vmin=, vmax= fix the range, essential when comparing panels
• center=0 puts the neutral colour at zero
• ax= which panel to draw in

Sequential, one hue light to dark, for magnitudes: temperature, rainfall, height.

Diverging, two hues around a neutral middle, for anomalies and differences. Always center=0.

Avoid rainbow maps.

• The bright bands create edges that are not in the data
• Lightness goes up and down, so it does not read as "more"
• Red-green readers cannot tell half the scale apart

Full list: matplotlib colormaps reference. Add _r to any name to reverse it.

import cartopy.crs as ccrs
import cartopy.feature as cfeature

fig = plt.figure(figsize=(9, 7))
ax = plt.axes(projection=ccrs.PlateCarree())

t_map.sel(lat=slice(0, -40), lon=slice(10, 52)).plot(
    ax=ax, transform=ccrs.PlateCarree(), cmap="RdYlBu_r",
    cbar_kwargs={"label": "Temperature (°C)"},
)
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
ax.gridlines(draw_labels=True)
plt.show()

cartopy gives matplotlib a map projection. Two things change:

• The axes gets a projection=
• The data gets a transform= saying which coordinates it is in. Lat/lon data is always PlateCarree(), whatever the projection.

Other projections: ccrs.Robinson() for the globe, ccrs.LambertConformal() for a region, ccrs.SouthPolarStereo() for Antarctica.

The first run downloads coastline data, so it needs internet.

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

pr_monthly.sel(month=2).plot(
    ax=axes[0], cmap="YlGnBu", vmin=0, vmax=15)
axes[0].set_title("February 2023")

pr_monthly.sel(month=3).plot(
    ax=axes[1], cmap="YlGnBu", vmin=0, vmax=15)
axes[1].set_title("March 2023")

plt.tight_layout()
plt.show()

• plt.subplots(rows, cols) returns the figure and a grid of axes
• Pass ax=axes[i] to say where to draw
• With an axes object the calls change name: ax.set_title, ax.set_xlabel, ax.grid
• Same vmin/vmax in every panel, or the comparison lies
• sharex=True when panels share a time axis

when = "2023-02-15T12"
region = dict(lat=slice(0, -40), lon=slice(10, 52))

gh500 = ds["gh"].sel(plev=50000, time=when, **region)
u850  = ds["u_2"].sel(plev=85000, time=when, **region)
v850  = ds["v_2"].sel(plev=85000, time=when, **region)

fig, ax = plt.subplots(figsize=(9, 7))
cs = gh500.plot.contour(ax=ax, levels=15, colors="black")
ax.clabel(cs, fontsize=8, fmt="%.0f")

ax.quiver(u850.lon[::3], u850.lat[::3],
          u850[::3, ::3], v850[::3, ::3])
plt.show()

• .plot.contour() draws lines, .plot.contourf() fills between them
• ax.clabel writes the value on each contour
• quiver(x, y, u, v) draws arrows. [::3] keeps every third grid point so the arrows do not pile up

The 500 hPa height field with 850 hPa winds is the classic synoptic chart. You now know how to make one from any reanalysis file.

pr_monthly.name = "pr"
pr_monthly.to_netcdf("pr_sa_monthly_2023.nc")

pr_daily.to_dataframe(name="pr_mm_day") \
        .to_csv("pr_sa_daily_2023.csv")

plt.savefig("figure.png", dpi=150,
            bbox_inches="tight")

Give the array a name and keep its units and long_name. Your future self, and your colleagues, will open the file with no idea what it is.

cds.climate.copernicus.eu holds ERA5 reanalysis, seasonal forecasts, CMIP6 projections, satellite records and more. Free, with an account.

One-time setup

1. Register and log in
2. Copy your Personal Access Token from your profile page
3. Save it in ~/.cdsapirc:

url: https://cds.climate.copernicus.eu/api
key: YOUR-PERSONAL-ACCESS-TOKEN

4. Accept the licence on each dataset's page
5. pip install cdsapi

Requesting data from Python

import cdsapi
client = cdsapi.Client()

client.retrieve(
    "reanalysis-era5-single-levels",
    {
        "product_type": "reanalysis",
        "variable": "2m_temperature",
        "year": "2023", "month": "02",
        "day": ["01", "02", "03"],
        "time": "12:00",
        "area": [0, 10, -40, 52],  # N, W, S, E
        "format": "netcdf",
    },
    "era5_t2m.nc",
)

Do not write requests by hand. Fill in the download form on the dataset page and click Show API request.

Everything today ran in JupyterLab on your own machine. The routine each time:

1. Open a terminal, conda activate climate
2. cd to the folder with your notebook and data
3. jupyter lab, it opens in the browser
4. File browser on the left, double-click a notebook
5. Kernel name top right should say Python 3
6. Finished: File → Shut Down, then close the tab

• Share a result: File → Save and Export Notebook As → HTML. Anyone can open it, no Python needed.
• Start clean: Kernel → Restart Kernel and Run All Cells. If it fails, the notebook depends on something you ran out of order.
• Add a package: conda install -c conda-forge scipy in the terminal, then restart the kernel.

• Read a CSV into pandas, summarise it, group by month and year
• Plot time series, anomalies and bars, and save the figure
• Open a NetCDF file with xarray and understand its dimensions, coordinates and attributes
• Select a level, a point, a region and a time range
• Convert units and reduce along a dimension
• Draw maps, pick a colour map on purpose, add coastlines, build panels
• Fetch ERA5 from the Climate Data Store

Slides and notebook: keep them, reuse the code. Every figure you make from now on starts from one of these cells.