Case Study 6.3 - Airshed PM2.5 Compliance

np.random.seed(103)
days = pd.date_range("2023-01-01", periods=365, freq="D")
doy = np.arange(365)
# synthesise daily PM2.5: baseline + seasonal swing + a late-year bump + noise,
# np.clip(..., 2, None) just keeps it physically positive
pm25 = np.clip(
    30
    + 12 * np.sin(2 * np.pi * (doy - 40) / 365)     # seasonal (drier months dustier)
    + (doy > 300) * 10                                 # boolean * 10 -> +10 only after day 300
    + np.random.normal(0, 8, 365),
    2, None,
)
air = pd.DataFrame({"date": days, "pm25": np.round(pm25, 1)})

daily_guideline = 50.0     # ug/m3, 24-hr limit
annual_standard = 25.0     # ug/m3, annual-mean limit
exceed_days = int((air["pm25"] > daily_guideline).sum())   # count of days over the daily limit
annual_mean = air["pm25"].mean()

print(f"Days above the 24-hr guideline ({daily_guideline}): {exceed_days} of 365")
print(f"Annual mean PM2.5: {annual_mean:.1f} ug/m3  (standard {annual_standard}) -> "
      f"{'FAIL' if annual_mean > annual_standard else 'PASS'}")   # inline if/else in an f-string
# .idxmax() -> index label of the worst day; .loc fetches that row's date
print(f"Worst day: {air.loc[air['pm25'].idxmax(), 'date'].date()} at {air['pm25'].max():.1f} ug/m3")

# a 30-day rolling mean makes the slow end-of-year rise visible under the daily noise
air["roll30"] = air["pm25"].rolling(30, min_periods=1).mean()
fig, ax = plt.subplots(figsize=(11, 4.5))
ax.plot(air["date"], air["pm25"], color="gray", alpha=0.5, label="Daily PM2.5")
ax.plot(air["date"], air["roll30"], color="darkgreen", label="30-day mean")
ax.axhline(daily_guideline, color="red", linestyle="--", label="24-hr guideline")
ax.set_ylabel("PM2.5 (ug/m3)"); ax.set_title("Fenceline PM2.5, 2023"); ax.legend()
plt.tight_layout(); plt.show()