GHCN Visualisations

This notebook produces the six Visualisations-section plots (Vis1-6) for the GHCND data source page: the raw GHCN station record on its own, then GHCN compared against ERA5.

Inputs (from ../data/, produced by GHCN_Download&Process.ipynb): - pr_GHCN_daily_Skukuza-SF000068296.csv, temp_GHCN_daily_Skukuza-SF000068296.csv - PRCPTOT_daily_ERA5_..._Skukuza.csv, tas_daily_..._Skukuza.csv, tasmax_daily_..._Skukuza.csv, tasmin_daily_..._Skukuza.csv

Outputs: Vis1 to Vis6 PNGs, described above each section below.

Step 1: Setup

Show the code

import os
import numpy as np
import pandas as pd
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 station and ERA5 data

The rainfall and temperature station files are loaded separately (since that’s how GHCN_Download&Process.ipynb saved them) and joined back into one station_data table, so the rest of this notebook can treat GHCN as a single dataset.

Show the code

pr_df = pd.read_csv('../data/pr_GHCN_daily_Skukuza-SF000068296.csv', index_col=0, parse_dates=True)
temp_df = pd.read_csv('../data/temp_GHCN_daily_Skukuza-SF000068296.csv', index_col=0, parse_dates=True)

station_data = pr_df.join(
    temp_df[['tasmax', 'TMAX_ATTRIBUTES', 'tasmin', 'TMIN_ATTRIBUTES', 'tas', 'TAVG_ATTRIBUTES']],
    how='outer'
)
station_data.index = pd.to_datetime(station_data.index)

daily_rainfall_df = pd.read_csv('../data/PRCPTOT_daily_ERA5_19590101-20241231_Skukuza.csv', parse_dates=True, index_col=0)
daily_temp_df = pd.read_csv('../data/tas_daily_ERA5_19590101-20241231_Skukuza.csv', parse_dates=True, index_col=0)
daily_tmin_df = pd.read_csv('../data/tasmin_daily_ERA5_19590101-20241231_Skukuza.csv', parse_dates=True, index_col=0)
daily_tmax_df = pd.read_csv('../data/tasmax_daily_ERA5_19590101-20241231_Skukuza.csv', parse_dates=True, index_col=0)

Step 3: Helper function used by Vis1 and Vis2

Both raw GHCN time series plots need to break the line wherever there’s a gap or a change between “good” quality data and GSOD-sourced data (a different processing source, shown in a lighter colour), so they share this one helper.

Show the code

def plot_segments(ax, dates, values, flagged_mask, base_color, flagged_color):
    """Plot a line in segments, starting a new segment whenever the data goes
    missing or switches between a flagged lower-confidence source (GSOD/SSOD)
    and the primary source."""
    start_idx = 0
    n = len(values)
    for i in range(1, n + 1):
        if i == n or np.isnan(values[i]) or flagged_mask[i] != flagged_mask[start_idx] or np.isnan(values[start_idx]):
            if not np.isnan(values[start_idx:i]).all():
                ax.plot(
                    dates[start_idx:i],
                    values[start_idx:i],
                    color=flagged_color if flagged_mask[start_idx] else base_color
                )
            start_idx = i

Vis1: Daily rainfall time series (GHCN only)

Plots the full raw rainfall record, highlighting the stretches sourced from GSOD in a lighter colour.

Show the code

var = 'pr'
attr_col = 'PRCP_ATTRIBUTES'

base_color = 'tab:blue'     # primary source data
flagged_color = 'lightblue' # GSOD ('S') or SSOD ('2') sourced data - dulled, lower-confidence

dates = station_data.index
values = station_data[var].values

# SFLAG is the last character of the attribute string (MFLAG,QFLAG,SFLAG).
# '2' (SSOD) is documented by NOAA as the direct successor to 'S' (GSOD), so both
# are treated the same way here - they share the same "use with caution" caveat.
sflag = station_data[attr_col].str.split(',').str[-1]
flagged_mask = sflag.isin(['S', '2']).values

plt.figure(figsize=(15, 5))
ax = plt.gca()
plot_segments(ax, dates, values, flagged_mask, base_color, flagged_color)

ax.set_xlabel('Date')
ax.set_ylabel('Daily Rainfall (mm)')
ax.set_title('Daily Rainfall (with GSOD/SSOD data dulled) for Skukuza, RSA')

plt.savefig('../images/Vis1_GHCN_daily_TimeSeries_Rainfall_Skukuza.png', dpi=300, bbox_inches='tight')
plt.show()

Vis2: Daily temperature time series (GHCN only)

Same idea as Vis1, but one subplot each for maximum, average, and minimum daily temperature.

