sciencehub-team/keeling-curvePublic

# /// script
# requires-python = ">=3.11"
# dependencies = ["marimo>=0.25.1", "numpy", "pandas", "matplotlib"]
# ///

import marimo

__generated_with = "0.25.1"
app = marimo.App(width="medium", app_title="The Keeling curve")

with app.setup:
    import numpy as np
    import pandas as pd

    DATA = "data/co2_mm_mlo.csv"


@app.function
def load_co2(path: str = DATA) -> pd.DataFrame:
    """Monthly mean CO₂ in the air at Mauna Loa, Hawaii, in ppm, from March 1958.

    Columns: year, month, decimal_date, average, deseasonalized, ndays, sdev,
    unc. NOAA fills months without measurements by interpolation and marks
    them with a negative ndays.
    """
    co2 = pd.read_csv(path, comment="#")
    return co2.rename(columns={"decimal date": "decimal_date"})


@app.function
def yearly_mean(path: str = DATA) -> pd.Series:
    """The deseasonalized concentration averaged over each complete year, in ppm, by year."""
    co2 = load_co2(path)
    complete = co2.groupby("year").filter(lambda year: len(year) == 12)
    return complete.groupby("year")["deseasonalized"].mean()


@app.function
def growth_rate(start_year: int, end_year: int, path: str = DATA) -> float:
    """Average growth of CO₂ between two complete years, in ppm a year.

    >>> round(growth_rate(1960, 1970), 2)
    0.88
    """
    yearly = yearly_mean(path)
    if start_year not in yearly.index or end_year not in yearly.index:
        raise ValueError(f"Both years need twelve months of data; there are {yearly.index.min()} to {yearly.index.max()}.")
    return float((yearly[end_year] - yearly[start_year]) / (end_year - start_year))


@app.function
def seasonal_cycle(path: str = DATA) -> pd.Series:
    """Each month's mean departure from the deseasonalized trend, in ppm: the breath of the northern forests."""
    co2 = load_co2(path)
    return (co2["average"] - co2["deseasonalized"]).groupby(co2["month"]).mean()


@app.cell
def _():
    import marimo as mo
    import matplotlib.pyplot as plt

    return mo, plt


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    # The Keeling curve

    In March 1958 Charles David Keeling began measuring carbon dioxide in the
    air on the slope of Mauna Loa, Hawaii, far from cities and forests. The
    record has run almost without a break since; NOAA's Global Monitoring
    Laboratory continues it today. It shows two things at once: a seasonal
    rhythm as the land plants of the northern hemisphere grow and decay, and a
    rise that has not paused in any year.

    This project is the data and the notebook together. The file
    `data/co2_mm_mlo.csv` is NOAA's monthly record as published; the functions
    below read it, and anyone can call them from Python with
    `sciencehub.load("sciencehub-team/keeling-curve")`.
    """)
    return


@app.cell
def _():
    co2 = load_co2()
    first_year, last_year = int(co2["year"].min()), int(co2["year"].max())
    return co2, first_year, last_year


@app.cell
def _(first_year, last_year, mo):
    years = mo.ui.range_slider(
        start=first_year,
        stop=last_year,
        step=1,
        value=[first_year, last_year],
        label="Years",
        full_width=True,
    )
    years
    return (years,)


@app.cell
def _(co2, mo, plt, years):
    _start, _end = years.value
    _shown = co2[(co2["year"] >= _start) & (co2["year"] <= _end)]
    fig, ax = plt.subplots(figsize=(7, 3.8))
    ax.plot(_shown["decimal_date"], _shown["average"], color="#0969da", lw=1, label="monthly mean")
    ax.plot(_shown["decimal_date"], _shown["deseasonalized"], color="#1f2328", lw=2, label="without the seasons")
    ax.set_xlabel("year")
    ax.set_ylabel("CO₂ (ppm)")
    ax.legend(frameon=False, loc="upper left")
    for _side in ("top", "right"):
        ax.spines[_side].set_visible(False)
    fig.tight_layout()

    _rise = _shown["deseasonalized"].iloc[-1] - _shown["deseasonalized"].iloc[0]
    mo.vstack(
        [
            fig,
            mo.md(
                f"From {_start} to {_end} the concentration rose by **{_rise:.1f} ppm**, "
                f"from {_shown['deseasonalized'].iloc[0]:.1f} to {_shown['deseasonalized'].iloc[-1]:.1f} ppm."
            ),
        ]
    )
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    ## The rise is speeding up

    The growth rate is the average yearly rise over a decade. In the 1960s the
    concentration grew by under one ppm a year; it now grows by more than two.
    """)
    return


@app.cell
def _(mo, plt):
    _yearly = yearly_mean()
    _decades = list(range(int(_yearly.index.min()) // 10 * 10 + 10, int(_yearly.index.max()), 10))
    _rates = [growth_rate(_d - 10, _d) if _d - 10 in _yearly.index else np.nan for _d in _decades]
    fig2, ax2 = plt.subplots(figsize=(7, 3))
    ax2.bar([f"{_d - 10}s" for _d in _decades], _rates, color="#1f2328")
    ax2.set_ylabel("ppm a year")
    for _side in ("top", "right"):
        ax2.spines[_side].set_visible(False)
    fig2.tight_layout()
    mo.vstack([fig2, mo.md(f"Over the last complete decade: **{_rates[-1]:.2f} ppm a year**.")])
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    ## The seasons

    Every year the concentration falls from May to October, while northern
    forests and crops take up carbon, and rises again through winter. The
    swing at Mauna Loa is about six ppm.
    """)
    return


@app.cell
def _(plt):
    _cycle = seasonal_cycle()
    fig3, ax3 = plt.subplots(figsize=(7, 2.8))
    ax3.bar(["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"], _cycle.values, color=["#9a6700" if _v > 0 else "#0969da" for _v in _cycle.values])
    ax3.axhline(0, color="#1f2328", lw=0.8)
    ax3.set_ylabel("ppm from the trend")
    for _side in ("top", "right"):
        ax3.spines[_side].set_visible(False)
    fig3.tight_layout()
    fig3
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""
    ## Data, limits and references

    **Data.** `data/co2_mm_mlo.csv`, the monthly means from NOAA's Global
    Monitoring Laboratory as published on 9 October 2026, unchanged. Values
    before May 1974 are from the Scripps Institution of Oceanography. Months
    without measurements are interpolated by NOAA (negative `ndays`). From
    December 2022 to July 2023, after the eruption of Mauna Loa, the
    measurements come from Maunakea, 21 miles to the north.

    **Limits.** One site measures the air over the central Pacific, which
    represents the well-mixed atmosphere well but not every region. The
    growth rates here are differences of yearly means, not the smoothed rates
    NOAA publishes.

    **Cite the data.** Lan, X., Tans, P. and Thoning, K. W.: Trends in
    globally-averaged CO₂ determined from NOAA Global Monitoring Laboratory
    measurements, https://gml.noaa.gov/ccgg/trends/; Keeling, C. D. et al.,
    Scripps Institution of Oceanography, https://scrippsco2.ucsd.edu/. NOAA
    makes the data freely available and asks for this credit.

    Keeling, C. D. (1960). The concentration and isotopic abundances of carbon
    dioxide in the atmosphere. *Tellus* 12(2), 200–203.
    """)
    return


if __name__ == "__main__":
    app.run()