"""Twenty months of the toll: traffic into Manhattan's core since the congestion toll began on 5 January 2025, set against thirty years of traffic before it — every number and figure in the article, with the working. py src/article_toll.py py src/article_toll.py --manifest-only # rewrite build.json after the viz script 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 ====================================================================================== "Is the traffic creeping back?" is the question the toll's anniversary raised. The data that answer it come from counting programmes that were never designed to be joined, so each is read on its own terms and none is spliced onto another. * **The toll's own count** (administrative). The MTA's toll gantries count the vehicles they detect entering the Congestion Relief Zone (Manhattan south of and including 60th Street, less the FDR Drive and the West Side Highway), from the first day of the toll, 5 January 2025, to the latest day published (data/raw/transit/mta_crz_entries_daily.csv, summed on the server by day, vehicle class, crossing and toll period). There is no count before the toll in this file, so what it can answer is 2026 against 2025. * **Thirty years of the cordon count** (administrative). NYMTC's Hub Bound reports count the vehicles entering Manhattan south of 60th Street on a typical fall business day, every year from 1993 and in 1990 and 1992 before that (Table 1A of the 2024 edition, with the 2024 count from Table 16). Road tubes at 60th Street and on the East River bridges over a two-week October window, toll transactions at the Hudson tunnels. The report's own long table prints two wrong totals, corrected below from its one-year tables, and has two known breaks (the 2019 Brooklyn Bridge count, the Queensboro Bridge from 2021). * **Speed in the core** (administrative). The city's performance indicator: the average speed of yellow taxis carrying passengers, 8 a.m. to 6 p.m. on weekdays, south of 60th Street, one figure a fiscal year (July to June), fiscal 2009-2012 and 2016-2026 from the Mayor's Management Report's open files; DOT's Mobility Reports give the calendar years 2010-2018 by the same taxi GPS. The MTA's own monthly taxi and ride-hail speeds by zone run from October 2019 to August 2025, and its bus speeds on segments inside the zone from 2023 to August 2026. * **Ride-hail** (administrative). The TLC's monthly aggregate report: trips a day and vehicles for yellow taxis (from 2010), green cabs (from August 2013) and the app companies (high-volume for-hire services, from 2015), citywide. * **Street space** (administrative). DOT's bike-route file, which keeps retired segments with their install and retirement dates, and its bus-lane file (lanes in effect now, with the year each began), both placed inside the MTA's zone polygon; the outdoor-dining files for the curb lane. * **Published figures**, each checked against the saved page it comes from: the governor's anniversary and September 2025 releases, the MTA's first-weeks figures, the toll schedule, the state surcharge, DOT's Mobility Report, the city's and consultants' ride-hail studies. ====================================================================================== 2. Definitions ====================================================================================== * **Entries** are vehicle crossings detected by the gantries: a vehicle entering twice counts twice, and a taxi entering empty counts. Tolled-zone entries leave out the FDR Drive and West Side Highway through-traffic, which the gantries also count ("excluded roadways"); the CBD count the MTA's releases quote includes them. * **Clean days** for the year-on-year comparison: not a federal holiday, Good Friday, the day after Thanksgiving or the week between Christmas and New Year (a calendar fixed in advance), and within a quarter of the median of its month and day type (the two 2026 snowstorms and their travel ban fall out by that rule). The second rule reads the counts themselves, so it is checked: from March, the window of every pooled comparison, it drops no day the calendar has not (results: crz.compare.screened_from_march), and every weekday kept in is reported as a sensitivity. Weekdays and weekends are compared separately. * **Vehicle classes** are the MTA's. Class 1, cars, pickups and vans, leaves out the taxis and for-hire cars billed through the per-trip charge; those form the "TLC Taxi/FHV" class, and a TLC vehicle not billed per trip is counted in its ordinary class (the entries file's column description, saved as reference/mta_crz_entries_dictionary.json). Neither class is "private cars". * **Hub Bound vehicles** are autos, taxis, vans and trucks entering the Hub in 24 hours on the count day; buses are counted apart. App-based for-hire cars are not separated from autos. Sectors are the crossings: 60th Street (with the FDR Drive and West Side Highway), Brooklyn, Queens, New Jersey. * **Bike lanes**: DOT facility class I (protected), II (painted) and III (shared or signed route), centreline miles on streets (off-street paths apart), on 31 December of each year, inside the zone polygon. ====================================================================================== 3. Uncertainty ====================================================================================== None of these series is a sample survey, so there is no sampling margin. What varies is the day: weather, events, a bridge closure. The 2026-against-2025 comparisons therefore carry a 90% matched-week bootstrap interval: 7-day blocks counted from 5 January in each year, the same block drawn from both years together (fixed seed), so a March week of 2025 is always set against the same March week of 2026 and the seasonal matching survives the resampling (an outside review of v0.9.1 pointed out that resampling each year's weeks separately did not). It describes how much the comparison moves with the weeks that happened to fall in the window. Comparisons across counting programmes (the October cordon count against October gantry days) are stated as estimates with their method differences. Outputs (all regenerable; nothing is edited by hand): output/articles/twenty_months_of_the_toll/results.json every number quoted in the prose output/articles/twenty_months_of_the_toll/charts.json cubes behind the bar charts in the questions output/articles/twenty_months_of_the_toll/monthly_cube.json monthly entries for the template chart output/articles/twenty_months_of_the_toll/data.csv gantry entries by day, class, crossing and period output/articles/twenty_months_of_the_toll/long_view.csv the long series, one row per source, series and year output/articles/twenty_months_of_the_toll/build.json version and checksums of the served files site/public/articles/twenty-months-of-the-toll/ the same, copied for the website """ from __future__ import annotations import argparse import hashlib import html as _html import json import re import shutil from pathlib import Path import numpy as np import pandas as pd from pandas.tseries.holiday import USFederalHolidayCalendar from _paths import P, OUTPUT, RAW, REFERENCE SLUG = "twenty-months-of-the-toll" OUT = OUTPUT / "articles" / "twenty_months_of_the_toll" SITE_OUT = P.root / "site" / "public" / "articles" / SLUG T = RAW / "transit" REPORTS = REFERENCE / "reports" HB = T / "nymtc_hub_bound_2024_tables" TOLL_START = pd.Timestamp("2025-01-05") SEED = 20260928 MONTHS = ["January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December"] MON3 = [m[:3] for m in MONTHS] CLASSES = [("all", "All vehicles"), ("cars", "Cars, pickups and vans"), ("tlc", "Taxis and for-hire cars on the per-trip charge"), ("trucks_single", "Single-unit trucks"), ("trucks_multi", "Multi-unit trucks"), ("buses", "Buses"), ("motorcycles", "Motorcycles")] CLASS_OF = {"1 - Cars, Pickups and Vans": "cars", "2 - Single-Unit Trucks": "trucks_single", "3 - Multi-Unit Trucks": "trucks_multi", "4 - Buses": "buses", "5 - Motorcycles": "motorcycles", "TLC Taxi/FHV": "tlc"} REGIONS = [("all", "Every crossing"), ("New Jersey", "From New Jersey"), ("Queens", "From Queens"), ("Brooklyn", "From Brooklyn"), ("East 60th St", "Across 60th Street, East Side"), ("West 60th St", "Across 60th Street, West Side")] EXCLUDED_REGIONS = ("FDR Drive", "West Side Highway") # through-traffic on the two highways: counted, not tolled PERIODS = [("all", "All day"), ("Peak", "Peak hours"), ("Overnight", "Overnight")] DAYTYPES = [("weekday", "Weekdays"), ("weekend", "Weekends")] # Gantry crossings grouped into the Hub Bound sectors (the cordon count includes the two highways). HUB_SECTOR_OF = {"East 60th St": "60th", "West 60th St": "60th", "FDR Drive": "60th", "West Side Highway": "60th", "Brooklyn": "brooklyn", "Queens": "queens", "New Jersey": "nj"} SECTORS = [("60th", "60th Street (with the FDR Drive and West Side Highway)"), ("brooklyn", "Brooklyn"), ("queens", "Queens"), ("nj", "New Jersey")] res: dict = {} def _r(x, nd=1): """Round for the results file; a missing value stays missing.""" if x is None: return None x = float(x) if x != x or abs(x) == float("inf"): return None return round(x, nd) def month_label(ym: str) -> str: y, m = ym.split("-"); return f"{MONTHS[int(m) - 1]} {y}" def pct(a, b) -> float: """Percentage change from a to b.""" return 100.0 * (float(b) / float(a) - 1.0) # ==================================================================================== # 0. Quoted figures: every number taken from a report is found in the saved copy first # ==================================================================================== _TEXT: dict = {} def _norm(s: str) -> str: return re.sub(r"\s+", " ", s.replace("’", "'").replace(" ", " ")) def page_texts(path: Path) -> list[str]: """Text of a saved source, one entry per PDF page (an HTML page is one entry).""" if path not in _TEXT: if path.suffix.lower() == ".pdf": import fitz doc = fitz.open(path) _TEXT[path] = [_norm(doc[i].get_text()) for i in range(doc.page_count)] else: raw = path.read_bytes() if raw[:2] == b"\x1f\x8b": import gzip raw = gzip.decompress(raw) t = raw.decode("utf-8", "replace") t = re.sub(r"(?s)<(script|style).*?", " ", t) _TEXT[path] = [_norm(_html.unescape(re.sub(r"<[^>]+>", " ", t)))] return _TEXT[path] def quote(source: str, fname: str, needle: str, page: int | None = None) -> dict: """Assert that `needle` appears in the saved source (on the given PDF page, 1-based) and return a short citation with its context.""" pages = page_texts(REPORTS / fname) where = [page - 1] if page else range(len(pages)) for i in where: k = pages[i].find(needle) if k >= 0: return {"source": source, "page": (i + 1) if fname.endswith(".pdf") else None, "quote": pages[i][max(0, k - 120): k + len(needle) + 120]} raise AssertionError(f"{source}: {needle!r} not found in {fname}" + (f" page {page}" if page else "")) def published() -> dict: """The published figures the article quotes or reconciles, each verified in its file.""" q = {} q["anniversary"] = [quote("governor_crz_anniversary_2026", "governor_crz_anniversary_2026.html", "27 million fewer vehicles entering the Congestion Relief Zone"), quote("governor_crz_anniversary_2026", "governor_crz_anniversary_2026.html", "over 73,000 fewer vehicles are entering the zone, an 11 percent reduction on average"), quote("governor_crz_anniversary_2026", "governor_crz_anniversary_2026.html", "has generated over $550 million in net revenue in its first year")] q["september"] = [quote("governor_crz_sep_2025", "governor_crz_sep_2025.html", "down by 12 percent"), quote("governor_crz_sep_2025", "governor_crz_sep_2025.html", "87,000 fewer vehicles enter the zone"), quote("governor_crz_sep_2025", "governor_crz_sep_2025.html", "17.6 million fewer vehicles have entered the zone compared to last year")] q["streetsblog"] = [quote("streetsblog_crz_first_weeks_2025", "streetsblog_crz_first_weeks_2025.html", "average of 538,955 vehicle entries"), quote("streetsblog_crz_first_weeks_2025", "streetsblog_crz_first_weeks_2025.html", "average of 556,381 vehicles"), quote("streetsblog_crz_first_weeks_2025", "streetsblog_crz_first_weeks_2025.html", "before the traffic toll started was 583,000"), quote("streetsblog_crz_first_weeks_2025", "streetsblog_crz_first_weeks_2025.html", "Hub Bound travel reports")] q["tolls"] = [quote("mta_crz_toll_phasing_2024", "mta_crz_toll_phasing_2024.html", "Automobiles will be charged $9 during peak periods and $2.25 overnight"), quote("mta_crz_toll_phasing_2024", "mta_crz_toll_phasing_2024.html", "$12 for automobiles during peak times and $3 overnight"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "from 5 a.m. to 9 p.m. on weekdays, and from 9 a.m. to 9 p.m. on weekends"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "the per-trip charge for high-volume for-hire vehicles is $1.50"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "the per-trip charge is $0.75"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "up to $3 for passenger vehicles"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "four tolled entries: Lincoln Tunnel, Holland Tunnel, Queens-Midtown Tunnel, and Hugh L. Carey"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "Instead of paying the daily toll, taxis and for-hire vehicles licensed with the NYC Taxi & Limousine Commission are eligible for a smaller per-trip charge"), quote("mta_crz_tolling_page_2026", "mta_crz_tolling_page_2026.html", "increase to $12 in 2028 and then $15 in 2031"), quote("nys_dtf_congestion_surcharge", "nys_dtf_congestion_surcharge.html", "$2.75 for each for-hire transportation trip"), quote("nys_dtf_congestion_surcharge", "nys_dtf_congestion_surcharge.html", "$2.50 for each trip when the transportation is provided by a medallion taxicab"), quote("nys_dtf_congestion_surcharge", "nys_dtf_congestion_surcharge.html", "south of and excluding 96th Street"), quote("nys_dtf_notice_n19_2", "nys_dtf_notice_n19_2.pdf", "12:01 a.m. on Saturday, February 2, 2019", 1), quote("tlc_commissioners_corner_2018_09", "tlc_commissioners_corner_2018_09.pdf", "On August 14, Mayor Bill de Blasio signed into law", 1)] q["ride_hail"] = [quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "Total taxi/TNC weekday mileage in the CBD increased by 36 percent from 2013 to 2017", 5), quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "The number of taxi/TNC vehicles in the CBD increased by 59 percent", 5), quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "Vehicles increased more rapidly than mileage due to slower traffic speeds", 5), quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "Taxi trips declined from 378,000 to 250,000 on an average weekday", 10), quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "TNC trips increased from virtually none to 202,000 trips per day", 10), quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "bringing average June speeds down to 6.8 mph from 8.2 mph in 2013", 10), quote("schaller_empty_seats_2017", "schaller_empty_seats_2017.pdf", "Uber trip volumes were about 1 percent of those for taxis", 27), quote("tlc_dot_fhv_study_2019", "tlc_dot_fhv_study_2019.pdf", "dropping from 6.1 mph in November 2010 to 4.3 mph in November 2018", 3), quote("tlc_dot_fhv_study_2019", "tlc_dot_fhv_study_2019.pdf", "In Manhattan, FHVs now make up nearly 30% of all traffic", 3), quote("tlc_dot_fhv_study_2019", "tlc_dot_fhv_study_2019.pdf", "Though not the only cause, the explosive growth of the for-hire vehicle (FHV) sector", 3), quote("tlc_dot_fhv_study_2019", "tlc_dot_fhv_study_2019.pdf", "to over 120,000 in 2019, is certainly an important factor", 3), quote("tlc_dot_fhv_study_2019", "tlc_dot_fhv_study_2019.pdf", "FHVs and taxis together make up roughly half of the traffic sampled in Manhattan", 22)] q["hub_bound"] = [quote("nymtc_hub_bound_2024_report", "nymtc_hub_bound_2024.pdf", "Beginning from 2021, vehicles that use the northern Queensboro bridge exit were excluded"), quote("nymtc_hub_bound_2024_report", "nymtc_hub_bound_2024.pdf", "In 2018, the Ed Koch Queensboro Bridge volume is noticeably lower than"), quote("nymtc_hub_bound_2024_report", "nymtc_hub_bound_2024.pdf", "Due to the loss of NYMTC's files and databases on Sept. 11, 2001, data for 1999 is largely absent", 10), quote("nymtc_hub_bound_2024_report", "nymtc_hub_bound_2024.pdf", "specifically on October 16, 2024"), quote("nymtc_hub_bound_2019_addendum", "nymtc_hub_bound_2019_addendum.pdf", "TOTAL 64,560 56,943", 4)] q["dot_speeds"] = [quote("nyc_dot_mobility_report_2019", "nyc_dot_mobility_report_2019.pdf", "are down by 22% since 2010", 7), quote("nyc_dot_mobility_report_2019", "nyc_dot_mobility_report_2019.pdf", "7.0 mph", 10)] q["mmr"] = [quote("nyc_mmr_fy2026_dot", "mmr_fy2026_dot.pdf", "Average vehicular travel speed in the Manhattan Central Business District 8.4 7.8 6.9 6.9 6.5", 7)] # Other speed evidence, each read from its saved copy (added after the outside review of v0.9). q["speed_other"] = [ quote("governor_crz_anniversary_2026", "governor_crz_anniversary_2026.html", "Weekday vehicle speeds in the CRZ were up 4 percent compared to 2024, with weekends seeing a 6.2 percent improvement"), quote("governor_crz_anniversary_2026", "governor_crz_anniversary_2026.html", "Within the CRZ, bus speeds increased 2.3 percent"), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "March 2025, Revised August 2026", 1), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "They have not been peer-reviewed", 1), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "The toll raised cordon-area speeds by 11%", 2), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "Speeds also rose outside the cordon, especially on roads frequently traversed by cordon-bound drivers", 2), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "The speed gains attenuate from 15% in the first four months to 10% thereafter", 4), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "The five control cities are Philadelphia, Boston, Chicago, Atlanta, and Baltimore", 9), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "The sample covers traffic conditions from September 2024 through June 2026", 10), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "the policy's causal impact can fully account for the change in raw average speeds on CBD segments", 13), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "under the assumptions of the GSC estimator", 13), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "The average speed on road segments in the NYC CBD increased from 7.1 mph in the four months preceding the policy's implementation to 7.7 mph in the eighteen months after implementation", 13), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "all analyses focus on data during priced hours (5am", 10), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "Because our data come from trips routed through Google Maps rather than a census of all cars, they are ill-suited for measuring changes in volumes", 5), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "Cody Cook Yale University", 2), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "Aboudy Kreidieh Google Research", 2), quote("nber_w33584_network_effects_2026", "nber_w33584_network_effects_2026.pdf", "Shoshana Vasserman Stanford University", 2), quote("tobin_yale_congestion_brief_2026", "tobin_yale_congestion_brief_2026.html", "Introducing congestion pricing made driving in the CBD about 11% faster"), quote("tobin_yale_congestion_brief_2026", "tobin_yale_congestion_brief_2026.html", "have persisted since the policy's launch"), ] q["baseline"] = [ quote("mta_crz_week_one_update_2025", "mta_crz_week_one_update_2025.pdf", "2022-2023-2024 actual/observed entries (October) 642,000", 7), quote("mta_crz_week_one_update_2025", "mta_crz_week_one_update_2025.pdf", "2022-2023-2024 entries (adjusted for January) 583,000", 7), ] # The MTA's First Evaluation Report (January 2026): how the baseline behind "11 percent" and # "27 million" was built, and that it covers the zone and the two highways together. ev = "mta_crz_first_year_report_2026.pdf" q["evaluation"] = [ quote("mta_crz_first_year_report_2026", ev, "FIRST EVALUATION REPORT January 2026", 1), quote("mta_crz_first_year_report_2026", ev, "for accuracy in comparison, certain sections of this chapter combine CRZ entries and excluded roadway entries", 14), quote("mta_crz_first_year_report_2026", ev, "monthly adjustment factors, accounting for these seasonal trends, were calculated to adjust October data", 15), quote("mta_crz_first_year_report_2026", ev, "These adjustment factors were calculated from inbound vehicle volumes on bridges and tunnels in 2022, 2023, and 2024", 15), quote("mta_crz_first_year_report_2026", ev, "applied to the approximately 640,000 October entries from the NYMTC Hub Bound Travel Data Report 2022 and 2023", 15), quote("mta_crz_first_year_report_2026", ev, "entries to the CRZ and excluded roadways are approximately 11 percent lower on average than the baseline", 16), quote("mta_crz_first_year_report_2026", ev, "totaling more than 21.5 million fewer entries between January 5, 2025, and October 31, 2025", 16), quote("mta_crz_first_year_report_2026", ev, "did not consistently differentiate between trips into the CRZ and those made solely on the excluded roadways", 17), ] p16 = page_texts(REPORTS / ev)[15] base_k = re.search(r"((?:\d{3}k ){11}\d{3}k)", p16) assert base_k, "the report's monthly baselines (Figure 2-2) were not found on page 16" baseline_k = [int(x[:-1]) for x in base_k.group(1).split()] table = {} for m in re.finditer(r"(January|February|March|April|May|June|July|August|September|October) 2025 ([\d,]+) ([\d,]+) -([\d,]+) -(\d+)% -([\d,]+)", page_texts(REPORTS / ev)[16]): table[MONTHS.index(m.group(1)) + 1] = (int(m.group(2).replace(",", "")), int(m.group(3).replace(",", ""))) assert len(table) == 10 and len(baseline_k) == 12, "Table 2-1 or Figure 2-2 of the evaluation report changed shape" q["purpose"] = [quote("fhwa_cbd_tolling_fonsi_2023", "fhwa_cbd_tolling_fonsi_2023.pdf", "The Project purpose is to reduce traffic congestion in the Manhattan CBD in a manner that will generate revenue for future transportation improvements", 9)] # What the vehicle classes are, from the entries file's own column description (an outside # review of v0.9.2): class 1 leaves out the taxis and for-hire cars billed per trip, and a # TLC vehicle not billed per trip is counted in its ordinary class. So the "TLC Taxi/FHV" # class is the per-trip class, not every TLC-plated vehicle, and "cars, pickups and vans" # are not private cars. dic = json.loads((REFERENCE / "mta_crz_entries_dictionary.json").read_text(encoding="utf-8")) vc = _norm(next(c["description"] for c in dic["columns"] if c["fieldName"] == "vehicle_class")) for needle in ("Not including TLC Taxi/ FHV vehicles, which are subject to a per-trip toll", "Taxis and FHVs not charged in a Per-Trip Charge Plan (PTCP) are tolled by vehicle class and appear in those categories"): assert needle in vc, f"mta_crz_entries_dictionary: {needle!r} not in the vehicle_class description" q["classes"] = [{"source": "mta_crz_entries_dictionary", "page": None, "quote": vc}] return {"quotes": q, "anniversary": {"fewer_total": 27_000_000, "fewer_per_day": 73_000, "pct": 11.0, "net_revenue_year_one": 550_000_000, "date": "2026-01-05"}, "september": {"fewer_per_day": 87_000, "pct": 12.0, "fewer_total_to_date": 17_600_000, "date": "2025-09-09"}, "streetsblog": {"cbd_week1": 538_955, "cbd_week2": 556_381, "baseline_jan": 583_000, "date": "2025-01-27"}, "tolls": {"car_peak": 9.0, "car_overnight": 2.25, "car_peak_2028": 12.0, "car_overnight_2028": 3.0, "car_peak_2031": 15.0, "hvfhv_trip": 1.50, "taxi_trip": 0.75, "crossing_credit_car": 3.0, "surcharge_fhv": 2.75, "surcharge_taxi": 2.50, "surcharge_pool": 0.75, "surcharge_start": "2019-02-02", "fhv_licence_pause": "2018-08-14"}, "ride_hail": {"cbd_miles_2013_2017_pct": 36, "cbd_vehicles_2013_2017_pct": 59, "taxi_trips_2013": 378_000, "taxi_trips_2017": 250_000, "tnc_trips_2017": 202_000, "june_speed_2013": 8.2, "june_speed_2017": 6.8, "midtown_nov2010": 6.1, "midtown_nov2018": 4.3, "fhv_share_manhattan": 30}, "brooklyn_bridge_2019_video_inbound": 64_560, "brooklyn_bridge_2019_video_outbound": 56_943, "mta_first_year_speeds": {"weekday_crz_pct": 4.0, "weekend_crz_pct": 6.2, "bus_crz_pct": 2.3}, "nber": {"att_pct": 11, "first_four_months_pct": 15, "thereafter_pct": 10, "sample": "September 2024 to June 2026", "raw_cbd_mph_before": 7.1, "raw_cbd_mph_after": 7.7, "hours": "priced hours, 5am-9pm weekdays and 9am-9pm weekends", "controls": ["Philadelphia", "Boston", "Chicago", "Atlanta", "Baltimore"], "revised": "August 2026"}, "mta_baseline": {"october_2022_2024": 642_000, "january_adjusted": 583_000}, "evaluation": {"october_hub_bound": 640_000, "pct_jan_oct": 11, "fewer_jan_oct": 21_500_000, "baseline_k": baseline_k, "table": table, "date": "2026-01"}} # ==================================================================================== # A. The toll's own count: gantry entries, 5 January 2025 to the latest date # ==================================================================================== def holidays() -> set: """Days on which traffic is not a normal working day: federal holidays, the week between Christmas and New Year, Good Friday and the day after Thanksgiving.""" cal = USFederalHolidayCalendar().holidays("2024-01-01", "2026-12-31") extra = [] for y in (2024, 2025, 2026): extra += [f"{y}-12-{dd}" for dd in range(24, 32)] + [f"{y}-01-01", f"{y}-07-04"] extra += ["2024-03-29", "2025-04-18", "2026-04-03"] # Good Friday extra += ["2024-11-29", "2025-11-28", "2026-11-27"] # the day after Thanksgiving return set(cal) | set(pd.to_datetime(extra)) def load_entries() -> pd.DataFrame: e = pd.read_csv(T / "mta_crz_entries_daily.csv"); e["d"] = pd.to_datetime(e["toll_date"]) e["cls"] = e["vehicle_class"].map(CLASS_OF) assert e["cls"].notna().all(), f"unknown vehicle class: {set(e.loc[e['cls'].isna(), 'vehicle_class'])}" known = {k for k, _ in REGIONS[1:]} | set(EXCLUDED_REGIONS) assert set(e["detection_region"]) <= known, f"unknown region: {set(e['detection_region']) - known}" return e def daily_totals(e: pd.DataFrame) -> pd.DataFrame: """One row per day: tolled-zone entries (CRZ), the untolled highway through-traffic, and their sum (the CBD count the MTA's releases use), with the flags that decide which days enter the year-on-year comparison.""" d = e.groupby("d")[["crz_entries", "excluded_roadway_entries"]].sum() d["cbd"] = d["crz_entries"] + d["excluded_roadway_entries"] d["weekday"] = d.index.dayofweek < 5; d["year"] = d.index.year; d["month"] = d.index.month d["holiday"] = d.index.isin(holidays()) med = d.groupby(["year", "month", "weekday"])["crz_entries"].transform("median") d["ratio_to_month_median"] = d["crz_entries"] / med d["disrupted"] = (d["ratio_to_month_median"] < 0.75) | (d["ratio_to_month_median"] > 1.25) d["clean"] = ~d["holiday"] & ~d["disrupted"] gaps = pd.date_range(d.index.min(), d.index.max()).difference(d.index) assert gaps.size == 0, f"days missing from the entries file: {list(gaps[:5])}" return d def same_window(idx: pd.DatetimeIndex, last: pd.Timestamp) -> np.ndarray: """Days comparable across 2025 and 2026: on or after 5 January and no later in the year than the latest day of 2026 in the file.""" doy_ok = (idx.month < last.month) | ((idx.month == last.month) & (idx.day <= last.day)) return ((idx.month > 1) | (idx.day >= TOLL_START.day)) & doy_ok & idx.year.isin([2025, 2026]) UNPAIRED_DAYS = [0] # days in matched-week blocks that only one year has (see boot_ratio) def boot_ratio(x25: pd.Series, x26: pd.Series, rng: np.random.Generator, n: int = 4000) -> tuple[float, float, float]: """Ratio of the 2026 mean to the 2025 mean, minus one, with a 90% matched-week bootstrap interval: the days are cut into 7-day blocks counted from 5 January in each year, and the same block is drawn from both years together, so the resampling keeps the seasonal matching the comparison is built on and every block holds one of each day of the week in both years. (ISO week numbers were tried first: they pair 5 January 2025, a Sunday, with no 2026 day and leave the last days of a window unmatched, which put some point estimates outside their own ranges.) The point estimate is computed on the same matched blocks as the resampling, so the two describe one estimator (an outside review of v0.9.2 noted that a block present in only one year had stayed in the point estimate); the days such blocks hold are tallied in UNPAIRED_DAYS and reported (crz.compare.unpaired_days). They arise where a block's few days inside a window are weekdays in one year and weekend days in the other: at month edges in the monthly comparisons, and possibly at the end of the window. In the file to 19 September 2026 there are seven, all in monthly comparisons; matching them moved the August weekday change from 4.0% to 3.7% (it had sat at the edge of its own range) and four weekend months by 0.2 to 0.7 points.""" def weeks(x: pd.Series) -> pd.DataFrame: k = ((x.index.dayofyear - TOLL_START.dayofyear) // 7).to_numpy() return pd.DataFrame({"k": k, "v": x.to_numpy()}).groupby("k")["v"].agg(["sum", "count"]) w25, w26 = weeks(x25), weeks(x26) common = w25.index.intersection(w26.index) a, b = w25.loc[common].to_numpy(), w26.loc[common].to_numpy() UNPAIRED_DAYS[0] += int(w25.drop(common)["count"].sum() + w26.drop(common)["count"].sum()) full = (b[:, 0].sum() / b[:, 1].sum()) / (a[:, 0].sum() / a[:, 1].sum()) - 1 i = rng.integers(0, len(common), (n, len(common))) m25 = a[i, 0].sum(1) / a[i, 1].sum(1); m26 = b[i, 0].sum(1) / b[i, 1].sum(1) lo, hi = np.percentile(m26 / m25 - 1, [5, 95]) return float(full), float(lo), float(hi) def _chg(f, lo, hi, a=None, b=None) -> dict: out = {"change_pct": _r(100 * f), "lo_pct": _r(100 * lo), "hi_pct": _r(100 * hi)} if a is not None: out.update({"mean_2025": _r(a.mean(), 0), "mean_2026": _r(b.mean(), 0), "diff_per_day": _r(b.mean() - a.mean(), 0), "days_2025": int(len(a)), "days_2026": int(len(b))}) return out def year_one(d: pd.DataFrame, e: pd.DataFrame) -> dict: y1 = d.loc["2025-01-05":"2026-01-04"] e1 = e[(e["d"] >= "2025-01-05") & (e["d"] <= "2026-01-04")] wd = e1[e1["d"].dt.dayofweek < 5] tlc_share = wd.loc[wd["cls"] == "tlc", "crz_entries"].sum() / wd["crz_entries"].sum() car_share = wd.loc[wd["cls"] == "cars", "crz_entries"].sum() / wd["crz_entries"].sum() return {"days": int(len(y1)), "crz_mean": _r(y1["crz_entries"].mean(), 0), "cbd_mean": _r(y1["cbd"].mean(), 0), "crz_weekday_mean": _r(y1.loc[y1.weekday, "crz_entries"].mean(), 0), "crz_weekend_mean": _r(y1.loc[~y1.weekday, "crz_entries"].mean(), 0), "crz_total": _r(y1["crz_entries"].sum(), 0), "cbd_total": _r(y1["cbd"].sum(), 0), "weekday_tlc_share_pct": _r(100 * tlc_share), "weekday_car_share_pct": _r(100 * car_share)} def reconcile(d: pd.DataFrame, pub: dict) -> dict: """What the entries file reproduces of the published figures, and what it cannot.""" y1 = d.loc["2025-01-05":"2026-01-04"]; a = pub["anniversary"]; s = pub["september"]; sb = pub["streetsblog"] w1 = d.loc["2025-01-06":"2025-01-10"]; w2 = d.loc["2025-01-13":"2025-01-17"] ja = d.loc["2025-01-05":"2025-08-31"] cbd = y1["cbd"].mean() return {"week1_cbd_file": _r(w1["cbd"].mean(), 0), "week1_cbd_published": sb["cbd_week1"], "week1_diff_pct": _r(pct(sb["cbd_week1"], w1["cbd"].mean()), 2), "week2_cbd_file": _r(w2["cbd"].mean(), 0), "week2_cbd_published": sb["cbd_week2"], "week2_diff_pct": _r(pct(sb["cbd_week2"], w2["cbd"].mean()), 2), "fewer_total_per_day": _r(a["fewer_total"] / len(y1), 0), # The 73,000 a day comes from a baseline for the zone and the two highways together, # so it is set against that count only. (v0.9.2 also set it against the tolled zone # alone, giving 12.9%; there is no zone-only baseline behind the 73,000, so that # pairing has no reading and was dropped after an outside review.) "implied_baseline_cbd": _r(cbd + a["fewer_per_day"], 0), "implied_pct_cbd": _r(100 * a["fewer_per_day"] / (cbd + a["fewer_per_day"])), "sep_days": int(len(ja)), "sep_total_per_day": _r(s["fewer_total_to_date"] / len(ja), 0), "sep_cbd_mean": _r(ja["cbd"].mean(), 0), **evaluation_check(d, pub["evaluation"], cbd + a["fewer_per_day"])} def evaluation_check(d: pd.DataFrame, ev: dict, implied_baseline: float) -> dict: """The MTA's First Evaluation Report against the open file: its monthly average entries (the zone and the two highways together, Table 2-1) recomputed from the gantry file, and its twelve monthly baselines (Figure 2-2) against the baseline the anniversary figures imply.""" months = {} for m, (entries, base) in ev["table"].items(): g = d[(d.index.year == 2025) & (d.index.month == m) & (d.index >= TOLL_START)] months[MON3[m - 1]] = {"mta_entries": entries, "file_entries": _r(g["cbd"].mean(), 0), "diff_pct": _r(pct(entries, g["cbd"].mean()), 2), "mta_baseline": base} base_mean = 1000 * float(np.mean(ev["baseline_k"])) return {"evaluation_months": months, "evaluation_max_abs_diff_pct": _r(max(abs(v["diff_pct"]) for v in months.values()), 2), "evaluation_baseline_mean": _r(base_mean, 0), "evaluation_baseline_min_k": min(ev["baseline_k"]), "evaluation_baseline_max_k": max(ev["baseline_k"]), "implied_vs_evaluation_baseline_pct": _r(pct(base_mean, implied_baseline), 2)} def compare_years(e: pd.DataFrame, d: pd.DataFrame, rng: np.random.Generator) -> dict: """Is traffic creeping back? 2026 against the same days of 2025, clean days only.""" UNPAIRED_DAYS[0] = 0 last = d.index.max() ok = same_window(d.index, last) & d["clean"].to_numpy() dd = d[ok] out = {"last_date": f"{last:%Y-%m-%d}", "last_label": f"{last.day} {MONTHS[last.month - 1]} {last.year}"} monthly = [] for m in range(1, last.month + 1): for wd in (True, False): g = dd[(dd["month"] == m) & (dd["weekday"] == wd)] a, b = g.loc[g["year"] == 2025, "crz_entries"], g.loc[g["year"] == 2026, "crz_entries"] if len(a) < 3 or len(b) < 3: continue monthly.append({"month": m, "daytype": "weekday" if wd else "weekend", **_chg(*boot_ratio(a, b, rng), a, b)}) out["monthly"] = monthly spring = dd[dd["month"] >= 3] # past the toll's first weeks and the 2026 snowstorms for wd, key in ((True, "weekday"), (False, "weekend")): g = spring[spring["weekday"] == wd] a, b = g.loc[g["year"] == 2025, "crz_entries"], g.loc[g["year"] == 2026, "crz_entries"] out[f"pooled_{key}"] = _chg(*boot_ratio(a, b, rng), a, b) wk = [r for r in monthly if r["daytype"] == "weekday"] out["weekday_months_all"] = len(wk) out["weekday_months_below"] = sum(r["change_pct"] < 0 for r in wk) out["weekday_months_interval_below_zero"] = sum(r["hi_pct"] < 0 for r in wk) days = spring.index[spring["weekday"]] e2 = e[e["d"].isin(days)] groups = {} for col, opts in (("cls", CLASSES[1:]), ("detection_region", REGIONS[1:]), ("time_period", PERIODS[1:])): for k, _ in opts: s = e2[e2[col] == k].groupby("d")["crz_entries"].sum(); s.index = pd.DatetimeIndex(s.index) a, b = s[s.index.year == 2025], s[s.index.year == 2026] groups[f"{col}:{k}"] = _chg(*boot_ratio(a, b, rng), a, b) ex = e2.groupby("d")["excluded_roadway_entries"].sum(); ex.index = pd.DatetimeIndex(ex.index) a, b = ex[ex.index.year == 2025], ex[ex.index.year == 2026] groups["excluded:highways"] = _chg(*boot_ratio(a, b, rng), a, b) out["groups"] = groups allw = d[same_window(d.index, last) & (d["month"] >= 3).to_numpy() & d["weekday"].to_numpy()] sw = spring[spring["weekday"]] out["sensitivity"] = { "all_weekdays_mean_pct": _r(pct(allw.loc[allw.year == 2025, "crz_entries"].mean(), allw.loc[allw.year == 2026, "crz_entries"].mean())), "clean_weekday_median_pct": _r(pct(sw.loc[sw.year == 2025, "crz_entries"].median(), sw.loc[sw.year == 2026, "crz_entries"].median())), "cbd_incl_highways_pct": _r(pct(sw.loc[sw.year == 2025, "cbd"].mean(), sw.loc[sw.year == 2026, "cbd"].mean())), "from_january_weekday_pct": _r(pct(dd.loc[(dd.year == 2025) & dd.weekday, "crz_entries"].mean(), dd.loc[(dd.year == 2026) & dd.weekday, "crz_entries"].mean())), } out["dropped_days"] = {"holiday": int(d.loc[same_window(d.index, last), "holiday"].sum()), "disrupted": [f"{x:%Y-%m-%d}" for x in d.index[same_window(d.index, last) & d["disrupted"].to_numpy()]]} # The traffic screen uses the counts themselves to drop a day, so a real change in behaviour # could in principle be dropped as a disruption (an outside review of v0.9.2). From March, # the window of every pooled comparison, it drops no day the fixed holiday calendar has not # already dropped: those comparisons rest on the calendar alone, and the screen matters only # for the January and February snowstorm days of the monthly figures. Listed so the prose # can say so and the test can hold it to that. spring_window = same_window(d.index, last) & (d["month"] >= 3).to_numpy() out["screened_from_march"] = [f"{x:%Y-%m-%d}" for x in d.index[spring_window & d["disrupted"].to_numpy() & ~d["holiday"].to_numpy()]] out["unpaired_days"] = UNPAIRED_DAYS[0] return out def monthly_cube(e: pd.DataFrame, d: pd.DataFrame) -> dict: """Average daily entries by calendar month for 2025 and 2026, for every combination of vehicle class, crossing, toll period and day type (clean days). The template chart's data.""" last = d.index.max() # The same window as the year-on-year comparison where the two years overlap: January from # the 5th in both years, and the latest month only to the latest day in both. 2025's months # after that have no 2026 counterpart yet and are kept whole. (v0.9.1 took every clean day of # each month, so January 2026 included the 2nd to 4th and September 2025 ran to the 30th; an # outside review caught it.) later_2025 = ((d.index.year == 2025) & (d.index.month > last.month)) keep = d["clean"].to_numpy() & (same_window(d.index, last) | later_2025) e = e[e["d"].isin(d.index[keep])].copy() e["year"] = e["d"].dt.year; e["month"] = e["d"].dt.month cells = {} for ck, _ in CLASSES: mc = np.ones(len(e), bool) if ck == "all" else (e["cls"] == ck).to_numpy() for rk, _ in REGIONS: mr = np.ones(len(e), bool) if rk == "all" else (e["detection_region"] == rk).to_numpy() for pk, _ in PERIODS: mp = np.ones(len(e), bool) if pk == "all" else (e["time_period"] == pk).to_numpy() per_day = e[mc & mr & mp].groupby(["d", "year", "month"])["crz_entries"].sum().reset_index() per_day["daytype"] = np.where(per_day["d"].dt.dayofweek < 5, "weekday", "weekend") means = per_day.groupby(["daytype", "year", "month"])["crz_entries"].mean() for dk, _ in DAYTYPES: cells[f"{ck}|{rk}|{pk}|{dk}"] = {str(y): {MON3[m - 1]: (_r(means[(dk, y, m)], 0) if (dk, y, m) in means.index else None) for m in range(1, 13)} for y in (2025, 2026)} return {"cells": cells, "last_date": f"{last:%Y-%m-%d}"} # ==================================================================================== # B. Thirty years of the cordon count (NYMTC Hub Bound) # ==================================================================================== def _lab(v) -> str: return re.sub(r"\s+", " ", str(v)).strip().upper() if v is not None else "" def _num(v): if isinstance(v, (int, float)) and not isinstance(v, bool): return float(v) return None def _year_cols(rows: list[tuple]) -> dict: for r in rows[:12]: cols = {} for j, v in enumerate(r): try: y = int(str(v).strip()) except (TypeError, ValueError): continue if 1900 <= y <= 2030: cols[y] = j if len(cols) >= 5: return cols raise AssertionError("no year header found") def hub_bound() -> dict: import openpyxl wb = openpyxl.load_workbook(HB / "DM_TDS_Hub_Bound_Travel_AppendixI_2024.xlsx", read_only=True, data_only=True) veh, persons, printed_total = {}, {}, {} sector_key = {"NORTH OF 60TH STREET": "60th", "BROOKLYN": "brooklyn", "QUEENS": "queens", "NEW JERSEY": "nj", "STATEN ISLAND": "si"} mode_key = {"AUTO, TAXI, VAN & TRUCK": "auto", "BUS": "bus", "SUBWAY": "subway", "RAILROAD": "rail", "PASSENGER FERRY": "ferry", "TRAMWAY": "tram", "TROLLEY/TRAM": "tram", "BICYCLE": "bicycle", "TOTAL": "total"} for k in range(1, 6): rows = list(wb[f"Table1A ({k})"].iter_rows(values_only=True)) ycol = _year_cols(rows); section = None for r in rows: lab = _lab(r[0]) if lab.startswith("PERSONS BY MODE"): section = "mode"; continue if lab.startswith("PERSONS BY SECTOR"): section = "psector"; continue if lab.startswith("MOTOR VEHICLES BY SECTOR"): section = "veh"; continue if section == "veh" and (lab in sector_key or lab == "TOTAL"): for y, c in ycol.items(): v = _num(r[c]) if c < len(r) else None if v is None: continue if lab == "TOTAL": printed_total[y] = v else: veh.setdefault(y, {})[sector_key[lab]] = v if section == "mode" and lab in mode_key: for y, c in ycol.items(): v = _num(r[c]) if c < len(r) else None if v is not None: persons.setdefault(y, {})[mode_key[lab]] = v # The one-year tables for the corrections. wb2 = openpyxl.load_workbook(HB / "DM_TDS_Hub_Bound_Travel_AppendixII_2024.xlsx", read_only=True, data_only=True) t16 = list(wb2["Table16"].iter_rows(values_only=True)) t16_sector, sec, t16_fac = {}, None, {} names16 = {"60TH STREET SECTOR": "60th", "BROOKLYN SECTOR": "brooklyn", "QUEENS SECTOR": "queens", "NEW JERSEY SECTOR": "nj"} for r in t16: lab = _lab(r[1]) if len(r) > 1 else "" if lab in names16: sec = names16[lab]; continue if lab == "SECTOR TOTAL" and sec: t16_sector[sec] = float(r[4]); sec = None; continue if lab == "TOTAL, ALL SECTORS": t16_total, t16_bus = float(r[4]), float(r[7]); continue if sec and _num(r[4]) is not None: t16_fac[lab] = float(r[4]) t23 = list(wb2["Table23A"].iter_rows(values_only=True)) tot23 = next(r for r in t23 if len(r) > 1 and _lab(r[1]) == "TOTAL") inbound = {2022: float(tot23[2]), 2023: float(tot23[5]), 2024: float(tot23[8])} t11 = list(wb["Table11"].iter_rows(values_only=True)); ycol11 = _year_cols(t11) bb = next(r for r in t11 if len(r) > 1 and _lab(r[1]) == "BROOKLYN BRIDGE") bb_twoway = {y: float(bb[c]) for y, c in ycol11.items() if _num(bb[c]) is not None} # Corrections to Table 1A, each checked against the report's own one-year tables. s23 = sum(veh[2023][k] for k in ("60th", "brooklyn", "queens", "nj")) assert abs(printed_total[2023] - 699.294) < 1e-6 and abs(s23 - 629.889) < 1e-6 and abs(s23 * 1000 - inbound[2023]) < 1, "Table 1A 2023 no longer as documented" assert abs(printed_total[2024] - 785.354) < 1e-6 and abs(t16_total - 618_970) < 1 and abs(t16_total - inbound[2024]) < 1, "Table 1A 2024 no longer as documented" assert abs(printed_total[2022] * 1000 - inbound[2022]) < 1, "Table 1A 2022 disagrees with Table 23A" veh[2024] = {k: t16_sector[k] / 1000 for k in ("60th", "brooklyn", "queens", "nj")} total = {} for y, s in veh.items(): total[y] = sum(s.get(k, 0.0) for k in ("60th", "brooklyn", "queens", "nj", "si")) if y not in (2023, 2024): assert abs(total[y] - printed_total[y]) < 2.5, f"Hub Bound {y}: sectors {total[y]} vs printed total {printed_total[y]}" years = sorted(y for y in total if y >= 1990) same = {y: sum(veh[y][k] for k in ("60th", "brooklyn", "nj")) for y in years} # The highest count with its data intact: NYMTC lost most of its 1999 files on 11 September # 2001 ("data for 1999 is largely absent", report p. 10), so 1999 is drawn but not used as the peak. peak = max((y for y in years if y <= 2019 and y != 1999), key=lambda y: total[y]) return { "years": years, "total": {y: _r(total[y] * 1000, 0) for y in years}, "sectors": {y: {k: _r(veh[y].get(k, 0.0) * 1000, 0) for k in ("60th", "brooklyn", "queens", "nj")} for y in years}, "same_basis": {y: _r(same[y] * 1000, 0) for y in years}, "persons": {y: {k: _r(v * 1000, 0) for k, v in persons[y].items()} for y in sorted(persons) if y >= 1990}, "peak_year": int(peak), "peak": _r(total[peak] * 1000, 0), # 2017 is the last count before the pandemic with no problem flagged by NYMTC (2018: a low # Queensboro Bridge count, its note (3); 2019: the Brooklyn Bridge tube count, the addendum). "change_peak_to_2017_pct": _r(pct(total[peak], total[2017])), "change_peak_to_2018_pct": _r(pct(total[peak], total[2018])), "change_2018_to_2024_pct": _r(pct(total[2018], total[2024])), "change_peak_to_2024_pct": _r(pct(total[peak], total[2024])), "same_basis_peak_to_2018_pct": _r(pct(same[peak], same[2018])), "same_basis_2018_to_2024_pct": _r(pct(same[2018], same[2024])), "same_basis_peak_to_2024_pct": _r(pct(same[peak], same[2024])), "table16_2024": {"sectors": {k: _r(v, 0) for k, v in t16_sector.items()}, "total": _r(t16_total, 0), "buses": _r(t16_bus, 0), "fdr": _r(t16_fac.get("FDR DRIVE"), 0), "wsh": _r(t16_fac.get("WEST SIDE HIGHWAY"), 0), "queensboro": _r(t16_fac.get("ED KOCH QUEENSBORO BRIDGE"), 0)}, "table23a_inbound": {y: _r(v, 0) for y, v in inbound.items()}, "table1a_printed_wrong": {"2023": printed_total[2023] * 1000, "2024": printed_total[2024] * 1000}, "brooklyn_bridge_twoway": {y: _r(v, 0) for y, v in bb_twoway.items()}, "persons_2019_2020": {"auto_pct": _r(pct(persons[2019]["auto"], persons[2020]["auto"])), "subway_pct": _r(pct(persons[2019]["subway"], persons[2020]["subway"]))}, "persons_2019_2023": {"auto_pct": _r(pct(persons[2019]["auto"], persons[2023]["auto"])), "subway_pct": _r(pct(persons[2019]["subway"], persons[2023]["subway"]))}, "y1999_partial": _r(total[1999] * 1000, 0), "y2000": _r(total[2000] * 1000, 0), "queens_2020_2021": {"2020": _r(veh[2020]["queens"] * 1000, 0), "2021": _r(veh[2021]["queens"] * 1000, 0)}, } def hub_vs_gantry(e: pd.DataFrame, hub: dict) -> dict: """The October cordon count set against October gantry days: the only way to set the toll's first year against the years before it on counts of the same crossings. Two counting programmes, so an estimate, not a before-and-after count. Buses are left out of both (the cordon count reports them apart); the two highways are kept in both.""" nb = e[e["cls"] != "buses"].copy() nb["tot"] = nb["crz_entries"] + nb["excluded_roadway_entries"] nb["sector"] = nb["detection_region"].map(HUB_SECTOR_OF) out = {} def window(y, m, dows): w = nb[(nb["d"].dt.year == y) & (nb["d"].dt.month == m) & nb["d"].dt.dayofweek.isin(dows) & ~nb["d"].isin(list(holidays()))] n = w["d"].nunique() return w.groupby("sector")["tot"].sum() / n, w.groupby("d")["tot"].sum().mean(), n sec, tot, n = window(2025, 10, [1, 2, 3]) # Tuesday to Thursday, the cordon count's kind of day h = hub["table16_2024"] out["october_2025_midweek"] = {"days": int(n), "total": _r(tot, 0), "sectors": {k: _r(sec[k], 0) for k, _ in SECTORS}} out["hub_2024"] = {"total": h["total"], "sectors": h["sectors"]} out["change_pct"] = _r(pct(h["total"], tot)) out["sector_change_pct"] = {k: _r(pct(h["sectors"][k], sec[k])) for k, _ in SECTORS} sens = {} for m in (9, 11): _, t, nn = window(2025, m, [1, 2, 3]); sens[f"month_{m}_midweek_pct"] = _r(pct(h["total"], t)) _, t, _ = window(2025, 10, [0, 1, 2, 3, 4]); sens["october_all_weekdays_pct"] = _r(pct(h["total"], t)) base3 = np.mean([hub["table23a_inbound"][y] for y in (2022, 2023, 2024)]) sens["vs_2022_2024_mean_pct"] = _r(pct(base3, tot)); sens["hub_2022_2024_mean"] = _r(base3, 0) out["sensitivity"] = sens return out # ==================================================================================== # C. The two MTA tunnels into the zone, 2005-2026 # ==================================================================================== def tunnels() -> dict: a = pd.read_csv(T / "mta_bt_daily_2005_2024.csv"); a["d"] = pd.to_datetime(a["collection_date"]) meta = json.loads((REFERENCE / "mta_bt_daily_dictionary.json").read_text(encoding="utf-8")) desc = next(c["description"] for c in meta["columns"] if c["fieldName"] == "plaza") assert "27 - Queens Midtown Tunnel" in desc and "28 - Hugh L. Carey Tunnel" in desc, "plaza codes changed" a = a[a["plaza"].isin([27, 28])] ann = a.groupby([a["d"].dt.year, "plaza"])["total"].sum().unstack() b = pd.read_csv(T / "mta_bt_crossings_daily_2019.csv"); b["d"] = pd.to_datetime(b["date"]) b = b[b["facility"].isin(["Queens Midtown Tunnel", "Hugh L. Carey Tunnel"])] last = b["d"].max() annb = b.groupby([b["d"].dt.year, "facility"])["traffic"].sum().unstack() to_mh = b[b["direction"].isin(["Westbound to Manhattan", "Northbound to Manhattan"])] assert set(to_mh["facility"]) == {"Queens Midtown Tunnel", "Hugh L. Carey Tunnel"}, "direction labels changed" def ytd(frame, y): f = frame[(frame["d"].dt.year == y) & ((frame["d"].dt.month < last.month) | ((frame["d"].dt.month == last.month) & (frame["d"].dt.day <= last.day)))] return float(f["traffic"].sum()) annual = {int(y): _r(ann.loc[y].sum(), 0) for y in ann.index} annual[2025] = _r(annb.loc[2025].sum(), 0) overlap = {int(y): _r(pct(ann.loc[y].sum(), annb.loc[y].sum()), 2) for y in range(2019, 2025)} return {"annual_both_directions": annual, "overlap_ebfx_vs_reconciled_pct": overlap, "last_date": f"{last:%Y-%m-%d}", "window": f"1 January to {last.day} {MONTHS[last.month - 1]}", "ytd_both": {y: _r(ytd(b, y), 0) for y in (2019, 2024, 2025, 2026)}, "ytd_to_manhattan": {y: _r(ytd(to_mh, y), 0) for y in (2019, 2024, 2025, 2026)}, "ytd_to_manhattan_2025_vs_2024_pct": _r(pct(ytd(to_mh, 2024), ytd(to_mh, 2025))), "ytd_to_manhattan_2026_vs_2025_pct": _r(pct(ytd(to_mh, 2025), ytd(to_mh, 2026))), "peak_year": int(max(annual, key=lambda y: annual[y] if y <= 2024 else 0)), "change_2019_2020_pct": _r(pct(annual[2019], annual[2020]))} # ==================================================================================== # D. Ride-hail: TLC monthly trips, 2010-2026 # ==================================================================================== TLC_CLASSES = [("yellow", "Yellow"), ("green", "Green"), ("hv", "FHV - High Volume")] def tlc() -> dict: t = pd.read_csv(T / "tlc_data_reports_monthly.csv", dtype=str) num = lambda s: pd.to_numeric(s.str.replace(",", "", regex=False).replace("-", None), errors="coerce") t["trips"] = num(t["Trips Per Day"]); t["veh"] = num(t["Unique Vehicles"]); t["ym"] = t["Month/Year"] out = {"monthly": {}, "annual": {}, "june": {}, "vehicles_june": {}} for k, lab in TLC_CLASSES: s = t[t["License Class"] == lab].set_index("ym").sort_index() assert len(s) > 100, f"TLC class {lab} missing" out["monthly"][k] = {ym: _r(v, 0) for ym, v in s["trips"].items()} yrs = s.groupby(s.index.str[:4])["trips"].agg(["mean", "size"]) out["annual"][k] = {y: {"mean": _r(r["mean"], 0), "months": int(r["size"])} for y, r in yrs.iterrows()} out["june"][k] = {ym[:4]: _r(v, 0) for ym, v in s["trips"].items() if ym.endswith("-06")} out["vehicles_june"][k] = {ym[:4]: _r(v, 0) for ym, v in s["veh"].items() if ym.endswith("-06")} latest = max(max(out["monthly"][k]) for k, _ in TLC_CLASSES) out["latest_month"] = latest y, h = out["june"]["yellow"], out["june"]["hv"] out["headline"] = {"yellow_june_2010": y["2010"], "yellow_june_2019": y["2019"], "yellow_june_2026": y["2026"], "hv_june_2015": h["2015"], "hv_june_2019": h["2019"], "hv_june_2026": h["2026"], "hv_vehicles_june_2015": out["vehicles_june"]["hv"]["2015"], "hv_vehicles_june_2019": out["vehicles_june"]["hv"]["2019"], "hv_vehicles_june_2026": out["vehicles_june"]["hv"]["2026"], "street_hail_june_2013": _r(y["2013"] + out["june"]["green"].get("2013", 0), 0)} return out # ==================================================================================== # E. Speed in the core # ==================================================================================== DOT_CY = {2010: 9.1, 2011: 8.8, 2012: 9.1, 2013: 8.5, 2014: 8.0, 2015: 7.4, 2016: 7.2, 2017: 7.1, 2018: 7.0} # Mobility Report 2019, pp. 7 and 10 def speeds() -> dict: # The city's indicator: fiscal 2009-2012 from the older file, 2016-2026 from the current one. old = pd.read_csv(T / "nyc_mmr_fy03_fy12_dot.csv", dtype=str) meta = json.loads((REFERENCE / "nyc_mmr_fy03_fy12_dictionary.json").read_text(encoding="utf-8")) fy_of = {c["fieldName"]: int(c["name"]) for c in meta["columns"] if re.fullmatch(r"_\d+", c["fieldName"])} row = old[old["indicator_name"].str.contains("travel speed - Manhattan Central Business District", regex=False)].iloc[0] fy = {fy_of[c]: float(row[c]) for c in fy_of if isinstance(row[c], str) and row[c].strip()} cur = pd.read_csv(T / "nyc_mmr_indicators_cbd_speed.csv") definition = cur["description"].dropna().iloc[0] assert "yellow taxis traveling with passengers" in definition and "8AM-6PM" in definition, "indicator definition changed" cur = cur[pd.to_datetime(cur["valuedate"]).dt.month == 6] for _, r in cur.iterrows(): v = r["acceptedvalueytd"] if pd.notna(r["acceptedvalueytd"]) and str(r["acceptedvalueytd"]) not in ("NA", "") else r["acceptedvalue"] if pd.notna(v) and str(v) not in ("NA", ""): fy[int(r["fiscalyear"])] = float(v) fy = dict(sorted(fy.items())) assert fy[2026] == 6.5 and fy[2025] == 6.9 and fy[2009] == 9.1, "MMR speed values moved" # The calendar-year speeds as the Mobility Report prints them: the chart on PDF page 10 labels every # year's CBD average ("9.1 mph" ... "7.0 mph"), and the indicator table on page 7 carries 2010-2017. chart = page_texts(REPORTS / "nyc_dot_mobility_report_2019.pdf")[9] table = page_texts(REPORTS / "nyc_dot_mobility_report_2019.pdf")[6] for y, v in DOT_CY.items(): assert f"{v:.1f} mph" in chart, f"Mobility Report chart lacks {v:.1f} mph ({y})" if y <= 2017: assert f"{v:.1f}" in table, f"Mobility Report table lacks {v:.1f} ({y})" # MTA taxi and ride-hail speeds by zone, October 2019 to August 2025. s = pd.read_csv(T / "mta_cbd_taxi_fhv_speeds.csv"); s["Month"] = pd.to_datetime(s["Month"]) p = s.pivot(index="Month", columns="Zone", values="Zonal Speed").sort_index(); last = p.index.max() m8 = p[p.index.month <= last.month]; ja = m8.groupby(m8.index.year).mean() mta = {"latest": f"{last:%Y-%m}", "monthly": {z: {f"{i:%Y-%m}": _r(v, 2) for i, v in p[z].items() if pd.notna(v)} for z in p.columns}, "janaug": {int(y): {z: _r(ja.at[y, z], 2) for z in p.columns} for y in ja.index if y >= 2020}, "cbd_2025_vs_2024_pct": _r(pct(ja.at[2024, "CBD"], ja.at[2025, "CBD"])), "greater_2025_vs_2024_pct": _r(pct(ja.at[2024, "Greater NYC"], ja.at[2025, "Greater NYC"])), "cbd_2024_vs_2023_pct": _r(pct(ja.at[2023, "CBD"], ja.at[2024, "CBD"])), "greater_2024_vs_2023_pct": _r(pct(ja.at[2023, "Greater NYC"], ja.at[2024, "Greater NYC"]))} # MTA bus speeds on segments inside and outside the zone, weekday local/limited/SBS. b = pd.read_csv(T / "mta_cbd_bus_speeds_monthly.csv"); b["month"] = pd.to_datetime(b["month"]) b = b[b["route_type"].isin(["Local", "Limited", "SBS"]) & (b["day_type"] == "Weekday")] blast = b["month"].max(); b = b[b["month"].dt.month <= blast.month] g = b.groupby([b["month"].dt.year, "cbd_relation", "time_period"])[["sum_mileage", "sum_time"]].sum() mph = (g["sum_mileage"] / g["sum_time"]).unstack(["cbd_relation", "time_period"]) bus = {"latest": f"{blast:%Y-%m}", "mph": {int(y): {f"{a}|{c}": _r(mph.at[y, (a, c)], 2) for a, c in mph.columns} for y in mph.index}, "cbd_peak_2025_vs_2024_pct": _r(pct(mph.at[2024, ("CBD", "Peak")], mph.at[2025, ("CBD", "Peak")])), "cbd_peak_2026_vs_2025_pct": _r(pct(mph.at[2025, ("CBD", "Peak")], mph.at[2026, ("CBD", "Peak")])), "cbd_peak_2026_vs_2024_pct": _r(pct(mph.at[2024, ("CBD", "Peak")], mph.at[2026, ("CBD", "Peak")])), "noncbd_peak_2025_vs_2024_pct": _r(pct(mph.at[2024, ("Non-CBD", "Peak")], mph.at[2025, ("Non-CBD", "Peak")])), "noncbd_peak_2026_vs_2025_pct": _r(pct(mph.at[2025, ("Non-CBD", "Peak")], mph.at[2026, ("Non-CBD", "Peak")]))} # Vehicle-miles in the zone, two probe models. v = pd.read_csv(T / "mta_cbd_vmt.csv"); v.columns = ["date", "area", "year", "month", "agps", "cvd", "index"] crz = v[v["area"] == "CRZ"].set_index(["year", "month"]).sort_index() def total_chg(col, a, bb, months=None): """Change in total vehicle-miles over the months both years have: each month's daily figure times its days, summed, so a long month weighs more than a short one. (v0.9.1 averaged the monthly percentage changes; an outside review asked for the total.)""" ta = tb = 0.0; used = 0 for m in range(1, 13): if months and m not in months: continue if (a, m) in crz.index and (bb, m) in crz.index and pd.notna(crz.at[(a, m), col]) and pd.notna(crz.at[(bb, m), col]): ta += crz.at[(a, m), col] * pd.Period(f"{a}-{m:02d}").days_in_month tb += crz.at[(bb, m), col] * pd.Period(f"{bb}-{m:02d}").days_in_month used += 1 return (_r(pct(ta, tb)) if used else None), used vmt = {} for col in ("agps", "cvd"): vmt[f"{col}_2025_vs_2024_pct"], vmt[f"{col}_2025_months"] = total_chg(col, 2024, 2025) vmt[f"{col}_2026_vs_2025_pct"], vmt[f"{col}_2026_months"] = total_chg(col, 2025, 2026, range(3, 13)) fy_rows = {y: v for y, v in fy.items()} return {"mmr_fy": fy_rows, "mmr_definition": definition, "dot_cy": DOT_CY, "mta": mta, "bus": bus, "vmt": vmt, "mmr_fy2026_vs_fy2025_pct": _r(pct(fy[2025], fy[2026])), "mmr_fy2026_vs_fy2019_pct": _r(pct(fy[2019], fy[2026])), "mmr_fy2012_vs_fy2019_pct": _r(pct(fy[2012], fy[2019])), "mmr_min_before_2026": _r(min(v for y, v in fy.items() if y < 2026)), "dot_2010_2018_pct": _r(pct(DOT_CY[2010], DOT_CY[2018]))} # ==================================================================================== # F. Street space inside the zone: bike lanes, bus lanes, the curb lane # ==================================================================================== def zone_polygons() -> list[np.ndarray]: g = pd.read_csv(T / "mta_cbd_geofence.csv") polys = [] for w in g["polygon"]: for ring in re.findall(r"\(\(([^()]+)\)", w) + re.findall(r"\(([^()]+)\)", w): pts = np.array([[float(a) for a in xy.split()] for xy in ring.split(",")]) if len(pts) >= 4: polys.append(pts) # Deduplicate rings found by both patterns. uniq, seen = [], set() for p in polys: k = (len(p), round(p[0, 0], 7), round(p[0, 1], 7)) if k not in seen: seen.add(k); uniq.append(p) return uniq def inside(px: np.ndarray, py: np.ndarray, polys: list[np.ndarray]) -> np.ndarray: """Even-odd ray casting; a point is in the zone if it is inside any polygon.""" res_ = np.zeros(px.shape, bool) for ring in polys: c = np.zeros(px.shape, bool) x1, y1, x2, y2 = ring[:-1, 0], ring[:-1, 1], ring[1:, 0], ring[1:, 1] for a, b, cc, dd in zip(x1, y1, x2, y2): cond = ((b > py) != (dd > py)) xint = a + (py - b) * (cc - a) / np.where(dd - b == 0, 1e-12, dd - b) c ^= cond & (px < xint) res_ |= c return res_ def _m_per_deg(lat: float) -> tuple[float, float]: p = np.radians(lat) return (111132.954 - 559.822 * np.cos(2 * p) + 1.175 * np.cos(4 * p), 111412.84 * np.cos(p) - 93.5 * np.cos(3 * p) + 0.118 * np.cos(5 * p)) def in_zone_miles(wkts: pd.Series, polys: list[np.ndarray], step_m: float = 10.0) -> np.ndarray: """Length inside the zone of each line geometry, in miles: every piece is cut into chunks of at most 10 m and a chunk counts if its midpoint is inside the zone.""" seg_id, mids_x, mids_y, lens = [], [], [], [] for i, w in enumerate(wkts): if not isinstance(w, str): continue for part in re.findall(r"\(([^()]+)\)", w): pts = np.array([[float(a) for a in xy.split()] for xy in part.split(",")]) if len(pts) < 2: continue my, mx = _m_per_deg(float(pts[:, 1].mean())) for (x0, y0), (x1, y1) in zip(pts[:-1], pts[1:]): L = float(np.hypot((x1 - x0) * mx, (y1 - y0) * my)) if L == 0: continue n = int(np.ceil(L / step_m)); t = (np.arange(n) + 0.5) / n mids_x.append(x0 + (x1 - x0) * t); mids_y.append(y0 + (y1 - y0) * t); lens.append(np.full(n, L / n)); seg_id.append(np.full(n, i)) out = np.zeros(len(wkts)) if not lens: return out mx_, my_, L_, sid = np.concatenate(mids_x), np.concatenate(mids_y), np.concatenate(lens), np.concatenate(seg_id) box = (mx_ > -74.03) & (mx_ < -73.95) & (my_ > 40.69) & (my_ < 40.785) ok = np.zeros(len(mx_), bool); ok[box] = inside(mx_[box], my_[box], polys) np.add.at(out, sid[ok], L_[ok]) return out / 1609.344 def streets() -> dict: polys = zone_polygons() assert len(polys) >= 1 import openpyxl # the dictionary's warning about retired records, quoted in the questions wb = openpyxl.load_workbook(REFERENCE / "nyc_dot_bike_routes_dictionary.xlsx", read_only=True) text = " ".join(str(x) for ws in wb for row in ws.iter_rows(values_only=True) for x in row if x) assert "retired bicycle facilities do not inherently indicate the re" in text, "bike dictionary wording changed" # Bike lanes: stock on 31 December of each year, one record per street segment on each date. b = pd.read_csv(T / "nyc_dot_bike_routes.csv", dtype=str) b = b[b["boro"] == "1"].copy() # Manhattan; the zone lies wholly inside it b["inst"] = pd.to_datetime(b["instdate"], format="%m/%d/%Y", errors="coerce") b["ret"] = pd.to_datetime(b["ret_date"], format="%m/%d/%Y", errors="coerce") b["mi"] = in_zone_miles(b["the_geom"], polys) b = b[b["mi"] > 0] kind = np.select([(b["facilitycl"] == "I") & (b["onoffst"] == "ON"), b["facilitycl"] == "I", b["facilitycl"] == "II", b["facilitycl"] == "III"], ["protected", "offstreet", "painted", "shared"], "other") b["kind"] = kind bike = {} for y in list(range(2000, 2026)) + ["current"]: D = pd.Timestamp("2099-12-31") if y == "current" else pd.Timestamp(f"{y}-12-31") live = b[(b["inst"] <= D) & (b["ret"].isna() | (b["ret"] > D))] live = live.sort_values(["inst", "ret"], na_position="last").groupby("segmentid").tail(1) bike[str(y)] = {k: _r(live.loc[live["kind"] == k, "mi"].sum(), 1) for k in ("protected", "painted", "shared", "offstreet")} D = pd.Timestamp("2007-12-31") # the same stock without one-record-per-segment raw07 = b[(b["inst"] <= D) & (b["ret"].isna() | (b["ret"] > D))] protected_2007_no_dedup = _r(raw07.loc[raw07["kind"] == "protected", "mi"].sum(), 1) # Bus lanes in effect now, by the year each began (lanes since removed are not in the file). L = pd.read_csv(T / "nyc_dot_bus_lanes.csv", dtype=str) L = L[L["Boro"].isin(["MAN", "BX, MN"])].copy() n_all = len(L) # The file repeats some lanes exactly (same geometry, direction, hours, type and year): count each once. L = L.drop_duplicates(subset=["the_geom", "Direction", "Hours", "Days", "Lane_Type", "Lane_Type1", "Year1"]) dup_removed = n_all - len(L) L["mi"] = in_zone_miles(L["the_geom"], polys) L = L[L["mi"] > 0]; L["y1"] = pd.to_numeric(L["Year1"], errors="coerce") bus = {str(y): _r(L.loc[(L["y1"] > 0) & (L["y1"] <= y), "mi"].sum(), 1) for y in range(1980, 2026)} busway = _r(L.loc[L["Lane_Type"] == "Busway", "mi"].sum(), 2) bus_undated = _r(L.loc[~(L["y1"] > 0), "mi"].sum(), 2) # The curb lane: pandemic outdoor dining (applications) and the permanent programme (licences). o = pd.read_csv(T / "nyc_open_restaurants_applications.csv", dtype=str) o = o[o["Approved for Roadway Seating"].str.lower() == "yes"].copy() o["lat"] = pd.to_numeric(o["Latitude"], errors="coerce"); o["lon"] = pd.to_numeric(o["Longitude"], errors="coerce") o = o.dropna(subset=["lat", "lon"]) o["zone"] = inside(o["lon"].to_numpy(), o["lat"].to_numpy(), polys) o["year"] = pd.to_datetime(o["Time of Submission"], format="%m/%d/%Y %I:%M:%S %p", errors="coerce").dt.year oz = o[o["zone"]] places = oz.drop_duplicates(subset=["Business Address"]) dn = pd.read_csv(T / "nyc_dining_out_licences.csv", dtype=str) dn["lat"] = pd.to_numeric(dn["Latitude"], errors="coerce"); dn["lon"] = pd.to_numeric(dn["Longitude"], errors="coerce") dn = dn.dropna(subset=["lat", "lon"]); dn["zone"] = inside(dn["lon"].to_numpy(), dn["lat"].to_numpy(), polys) dz = dn[dn["zone"]] return {"bike": bike, "bus_lane_miles": bus, "busway_miles": busway, "bus_undated_miles": bus_undated, "bus_duplicates_removed_manhattan": int(dup_removed), "bus_shared_lane_miles_2025": _r(L.loc[(L["Lane_Type1"] == "Shared Lane") & (L["y1"] > 0) & (L["y1"] <= 2025), "mi"].sum(), 1), "bike_protected_2007": bike["2007"]["protected"], "bike_protected_2007_no_dedup": protected_2007_no_dedup, "bike_protected_current": bike["current"]["protected"], "bike_protected_2019": bike["2019"]["protected"], "bike_painted_current": bike["current"]["painted"], "dining": {"roadway_applications_zone_by_year": {int(k): int(v) for k, v in oz["year"].value_counts().sort_index().items()}, "roadway_places_zone": int(len(places)), "dining_out_roadway_zone": int((dz["License Type"] == "Roadway").sum()), "dining_out_sidewalk_zone": int((dz["License Type"] == "Sidewalk").sum()), "dining_out_latest_issue": str(pd.to_datetime(dn["License Issue Date"], format="%m/%d/%Y %I:%M:%S %p", errors="coerce").max().date())}} # ==================================================================================== # G. Bar-chart cubes for the questions, the companion data, the manifest # ==================================================================================== def chart_cubes(crz: dict, hub: dict, tun: dict) -> dict: # Each bar carries its matched-week range as it is, low and high ends (change_lo, change_hi), # not a symmetric half-width: v0.9.2 drew the wider half on both sides and the hover read # "± x", which looks like a sampling margin (an outside review). g = crz["compare"]["groups"]; charts = {} series_all, cells = [], {} rng_of = lambda r: {"n": 1, "change": r["change_pct"], "change_lo": r["lo_pct"], "change_hi": r["hi_pct"]} for cut, col, opts in (("cls", "cls", CLASSES[1:]), ("region", "detection_region", REGIONS[1:]), ("period", "time_period", PERIODS[1:])): row = {} for k, lab in opts: row[f"{cut}:{k}"] = rng_of(g[f"{col}:{k}"]) series_all.append([f"{cut}:{k}", lab]) if cut == "region": row["region:highways"] = rng_of(g["excluded:highways"]) series_all.append(["region:highways", "FDR Drive and West Side Highway (not tolled)"]) cells[cut] = row last = crz["compare"]["last_label"] charts["groups"] = {"title": "Which entries fell in 2026, and which rose", "subtitle": f"Average weekday entries into the zone, 2026 against the comparable period of 2025, March to {last}", "note": "Weekdays from March, holidays left out. The tick is the 90% matched-week range: the same week of the year is drawn from both years together, and the range shows how much the change moves with the weeks that fell in the window. It is not a sampling margin, since the file counts the vehicles the tolling system recorded rather than a sample, and it is not an estimate of the toll's effect. Cars, pickups and vans are the MTA's class 1, which leaves out the taxis and for-hire cars billed through the per-trip charge; a taxi or for-hire car not billed that way is counted in its ordinary class.", "controls": [{"key": "cut", "label": "Group by", "options": [["cls", "Vehicle type"], ["region", "Crossing"], ["period", "Toll period"]]}], "metrics": [{"key": "change", "label": "Change on 2025", "format": "pct", "axis": "Change in average weekday entries, 2026 against 2025 (%)"}], "series": series_all, "series_by_cell": True, "cells": cells} rows = {} for r in crz["compare"]["monthly"]: k = f"{r['month']:02d}" rows.setdefault(r["daytype"], {})[k] = rng_of(r) months = sorted({f"{r['month']:02d}" for r in crz["compare"]["monthly"]}) charts["months"] = {"title": "Every month of 2026 so far has had fewer weekday entries than 2025", "subtitle": f"Average daily entries into the zone, 2026 against the comparable period of 2025, January to {last}", "note": "Holidays left out, and in January and February 2026 the two snowstorms; January starts on the 5th in both years, the toll's first day. The tick is the 90% matched-week range; a month holds only three to five weeks, so the ranges are rough.", "controls": [{"key": "daytype", "label": "Days", "options": [["weekday", "Weekdays"], ["weekend", "Weekends"]]}], "metrics": [{"key": "change", "label": "Change on 2025", "format": "pct", "axis": "Change in average daily entries, 2026 against 2025 (%)"}], "series": [[m, MONTHS[int(m) - 1]] for m in months], "cells": rows} yrs = sorted(int(y) for y in tun["annual_both_directions"]) charts["tunnels"] = {"title": "The two MTA tunnels into the zone, year by year", "subtitle": f"Vehicle crossings a year, both directions, Queens-Midtown and Hugh L. Carey tunnels, 2005 to 2025", "note": "Reconciled toll transactions to 2024; 2025 from the MTA's hourly crossings file, which counts every vehicle whether or not it pays and reads 0.2 to 0.6% higher where the two overlap.", "controls": [], "metrics": [{"key": "total", "label": "Crossings a year", "format": "count", "axis": "Vehicle crossings a year, both directions"}], "series": [[str(y), str(y)] for y in yrs], "cells": {"": {str(y): {"n": 1, "total": tun["annual_both_directions"][y]} for y in yrs}}} return {"charts": charts, "min_records": 1, "source": "MTA Congestion Relief Zone vehicle entries; MTA Bridges and Tunnels traffic (data.ny.gov)."} def fig_change(spec: dict, path: Path) -> None: """Static twin of a change chart, for readers without scripts: the default view as horizontal bars, each with its matched-week range drawn from its low to its high end and the change to one decimal (the shared bar-chart twin draws a symmetric tick and whole percentages, which hid the difference between, say, 3.4% and 3.0%).""" from article_migration import _mpl, _save, _style, BLUE, INK, TICK, FIG_W state = "|".join(str(c.get("default", c["options"][0][0])) for c in spec.get("controls", [])) row = spec["cells"].get(state) or spec["cells"].get("") or {} series = [s for s in spec["series"] if s[0] in row] v = np.array([row[s[0]]["change"] for s in series], float) lo = np.array([row[s[0]]["change_lo"] for s in series], float); hi = np.array([row[s[0]]["change_hi"] for s in series], float) plt = _mpl(); y = np.arange(len(series))[::-1] fig, ax = plt.subplots(figsize=(FIG_W, 1.4 + 0.34 * len(series)), dpi=100); fig.patch.set_alpha(0) ax.barh(y, v, height=0.7, color=BLUE) ax.errorbar(v, y, xerr=[v - lo, hi - v], fmt="none", ecolor=TICK, elinewidth=1.1) top = float(max(np.abs(lo).max(), np.abs(hi).max())) for yi, vi, l, h in zip(y, v, lo, hi): lab = f"{'+' if vi > 0 else '−' if vi < 0 else ''}{abs(vi):.1f}%" ax.text(h + top * 0.03 if vi >= 0 else l - top * 0.03, yi, lab, va="center", ha="left" if vi >= 0 else "right", fontsize=8, color=INK) ax.set_yticks(y); ax.set_yticklabels([s[1] for s in series], fontsize=9) ax.set_xlim(-top * 1.35, top * 1.35); ax.axvline(0, color=TICK, linewidth=0.8) ax.set_xlabel(spec["metrics"][0]["axis"], fontsize=9, color=INK) _style(ax, spec["title"], spec["subtitle"]) fig.tight_layout(); _save(fig, path); plt.close(fig) def companion(e: pd.DataFrame) -> pd.DataFrame: c = e[["toll_date", "vehicle_class", "detection_region", "time_period", "crz_entries", "excluded_roadway_entries"]].copy() c["toll_date"] = pd.to_datetime(c["toll_date"]).dt.strftime("%Y-%m-%d") return c.sort_values(["toll_date", "vehicle_class", "detection_region", "time_period"]).reset_index(drop=True) def long_view(hub: dict, sp: dict, tl: dict, st: dict, tun: dict) -> pd.DataFrame: rows = [] for y in hub["years"]: rows.append(("NYMTC Hub Bound", "vehicles entering, all sectors", y, hub["total"][y])) for k, _ in SECTORS: rows.append(("NYMTC Hub Bound", f"vehicles entering, {k}", y, hub["sectors"][y][k])) for y, p in hub["persons"].items(): for k in ("auto", "subway", "bus", "rail", "total"): if k in p: rows.append(("NYMTC Hub Bound", f"persons entering, {k}", y, p[k])) for y, v in sp["mmr_fy"].items(): rows.append(("MMR", "CBD taxi speed, fiscal year (mph)", y, v)) for y, v in sp["dot_cy"].items(): rows.append(("DOT Mobility Report 2019", "CBD taxi speed, calendar year (mph)", y, v)) for k, _ in TLC_CLASSES: for y, r in tl["annual"][k].items(): rows.append(("TLC", f"trips per day, {k} (mean of months)", int(y), r["mean"])) for y, r in st["bike"].items(): for k, v in r.items(): rows.append(("NYC DOT bike routes", f"bike lane centreline miles in the zone, {k}", y, v)) for y, v in st["bus_lane_miles"].items(): rows.append(("NYC DOT bus lanes", "bus lane miles in the zone (lanes in effect now, by year begun)", int(y), v)) for y, v in tun["annual_both_directions"].items(): rows.append(("MTA B&T", "Queens-Midtown and Hugh L. Carey crossings, both directions", int(y), v)) return pd.DataFrame(rows, columns=["source", "series", "year", "value"]) def write_build_manifest() -> None: mdx = P.root / "site" / "src" / "content" / "articles" / f"{SLUG}.mdx" m = re.search(r'^version:\s*"([^"]+)"', mdx.read_text(encoding="utf-8"), re.M) if mdx.exists() else None files = sorted(f for f in SITE_OUT.glob("*") if f.is_file() and f.name != "build.json" and f.suffix in (".svg", ".csv", ".json", ".py", ".gz")) manifest = {"article": SLUG, "article_version": m.group(1) if m else None, "seed": SEED, "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "files": {f.name: hashlib.sha256(f.read_bytes()).hexdigest() for f in files}, "note": "Every file listed comes from one build: analysis.py (checksum script_sha256), then the template-chart step. The version box shows " "article_version; if the two disagree the page and the files are from different builds."} text = json.dumps(manifest, indent=1) (OUT / "build.json").write_text(text, encoding="utf-8"); (SITE_OUT / "build.json").write_text(text, encoding="utf-8") def main(argv=None) -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--seed", type=int, default=SEED) ap.add_argument("--manifest-only", action="store_true") args = ap.parse_args(argv) if args.manifest_only: write_build_manifest(); print(f"build manifest -> {SITE_OUT / 'build.json'}"); return 0 rng = np.random.default_rng(args.seed); OUT.mkdir(parents=True, exist_ok=True); SITE_OUT.mkdir(parents=True, exist_ok=True) pub = published() e = load_entries(); d = daily_totals(e) res["published"] = pub res["crz"] = {"first_date": f"{d.index.min():%Y-%m-%d}", "last_date": f"{d.index.max():%Y-%m-%d}", "days": int(len(d)), "year_one": year_one(d, e), "reconcile": reconcile(d, pub), "compare": compare_years(e, d, rng)} res["hub"] = hub_bound() res["hub_vs_gantry"] = hub_vs_gantry(e, res["hub"]) res["tunnels"] = tunnels() res["tlc"] = tlc() res["speeds"] = speeds() res["streets"] = streets() # The Reproduce box: the published entry counts the open file can be checked against. rc = res["crz"]["reconcile"] res["reproduction"] = [ {"statistic": "Vehicles entering the CBD, average of 6-10 January 2025 (MTA, via Streetsblog)", "published": rc["week1_cbd_published"], "reproduced": rc["week1_cbd_file"], "within_margin": abs(rc["week1_diff_pct"]) <= 1.0, "kind": "count"}, {"statistic": "Vehicles entering the CBD, average of 13-17 January 2025 (MTA, via Streetsblog)", "published": rc["week2_cbd_published"], "reproduced": rc["week2_cbd_file"], "within_margin": abs(rc["week2_diff_pct"]) <= 1.0, "kind": "count"}, {"statistic": "Vehicles entering Manhattan south of 60th Street, fall 2022 (NYMTC Table 1A against Table 23A)", "published": res["hub"]["table23a_inbound"][2022], "reproduced": res["hub"]["total"][2022], "within_margin": abs(res["hub"]["total"][2022] - res["hub"]["table23a_inbound"][2022]) <= 1, "kind": "count"}, ] res["reproduction_source"] = "the MTA's first published entry counts and NYMTC's own tables" res["reproduction_note"] = (f"Tolerance: 1% for the MTA's weekly averages, one vehicle for NYMTC's table. The first week reproduces " f"{abs(rc['week1_diff_pct']):.1f}% {'below' if rc['week1_diff_pct'] < 0 else 'above'} the published figure (a revised file " f"or different days), the second within {max(abs(rc['week2_diff_pct']), 0.1):.1f}%; NYMTC's long table agrees with its " "one-year table for 2022, and its misprinted 2023 and 2024 totals are replaced by the one-year tables' figures.") res["as_of"] = res["crz"]["last_date"] cube = monthly_cube(e, d) charts = chart_cubes(res["crz"], res["hub"], res["tunnels"]) (OUT / "results.json").write_text(json.dumps(res, indent=1, default=float), encoding="utf-8") (OUT / "monthly_cube.json").write_text(json.dumps(cube, separators=(",", ":")), encoding="utf-8") (OUT / "charts.json").write_text(json.dumps(charts, separators=(",", ":"), allow_nan=False), encoding="utf-8") from article_parents import fig_from_cube # static twins of the bar charts, for readers without scripts for n, cid in ((6, "groups"), (7, "months"), (9, "tunnels")): (fig_change if cid in ("groups", "months") else fig_from_cube)(charts["charts"][cid], OUT / f"fig{n}_{cid}.svg") shutil.copy2(OUT / f"fig{n}_{cid}.svg", SITE_OUT / f"fig{n}_{cid}.svg") companion(e).to_csv(OUT / "data.csv", index=False) long_view(res["hub"], res["speeds"], res["tlc"], res["streets"], res["tunnels"]).to_csv(OUT / "long_view.csv", index=False) for f in ("results.json", "charts.json", "data.csv", "long_view.csv"): shutil.copy2(OUT / f, SITE_OUT / f) shutil.copy2(Path(__file__), SITE_OUT / "analysis.py") write_build_manifest() c = res["crz"]["compare"] print(f"entries to {c['last_date']}: weekday {c['pooled_weekday']['change_pct']}% ({c['pooled_weekday']['lo_pct']} to {c['pooled_weekday']['hi_pct']})") print(f"hub peak {res['hub']['peak_year']} {res['hub']['peak']:,.0f}; 2018 {res['hub']['total'][2018]:,.0f}; 2024 {res['hub']['total'][2024]:,.0f}; gantry Oct 2025 {res['hub_vs_gantry']['change_pct']}%") print(f"outputs -> {OUT.relative_to(P.root)} and {SITE_OUT.relative_to(P.root)}") return 0 if __name__ == "__main__": raise SystemExit(main())