Show the code

temp_vars = ['tasmax', 'tas', 'tasmin']
attr_cols = ['TMAX_ATTRIBUTES', 'TAVG_ATTRIBUTES', 'TMIN_ATTRIBUTES']
labels = ['Daily Maximum Temperature', 'Daily Average Temperature', 'Daily Minimum Temperature']

base_colors = ['tab:red', 'tab:orange', 'y']       # primary source data
flagged_colors = ['lightcoral', 'moccasin', 'khaki']  # GSOD ('S') or SSOD ('2') sourced data - dulled

fig, axes = plt.subplots(len(temp_vars), 1, figsize=(15, 10), sharex=True)
dates = station_data.index

for ax, var, attr_col, base_color, flagged_color, label in zip(axes, temp_vars, attr_cols, base_colors, flagged_colors, labels):
    values = station_data[var].values

    # SFLAG is the last character of the attribute string (MFLAG,QFLAG,SFLAG).
    # '2' (SSOD) is NOAA's documented successor to 'S' (GSOD), so both are flagged the same way.
    sflag = station_data[attr_col].str.split(',').str[-1]
    flagged_mask = sflag.isin(['S', '2']).values

    plot_segments(ax, dates, values, flagged_mask, base_color, flagged_color)
    ax.set_ylabel('Temperature (°C)')
    ax.set_title(f'{label} (with GSOD/SSOD data dulled) for Skukuza, RSA')

axes[-1].set_xlabel('Date')
plt.tight_layout()

plt.savefig('../images/Vis2_GHCN_daily_TimeSeries_Temperatures_Skukuza.png', dpi=300, bbox_inches='tight')
plt.show()

Vis3: Annual rainfall - GHCN vs ERA5

Compares yearly rainfall totals from the station against ERA5, with a secondary line showing what percentage of each GHCN year is made up of valid (non-missing) days - low validity years should be read with more caution.

Show the code

ghcn_pr_year = station_data['pr'].resample('YE').sum()
era5_pr_year = daily_rainfall_df['pr'].resample('YE').sum()

# Work out what fraction of each year actually has valid GHCN data
ghcn_valid_counts = station_data['pr'].groupby(station_data.index.year).count()
ghcn_days_per_year = pd.Series({
    y: 366 if pd.Timestamp(f'{y}-12-31').is_leap_year else 365
    for y in ghcn_valid_counts.index
})
ghcn_valid_frac = (ghcn_valid_counts / ghcn_days_per_year) * 100

width = 0.4
fig, ax1 = plt.subplots(figsize=(12, 6))

ghcn_years = ghcn_pr_year.index.year
era5_years = era5_pr_year.index.year

# Treat missing years as zero just for plotting, so bars aren't skipped
ghcn_plot = ghcn_pr_year.fillna(0)
era5_plot = era5_pr_year.fillna(0)

ax1.bar(ghcn_years - width / 2, ghcn_plot.values, width=width, label='Station Rainfall (GHCN)', color='tab:blue')
ax1.bar(era5_years + width / 2, era5_plot.values, width=width, label='ERA5 Rainfall', color='tab:orange')
ax1.set_xlabel('Year')
ax1.set_ylabel('Annual Rainfall (mm)')
ax1.grid(True, linestyle='--', alpha=0.5)

# Secondary axis: % of the year that has valid GHCN data
ax2 = ax1.twinx()
ax2.plot(ghcn_years - width / 2, ghcn_valid_frac.values, label='% Valid (GHCN)', color='black', marker='.', linestyle='--')
ax2.set_ylabel('% Valid', color='black')
ax2.set_ylim(0, 100)

lines_labels = [ax.get_legend_handles_labels() for ax in [ax1, ax2]]
lines, labels = [sum(lol, []) for lol in zip(*lines_labels)]
ax1.legend(lines, labels, loc='upper right')

plt.title('Annual Rainfall Comparison at Skukuza, South Africa')
plt.tight_layout()

plt.savefig('../images/Vis3_GHCN-ERA5_year_TimeSeries_Rainfall_Skukuza.png', dpi=300, bbox_inches='tight')
plt.show()

Vis4: Annual mean temperature - GHCN vs ERA5

Same comparison idea as Vis3, but for temperature (max/average/min), each with its own subplot and its own validity line.

Show the code

min_frac = 0  # minimum fraction of valid daily data required to keep a yearly value

