import os, json, glob
from pathlib import Path
import numpy as np, pandas as pd, duckdb
from scipy.spatial import cKDTree
from transbordo_utils import novo_mapa, resolver_cache

START_DATE, END_DATE = "2026-01-01", "2026-01-31"
os.makedirs("/tmp/duckdb_spill", exist_ok=True)
PARQUET    = resolver_cache(f"usuarios_limpo_{START_DATE}_{END_DATE}.parquet")
PONTOS_DIR = resolver_cache("itinerarios/pontos")
GAP_MAX_KM = 0.5      # confiança: próximo embarque a <=500m de uma parada da linha
DESLOC_MIN_KM = 0.3   # deslocamento mínimo p/ contar como O-D real
print("parquet:", PARQUET); print("pontos:", PONTOS_DIR)

parquet: cache_urbs/usuarios_limpo_2026-01-01_2026-01-31.parquet
pontos: cache_urbs/itinerarios/pontos

# --- Paradas por linha -> KDTree ---
from transbordo_utils import haversine_km as hav
stops={}
for fp in glob.glob(f"{PONTOS_DIR}/*.json"):
    cod=os.path.basename(fp)[:-5]
    try: d=json.load(open(fp))
    except: continue
    pts=[]
    for s in d:
        try: pts.append((float(str(s["LAT"]).replace(",",".")), float(str(s["LON"]).replace(",","."))))
        except: pass
    if pts:
        arr=np.array(pts); stops[cod]=(cKDTree(arr), arr)
print(f"linhas com paradas: {len(stops)}")

linhas com paradas: 308

# --- Viagens a resolver (não-última de jornadas >=2 usos) via DuckDB ---
con=duckdb.connect(config={"memory_limit":"500MB","temp_directory":"/tmp/duckdb_spill","threads":"2"})
df=con.execute(f"""
WITH base AS (SELECT NUMEROCARTAO, CAST(DATA_HORA_UTILIZACAO AS DATE) dia, DATA_HORA_UTILIZACAO dh,
   CODLINHA, LATITUDE lat, LONGITUDE lon FROM read_parquet('{PARQUET}')
   WHERE LATITUDE IS NOT NULL AND LONGITUDE IS NOT NULL),
j AS (SELECT *, dense_rank() OVER (ORDER BY NUMEROCARTAO,dia) jid,
   COUNT(*) OVER (PARTITION BY NUMEROCARTAO,dia) n,
   LEAD(lat) OVER w nlat, LEAD(lon) OVER w nlon,
   FIRST_VALUE(lat) OVER w olat, FIRST_VALUE(lon) OVER w olon,
   FIRST_VALUE(CODLINHA) OVER w olinha
   FROM base WINDOW w AS (PARTITION BY NUMEROCARTAO,dia ORDER BY dh))
SELECT CAST(jid AS INTEGER) jid, CODLINHA,
   CAST(nlat AS FLOAT) nlat, CAST(nlon AS FLOAT) nlon,
   CAST(olat AS FLOAT) olat, CAST(olon AS FLOAT) olon, olinha
FROM j WHERE n>=2 AND nlat IS NOT NULL
""").df()
con.close()
print(f"trips a resolver: {len(df):,} | mem {df.memory_usage(deep=True).sum()/1e6:.0f} MB")

trips a resolver: 217,515 | mem 9 MB

# --- Desembarque estimado = parada da linha i mais próxima do PRÓXIMO embarque ---
df["alat"]=np.nan; df["alon"]=np.nan; df["gap_km"]=np.nan
for cod,grp in df.groupby("CODLINHA"):
    if cod not in stops: continue
    tree,arr=stops[cod]
    q=grp[["nlat","nlon"]].to_numpy()
    dist,idx=tree.query(q); al=arr[idx]
    df.loc[grp.index,"alat"]=al[:,0]; df.loc[grp.index,"alon"]=al[:,1]
    df.loc[grp.index,"gap_km"]=hav(q[:,0],q[:,1],al[:,0],al[:,1])
df["dist_od_km"]=hav(df["olat"].values,df["olon"].values,df["alat"].values,df["alon"].values)
conf = df[df["gap_km"]<=GAP_MAX_KM]
print(f"trips confiáveis (gap<={GAP_MAX_KM}km): {len(conf):,} ({100*len(conf)/len(df):.1f}%)")
print(f"gap mediano: {df['gap_km'].median()*1000:.0f} m")

trips confiáveis (gap<=0.5km): 177,732 (81.7%)
gap mediano: 17 m

