"""The hour-plus commute — every number and figure in the article, with the working. py src/article_commute.py This file is the article's notebook. The article keeps to what was found; the how and the why live here, next to the code that computes them. ====================================================================================== 1. The question, and the data that answer it ====================================================================================== A quarter of New York City's commuters spend an hour or more getting to work, one way. Who are they, where do they live, how do they travel, what do they earn, and has it always been like this? Three sources, and the article says which one each number comes from. * **The published journey-to-work tables.** Census tables B08303 (workers by travel time), B08013 (aggregate travel time, which divided by workers gives the mean) and B08301 (workers by means of transportation), for New York City and the United States from the 2024 one-year ACS, and for the city and its five counties from the 2020-2024 five-year ACS (data/reference/tables/acs_commute_*.json). These are the numbers a reader could find on their own; the article reproduces them before going anywhere they do not. * **The census microdata.** The 2020-2024 five-year PUMS for the study area with its 80 person replicate weights is the main file: every count, share and margin labelled "census" is a weighted tabulation of NYC resident workers in it. The 2024 one-year file is the sensitivity check. For the thirty-year view the one-year files for 2005-2024 (data/interim/acs1/), the Census 2000 5% PUMS and the Census 1990 5% PUMS (data/interim/acs1/commute_1990.csv.gz, commute_2000.csv.gz; see extract_commute_history.py) give the same statistics year by year. * **MTA ridership.** Daily subway, bus and bridge-and-tunnel counts through the first days of September 2026, administrative data from data.ny.gov. The census runs through 2024; the ridership series is how the article says where things stand now. It is not a survey and carries no margin, and it counts trips, not commuters. ====================================================================================== 2. Universe and definitions ====================================================================================== * **Workers**: residents of the five boroughs, 16 and over, employed and at work in the survey reference week, civilian or armed forces, so with a means of transportation recorded (JWTRNS not blank). This is the universe of table B08301. * **Commuters**: workers who did not work from home (JWTRNS != 11), the universe of B08303 and B08013, with a travel time (JWMNP). "Commute" means one-way door-to-door minutes for the usual trip, as the respondent reports it. * **Hour-plus**: 60 minutes or more; **ninety-plus**: 90 or more (B08303's own top band). Travel time is self-reported and heaps at 30, 45, 60 and 90. * **Mode**: the census's single "usually" mode, the one used for most of the distance. Multi-leg trips are one mode. * **Earnings**: the worker's own earnings (PERNP) in 2024 dollars, in fixed bands. * **Destination**: place of work, coarse. Manhattan; the worker's own borough; another borough; outside the city. 2020-vintage place-of-work areas: Manhattan 04100, Bronx 04200, Brooklyn 04300, Queens 04400, Staten Island 04500. * **Race/ethnicity**: Hispanic of any race; otherwise White, Black, Asian, Other. * **Occupation groups**: the census occupation codes collapsed into the Bureau's own major groups. * **Leaves before 6 a.m.**: departure-time code 30 or below (5:55-5:59 a.m. is 30). ====================================================================================== 3. Margins ====================================================================================== PUMS ships 80 replicate weights per person. For any statistic θ computed with the full weight, the same statistic is recomputed with each replicate weight r, and Var(θ) = (4 / 80) · Σ_r (θ_r − θ)² Margins are 90% (1.645 SE), the Bureau's convention. Medians are weighted medians without a margin. The 1990 and 2000 census points are 5% samples with no replicate weights and are reported without a margin; a 5% sample of the city's workers is about 150,000 records and the sampling error on a citywide share is well under a point. ====================================================================================== 4. The thirty-year series, and what changed underneath it ====================================================================================== Travel time, means of transportation and place of work have been asked the same way since the 1990 census, so the statistics line up. What does not line up is the geography and the sample: 1990 and 2000 are decennial long-form samples (5% of housing units), 2005-2024 are annual surveys of about 1% with different weighting, and PUMA vintages change in 2012 and 2022 (the five boroughs are always identified, so the city total is unaffected). Means-of-transportation codes were harmonised to the 2019 scheme (extract_commute_history.py). "Worked from home" is a mode in every year. Dollars are not used in the series; earnings bands appear only in the five-year cross-section. Outputs (all regenerable; nothing is edited by hand): output/articles/hour_plus_commute/results.json every number quoted in the prose output/articles/hour_plus_commute/pums_direct.csv the census-side estimates with margins output/articles/hour_plus_commute/reproduction.csv the published-table reproduction output/articles/hour_plus_commute/charts.json cubes behind the interactive charts output/articles/hour_plus_commute/fig*.svg static twins of the same figures output/articles/hour_plus_commute/data.csv companion dataset (re-randomised ids) site/public/articles/hour-plus-commute/ the same, copied for the website """ from __future__ import annotations import argparse import json import shutil from pathlib import Path import numpy as np import pandas as pd from _paths import P, INTERIM, OUTPUT, RAW, REFERENCE import article_commute_charts as CH import article_commute_model as MODEL SLUG = "hour-plus-commute" OUT = OUTPUT / "articles" / "hour_plus_commute" SITE_OUT = P.root / "site" / "public" / "articles" / SLUG Z90 = 1.645 LONG = 60 VLONG = 90 MIN_RECORDS = 30 POW_2020 = {"Manhattan": 4100, "Bronx": 4200, "Brooklyn": 4300, "Queens": 4400, "Staten Island": 4500} BOROUGHS = ["Manhattan", "Brooklyn", "Queens", "Bronx", "Staten Island"] MODE_LABEL = {1: "car", 2: "bus", 3: "subway", 4: "commuter rail", 5: "light rail", 6: "ferry", 7: "taxi", 8: "motorcycle", 9: "bicycle", 10: "walked", 11: "worked from home", 12: "other"} # Mode groups used by the charts: car includes motorcycle; taxi, ferry, light rail and # "other" are grouped as other. def mode_group(code: pd.Series) -> np.ndarray: c = code.to_numpy(dtype=float) return np.select([c == 3, c == 2, np.isin(c, [1, 8]), c == 4, c == 10, c == 9, c == 11], ["subway", "bus", "car", "commuter rail", "walked", "bicycle", "home"], "other") MODE_GROUPS = [("subway", "Subway"), ("bus", "Bus"), ("car", "Car, truck or van"), ("commuter rail", "Commuter rail"), ("walked", "Walked"), ("bicycle", "Bicycle"), ("other", "Taxi, ferry and other")] EARN_BANDS = [("lt30k", "Under $30k", -np.inf, 30_000), ("30k_60k", "$30k–60k", 30_000, 60_000), ("60k_100k", "$60k–100k", 60_000, 100_000), ("100k_150k", "$100k–150k", 100_000, 150_000), ("150k_plus", "$150k and over", 150_000, np.inf)] TIME_BANDS = [("lt15", "Under 15 min", 0, 15), ("15_29", "15–29 min", 15, 30), ("30_44", "30–44 min", 30, 45), ("45_59", "45–59 min", 45, 60), ("60_89", "60–89 min", 60, 90), ("90_plus", "90 min or more", 90, 10_000)] RACES = ["White", "Black", "Hispanic", "Asian", "Other"] DESTS = [("manhattan", "Manhattan"), ("own", "Own borough"), ("other_boro", "Another borough"), ("outside", "Outside the city")] OCC_GROUPS = [(10, 440, "Management"), (500, 750, "Business and finance"), (800, 960, "Business and finance"), (1005, 1240, "Computer and math"), (1305, 1560, "Architecture and engineering"), (1600, 1980, "Science"), (2001, 2060, "Community and social service"), (2100, 2180, "Legal"), (2205, 2555, "Education"), (2600, 2920, "Arts, media and entertainment"), (3000, 3550, "Health practitioners"), (3601, 3655, "Health support"), (3700, 3960, "Protective service"), (4000, 4160, "Food service"), (4200, 4255, "Cleaning and maintenance"), (4330, 4655, "Personal care and service"), (4700, 4965, "Sales"), (5000, 5940, "Office and administrative"), (6005, 6130, "Farming"), (6200, 6765, "Construction"), (6800, 6950, "Extraction"), (7000, 7640, "Installation and repair"), (7700, 8990, "Production"), (9005, 9760, "Transportation and moving"), (9800, 9830, "Military")] INK = "#6f6e69"; GRID = "#d5d4cd"; BLUE = "#2a78d6"; ORANGE = "#eb6834"; TICK = "#1a1a1a"; FIG_W = 6.4 HISTORY_YEARS = [1990, 2000] + list(range(2005, 2025)) COMMUTE_TABLES_1YR = REFERENCE / "tables" / "acs_commute_nyc_us_2024_1yr.json" COMMUTE_TABLES_5YR = REFERENCE / "tables" / "acs_commute_nyc_2024_5yr.json" MTA_CURRENT = RAW / "transit" / "mta_daily_ridership.csv" MTA_FROZEN = RAW / "transit" / "mta_daily_ridership_2020_2025.csv" def occ_group(s: pd.Series) -> pd.Series: c = pd.to_numeric(s, errors="coerce"); out = pd.Series("Other", index=s.index) for lo, hi, lab in OCC_GROUPS: out[(c >= lo) & (c <= hi)] = lab return out def race_label(rac1p, hisp) -> np.ndarray: r = pd.to_numeric(rac1p, errors="coerce").to_numpy(); h = pd.to_numeric(hisp, errors="coerce").to_numpy() return np.where(h > 1, "Hispanic", np.where(r == 1, "White", np.where(r == 2, "Black", np.where(r == 6, "Asian", "Other")))) def earn_band(e: pd.Series) -> np.ndarray: v = e.to_numpy(dtype=float) return np.select([v < 30_000, v < 60_000, v < 100_000, v < 150_000], ["lt30k", "30k_60k", "60k_100k", "100k_150k"], "150k_plus") # ------------------------------------------------------------------ census side def load_pums() -> tuple[pd.DataFrame, np.ndarray]: """NYC resident workers in the five-year file, with replicate weights and the article's derived variables. Only rows with a means of transportation (workers at work) are kept, so the replicate matrix is cut the same way.""" cols = ["SERIALNO", "PUMA", "STATE", "PWGTP", "AGEP", "SEX", "ESR", "JWTRNS", "JWMNP", "JWDP", "POWSP", "POWPUMA", "PERNP", "WKHP", "RAC1P", "HISP", "SCHL", "OCCP", "NATIVITY", "COW", "JWRIP"] p = pd.read_csv(INTERIM / "metro_person.csv.gz", usecols=cols, dtype={"STATE": str, "PUMA": str, "OCCP": str}) rw = np.load(INTERIM / "metro_person_repwts.npy", mmap_mode="r") assert len(rw) == len(p) geo = pd.read_csv(REFERENCE / "geography_pumas.csv", dtype={"state_fips": str, "puma": str}) geo["puma_geoid"] = geo["state_fips"].str.zfill(2) + geo["puma"].str.zfill(5) p["puma_geoid"] = p["STATE"].str.zfill(2) + p["PUMA"].str.zfill(5) p = p.merge(geo[["puma_geoid", "in_nyc", "borough"]], on="puma_geoid", how="left") keep = ((p["in_nyc"] == 1) & p["JWTRNS"].notna()).to_numpy() h = p[keep].reset_index(drop=True); rw = np.asarray(rw[keep]) h = derive(h) return h, rw def derive(h: pd.DataFrame) -> pd.DataFrame: h["mode"] = mode_group(h["JWTRNS"]) h["commuter"] = (h["JWTRNS"] != 11) & h["JWMNP"].notna() h["minutes"] = h["JWMNP"].astype(float) h["long"] = h["commuter"] & (h["minutes"] >= LONG) h["vlong"] = h["commuter"] & (h["minutes"] >= VLONG) h["earn_band"] = earn_band(h["PERNP"].fillna(0)) h["race"] = race_label(h["RAC1P"], h["HISP"]) h["foreign_born"] = (h["NATIVITY"] == 2) h["ba_plus"] = (h["SCHL"] >= 21) h["occ"] = occ_group(h["OCCP"]) pow_state = pd.to_numeric(h["POWSP"], errors="coerce"); pow_puma = pd.to_numeric(h["POWPUMA"], errors="coerce") own_code = h["borough"].map(POW_2020) in_city = (pow_state == 36) & pow_puma.isin(list(POW_2020.values())) h["dest"] = np.select([(pow_state == 36) & (pow_puma == POW_2020["Manhattan"]) & (h["borough"] != "Manhattan"), in_city & (pow_puma == own_code), in_city], ["manhattan", "own", "other_boro"], "outside") h["to_manhattan"] = (pow_state == 36) & (pow_puma == POW_2020["Manhattan"]) h["early"] = pd.to_numeric(h["JWDP"], errors="coerce") <= 30 h["female"] = h["SEX"] == 2 h["full_time"] = h["WKHP"] >= 35 h["time_band"] = pd.cut(h["minutes"], [0, 15, 30, 45, 60, 90, 10_000], right=False, labels=[k for k, *_ in TIME_BANDS]).astype(object) return h class Census: """Weighted estimates with successive-difference-replication margins (90%).""" def __init__(self, h: pd.DataFrame, rw: np.ndarray): self.h, self.rw, self.w = h, rw, h["PWGTP"].to_numpy(dtype=float) def share(self, mask, num) -> tuple[float, float]: mask = np.asarray(mask); num = np.asarray(num) w, rw = self.w[mask], self.rw[mask]; x = num[mask].astype(float) if w.sum() == 0: return float("nan"), float("nan") full = (w * x).sum() / w.sum(); reps = (rw * x[:, None]).sum(0) / rw.sum(0) return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def count(self, mask) -> tuple[float, float]: mask = np.asarray(mask); full = self.w[mask].sum(); reps = self.rw[mask].sum(0) return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def mean(self, mask, col: str) -> tuple[float, float]: mask = np.asarray(mask); v = self.h.loc[mask, col].to_numpy(dtype=float); w = self.w[mask]; rw = self.rw[mask] keep = ~np.isnan(v); v, w, rw = v[keep], w[keep], rw[keep] if w.sum() == 0: return float("nan"), float("nan") full = (w * v).sum() / w.sum(); reps = (rw * v[:, None]).sum(0) / rw.sum(0) return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def median(self, mask, col: str) -> float: mask = np.asarray(mask); v = self.h.loc[mask, col].to_numpy(dtype=float); w = self.w[mask] keep = ~np.isnan(v); v, w = v[keep], w[keep] if w.sum() == 0: return float("nan") o = np.argsort(v); cw = np.cumsum(w[o]); return float(v[o][np.searchsorted(cw, cw[-1] / 2)]) def weighted_median(v, w) -> float: v = np.asarray(v, float); w = np.asarray(w, float); keep = ~np.isnan(v); v, w = v[keep], w[keep] if w.sum() == 0: return float("nan") o = np.argsort(v); cw = np.cumsum(w[o]); return float(v[o][np.searchsorted(cw, cw[-1] / 2)]) # ------------------------------------------------------------------ published tables def published_tables() -> dict: """B08303, B08013 and B08301 as published, for the article's anchors and the reproduction.""" out = {} for key, path in (("one_year_2024", COMMUTE_TABLES_1YR), ("five_year_2020_2024", COMMUTE_TABLES_5YR)): d = json.loads(path.read_text(encoding="utf-8")) geos = {"01000US": "United States", "16000US3651000": "New York City", "05000US36005": "Bronx", "05000US36047": "Brooklyn", "05000US36061": "Manhattan", "05000US36081": "Queens", "05000US36085": "Staten Island"} block = {"release": d["release"]["name"], "table_note": "B08303 workers 16+ who did not work from home, by travel time; " "B08013 aggregate travel time (minutes); B08301 workers 16+ by means of transportation."} for g, label in geos.items(): if g not in d["data"]: continue b = d["data"][g]; t = b["B08303"]["estimate"]; te = b["B08303"]["error"]; a = b["B08013"]["estimate"]; ae = b["B08013"]["error"] m = b["B08301"]["estimate"] total, total_moe = t["B08303001"], te["B08303001"] long_ = t["B08303012"] + t["B08303013"]; long_moe = float(np.sqrt(te["B08303012"] ** 2 + te["B08303013"] ** 2)) row = {"commuters": total, "commuters_moe90": total_moe, "long_commuters": long_, "long_commuters_moe90": long_moe, "long_share": long_ / total, "vlong_share": t["B08303013"] / total, "vlong_commuters": t["B08303013"], "mean_minutes": a["B08013001"] / total, "aggregate_minutes_moe90": ae["B08013001"], "workers": m["B08301001"], "worked_from_home": m["B08301021"], "wfh_share": m["B08301021"] / m["B08301001"], "subway_share": m["B08301013"] / m["B08301001"], "bus_share": m["B08301011"] / m["B08301001"], "car_share": m["B08301002"] / m["B08301001"], "walk_share": m["B08301019"] / m["B08301001"], "bands": {k: t[f"B08303{i:03d}"] for i, k in zip(range(2, 14), ["lt5", "5_9", "10_14", "15_19", "20_24", "25_29", "30_34", "35_39", "40_44", "45_59", "60_89", "90_plus"])}} block[label] = row out[key] = block return out # ------------------------------------------------------------------ estimates def census_estimates(cs: Census) -> tuple[dict, pd.DataFrame]: h = cs.h; rows: list[dict] = [] def put(key, label, value, moe, n, kind): rows.append({"key": key, "label": label, "value": value, "moe90": moe, "n_records": int(n), "kind": kind}) workers = np.ones(len(h), dtype=bool) cm = h["commuter"].to_numpy(); long_ = h["long"].to_numpy(); vlong = h["vlong"].to_numpy() boro = h["borough"].to_numpy(); mode = h["mode"].to_numpy(); eb = h["earn_band"].to_numpy(); race = h["race"].to_numpy() dest = h["dest"].to_numpy(); occ = h["occ"].to_numpy() res: dict = {} def block(mask, with_profile=False): n = int(mask.sum()) c, cmo = cs.count(mask); p, pm = cs.share(mask, long_); v, vm = cs.share(mask, vlong); mu, mum = cs.mean(mask, "minutes") lc, lcm = cs.count(mask & long_) d = {"commuters": c, "commuters_moe90": cmo, "long_share": p, "long_share_moe90": pm, "vlong_share": v, "vlong_share_moe90": vm, "long_commuters": lc, "long_commuters_moe90": lcm, "mean_minutes": mu, "mean_minutes_moe90": mum, "median_minutes": cs.median(mask, "minutes"), "n_records": n} if with_profile: d["to_manhattan_share"] = cs.share(mask, h["to_manhattan"].to_numpy())[0] d["median_earnings"] = cs.median(mask, "PERNP") return d # Workers and commuters. wc, wcm = cs.count(workers); res["workers"] = {"value": wc, "moe90": wcm, "n_records": int(workers.sum())} wfh, wfhm = cs.share(workers, (mode == "home")); res["wfh_share"] = {"value": wfh, "moe90": wfhm} put("workers", "NYC resident workers at work", wc, wcm, workers.sum(), "count") put("wfh_share", "Share of workers who worked from home", wfh, wfhm, workers.sum(), "share") res["all"] = block(cm, True) put("commuters", "NYC resident commuters (did not work from home)", res["all"]["commuters"], res["all"]["commuters_moe90"], cm.sum(), "count") put("long_share", "Share of commuters travelling 60+ minutes", res["all"]["long_share"], res["all"]["long_share_moe90"], cm.sum(), "share") put("vlong_share", "Share of commuters travelling 90+ minutes", res["all"]["vlong_share"], res["all"]["vlong_share_moe90"], cm.sum(), "share") put("mean_minutes", "Mean one-way travel time (minutes)", res["all"]["mean_minutes"], res["all"]["mean_minutes_moe90"], cm.sum(), "mean") # Travel-time bands. res["bands"] = [] tb = h["time_band"].to_numpy() for key, label, lo, hi in TIME_BANDS: mk = cm & (tb == key); s, sm = cs.share(cm, mk); c, cmo = cs.count(mk) res["bands"].append({"band": key, "label": label, "share": s, "share_moe90": sm, "commuters": c, "commuters_moe90": cmo, "n_records": int(mk.sum())}) put(f"band_{key}", f"Share of commuters, {label}", s, sm, cm.sum(), "share") # Borough, mode, borough x mode. res["boroughs"] = {b: block(cm & (boro == b), True) for b in BOROUGHS} for b in BOROUGHS: put(f"boro_{b}", f"60+ minute share, {b} residents", res["boroughs"][b]["long_share"], res["boroughs"][b]["long_share_moe90"], (cm & (boro == b)).sum(), "share") res["modes"] = {} for key, label in MODE_GROUPS: mk = cm & (mode == key); d = block(mk); d["label"] = label; d["share_of_commuters"] = cs.share(cm, mk)[0] res["modes"][key] = d put(f"mode_{key}", f"60+ minute share, {label}", d["long_share"], d["long_share_moe90"], mk.sum(), "share") res["borough_x_mode"] = {b: {k: block(cm & (boro == b) & (mode == k)) for k in ("subway", "bus", "car")} for b in BOROUGHS} res["mode_share_workers"] = {key: cs.share(workers, mode == key)[0] for key, _ in MODE_GROUPS + [("home", "home")]} # Earnings. res["earnings"] = [] for key, label, lo, hi in EARN_BANDS: mk = cm & (eb == key); d = block(mk, True); d.update({"band": key, "label": label, "share_of_commuters": cs.share(cm, mk)[0], "subway_share": cs.share(mk, mode == "subway")[0], "car_share": cs.share(mk, mode == "car")[0], "bus_share": cs.share(mk, mode == "bus")[0], "walk_share": cs.share(mk, mode == "walked")[0]}) res["earnings"].append(d) put(f"earn_{key}", f"60+ minute share, own earnings {label}", d["long_share"], d["long_share_moe90"], mk.sum(), "share") # Earnings within mode: does money buy a shorter subway ride? res["earnings_x_mode"] = {k: [{"band": key, "long_share": cs.share(cm & (eb == key) & (mode == k), long_)[0], "long_share_moe90": cs.share(cm & (eb == key) & (mode == k), long_)[1], "median_minutes": cs.median(cm & (eb == key) & (mode == k), "minutes"), "n_records": int((cm & (eb == key) & (mode == k)).sum())} for key, *_ in EARN_BANDS] for k in ("subway", "bus", "car")} res["earnings_x_borough"] = {b: [{"band": key, "long_share": cs.share(cm & (eb == key) & (boro == b), long_)[0], "long_share_moe90": cs.share(cm & (eb == key) & (boro == b), long_)[1], "n_records": int((cm & (eb == key) & (boro == b)).sum())} for key, *_ in EARN_BANDS] for b in BOROUGHS} # Race. res["race"] = {} for r in RACES: mk = cm & (race == r); d = block(mk, True); d["share_of_commuters"] = cs.share(cm, mk)[0] d["subway_share"] = cs.share(mk, mode == "subway")[0]; d["bus_share"] = cs.share(mk, mode == "bus")[0]; d["car_share"] = cs.share(mk, mode == "car")[0] d["outer_share"] = cs.share(mk, boro != "Manhattan")[0] res["race"][r] = d put(f"race_{r}", f"60+ minute share, {r} commuters", d["long_share"], d["long_share_moe90"], mk.sum(), "share") # Race within borough and mode: how much of the gap is geography? res["race_x_borough"] = {b: {r: {"long_share": cs.share(cm & (boro == b) & (race == r), long_)[0], "long_share_moe90": cs.share(cm & (boro == b) & (race == r), long_)[1], "n_records": int((cm & (boro == b) & (race == r)).sum())} for r in RACES} for b in BOROUGHS} res["race_x_mode"] = {k: {r: {"long_share": cs.share(cm & (mode == k) & (race == r), long_)[0], "long_share_moe90": cs.share(cm & (mode == k) & (race == r), long_)[1], "n_records": int((cm & (mode == k) & (race == r)).sum())} for r in RACES} for k in ("subway", "bus", "car")} # Occupation groups. res["occupations"] = [] for lab in sorted(set(occ)): mk = cm & (occ == lab) if mk.sum() < 500: continue d = block(mk, True); d["label"] = lab; res["occupations"].append(d) res["occupations"].sort(key=lambda d: -d["long_share"]) # Destination. res["destination"] = {} for key, label in DESTS: mk = cm & (dest == key); d = block(mk); d["label"] = label; d["share_of_commuters"] = cs.share(cm, mk)[0] d["share_moe90"] = cs.share(cm, mk)[1]; res["destination"][key] = d res["destination_by_borough"] = {b: {key: cs.share(cm & (boro == b), dest == key)[0] for key, _ in DESTS} for b in BOROUGHS} res["destination_x_borough_long"] = {b: {key: {"long_share": cs.share(cm & (boro == b) & (dest == key), long_)[0], "long_share_moe90": cs.share(cm & (boro == b) & (dest == key), long_)[1], "median_minutes": cs.median(cm & (boro == b) & (dest == key), "minutes"), "n_records": int((cm & (boro == b) & (dest == key)).sum())} for key, _ in DESTS} for b in BOROUGHS} # The long commuters as a group, against the short ones. short = cm & (h["minutes"].to_numpy() < 30) traits = [("outer", "Lives outside Manhattan", boro != "Manhattan"), ("subway", "Takes the subway", mode == "subway"), ("bus", "Takes the bus", mode == "bus"), ("car", "Drives", mode == "car"), ("to_manhattan", "Works in Manhattan", h["to_manhattan"].to_numpy()), ("black_hispanic", "Black or Hispanic", np.isin(race, ["Black", "Hispanic"])), ("foreign_born", "Foreign-born", h["foreign_born"].to_numpy()), ("no_degree", "No bachelor's degree", ~h["ba_plus"].to_numpy()), ("under_60k", "Earns under $60k", np.isin(eb, ["lt30k", "30k_60k"])), ("over_100k", "Earns $100k or more", np.isin(eb, ["100k_150k", "150k_plus"])), ("early", "Leaves before 6 a.m.", h["early"].to_numpy()), ("female", "Women", h["female"].to_numpy()), ("full_time", "Works 35+ hours", h["full_time"].to_numpy())] res["profile"] = [] for key, label, num in traits: pl, plm = cs.share(long_, num); ps, psm = cs.share(short, num); pa, _ = cs.share(cm, num) res["profile"].append({"trait": key, "label": label, "long": pl, "long_moe90": plm, "short": ps, "short_moe90": psm, "all": pa, "distinguishable": bool(abs(pl - ps) > plm + psm)}) res["long_group"] = {"median_earnings": cs.median(long_, "PERNP"), "short_median_earnings": cs.median(short, "PERNP"), "all_median_earnings": cs.median(cm, "PERNP"), "median_age": cs.median(long_, "AGEP"), "share_of_commuters": res["all"]["long_share"], "n_records": int(long_.sum()), "borough_shares": {b: cs.share(long_, boro == b)[0] for b in BOROUGHS}, "mode_shares": {k: cs.share(long_, mode == k)[0] for k, _ in MODE_GROUPS}, "race_shares": {r: cs.share(long_, race == r)[0] for r in RACES}, "all_race_shares": {r: cs.share(cm, race == r)[0] for r in RACES}, "dest_shares": {k: cs.share(long_, dest == k)[0] for k, _ in DESTS}} # Hours on the road: a 60-minute one-way commute, ten trips a week. res["time_budget"] = {"long_weekly_hours": 2 * 60 * 5 / 60, "median_weekly_hours": 2 * res["all"]["median_minutes"] * 5 / 60, "note": "two trips a day, five days a week, at the stated one-way time"} return res, pd.DataFrame(rows) # ------------------------------------------------------------------ the thirty-year series def _boro_of_puma_2010(code: str): c = int(code) if 3701 <= c <= 3710: return "Bronx" if 3801 <= c <= 3810: return "Manhattan" if 3901 <= c <= 3903: return "Staten Island" if 4001 <= c <= 4018: return "Brooklyn" if 4101 <= c <= 4114: return "Queens" return None def load_history_year(year: int) -> tuple[pd.DataFrame, np.ndarray | None]: """One year in the canonical columns: weight, borough, esr, mode (JWTRNS coding), minutes, pow_manhattan, plus replicate weights where the year has them.""" if year <= 2014: d = pd.read_csv(INTERIM / "acs1" / f"commute_{year}.csv.gz") rp = INTERIM / "acs1" / f"commute_{year}_repwts.npy" rw = np.load(rp) if rp.exists() else None d = d[d["mode"].notna()].reset_index(drop=True) if rw is None else d if rw is not None: keep = d["mode"].notna().to_numpy(); d = d[keep].reset_index(drop=True); rw = rw[keep] return d, rw cols = ["PUMA", "STATE", "PWGTP", "ESR", "JWTRNS", "JWMNP", "POWSP", "POWPUMA"] p = pd.read_csv(INTERIM / "acs1" / f"metro_person_{year}.csv.gz", usecols=cols, dtype={"STATE": str, "PUMA": str}) rw = np.load(INTERIM / "acs1" / f"metro_person_{year}_repwts.npy", mmap_mode="r") assert len(rw) == len(p) p["PUMA"] = p["PUMA"].str.zfill(5) if year <= 2021: boro = p["PUMA"].where(p["STATE"].str.zfill(2) == "36").map(lambda c: _boro_of_puma_2010(c) if isinstance(c, str) else None) man = 3800 else: geo = pd.read_csv(REFERENCE / "geography_pumas.csv", dtype={"state_fips": str, "puma": str}) geo = geo[geo["in_nyc"] == 1]; lut = dict(zip(geo["state_fips"].str.zfill(2) + geo["puma"].str.zfill(5), geo["borough"])) boro = (p["STATE"].str.zfill(2) + p["PUMA"]).map(lut); man = 4100 keep = (boro.notna() & p["JWTRNS"].notna()).to_numpy() d = p[keep].reset_index(drop=True); rw = np.asarray(rw[keep]) pow_state = pd.to_numeric(d["POWSP"], errors="coerce"); pow_puma = pd.to_numeric(d["POWPUMA"], errors="coerce") out = pd.DataFrame({"year": year, "borough": boro[keep].to_numpy(), "weight": d["PWGTP"], "esr": d["ESR"], "mode": d["JWTRNS"], "minutes": d["JWMNP"], "pow_manhattan": ((pow_state == 36) & (pow_puma == man)).astype(int)}) return out, rw def history() -> dict: """The same statistics for every year from 1990: workers, commuters, the share working from home, mode shares, the 60+ and 90+ shares, mean and median minutes, by borough and by mode group. Margins where the year has replicate weights.""" out = {"years": [], "note": "1990 and 2000 are 5% decennial samples without replicate weights (no margin); " "2005-2024 are one-year ACS files with 80 replicate weights (90% margins). The 2020 one-year file " "carries the Bureau's experimental weights and is shown with that label."} for year in HISTORY_YEARS: try: d, rw = load_history_year(year) except FileNotFoundError: print(f" history: {year} not on disk, skipped"); continue w = d["weight"].to_numpy(dtype=float); mode = d["mode"].to_numpy(dtype=float); mins = d["minutes"].to_numpy(dtype=float) mg = mode_group(d["mode"]); boro = d["borough"].to_numpy() cm = (mode != 11) & ~np.isnan(mins); long_ = cm & (mins >= LONG); vlong = cm & (mins >= VLONG) def est(mask, num=None, col=None): """share (num within mask), or mean of col within mask, with a margin if rw.""" ww = w[mask] if ww.sum() == 0: return float("nan"), float("nan") if num is not None: x = num[mask].astype(float); full = (ww * x).sum() / ww.sum() if rw is None: return full, float("nan") r = rw[mask]; reps = (r * x[:, None]).sum(0) / r.sum(0) else: x = col[mask]; k = ~np.isnan(x); full = (ww[k] * x[k]).sum() / ww[k].sum() if rw is None: return full, float("nan") r = rw[mask][k]; reps = (r * x[k][:, None]).sum(0) / r.sum(0) return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) def cnt(mask): full = w[mask].sum() if rw is None: return float(full), float("nan") reps = rw[mask].sum(0); return float(full), float(Z90 * np.sqrt(4 / 80 * ((reps - full) ** 2).sum())) row = {"year": year, "source": "Census 5% PUMS" if year <= 2000 else "ACS 1-year PUMS", "n_records": int(cm.sum()), "has_margin": rw is not None, "experimental_weights": year == 2020} row["workers"], row["workers_moe90"] = cnt(np.ones(len(d), bool)) row["commuters"], row["commuters_moe90"] = cnt(cm) row["wfh_share"], row["wfh_share_moe90"] = est(np.ones(len(d), bool), mode == 11) for k in ("subway", "bus", "car", "walked", "commuter rail", "bicycle"): row[f"{k}_share"] = est(cm, mg == k)[0] row["long_share"], row["long_share_moe90"] = est(cm, long_) row["vlong_share"], row["vlong_share_moe90"] = est(cm, vlong) row["mean_minutes"], row["mean_minutes_moe90"] = est(cm, col=mins) row["median_minutes"] = weighted_median(mins[cm], w[cm]) row["to_manhattan_share"] = est(cm, d["pow_manhattan"].to_numpy() == 1)[0] row["boroughs"] = {} for b in BOROUGHS: mk = cm & (boro == b); p_, pm = est(mk, long_); mu, mum = est(mk, col=mins) row["boroughs"][b] = {"long_share": p_, "long_share_moe90": pm, "mean_minutes": mu, "mean_minutes_moe90": mum, "commuters": cnt(mk)[0], "n_records": int(mk.sum())} row["modes"] = {} for k in ("subway", "bus", "car", "walked"): mk = cm & (mg == k); p_, pm = est(mk, long_); mu, mum = est(mk, col=mins) row["modes"][k] = {"long_share": p_, "long_share_moe90": pm, "mean_minutes": mu, "mean_minutes_moe90": mum, "commuters": cnt(mk)[0], "n_records": int(mk.sum())} row["borough_x_mode"] = {b: {k: {"long_share": est(cm & (boro == b) & (mg == k), long_)[0], "long_share_moe90": est(cm & (boro == b) & (mg == k), long_)[1], "mean_minutes": est(cm & (boro == b) & (mg == k), col=mins)[0], "mean_minutes_moe90": est(cm & (boro == b) & (mg == k), col=mins)[1], "commuters": cnt(cm & (boro == b) & (mg == k))[0], "n_records": int((cm & (boro == b) & (mg == k)).sum())} for k in ("subway", "bus", "car", "walked")} for b in BOROUGHS} out["years"].append(row) print(f" {year}: commuters {row['commuters']:,.0f}, 60+ {row['long_share']:.1%}, mean {row['mean_minutes']:.1f}, " f"wfh {row['wfh_share']:.1%}, subway {row['subway_share']:.1%}") return out # ------------------------------------------------------------------ MTA ridership, 2020-2026 def mta_series() -> dict: """Average daily counts by year and day type, set against the MTA's own 2019 baseline recovered from the frozen 2020-2025 file (count / share-of-comparable-day).""" cur = pd.read_csv(MTA_CURRENT); cur["Date"] = pd.to_datetime(cur["Date"]) old = pd.read_csv(MTA_FROZEN); old["Date"] = pd.to_datetime(old["Date"]) modes = {"Subway": ("Subways: Total Estimated Ridership", "Subways: % of Comparable Pre-Pandemic Day"), "Bus": ("Buses: Total Estimated Ridership", "Buses: % of Comparable Pre-Pandemic Day"), "BT": ("Bridges and Tunnels: Total Traffic", "Bridges and Tunnels: % of Comparable Pre-Pandemic Day"), "LIRR": ("LIRR: Total Estimated Ridership", "LIRR: % of Comparable Pre-Pandemic Day"), "MNR": ("Metro-North: Total Estimated Ridership", "Metro-North: % of Comparable Pre-Pandemic Day")} old["weekday"] = old["Date"].dt.dayofweek < 5 base = {} for m, (c, pc) in modes.items(): o = old[(old["Date"].dt.year == 2024) & (old[pc] > 0)] implied = o[c] / o[pc] base[m] = {"weekday": float(implied[o["weekday"]].mean()), "weekend": float(implied[~o["weekday"]].mean())} cur["weekday"] = cur["Date"].dt.dayofweek < 5; cur["year"] = cur["Date"].dt.year; cur["month"] = cur["Date"].dt.to_period("M").astype(str) last = cur["Date"].max() out = {"baseline_note": "2019 baselines are the MTA's own comparable-day figures, recovered from the 2020-2025 file as " "count divided by its percentage-of-comparable-pre-pandemic-day, averaged over 2024 weekdays and weekends.", "last_date": last.strftime("%Y-%m-%d"), "baseline": base, "years": [], "months": []} for (y, wd), g in cur[cur["Mode"].isin(modes)].groupby(["year", "weekday"]): row = {"year": int(y), "day_type": "weekday" if wd else "weekend", "days": int(g["Date"].nunique())} for m in modes: v = g.loc[g["Mode"] == m, "Count"].mean(); row[m] = float(v); row[f"{m}_vs_2019"] = float(v / base[m]["weekday" if wd else "weekend"]) out["years"].append(row) for (mo, wd), g in cur[(cur["Date"] >= "2024-01-01") & cur["Mode"].isin(modes)].groupby(["month", "weekday"]): row = {"month": mo, "day_type": "weekday" if wd else "weekend"} for m in ("Subway", "Bus", "BT"): v = g.loc[g["Mode"] == m, "Count"].mean(); row[m] = float(v); row[f"{m}_vs_2019"] = float(v / base[m]["weekday" if wd else "weekend"]) out["months"].append(row) crz = cur[cur["Mode"].isin(["CBD Entries", "CRZ Entries"]) & cur["weekday"]].groupby(["year", "Mode"])["Count"].mean() out["zone_entries_weekday"] = {f"{y}_{m}": float(v) for (y, m), v in crz.items()} return out # ------------------------------------------------------------------ one-year check def one_year_2024() -> dict: p = pd.read_csv(INTERIM / "acs1" / "metro_person_2024.csv.gz", usecols=["PUMA", "STATE", "PWGTP", "ESR", "JWTRNS", "JWMNP", "POWSP", "POWPUMA", "PERNP", "RAC1P", "HISP", "JWDP", "SEX", "NATIVITY", "SCHL", "OCCP", "WKHP", "AGEP"], dtype={"STATE": str, "PUMA": str, "OCCP": str}) rw = np.load(INTERIM / "acs1" / "metro_person_2024_repwts.npy", mmap_mode="r") geo = pd.read_csv(REFERENCE / "geography_pumas.csv", dtype={"state_fips": str, "puma": str}) geo["puma_geoid"] = geo["state_fips"].str.zfill(2) + geo["puma"].str.zfill(5) p["puma_geoid"] = p["STATE"].str.zfill(2) + p["PUMA"].str.zfill(5) p = p.merge(geo[["puma_geoid", "in_nyc", "borough"]], on="puma_geoid", how="left") keep = ((p["in_nyc"] == 1) & p["JWTRNS"].notna()).to_numpy() h = derive(p[keep].reset_index(drop=True)); rw = np.asarray(rw[keep]); cs = Census(h, rw) cm = h["commuter"].to_numpy(); long_ = h["long"].to_numpy(); vlong = h["vlong"].to_numpy(); boro = h["borough"].to_numpy(); mode = h["mode"].to_numpy() out = {"year": 2024, "n_records": int(cm.sum())} out["workers"], out["workers_moe90"] = cs.count(np.ones(len(h), bool)) out["commuters"], out["commuters_moe90"] = cs.count(cm) out["wfh_share"], out["wfh_share_moe90"] = cs.share(np.ones(len(h), bool), mode == "home") out["long_share"], out["long_share_moe90"] = cs.share(cm, long_) out["vlong_share"], out["vlong_share_moe90"] = cs.share(cm, vlong) out["long_commuters"], out["long_commuters_moe90"] = cs.count(cm & long_) out["mean_minutes"], out["mean_minutes_moe90"] = cs.mean(cm, "minutes") out["median_minutes"] = cs.median(cm, "minutes") out["boroughs"] = {b: {"long_share": cs.share(cm & (boro == b), long_)[0], "long_share_moe90": cs.share(cm & (boro == b), long_)[1], "mean_minutes": cs.mean(cm & (boro == b), "minutes")[0]} for b in BOROUGHS} eb = h["earn_band"].to_numpy(); race = h["race"].to_numpy() out["earnings"] = {k: {"long_share": cs.share(cm & (eb == k), long_)[0], "long_share_moe90": cs.share(cm & (eb == k), long_)[1]} for k, *_ in EARN_BANDS} out["race"] = {r: {"long_share": cs.share(cm & (race == r), long_)[0], "long_share_moe90": cs.share(cm & (race == r), long_)[1]} for r in RACES} out["modes"] = {k: {"long_share": cs.share(cm & (mode == k), long_)[0], "long_share_moe90": cs.share(cm & (mode == k), long_)[1]} for k in ("subway", "bus", "car")} return out # ------------------------------------------------------------------ reproduction def reproduction(pub: dict, census: dict, one: dict) -> pd.DataFrame: rows = [] p1 = pub["one_year_2024"]["New York City"]; p5 = pub["five_year_2020_2024"] def add(label, published, pub_moe, ours, ours_moe, kind): rows.append({"statistic": label, "published": published, "published_moe90": pub_moe, "reproduced": ours, "reproduced_moe90": ours_moe, "difference": ours - published, "within_margin": bool(abs(ours - published) <= max(pub_moe if pub_moe == pub_moe else 0, ours_moe if ours_moe == ours_moe else 0) or abs(ours - published) <= 0.002), "kind": kind}) add("Commuters, NYC, 2024 one-year", p1["commuters"], p1["commuters_moe90"], one["commuters"], one["commuters_moe90"], "count") add("Share travelling 60+ minutes, NYC, 2024 one-year", p1["long_share"], np.nan, one["long_share"], one["long_share_moe90"], "share") add("Share travelling 90+ minutes, NYC, 2024 one-year", p1["vlong_share"], np.nan, one["vlong_share"], one["vlong_share_moe90"], "share") add("Mean travel time (minutes), NYC, 2024 one-year", p1["mean_minutes"], np.nan, one["mean_minutes"], one["mean_minutes_moe90"], "mean") add("Share of workers working from home, NYC, 2024 one-year", p1["wfh_share"], np.nan, one["wfh_share"], one["wfh_share_moe90"], "share") c5 = p5["New York City"] add("Commuters, NYC, 2020-2024 five-year", c5["commuters"], c5["commuters_moe90"], census["all"]["commuters"], census["all"]["commuters_moe90"], "count") add("Share travelling 60+ minutes, NYC, five-year", c5["long_share"], np.nan, census["all"]["long_share"], census["all"]["long_share_moe90"], "share") add("Mean travel time (minutes), NYC, five-year", c5["mean_minutes"], np.nan, census["all"]["mean_minutes"], census["all"]["mean_minutes_moe90"], "mean") for b in BOROUGHS: add(f"Share travelling 60+ minutes, {b}, five-year", p5[b]["long_share"], np.nan, census["boroughs"][b]["long_share"], census["boroughs"][b]["long_share_moe90"], "share") return pd.DataFrame(rows) # ------------------------------------------------------------------ figures def _mpl(): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"font.family": "sans-serif", "font.size": 10, "axes.edgecolor": GRID, "axes.labelcolor": INK, "xtick.color": INK, "ytick.color": INK, "text.color": INK, "svg.fonttype": "none"}) return plt def _style(ax, title, sub): ax.set_title(title, loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate(sub, xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") for s in ("top", "right", "left"): ax.spines[s].set_visible(False) ax.grid(axis="x", color=GRID, linewidth=0.8); ax.set_axisbelow(True); ax.tick_params(length=0) def _save(fig, path: Path): path.parent.mkdir(parents=True, exist_ok=True) fig.savefig(path, format="svg", bbox_inches="tight"); (path.parent / "pdf").mkdir(exist_ok=True) fig.savefig(path.parent / "pdf" / path.with_suffix(".pdf").name, format="pdf", bbox_inches="tight") def _hbar(rows, path, title, sub, xmax=None, fmt="{:.0f}%", xlab_ticks=None, height=None): plt = _mpl(); n = len(rows) fig, ax = plt.subplots(figsize=(FIG_W, height or (1.0 + 0.36 * n)), dpi=100); fig.patch.set_alpha(0) y = np.arange(n)[::-1]; v = np.array([r[1] for r in rows]); m = np.array([r[2] if r[2] == r[2] else 0 for r in rows]) ax.barh(y, v, color=BLUE, height=0.6); ax.errorbar(v, y, xerr=m, fmt="none", ecolor=TICK, elinewidth=1.5) for yi, vi in zip(y, v): ax.text(vi + (xmax or max(v)) * 0.012, yi, fmt.format(vi), va="center", fontsize=10, color=INK) ax.set_yticks(y); ax.set_yticklabels([r[0] for r in rows]); ax.set_xlim(0, xmax or max(v + m) * 1.25) if xlab_ticks: ax.set_xticks(xlab_ticks[0]); ax.set_xticklabels(xlab_ticks[1]) _style(ax, title, sub); fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_bands(c, path): rows = [(b["label"], b["share"] * 100, b["share_moe90"] * 100) for b in c["bands"]] _hbar(rows, path, "How long the trip to work takes", "Share of NYC resident commuters by one-way travel time (census)", xmax=40, xlab_ticks=([0, 10, 20, 30, 40], ["0%", "10%", "20%", "30%", "40%"])) def fig_boroughs(c, path): plt = _mpl(); modes = ("subway", "bus", "car"); cols = [BLUE, ORANGE, "#4a9c6d"] fig, ax = plt.subplots(figsize=(FIG_W, 3.6), dpi=100); fig.patch.set_alpha(0) y = np.arange(len(BOROUGHS))[::-1]; hgt = 0.24 for j, (k, col) in enumerate(zip(modes, cols)): v = np.array([c["borough_x_mode"][b][k]["long_share"] * 100 for b in BOROUGHS]); m = np.array([c["borough_x_mode"][b][k]["long_share_moe90"] * 100 for b in BOROUGHS]) yy = y + (1 - j) * (hgt + 0.02) ax.barh(yy, v, height=hgt, color=col, label=MODE_GROUPS[[x[0] for x in MODE_GROUPS].index(k)][1]) ax.errorbar(v, yy, xerr=m, fmt="none", ecolor=TICK, elinewidth=1.2) for yi, vi in zip(yy, v): ax.text(vi + 1, yi, f"{vi:.0f}%", va="center", fontsize=8.5, color=INK) ax.set_yticks(y); ax.set_yticklabels(BOROUGHS); ax.set_xlim(0, 100); ax.set_xticks([0, 25, 50, 75, 100]); ax.set_xticklabels(["0%", "25%", "50%", "75%", "100%"]) _style(ax, "Where the hour-plus commute lives", "Share of commuters travelling 60 minutes or more, by borough and mode (census)") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=3, frameon=False, fontsize=9, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_modes(c, path): rows = [(c["modes"][k]["label"], c["modes"][k]["long_share"] * 100, c["modes"][k]["long_share_moe90"] * 100) for k, _ in MODE_GROUPS] _hbar(rows, path, "How they travel, and how long it takes", "Share of NYC commuters travelling 60+ minutes, by usual mode (census)", xmax=80, xlab_ticks=([0, 20, 40, 60, 80], ["0%", "20%", "40%", "60%", "80%"])) def fig_earnings(c, path): rows = [(e["label"], e["long_share"] * 100, e["long_share_moe90"] * 100) for e in c["earnings"]] _hbar(rows, path, "Money buys a shorter commute only at the top", "Share of NYC commuters travelling 60+ minutes, by own earnings (census)", xmax=40, xlab_ticks=([0, 10, 20, 30, 40], ["0%", "10%", "20%", "30%", "40%"])) def fig_race(c, path): rows = [(r, c["race"][r]["long_share"] * 100, c["race"][r]["long_share_moe90"] * 100) for r in RACES] _hbar(rows, path, "The hour-plus commute by race and ethnicity", "Share of NYC commuters travelling 60+ minutes (census)", xmax=40, xlab_ticks=([0, 10, 20, 30, 40], ["0%", "10%", "20%", "30%", "40%"])) def fig_occupations(c, path): rows = [(o["label"], o["long_share"] * 100, o["long_share_moe90"] * 100) for o in c["occupations"]] _hbar(rows, path, "Who spends the hour", "Share of NYC commuters travelling 60+ minutes, by occupation group (census)", xmax=45, xlab_ticks=([0, 10, 20, 30, 40], ["0%", "10%", "20%", "30%", "40%"])) def fig_destination(c, path): rows = [(c["destination"][k]["label"], c["destination"][k]["long_share"] * 100, c["destination"][k]["long_share_moe90"] * 100) for k, _ in DESTS] _hbar(rows, path, "Where the trip goes", "Share of NYC commuters travelling 60+ minutes, by place of work (census)", xmax=40, xlab_ticks=([0, 10, 20, 30, 40], ["0%", "10%", "20%", "30%", "40%"])) def fig_history(hist, path): plt = _mpl(); ys = hist["years"] fig, ax = plt.subplots(figsize=(FIG_W, 3.4), dpi=100); fig.patch.set_alpha(0) x = [r["year"] for r in ys]; v = [r["long_share"] * 100 for r in ys]; m = [(r["long_share_moe90"] or 0) * 100 if r["has_margin"] else 0 for r in ys] ax.plot(x, v, color=BLUE, linewidth=2, marker="o", markersize=4) ax.fill_between(x, np.array(v) - np.array(m), np.array(v) + np.array(m), color=BLUE, alpha=0.15, linewidth=0) ax.set_ylim(0, 35); ax.set_yticks([0, 10, 20, 30]); ax.set_yticklabels(["0%", "10%", "20%", "30%"]) ax.set_xticks([1990, 2000, 2010, 2020]); ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for s in ("top", "right"): ax.spines[s].set_visible(False) ax.set_title("Thirty-six years of the hour-plus commute", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate("Share of NYC commuters travelling 60+ minutes; 1990 and 2000 census, then annual ACS (band: 90% margin)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") fig.tight_layout(); _save(fig, path); plt.close(fig) def fig_transit(mta, path): plt = _mpl(); rows = [r for r in mta["years"] if r["day_type"] == "weekday" and r["year"] >= 2020] fig, ax = plt.subplots(figsize=(FIG_W, 3.2), dpi=100); fig.patch.set_alpha(0) x = np.arange(len(rows)); w = 0.26 for j, (m, lab, col) in enumerate((("Subway", "Subway", BLUE), ("Bus", "Bus", ORANGE), ("BT", "Bridges and tunnels", "#4a9c6d"))): v = [r[f"{m}_vs_2019"] * 100 for r in rows] ax.bar(x + (j - 1) * w, v, width=w, color=col, label=lab) ax.axhline(100, color=TICK, linewidth=0.8, linestyle="--") ax.set_xticks(x); ax.set_xticklabels([str(r["year"]) + ("\n(to Sep)" if r["year"] == 2026 else "") for r in rows], fontsize=9) ax.set_ylim(0, 115); ax.set_yticks([0, 25, 50, 75, 100]); ax.set_yticklabels(["0%", "25%", "50%", "75%", "100%"]) ax.grid(axis="y", color=GRID, linewidth=0.8); ax.grid(axis="x", visible=False) for s in ("top", "right"): ax.spines[s].set_visible(False) ax.set_title("Where transit stands after the census stops", loc="left", fontsize=12, color="#1a1a1a", pad=22) ax.annotate("Average weekday count as a share of the MTA's 2019 comparable-day baseline (administrative data)", xy=(0, 1), xycoords="axes fraction", xytext=(0, 6), textcoords="offset points", fontsize=9.5, color=INK, va="bottom", ha="left") ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.16), ncol=3, frameon=False, fontsize=9, labelcolor=INK) fig.tight_layout(); _save(fig, path); plt.close(fig) # ------------------------------------------------------------------- companion data def companion_dataset(h: pd.DataFrame, rng: np.random.Generator) -> pd.DataFrame: """The census records the article is computed from: one row per NYC resident worker at work in the five-year file, with the survey weight and the derived variables. Public-use microdata, already disclosure-protected; the serial number is dropped and rows are shuffled.""" df = pd.DataFrame({ "weight": h["PWGTP"].astype(int), "borough": h["borough"], "worked_from_home": (h["mode"] == "home"), "mode": h["JWTRNS"].map(MODE_LABEL), "mode_group": h["mode"], "one_way_minutes": h["minutes"], "hour_plus": h["long"], "ninety_plus": h["vlong"], "place_of_work": h["dest"].map(dict(DESTS)), "works_in_manhattan": h["to_manhattan"], "own_earnings_2024": h["PERNP"].round(0), "earnings_band": h["earn_band"].map({k: l for k, l, *_ in EARN_BANDS}), "race_ethnicity": h["race"], "foreign_born": h["foreign_born"], "bachelors_or_higher": h["ba_plus"], "occupation_group": h["occ"], "usual_hours": h["WKHP"], "leaves_before_6am": h["early"], "departure_code": h["JWDP"], "age": h["AGEP"], "sex": h["SEX"].map({1: "male", 2: "female"}), }) df.insert(0, "row_id", [f"hpc-{i:06d}" for i in rng.permutation(len(df)) + 1]) return df.sort_values("row_id").reset_index(drop=True) # ------------------------------------------------------------------------- main def main(argv=None) -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--seed", type=int, default=20240101) ap.add_argument("--skip-history", action="store_true") args = ap.parse_args(argv) rng = np.random.default_rng(args.seed); OUT.mkdir(parents=True, exist_ok=True) print("published tables ...") pub = published_tables() print("census side: five-year PUMS with replicate weights ...") h, rw = load_pums(); cs = Census(h, rw) census, pums_rows = census_estimates(cs) print(f" commuters {census['all']['commuters']:,.0f}; 60+ {census['all']['long_share']:.1%} ± {census['all']['long_share_moe90']:.1%}; " f"mean {census['all']['mean_minutes']:.1f} min") print("holding everything else fixed: weighted logistic model with 80 replicate refits ...") census["adjusted"] = MODEL.adjusted(cs) for key in ("race", "earn_band", "mode", "borough"): f = census["adjusted"]["factors"][key] print(" " + key + ": " + "; ".join(f"{r['label']} {r['ame_points']:+.1f} (raw {r['raw_gap_points']:+.1f})" for r in f["rows"])) print("2024 one-year check ...") one = one_year_2024() print(f" 1-year 60+ {one['long_share']:.1%} ± {one['long_share_moe90']:.1%}; mean {one['mean_minutes']:.1f}; published mean {pub['one_year_2024']['New York City']['mean_minutes']:.1f}") print("reproduction ...") repro = reproduction(pub, census, one) print(repro[["statistic", "published", "reproduced", "within_margin"]].to_string(index=False)) hist = {"years": []} if args.skip_history else history() print("MTA ridership ...") mta = mta_series() for r in mta["years"]: if r["day_type"] == "weekday": print(f" {r['year']} weekday: subway {r['Subway_vs_2019']:.0%}, bus {r['Bus_vs_2019']:.0%}, bridges {r['BT_vs_2019']:.0%}") res = {"article": SLUG, "as_of": "2026-09-06", "definitions": {"long_minutes": LONG, "vlong_minutes": VLONG, "universe": "NYC resident workers 16+ at work in the reference week; commuters exclude those who worked from home", "dollar_year": 2024, "census_margin": "90% (successive difference replication, 80 replicate weights)", "modes": dict(MODE_GROUPS), "earnings_bands": {k: l for k, l, *_ in EARN_BANDS}, "destinations": dict(DESTS)}, "published": pub, "census": census, "one_year_2024": one, "history": hist, "mta": mta, "reproduction": repro.to_dict(orient="records"), "seed": args.seed} print("figures ...") fig_bands(census, OUT / "fig1_travel_time_bands.svg") fig_boroughs(census, OUT / "fig2_borough_by_mode.svg") fig_modes(census, OUT / "fig2b_modes.svg") fig_earnings(census, OUT / "fig3_earnings.svg") fig_race(census, OUT / "fig4_race.svg") fig_occupations(census, OUT / "fig5_occupations.svg") fig_destination(census, OUT / "fig6_destination.svg") if hist["years"]: fig_history(hist, OUT / "fig7_history.svg") fig_transit(mta, OUT / "fig8_transit_2026.svg") print("interactive chart cubes ...") cs.adjusted = census["adjusted"] charts = CH.build(cs, hist, mta) CH.write(OUT / "charts.json", charts) print(f" {len(charts['charts'])} charts, {sum(len(c['cells']) for c in charts['charts'].values())} filter combinations, {(OUT / 'charts.json').stat().st_size / 1024:.0f} KB") companion_dataset(h, rng).to_csv(OUT / "data.csv", index=False) pums_rows.to_csv(OUT / "pums_direct.csv", index=False) repro.to_csv(OUT / "reproduction.csv", index=False) (OUT / "results.json").write_text(json.dumps(res, indent=1, default=float), encoding="utf-8") SITE_OUT.mkdir(parents=True, exist_ok=True) for f in OUT.glob("*"): if f.suffix in (".svg", ".csv", ".json"): shutil.copy2(f, SITE_OUT / f.name) shutil.copy2(Path(__file__), SITE_OUT / "analysis.py") print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}") return 0 if __name__ == "__main__": raise SystemExit(main())