Plot the global land surface temperature anomaly from Berkeley Earth

plot-temperature-anomaly.py

#!/usr/bin/env python
# coding: utf-8
import matplotlib.pyplot as plt
import matplotlib as mpl
import pandas as pd
import numpy as np
import cmocean
import os
import scprep


# download data
os.system(
    "curl http://berkeleyearth.lbl.gov/auto/Global/Complete_TAVG_complete.txt | tail -n +35 > temp-anomaly_berkeley-project.csv"
)
# remove leading spaces
os.system('sed -i "s/^\\s*//" temp-anomaly_berkeley-project.csv')
# replace spaces with commas to convert to csv
os.system('sed -i "s/\\s\\s*/,/g" temp-anomaly_berkeley-project.csv')


df = pd.read_csv("temp-anomaly_berkeley-project.csv")
# name the columns ourselves as the native headers were multilevel
df.columns = [
    "year",
    "month",
    "anomaly",
    "uncertainty",
    "anomaly_year",
    "uncertainty_year",
    "anomaly_5year",
    "uncertainty_5year",
    "anomaly_10year",
    "uncertainty_10year",
    "anomaly_20year",
    "uncertainty_20year",
]

df["time"] = df["year"] + df["month"] / 12

# if we're missing the 20-year average, modern day uncertainty is good enough to take 5 or 10
for measure in ["anomaly", "uncertainty"]:
    for year in ["10", "5"]:
        df["{}_20year".format(measure)] = np.where(
            np.isfinite(df["{}_20year".format(measure)]) | (df["year"] < 2000),
            df["{}_20year".format(measure)],
            df["{}_{}year".format(measure, year)],
        )

# subset only important columns
df = df[["time", "anomaly_20year", "uncertainty_20year"]]
# remove years when we don't have reliable data (pre-1770)
df = df.dropna().reset_index(drop=True)


# define a color indicator between zero and one
c = df["anomaly_20year"] - np.min(df["anomaly_20year"])
c = c / np.max(c)
# rescale to between 0.2 and 0.8 for better contrast
c = c * 0.6 + 0.2
# use cmocean's beautiful colormaps
c = cmocean.cm.balance(c)


with plt.style.context("dark_background"):
    fig, ax = plt.subplots(figsize=(16, 9), frameon=False)
    # change color every 5 years
    step = 5
    for time in range(0, len(df) - 1, step):
        plot_df = df[time : time + step + 1]
        # use the mean color for this time period
        plot_c = np.mean(c[time : time + step + 1], axis=0)
        # plot the temperature anomaly
        ax.plot(plot_df["time"], plot_df["anomaly_20year"], c=plot_c, linewidth=3)
    # plot the uncertainty all in one color
    ax.fill_between(
        df["time"],
        df["anomaly_20year"] - df["uncertainty_20year"],
        df["anomaly_20year"] + df["uncertainty_20year"],
        color=np.mean(c, axis=0),
        linewidth=0.0,
        alpha=0.4,
    )
    # tighten plot range
    ax.set_xlim(np.min(df["time"]), np.max(df["time"]))
    # remove axes
    ax.set_xticks([])
    ax.set_yticks([])
    for spine in ax.spines:
        ax.spines[spine].set_visible(False)
    # label every 20 years
    for year in range(1770, 2020, 20):
        # place the label below the lowest point within a 5-year radius
        text_y = np.min(
            df["anomaly_20year"].loc[(df["time"] < year + 5) & (df["time"] > year - 5)]
            - 0.15
        )
        # color the label the same as the surrounding decade
        text_c = np.mean(
            c[df.index[(df["time"] < year + 5) & (df["time"] > year - 5)].to_numpy()],
            axis=0,
        )
        ax.text(
            year,
            text_y,
            str(year),
            horizontalalignment="center",
            verticalalignment="center",
            fontdict={"size": 20, "color": text_c, "weight": "bold"},
        )

    # set the background to a black-red gradient
    black = np.array([0, 0, 0, 1])
    # choose a color 95% along the balance colormap
    darkred = cmocean.cm.balance(np.linspace(0, 1, 21))[-2]
    # create a linear colormap between them
    cmap = scprep.plot.tools.create_colormap([black, darkred])
    # generate and plot as an image
    background = cmap(np.linspace(0, 1, 200)**2)[None,:]
    ax.imshow(background, aspect='auto', extent=(*ax.get_xlim(), *ax.get_ylim()), origin='lower')

    # save
    fig.tight_layout()
    fig.savefig(
        "temp-anomaly.png", dpi=300, frameon=False, bbox_inches="tight", pad_inches=0
    )