# Exposición de edificios, viviendas y personas a la propagación simulada (sub-estudio de exposición censal).
# Insumos: simulación (.npz), huellas de edificios Overture (dem/overture/edificios.geojson) y manzanas del Censo 2017
# (dem/censo/manzanas.geojson; misma fuente pública que usa Hirely Broker).
# Supuestos: (1) un edificio está expuesto si el espesor máximo simulado a ≤ 8 m de su huella supera el umbral;
# (2) viviendas y personas se distribuyen uniformemente en el área de cada manzana: la manzana aporta en proporción a la fracción
# de su área que supera cada umbral. Salida: JSON por escenario y GeoJSON de manzanas con su exposición.
import sys, json, numpy as np, shapely
from shapely.geometry import shape
from shapely.ops import transform as stf
from pyproj import Transformer
from scipy.spatial import cKDTree
tf = Transformer.from_crs('EPSG:4326', 'EPSG:32719', always_xy=True).transform
UMB = (0.05, 0.3, 0.5, 1.0)
def calcular(f):
    r = np.load(f); xy = np.c_[r['cx'], r['cy']]; tree = cKDTree(xy); h = r['hmax']; A = r['area']
    x0, x1, y0, y1 = r['cx'].min(), r['cx'].max(), r['cy'].min(), r['cy'].max()
    out = {'corrida': f.split('/')[-1].replace('.npz', ''), 'params': json.loads(str(r['params']))}
    B = [stf(tf, shape(g['geometry'])) for g in json.load(open('dem/overture/edificios.geojson'))['features']]
    area = np.array([b.area for b in B]); cen = np.array([[b.centroid.x, b.centroid.y] for b in B])
    dentro = (cen[:, 0] > x0) & (cen[:, 0] < x1) & (cen[:, 1] > y0) & (cen[:, 1] < y1)
    rad = np.sqrt(area / np.pi) + 8.0; hb = np.zeros(len(B))
    for i in np.nonzero(dentro)[0]:
        ii = tree.query_ball_point(cen[i], rad[i]); hb[i] = h[ii].max() if ii else 0.0
    out['edificios_dominio'] = int(dentro.sum()); out['edificios'] = {f'>{u}': int(((hb > u) & dentro).sum()) for u in UMB}
    casas = dentro & (area >= 40) & (area <= 400)
    out['huella_m2'] = {'media_todas': round(float(area[casas].mean()), 1), 'mediana_todas': round(float(np.median(area[casas])), 1),
                        'media_expuestas': round(float(area[casas & (hb > 0.05)].mean()), 1) if (casas & (hb > 0.05)).any() else None}
    M = json.load(open('dem/censo/manzanas.geojson'))['features']; viv = {u: 0.0 for u in UMB}; per = {u: 0.0 for u in UMB}; feats = []
    for m in M:
        g = stf(tf, shape(m['geometry']))
        if g.area <= 0: continue
        bx0, by0, bx1, by1 = g.bounds
        if bx1 < x0 or bx0 > x1 or by1 < y0 or by0 > y1: continue
        idx = np.nonzero((xy[:, 0] > bx0) & (xy[:, 0] < bx1) & (xy[:, 1] > by0) & (xy[:, 1] < by1))[0]
        if not len(idx): continue
        ins = idx[shapely.contains_xy(g, xy[idx, 0], xy[idx, 1])]
        if not len(ins): continue
        At = A[ins].sum(); V = m['properties']['TOTAL_VIVI'] or 0; P = m['properties']['TOTAL_PERS'] or 0; fr = {}
        for u in UMB:
            fr[u] = float(A[ins][h[ins] > u].sum() / At); viv[u] += fr[u] * V; per[u] += fr[u] * P
        if fr[0.05] > 0:
            feats.append({'type': 'Feature', 'properties': {'manzent': m['properties']['MANZENT_I'], 'comuna': m['properties']['NOM_COMUNA'], 'viviendas': V, 'personas': P,
                          'frac_05': round(fr[0.05], 3), 'frac_50': round(fr[0.5], 3), 'viv_expuestas': round(fr[0.05] * V, 1)}, 'geometry': m['geometry']})
    out['viviendas'] = {f'>{u}': round(v) for u, v in viv.items()}; out['personas'] = {f'>{u}': round(v) for u, v in per.items()}
    out['manzanas_alcanzadas'] = len(feats)
    json.dump({'type': 'FeatureCollection', 'features': feats}, open(f'res/exposicion_{out["corrida"]}.geojson', 'w'))
    return out
res = [calcular(f) for f in sys.argv[1:]]
json.dump(res, open('res/exposicion.json', 'w'), indent=1, ensure_ascii=False)
for o in res: print(o['corrida'], 'edif', o['edificios'], 'viv', o['viviendas'], 'pers', o['personas'], 'manz', o['manzanas_alcanzadas'], o['huella_m2'])