Temperature plot script

weather.py

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests

# --- pick fewer years here ---
YEARS = [2016, 2017, 2018]          # e.g. [2024, 2025] for 2 years
OUT_PNG = f"us_county_tmax_tmin_hist_CONUS_{min(YEARS)}-{max(YEARS)}.png"

BASE = "https://www.ncei.noaa.gov/pub/data/daily-grids/v1-0-0/averages"
CACHE_DIR = "./nclimgrid_county_cache"
os.makedirs(CACHE_DIR, exist_ok=True)

def url_for(variable, y, m):
    yyyymm = f"{y}{m:02d}"
    return f"{BASE}/{y}/{variable}-{yyyymm}-cty-scaled.csv"

def local_path(variable, y, m):
    yyyymm = f"{y}{m:02d}"
    return os.path.join(CACHE_DIR, f"{variable}-{yyyymm}-cty-scaled.csv")

# 1°F bins (adjust if you want)
bins = np.arange(-60, 131, 1)
hist_tmax = np.zeros(len(bins) - 1, dtype=np.int64)
hist_tmin = np.zeros(len(bins) - 1, dtype=np.int64)

stats = {
    "tmax": {
        "global_min": np.inf,
        "global_max": -np.inf,
        "count_below_0": 0,
        "count_above_100": 0,
        "n_vals": 0,
    },
    "tmin": {
        "global_min": np.inf,
        "global_max": -np.inf,
        "count_below_0": 0,
        "count_above_100": 0,
        "n_vals": 0,
    },
}

for variable, hist in [("tmax", hist_tmax), ("tmin", hist_tmin)]:
    for y in YEARS:
        for m in range(1, 13):
            fp = local_path(variable, y, m)
            if not os.path.exists(fp):
                r = requests.get(url_for(variable, y, m), timeout=120)
                r.raise_for_status()
                with open(fp, "wb") as f:
                    f.write(r.content)

            # Daily values live in columns 6..36 (31 day columns). Missing days use sentinels.
            df = pd.read_csv(
                fp,
                header=None,
                usecols=range(6, 37),
                na_values=[-999.99, -999.9, -9999, -9999.0, ""],
                engine="c",
            )

            vals_c = df.to_numpy(dtype=np.float32).ravel()
            vals_c = vals_c[np.isfinite(vals_c)]
            if vals_c.size == 0:
                continue

            vals_f = vals_c * 9 / 5 + 32  # °C -> °F

            stats[variable]["n_vals"] += vals_f.size
            stats[variable]["global_min"] = min(stats[variable]["global_min"], float(vals_f.min()))
            stats[variable]["global_max"] = max(stats[variable]["global_max"], float(vals_f.max()))
            stats[variable]["count_below_0"] += int((vals_f < 0).sum())
            stats[variable]["count_above_100"] += int((vals_f > 100).sum())

            hist += np.histogram(vals_f, bins=bins)[0]

# --- “Reddit-ready” plot (no explicit color picking) ---
plt.rcParams.update({
    "savefig.dpi": 300,
    "font.size": 12,
    "axes.titlesize": 18,
    "axes.labelsize": 12,
})

fig, ax = plt.subplots(figsize=(12, 7))
bin_centers = (bins[:-1] + bins[1:]) / 2
ax.hist(
    [bin_centers, bin_centers],
    bins=bins,
    weights=[hist_tmin, hist_tmax],
    histtype="barstacked",
    color=["purple", "#ff7f0e"],
    label=["Daily low (TMIN)", "Daily high (TMAX)"],
    rwidth=1.0,
)

ax.set_yscale("linear")
ax.set_xlim(-25, 125)
ax.set_xlabel("County area-average daily temperature (°F)")
ax.set_ylabel("County-day count")

period = f"Jan {min(YEARS)}–Dec {max(YEARS)}"
ax.set_title(
    f"How hot and cold do U.S. counties get? Daily highs and lows across CONUS\n"
    f"{period}  •  TMAX n = {stats['tmax']['n_vals']:,}, TMIN n = {stats['tmin']['n_vals']:,} county-days"
)

ax.axvline(32, linestyle="--", linewidth=1)
ax.axvline(100, linestyle="--", linewidth=1)
ax.legend()

ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.grid(True, axis="y", which="both", linewidth=0.6, alpha=0.25)

fig.text(
    0.01, 0.01,
    "Source: NOAA NCEI nClimGrid-Daily v1.0.0 (county area-averages, CONUS), variables=TMAX/TMIN.\n"
    "Notes: County area-average daily highs/lows (not station extremes). Alaska & Hawaii not included.",
    ha="left", va="bottom", fontsize=9
)

fig.tight_layout(rect=[0, 0.05, 1, 1])
fig.savefig(OUT_PNG)
print("Saved:", OUT_PNG)