# --- O-D por jornada: destino = desembarque estimado mais distante da origem ---
# filtro de caixa de Curitiba (remove coords-lixo, ex. 0,0 / outliers)
def _emcwb(la, lo):
    return (la.between(-25.75,-25.25)) & (lo.between(-49.55,-49.10))
conf = conf[_emcwb(conf["olat"],conf["olon"]) & _emcwb(conf["alat"],conf["alon"])]
dest = conf.loc[conf.groupby("jid")["dist_od_km"].idxmax()].copy()
od_real = dest[dest["dist_od_km"]>=DESLOC_MIN_KM]   # com deslocamento real
print(f"jornadas (com destino confiável): {len(dest):,}")
print(f"  com deslocamento real (>={DESLOC_MIN_KM}km): {len(od_real):,} ({100*len(od_real)/len(dest):.1f}%)")
print(f"  no mesmo local (<{DESLOC_MIN_KM}km): {len(dest)-len(od_real):,}")
print("\nDistância O-D estimada (km) — jornadas com deslocamento real:")
print(od_real["dist_od_km"].describe(percentiles=[.5,.75,.9])[["mean","50%","75%","90%","max"]].round(2).to_string())
print("\n>> Antes (04, ponto mais distante de EMBARQUE): mediana ~2,8 km.")
print(">> Agora estimamos o DESEMBARQUE (parada de descida), não só onde reembarcou.")
out = Path(PARQUET).parent / f"od_desembarque_{START_DATE}_{END_DATE}.parquet"
od_real[["jid","olinha","olat","olon","CODLINHA","alat","alon","dist_od_km","gap_km"]].to_parquet(out, index=False)
print(f"\nsalvo: {out}")

jornadas (com destino confiável): 141,961
  com deslocamento real (>=0.3km): 93,582 (65.9%)
  no mesmo local (<0.3km): 48,379

Distância O-D estimada (km) — jornadas com deslocamento real:
mean     4.05
50%      3.60
75%      5.91
90%      7.64
max     27.96

>> Antes (04, ponto mais distante de EMBARQUE): mediana ~2,8 km.
>> Agora estimamos o DESEMBARQUE (parada de descida), não só onde reembarcou.

salvo: cache_urbs/od_desembarque_2026-01-01_2026-01-31.parquet

# --- Mapa O-D estimado (origem -> desembarque), amostra ---
import folium
am = od_real.sample(n=min(600,len(od_real)), random_state=42)
mp = novo_mapa(location=[-25.45,-49.29], zoom_start=12)
for _,r in am.iterrows():
    folium.PolyLine([[r.olat,r.olon],[r.alat,r.alon]], color="#8e44ad", weight=2, opacity=0.5,
        popup=f"{r.olinha} → desembarque (linha {r.CODLINHA})<br>{r.dist_od_km:.1f} km").add_to(mp)
    folium.CircleMarker([r.olat,r.olon], radius=3, color="#2980b9", fill=True).add_to(mp)
    folium.CircleMarker([r.alat,r.alon], radius=3, color="#8e44ad", fill=True).add_to(mp)
mp.save("mapa_od_desembarque.html")
print(f"amostra: {len(am)} | azul=origem (embarque), roxo=desembarque estimado | salvo: mapa_od_desembarque.html")
mp

amostra: 600 | azul=origem (embarque), roxo=desembarque estimado | salvo: mapa_od_desembarque.html

import os, json
from pathlib import Path
import numpy as np, pandas as pd, geopandas as gpd
import matplotlib.pyplot as plt, folium
from transbordo_utils import resolver_cache, carregar_pares_oficiais

START_DATE, END_DATE = "2026-01-01", "2026-01-31"
bairros = gpd.read_file(resolver_cache("bairros_curitiba.geojson"))[["nome","nm_regional","geometry"]].to_crs(4326)
od = pd.read_parquet(resolver_cache(f"od_desembarque_{START_DATE}_{END_DATE}.parquet"))
pares = carregar_pares_oficiais(resolver_cache("matriz_integracao.json"))
print(f"bairros: {len(bairros)} | jornadas O-D: {len(od):,}")

bairros: 75 | jornadas O-D: 93,582

# --- Spatial join: cada O e D cai no bairro OFICIAL (sem ruído de nome) ---
def bairro_de(la, lo):
    g=gpd.GeoDataFrame(geometry=gpd.points_from_xy(lo, la), crs=4326)
    return gpd.sjoin(g, bairros, how="left", predicate="within")["nome"].values
