""" New York plumbing complaints against the thermometer, 2021 to 2026. Rebuilds the data set and the figures behind the SimpleKPI post https://www.simplekpi.com/blog/when-do-pipes-freeze-new-york-data Sources (both public, no key needed): - NYC Open Data, 311 Service Requests from 2010 to Present (erm2-nwe9), complaint type "PLUMBING" (all routed to HPD). - NOAA NCEI daily summaries (GHCN-Daily), Central Park, station USW00094728. Cite: Menne et al. (2012), Global Historical Climatology Network - Daily (GHCN-Daily), Version 3. "Normal" for a day is the median complaint count on the same weekday at 1 to 4 weeks before and after it (at least 5 of those 8 days must exist). Run: python nyc-plumbing-complaints-vs-temperature-analysis.py Needs Python 3.9+ and nothing outside the standard library. Records added to 311 after 27 September 2026 can shift recent daily counts slightly. Licence: CC BY 4.0, SimpleKPI. Provided as is, without warranty. """ import collections import csv import datetime as dt import json import statistics as st import urllib.parse import urllib.request START, END = "2021-01-01", "2026-09-22" OUT_CSV = "nyc-plumbing-complaints-vs-temperature-2021-2026.csv" TYPES = ["WATER SUPPLY", "BASIN/SINK", "BATHTUB/SHOWER", "TOILET", "RADIATOR", "STEAM PIPE/RISER"] def get_json(url): with urllib.request.urlopen(url, timeout=300) as r: return json.load(r) def fetch_complaints(): """Daily plumbing complaint counts by descriptor, from NYC Open Data.""" q = { "$select": "date_trunc_ymd(created_date) as d, descriptor, count(*) as n", "$where": f"complaint_type='PLUMBING' and created_date >= '{START}T00:00:00'", "$group": "d, descriptor", "$limit": "100000", } rows = get_json("https://data.cityofnewyork.us/resource/erm2-nwe9.json?" + urllib.parse.urlencode(q)) by_day = collections.defaultdict(collections.Counter) for r in rows: by_day[r["d"][:10]][r.get("descriptor") or "OTHER"] += int(r["n"]) return by_day def fetch_weather(): """Daily low, high, rain and snow at Central Park, from NOAA NCEI.""" q = { "dataset": "daily-summaries", "stations": "USW00094728", "dataTypes": "TMAX,TMIN,PRCP,SNOW", "startDate": START, "endDate": END, "units": "standard", "format": "json", } rows = get_json("https://www.ncei.noaa.gov/access/services/data/v1?" + urllib.parse.urlencode(q)) return {r["DATE"]: r for r in rows if r.get("TMIN") not in (None, "") and r.get("TMAX") not in (None, "")} def build_rows(by_day, weather): total = {d: sum(c.values()) for d, c in by_day.items()} days = sorted(d for d in total if START <= d <= END and d in weather) def normal(d): x = dt.date.fromisoformat(d) near = [(x + dt.timedelta(days=k)).isoformat() for k in (-28, -21, -14, -7, 7, 14, 21, 28)] vals = [total[k] for k in near if k in total and k <= END] return st.median(vals) if len(vals) >= 5 else None rows = [] for d in days: b = normal(d) w = weather[d] rows.append({"d": d, "n": total[d], "b": b, "u": total[d] / b if b else None, "tmin": int(w["TMIN"]), "tmax": int(w["TMAX"]), "prcp": w.get("PRCP"), "snow": w.get("SNOW"), "desc": by_day[d]}) return rows def write_csv(rows): with open(OUT_CSV, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["date", "plumbing_complaints", "normal_for_weekday", "ratio_to_normal", "low_f", "high_f", "precip_in", "snow_in", "water_supply", "basin_sink", "bathtub_shower", "toilet", "radiator", "steam_pipe_riser", "other"]) for r in rows: c = r["desc"] other = r["n"] - sum(c[k] for k in TYPES) w.writerow([r["d"], r["n"], "" if r["b"] is None else r["b"], "" if r["u"] is None else round(r["u"], 3), r["tmin"], r["tmax"], r["prcp"], r["snow"]] + [c[k] for k in TYPES] + [other]) def band(rows, key, lo, hi): u = [r["u"] for r in rows if lo <= r[key] < hi] return f"{len(u):5d} days median {st.median(u):.2f}x 30%+ above normal on {sum(x > 1.3 for x in u) / len(u):.0%}" def report(rows): a = [r for r in rows if r["u"] is not None] print(f"{sum(r['n'] for r in rows):,} complaints over {len(rows):,} days; {len(a):,} days with a full normal") print("\nBy the day's low (F):") for lab, lo, hi in [("35+", 35, 200), ("30-34", 30, 35), ("25-29", 25, 30), ("20-24", 20, 25), ("15-19", 15, 20), ("10-14", 10, 15), ("<10", -50, 10)]: print(f" {lab:>6} {band(a, 'tmin', lo, hi)}") print("\nBy the day's high (F):") for lab, lo, hi in [("<32", -50, 32), ("32-49", 32, 50), ("50-69", 50, 70), ("70-79", 70, 80), ("80-89", 80, 90), ("90-94", 90, 95), ("95+", 95, 200)]: print(f" {lab:>6} {band(a, 'tmax', lo, hi)}") # Cold snaps: the low falls below 15F after a day with a low at or above it. idx = {r["d"]: r for r in a} low = {r["d"]: r["tmin"] for r in rows} prev = lambda d: (dt.date.fromisoformat(d) - dt.timedelta(days=1)).isoformat() starts = [r["d"] for r in a if r["tmin"] < 15 and not (prev(r["d"]) in low and low[prev(r["d"])] < 15)] print(f"\n{len(starts)} cold snaps; median ratio to normal by day:") for k in range(-1, 5): u = [idx[x]["u"] for x in ((dt.date.fromisoformat(s) + dt.timedelta(days=k)).isoformat() for s in starts) if x in idx] print(f" day {k:+d}: {st.median(u):.2f}x") print("\nShare of complaints about water supply:") groups = {"normal (low 35F+, high <90F)": [r for r in rows if r["tmin"] >= 35 and r["tmax"] < 90], "freeze (low <15F)": [r for r in rows if r["tmin"] < 15], "heat (high 95F+)": [r for r in rows if r["tmax"] >= 95]} for name, g in groups.items(): ws = sum(r["desc"]["WATER SUPPLY"] for r in g) print(f" {name:30s} {ws / sum(r['n'] for r in g):.1%} ({ws / len(g):.0f} a day)") if __name__ == "__main__": rows = build_rows(fetch_complaints(), fetch_weather()) write_csv(rows) print(f"Wrote {OUT_CSV}\n") report(rows)