"""Reproduce Texas named-highway fatal-crash counts from official FARS annual ZIPs.
Python 3 stdlib only. Usage: python analyze_roads.py INPUT_ZIP_DIR OUTPUT_DIR
Download FARS2022NationalCSV.zip through FARS2024NationalCSV.zip from
https://static.nhtsa.gov/nhtsa/downloads/FARS/{year}/National/FARS{year}NationalCSV.zip
No traffic exposure data: output is crash burden, NOT individual travel risk.
"""
import csv, io, json, re, zipfile, pathlib, sys, collections, hashlib
YEARS=(2022,2023,2024)
# FARS uses SR for state highways. Preserve raw values in the audit.
PREFIX=re.compile(r'^(I|US|SR)-0*(\d+)(?=$|[^0-9])')
EXCLUDE=re.compile(r'\b(?:RAMP|FRONTAGE|FRNTG|FRONT|SERVICE|SERV|SVRD|SV|FEEDER|ACCESS|CONNECTOR|BUSINESS|BUS|BU|SPUR)\b')
def load(zpath):
    with zipfile.ZipFile(zpath) as z:
        name=next(n for n in z.namelist() if n.lower().endswith('accident.csv'))
        raw=z.read(name)
    try: text=raw.decode('utf-8-sig')
    except UnicodeDecodeError: text=raw.decode('cp1252')
    return list(csv.DictReader(io.StringIO(text)))
def classify(r):
    road=r['TWAY_ID'].strip().upper()
    if r['ROUTE'] not in ('1','2','3'): return None,'not_interstate_us_or_state_system'
    if EXCLUDE.search(road) or r['RELJCT2']=='20': return None,'explicit_ramp_frontage_business_or_connector'
    m=PREFIX.match(road)
    if not m: return None,'unmatched_primary_road_name'
    if {'I':'1','US':'2','SR':'3'}[m[1]]!=r['ROUTE']: return None,'route_system_name_disagreement'
    return ('SH' if m[1]=='SR' else m[1])+'-'+str(int(m[2])),'included'
def run(src,out):
    out.mkdir(parents=True,exist_ok=True)
    totals=[];matched=[];audit=[];seen=set();sources=[];counters=collections.Counter()
    for year in YEARS:
        p=src/f'FARS{year}NationalCSV.zip';rows=[r for r in load(p) if r['STATE']=='48'];assert rows
        sources.append({'year':year,'url':f'https://static.nhtsa.gov/nhtsa/downloads/FARS/{year}/National/FARS{year}NationalCSV.zip','sha256':hashlib.sha256(p.read_bytes()).hexdigest()})
        totals.append({'year':year,'fatalCrashes':len(rows),'deaths':sum(int(r['FATALS']) for r in rows)})
        for r in rows:
            key=(year,r['ST_CASE']);assert key not in seen;seen.add(key)
            assert int(r['YEAR'])==year and int(r['FATALS'])>0
            route,reason=classify(r);counters[reason]+=1
            audit.append({'year':year,'case':r['ST_CASE'],'primaryRoad':r['TWAY_ID'],'secondaryRoad':r['TWAY_ID2'],'routeSystem':r['ROUTE'],'relationToJunction':r['RELJCT2'],'group':route or '', 'status':reason})
            if not route:continue
            lat,lon=float(r['LATITUDE']),float(r['LONGITUD'])
            valid=25.80<=lat<=36.55 and -106.70<=lon<=-93.45
            matched.append({'year':year,'case':r['ST_CASE'],'route':route,'road':r['TWAY_ID'],'county':re.sub(r'\s*\(\d+\)$','',r['COUNTYNAME']).title(),'city':r['CITYNAME'].title() if r['CITY'] not in ('0','9997','9998','9999') else None,'deaths':int(r['FATALS']),'lat':round(lat,5) if valid else None,'lon':round(lon,5) if valid else None})
    groups=[];route_year=[];county_rows=[]
    for route in sorted(set(r['route'] for r in matched)):
        rs=[r for r in matched if r['route']==route]
        annual=[{'year':y,'fatalCrashes':sum(r['year']==y for r in rs),'deaths':sum(r['deaths'] for r in rs if r['year']==y)} for y in YEARS]
        counties=[]
        for county in sorted(set(r['county'] for r in rs)):
            cs=[r for r in rs if r['county']==county];counties.append({'county':county,'fatalCrashes':len(cs),'deaths':sum(r['deaths'] for r in cs)})
        counties.sort(key=lambda x:(-x['fatalCrashes'],x['county']))
        groups.append({'route':route,'fatalCrashes':len(rs),'deaths':sum(r['deaths'] for r in rs),'mappedCrashes':sum(r['lat'] is not None for r in rs),'annual':annual,'counties':counties,'points':[[r['lon'],r['lat'],r['year'],r['deaths'],r['county']] for r in rs if r['lat'] is not None]})
        route_year.extend({'route':route,**a} for a in annual);county_rows.extend({'route':route,**a} for a in counties)
    groups.sort(key=lambda x:(-x['fatalCrashes'],x['route']))
    result={'years':list(YEARS),'source':'NHTSA FARS','retrieved':'2026-09-18','scope':'Texas crashes with at least one death within 30 days; primary road named as an Interstate, U.S. or state highway matching its coded route system. Explicit ramp/frontage/business/connector/spur names and entrance/exit-ramp locations excluded. Branch suffixes grouped under the parent route number; secondary road IDs not used. Counts are not risk per mile or trip.','sources':sources,'statewide':totals,'includedCrashes':len(matched),'includedDeaths':sum(r['deaths'] for r in matched),'mappedCrashes':sum(r['lat'] is not None for r in matched),'classificationCounts':dict(counters),'groups':groups}
    assert sum(g['fatalCrashes'] for g in groups)==len(matched)
    assert sum(counters.values())==sum(r['fatalCrashes'] for r in totals)
    (out/'roads.json').write_text(json.dumps(result,separators=(',',':')))
    for filename,rows in [('route-year.csv',route_year),('route-county.csv',county_rows),('classification-audit.csv',audit)]:
        with (out/filename).open('w',newline='') as f:
            w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
    (out/'analysis-summary.json').write_text(json.dumps({k:v for k,v in result.items() if k!='groups'}|{'routeGroups':len(groups),'top15':[{k:v for k,v in g.items() if k not in ('points','counties')} for g in groups[:15]]},indent=2))
    print(json.dumps({'totalTexasCrashes':sum(t['fatalCrashes'] for t in totals),'includedCrashes':len(matched),'routeGroups':len(groups),'classification':dict(counters),'top10':[{k:g[k] for k in ['route','fatalCrashes','deaths']} for g in groups[:10]]},indent=2))
if __name__=='__main__':run(pathlib.Path(sys.argv[1]),pathlib.Path(sys.argv[2]))