od["b_o"]=bairro_de(od["olat"].values, od["olon"].values)
od["b_d"]=bairro_de(od["alat"].values, od["alon"].values)
od["tem_integ"]=[(str(o),str(d)) in pares for o,d in zip(od["olinha"],od["CODLINHA"])]
od=od.dropna(subset=["b_o","b_d"])
odd=od[od["b_o"]!=od["b_d"]].copy()   # com troca de bairro
print(f"jornadas em Curitiba: {len(od):,} | com troca de bairro: {len(odd):,}")
print(f"% com integração oficial no par de linhas: {100*odd['tem_integ'].mean():.1f}%")

jornadas em Curitiba: 92,580 | com troca de bairro: 72,119
% com integração oficial no par de linhas: 8.8%

# --- Agrega por bairro: origens, destinos, % sem integração (na chegada) ---
orig = odd.groupby("b_o").size().rename("origens")
dest = odd.groupby("b_d").size().rename("destinos")
semint = odd.groupby("b_d")["tem_integ"].apply(lambda s: 100*(1-s.mean())).rename("pct_sem_integ")
bg = bairros.merge(orig, left_on="nome", right_index=True, how="left")\
            .merge(dest, left_on="nome", right_index=True, how="left")\
            .merge(semint, left_on="nome", right_index=True, how="left")
for c in ["origens","destinos"]: bg[c]=bg[c].fillna(0)
bg["cent"]=bg.geometry.to_crs(31982).centroid.to_crs(4326)  # centróide em CRS projetado
print("agregado por bairro pronto")

agregado por bairro pronto

# --- Choropleths: ORIGENS, DESTINOS e % SEM INTEGRAÇÃO (na chegada) ---
fig, ax = plt.subplots(1, 3, figsize=(20,8))
bg.plot(column="origens", cmap="Blues", legend=True, ax=ax[0], edgecolor="white", linewidth=.3,
        legend_kwds={"shrink":.6})
ax[0].set_title("Origens por bairro (de onde saem)", fontsize=13)
bg.plot(column="destinos", cmap="OrRd", legend=True, ax=ax[1], edgecolor="white", linewidth=.3,
        legend_kwds={"shrink":.6})
ax[1].set_title("Destinos por bairro (onde chegam)", fontsize=13)
bg.plot(column="pct_sem_integ", cmap="RdYlGn_r", legend=True, ax=ax[2], edgecolor="white", linewidth=.3,
        vmin=0, vmax=100, legend_kwds={"shrink":.6}, missing_kwds={"color":"#eee"})
ax[2].set_title("% das chegadas SEM integração oficial", fontsize=13)
for a in ax: a.axis("off")
# rótulos dos top destinos
for _,r in bg.sort_values("destinos",ascending=False).head(8).iterrows():
    ax[1].annotate(r["nome"], (r["cent"].x, r["cent"].y), fontsize=7, ha="center")
plt.tight_layout(); plt.show()

# --- Mapa de corredores top-40 (linhas entre centróides; cor=integração, espessura=volume) ---
fl=(odd.groupby(["b_o","b_d"]).agg(vol=("jid","size"), integ=("tem_integ","mean")).reset_index()
      .sort_values("vol",ascending=False))
cmap_c={r["nome"]:(r["cent"].y, r["cent"].x) for _,r in bg.iterrows()}
mp=novo_mapa(location=[-25.45,-49.29], zoom_start=11)
folium.GeoJson(bairros, style_function=lambda x:{"fillOpacity":0,"color":"#bbb","weight":.5}).add_to(mp)
vmax=fl["vol"].max()
for r in fl.head(40).itertuples():
    if r.b_o not in cmap_c or r.b_d not in cmap_c: continue
    o=cmap_c[r.b_o]; d=cmap_c[r.b_d]
    cor="#27ae60" if r.integ>=0.5 else "#c0392b"
    folium.PolyLine([o,d], color=cor, weight=1+8*r.vol/vmax, opacity=0.65,
        popup=f"{r.b_o} → {r.b_d}<br>{r.vol} jornadas · {100*r.integ:.0f}% c/ integração").add_to(mp)
for _,r in bg.sort_values("destinos",ascending=False).head(10).iterrows():
    folium.CircleMarker([r["cent"].y, r["cent"].x], radius=3+10*r["destinos"]/bg["destinos"].max(),
        color="#2c3e50", fill=True, fill_opacity=.7, popup=f'{r["nome"]}: {int(r["destinos"])} chegadas').add_to(mp)