def compute_annual_mean(df, var_name, min_frac):
    """Return the annual mean of a variable, plus what fraction of each year was valid.
    Years below min_frac valid data are set to NaN rather than silently averaged."""
    valid_counts = df[var_name].groupby(df.index.year).count()
    annual_mean = df[var_name].resample('YE').mean()
    days_per_year = pd.Series({
        y: 366 if pd.Timestamp(f'{y}-12-31').is_leap_year else 365
        for y in valid_counts.index
    })
    frac_valid = valid_counts / days_per_year

    annual_mean_filtered = annual_mean.copy()
    annual_mean_filtered[:] = np.nan
    for y, val in zip(annual_mean.index.year, annual_mean.values):
        if frac_valid.get(y, 0) >= min_frac:
            annual_mean_filtered.loc[str(y)] = val
    return annual_mean_filtered, frac_valid

# GHCN annual means
ghcn_tmax_year, ghcn_tmax_valid = compute_annual_mean(station_data, 'tasmax', min_frac)
ghcn_tmean_year, ghcn_tmean_valid = compute_annual_mean(station_data, 'tas', min_frac)
ghcn_tmin_year, ghcn_tmin_valid = compute_annual_mean(station_data, 'tasmin', min_frac)

# ERA5 annual means
era5_tmax_year, _ = compute_annual_mean(daily_tmax_df, 'tasmax', min_frac)
era5_tmean_year, _ = compute_annual_mean(daily_temp_df, 'tas', min_frac)
era5_tmin_year, _ = compute_annual_mean(daily_tmin_df, 'tasmin', min_frac)

# Restrict everything to 1960 onward, matching the temperature record's start
start_year = 1960
ghcn_tmax_year = ghcn_tmax_year[ghcn_tmax_year.index.year >= start_year]
ghcn_tmean_year = ghcn_tmean_year[ghcn_tmean_year.index.year >= start_year]
ghcn_tmin_year = ghcn_tmin_year[ghcn_tmin_year.index.year >= start_year]

ghcn_tmax_valid = ghcn_tmax_valid[ghcn_tmax_valid.index >= start_year]
ghcn_tmean_valid = ghcn_tmean_valid[ghcn_tmean_valid.index >= start_year]
ghcn_tmin_valid = ghcn_tmin_valid[ghcn_tmin_valid.index >= start_year]

era5_tmax_year = era5_tmax_year[era5_tmax_year.index.year >= start_year]
era5_tmean_year = era5_tmean_year[era5_tmean_year.index.year >= start_year]
era5_tmin_year = era5_tmin_year[era5_tmin_year.index.year >= start_year]

fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)

plot_vars = [
    ('tasmax', ghcn_tmax_year, ghcn_tmax_valid, era5_tmax_year, 'Daily Maximum Temperature (°C)', 'tab:red', 'lightcoral'),
    ('tas', ghcn_tmean_year, ghcn_tmean_valid, era5_tmean_year, 'Daily Average Temperature (°C)', 'tab:orange', 'moccasin'),
    ('tasmin', ghcn_tmin_year, ghcn_tmin_valid, era5_tmin_year, 'Daily Minimum Temperature (°C)', 'y', 'khaki')
]

for ax, (var, ghcn_data, ghcn_valid, era5_data, label, base_color, gsod_color) in zip(axes, plot_vars):
    ax.plot(ghcn_data.index.year, ghcn_data.values, color=base_color, linewidth=2.2, label='Station (GHCN)')
    ax.plot(era5_data.index.year, era5_data.values, color=base_color, linewidth=2.2, linestyle='--', label='ERA5')
    ax.set_ylabel('Temperature (°C)')
    ax.set_title(f'Annual Mean {label}')
    ax.grid(True, linestyle='--', alpha=0.4)

    ax2 = ax.twinx()
    ax2.plot(ghcn_valid.index + 0.2, ghcn_valid.values * 100, color='black', marker='.', linestyle='--', label='% Valid (GHCN)')
    ax2.set_ylabel('% Valid', color='black')
    ax2.set_ylim(0, 105)
    ax2.tick_params(axis='y', labelcolor='black')

    lines_labels = [ax.get_legend_handles_labels() for ax in [ax, ax2]]
    lines, labels = [sum(lol, []) for lol in zip(*lines_labels)]
    ax.legend(lines, labels, loc='lower left')

axes[-1].set_xlabel('Year')
axes[-1].set_xlim(1959, max(era5_tmax_year.index.year))

fig.suptitle('Annual Mean Daily Temperature Comparison at Skukuza, South Africa', fontsize=14, y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])

plt.savefig('../images/Vis4_GHCN-ERA5_year_TimeSeries_Temperatures_Skukuza.png', dpi=300, bbox_inches='tight')
plt.show()

Vis5: Monthly rainfall climatology - GHCN vs ERA5

Averages rainfall by calendar month over 1960-2000 to show the typical seasonal pattern, with a secondary line showing the average number of rain days (≥1 mm) per month.

