← Back to browse

GHCN Plots

This notebook produces the “What the Data Looks Like” plots for the GHCND data source page: four maps showing where GHCND stations are located across Southern Africa, sized by record length and coloured by how recent each station’s data is.

Inputs: ../data/Fig{1,2,3,4}_ghcnd_{PRCP,TAVG,TMAX,TMIN}_stations_southern-africa.csv — produced by GHCN_Download&Process.ipynb.

Outputs: Fig1_ghcnd_PRCP_stations_size.png, Fig2_ghcnd_TAVG_stations_size.png, Fig3_ghcnd_TMAX_stations_size.png, Fig4_ghcnd_TMIN_stations_size.png, saved to ../images/.

Step 1: Setup

Show the code 
import os
import pandas as pd
import numpy as np

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

# Make sure the output folder exists before we try to save anything into it
os.makedirs('../images', exist_ok=True)

Step 2: Load the per-variable station tables

GHCN_Download&Process.ipynb filters the inventory to Southern Africa and splits it into one table per variable. We just load those back in here.

Show the code 
# Figure numbers match the order these four maps appear on the GHCND page,
# and pair each map with its data file the same way Vis1-Vis6 do elsewhere.
fig_numbers = {'PRCP': 1, 'TAVG': 2, 'TMAX': 3, 'TMIN': 4}

var_dfs = {
    var: pd.read_csv(f'../data/Fig{fig_num}_ghcnd_{var}_stations_southern-africa.csv')
    for var, fig_num in fig_numbers.items()
}

Step 3: Plot a station map for each variable

Each map shows: - Marker size — how many years of record the station has (LASTYEAR - FIRSTYEAR) - Marker colour — the most recent year the station reported data, on a shared colour scale across all four maps so they’re easy to compare

Show the code 
# Use the same colour scale across all four maps
all_last_years = np.concatenate([df['LASTYEAR'] for df in var_dfs.values()])
vmin = 1950
vmax = int(np.max(all_last_years))

def plot_station_map(var_df, variable_name, vmin, vmax):
    """Plot station locations for one variable, sized by record length and
    coloured by the last year of available data."""

    # Work out how many years of record each station has
    var_df = var_df.copy()
    var_df['RECORD_LENGTH'] = var_df['LASTYEAR'] - var_df['FIRSTYEAR']

    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

    # Basic map background
    ax.add_feature(cfeature.LAND, facecolor='lightgray')
    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
    ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
    ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)

    # Scale marker size between 10 and 100 based on record length
    min_size, max_size = 10, 100
    sizes = np.interp(
        var_df['RECORD_LENGTH'],
        (var_df['RECORD_LENGTH'].min(), var_df['RECORD_LENGTH'].max()),
        (min_size, max_size)
    )

    # Plot each station: position = location, size = record length, colour = last year of data
    sc = ax.scatter(
        x=var_df['Lon'],
        y=var_df['Lat'],
        c=var_df['LASTYEAR'],
        s=sizes,
        cmap='viridis',
        vmin=vmin,  # shared colour range across all four maps
        vmax=vmax,
        edgecolor='black',
        linewidth=0.5,
        alpha=0.8,
        transform=ccrs.PlateCarree()
    )

    cbar = plt.colorbar(sc, pad=0.02, extend='min')
    cbar.set_label('Last Year of Data')

    # Add a size legend showing what different record lengths look like
    for length in [10, 30, 50, 70]:
        ax.scatter(
            [], [], c='gray',
            s=np.interp(length,
                        (var_df['RECORD_LENGTH'].min(), var_df['RECORD_LENGTH'].max()),
                        (min_size, max_size)),
            label=f'{length} years'
        )
    ax.legend(title='Record Length', loc='lower right')

    ax.set_extent([10, 42, -35, 0])  # Southern Africa
    gl = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), linewidth=0.5, color='gray', alpha=0.7, linestyle='--')
    gl.top_labels = False
    gl.right_labels = False
    gl.xlabel_style = {'size': 10}
    gl.ylabel_style = {'size': 10}

    plt.title(
        f'GHCNd {variable_name} Stations Locations\nSize = Record Length | Color = Last Year\nTotal: {len(var_df)} stations',
        fontsize=12, pad=20
    )

    fig_num = fig_numbers[variable_name]
    plt.savefig(f'../images/Fig{fig_num}_ghcnd_{variable_name}_stations_size.png', dpi=300, bbox_inches='tight')
    plt.show()
    plt.close()

for var_name, var_data in var_dfs.items():
    plot_station_map(var_data, var_name, vmin, vmax)
Back to top