mp.save("mapa_corredores_bairros.html")
print("salvo: mapa_corredores_bairros.html (verde=tem integração, vermelho=sem)")
mp

salvo: mapa_corredores_bairros.html (verde=tem integração, vermelho=sem)

# --- Matriz bairro->bairro (top 15 bairros por fluxo total) ---
topb = pd.concat([odd["b_o"],odd["b_d"]]).value_counts().head(15).index.tolist()
sub=odd[odd["b_o"].isin(topb) & odd["b_d"].isin(topb)]
mat=sub.pivot_table(index="b_o", columns="b_d", values="jid", aggfunc="size", fill_value=0).reindex(index=topb, columns=topb, fill_value=0)
fig,ax=plt.subplots(figsize=(11,9))
im=ax.imshow(mat.values, cmap="magma_r")
ax.set_xticks(range(len(topb))); ax.set_xticklabels(topb, rotation=90, fontsize=8)
ax.set_yticks(range(len(topb))); ax.set_yticklabels(topb, fontsize=8)
ax.set_xlabel("DESTINO"); ax.set_ylabel("ORIGEM"); ax.set_title("Matriz O-D bairro→bairro (top 15)")
plt.colorbar(im, shrink=.7, label="jornadas"); plt.tight_layout(); plt.show()

# --- Top 15 corredores SEM integração oficial (candidatos) ---
sem=fl[fl["integ"]<0.5].head(15)[::-1]
plt.figure(figsize=(10,6))
plt.barh([f"{r.b_o} → {r.b_d}" for r in sem.itertuples()], sem["vol"], color="#c0392b")
plt.title("Top 15 corredores bairro→bairro SEM integração oficial"); plt.xlabel("jornadas")
plt.tight_layout(); plt.show()
fl.to_csv("fluxos_bairros_oficial.csv", index=False)
print("exportado: fluxos_bairros_oficial.csv")

exportado: fluxos_bairros_oficial.csv

# --- OD Map (small multiples): para cada bairro de ORIGEM, onde os trips vão ---
top_orig = odd["b_o"].value_counts().head(9).index.tolist()
fig, axes = plt.subplots(3,3, figsize=(18,16))
for ax, orig in zip(axes.ravel(), top_orig):
    sub = odd[odd["b_o"]==orig].groupby("b_d").size().rename("v")
    tmp = bairros.merge(sub, left_on="nome", right_index=True, how="left")
    tmp["v"] = tmp["v"].fillna(0)
    tmp.plot(column="v", cmap="OrRd", ax=ax, edgecolor="white", linewidth=.2)
    bairros[bairros["nome"]==orig].plot(ax=ax, facecolor="none", edgecolor="#2980b9", linewidth=2)
    ax.set_title(f"Origem: {orig}", fontsize=10); ax.axis("off")
plt.suptitle("OD Map (small multiples) — destinos de cada bairro de origem (azul = origem)", fontsize=15)
plt.tight_layout(); plt.show()
print("Técnica Wood/Dykes/Slingsby: cada painel isola UMA origem -> sem oclusão de linhas.")

Técnica Wood/Dykes/Slingsby: cada painel isola UMA origem -> sem oclusão de linhas.

# --- Export para flowmap.blue (interativo, sem código) ---
loc = bg[["nome","cent"]].copy()
loc["lat"]=loc["cent"].y; loc["lon"]=loc["cent"].x
loc = loc.rename(columns={"nome":"id"}); loc["name"]=loc["id"]
loc[["id","name","lat","lon"]].to_csv("flowmap_locations.csv", index=False)
flows = (odd.groupby(["b_o","b_d"]).size().reset_index(name="count")
            .rename(columns={"b_o":"origin","b_d":"dest"}))
flows.to_csv("flowmap_flows.csv", index=False)
print(f"flowmap_locations.csv ({len(loc)} bairros) e flowmap_flows.csv ({len(flows)} fluxos)")
print("Use em https://www.flowmap.blue/ : crie um Google Sheet com abas 'locations' (id,name,lat,lon)")
print("e 'flows' (origin,dest,count) colando esses CSVs — gera um mapa de fluxo interativo.")

flowmap_locations.csv (75 bairros) e flowmap_flows.csv (2149 fluxos)
Use em https://www.flowmap.blue/ : crie um Google Sheet com abas 'locations' (id,name,lat,lon)
e 'flows' (origin,dest,count) colando esses CSVs — gera um mapa de fluxo interativo.

07 — Origem-Destino e visualização por bairro¶

1. Desembarque estimado (trip-chaining)¶

2. Visualização O-D por bairro (oficial)¶