Show the code

ghcn_rain = station_data.loc['1960':'2000', 'pr']
era5_rain = daily_rainfall_df.loc['1960':'2000', 'pr']

ghcn_monthly = ghcn_rain.resample('ME').sum()
era5_monthly = era5_rain.resample('ME').sum()

ghcn_raindays = ghcn_rain.resample('ME').apply(lambda x: (x >= 1).sum())
era5_raindays = era5_rain.resample('ME').apply(lambda x: (x >= 1).sum())

# Average each calendar month across all years to get the typical seasonal pattern
ghcn_clim = ghcn_monthly.groupby(ghcn_monthly.index.month).mean()
era5_clim = era5_monthly.groupby(era5_monthly.index.month).mean()
ghcn_rd_clim = ghcn_raindays.groupby(ghcn_raindays.index.month).mean()
era5_rd_clim = era5_raindays.groupby(era5_raindays.index.month).mean()

months = np.arange(1, 13)
month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']

fig, ax1 = plt.subplots(figsize=(12, 6))
width = 0.35

ax1.bar(months - width / 2, ghcn_clim.values, width, label='GHCN Rainfall', color='tab:blue')
ax1.bar(months + width / 2, era5_clim.values, width, label='ERA5 Rainfall', color='tab:orange')
ax1.set_ylabel('Mean Monthly Rainfall (mm)')
ax1.set_xticks(months)
ax1.set_xticklabels(month_labels)
ax1.set_xlabel('Month')
ax1.grid(True, linestyle='--', alpha=0.4)

ax2 = ax1.twinx()
ax2.plot(months, ghcn_rd_clim.values, color='navy', marker='o', linestyle='-', label='GHCN Rain Days')
ax2.plot(months, era5_rd_clim.values, color='chocolate', marker='o', linestyle='--', label='ERA5 Rain Days')
ax2.set_ylabel('Mean Monthly Rain Days (≥1 mm)')

lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines + lines2, labels + labels2, loc='upper center')

plt.title('Monthly climatology of rainfall (1960-2000) over Skukuza, South Africa')
plt.tight_layout()

plt.savefig('../images/Vis5_GHCN-ERA5_monthly_Climatology_Rainfall_Skukuza.png', dpi=300, bbox_inches='tight')
plt.show()

Vis6: Monthly temperature climatology - GHCN vs ERA5

Same seasonal-average idea as Vis5, but for maximum and minimum temperature.

Show the code

ghcn_tmax = station_data.loc['1960':'2000', 'tasmax']
ghcn_tmin = station_data.loc['1960':'2000', 'tasmin']
era5_tmax = daily_tmax_df.loc['1960':'2000', 'tasmax']
era5_tmin = daily_tmin_df.loc['1960':'2000', 'tasmin']

ghcn_tmax_month = ghcn_tmax.resample('ME').mean()
ghcn_tmin_month = ghcn_tmin.resample('ME').mean()
era5_tmax_month = era5_tmax.resample('ME').mean()
era5_tmin_month = era5_tmin.resample('ME').mean()

ghcn_tmax_clim = ghcn_tmax_month.groupby(ghcn_tmax_month.index.month).mean()
ghcn_tmin_clim = ghcn_tmin_month.groupby(ghcn_tmin_month.index.month).mean()
era5_tmax_clim = era5_tmax_month.groupby(era5_tmax_month.index.month).mean()
era5_tmin_clim = era5_tmin_month.groupby(era5_tmin_month.index.month).mean()

months = np.arange(1, 13)
month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']

fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(months, ghcn_tmax_clim, color='tab:red', marker='o', linestyle='-', label='GHCN Tmax')
ax.plot(months, ghcn_tmin_clim, color='y', marker='o', linestyle='-', label='GHCN Tmin')
ax.plot(months, era5_tmax_clim, color='tab:red', marker='^', linestyle='--', label='ERA5 Tmax')
ax.plot(months, era5_tmin_clim, color='y', marker='^', linestyle='--', label='ERA5 Tmin')

ax.set_xticks(months)
ax.set_xticklabels(month_labels)
ax.set_xlabel('Month')
ax.set_ylabel('Temperature (°C)')
ax.set_title('Monthly climatology of daily maximum and minimum temperature (1960-2000) over Skukuza, South Africa')
ax.grid(True, linestyle='--', alpha=0.4)
ax.legend(loc='lower right')
plt.tight_layout()

plt.savefig('../images/Vis6_GHCN-ERA5_monthly_Climatology_Temperatures_Skukuza.png', dpi=300, bbox_inches='tight')
plt.show()