# Exposición y daño esperado (sub-estudio censal, v2). Para cada simulación:
#  · edificios expuestos (huellas Overture) por umbral h* y radio de búsqueda r (sensibilidad h* ∈ {0,02;0,05;0,10;0,30} m, r ∈ {4;8;16} m)
#  · viviendas/personas por dos métodos: (A) ponderación areal por manzana; (B) dasimétrico: viviendas de cada manzana repartidas
#    entre sus edificios en proporción a la huella (solo huellas 30–2.000 m²), expuestas si su edificio lo está
#  · daño esperado: curvas de vulnerabilidad de Fuchs et al. (2007) V = 0,11h² − 0,02h (h < 2,5 m, V≥0) y
#    Papathoma-Köhle et al. (2012) V = 1 − exp(−0,270((h+1,287)/1,287 − 1)^2,974), con h = espesor máximo junto al edificio
#    (cota superior del depósito). Suma de V = "edificios equivalentes destruidos"; conteo con V ≥ 0,2 y ≥ 0,5.
import sys, json, numpy as np, shapely
from shapely.geometry import shape
from shapely.ops import transform as stf
from shapely.strtree import STRtree
from pyproj import Transformer
from scipy.spatial import cKDTree
tf = Transformer.from_crs('EPSG:4326', 'EPSG:32719', always_xy=True).transform
FU = lambda h: np.clip(0.11 * np.minimum(h, 2.5) ** 2 - 0.02 * np.minimum(h, 2.5), 0, 1)
PK = lambda h: 1 - np.exp(-0.270 * (((np.maximum(h, 0) + 1.287) / 1.287 - 1) ** 2.974))
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])
M = json.load(open('dem/censo/manzanas.geojson'))['features']; MG = [stf(tf, shape(m['geometry'])) for m in M]
V = np.array([m['properties']['TOTAL_VIVI'] or 0 for m in M], float); Pp = np.array([m['properties']['TOTAL_PERS'] or 0 for m in M], float)
# asignación dasimétrica: edificio → manzana (centroide dentro), peso = huella
res_ok = (area >= 30) & (area <= 2000); bm = np.full(len(B), -1)
tree_m = STRtree(MG)
for i in np.nonzero(res_ok)[0]:
    for j in tree_m.query(shapely.Point(cen[i]), predicate='within'): bm[i] = j; break
w = np.zeros(len(B)); tot = np.zeros(len(M))
np.add.at(tot, bm[bm >= 0], area[bm >= 0])
w[bm >= 0] = area[bm >= 0] / tot[bm[bm >= 0]]
viv_b = np.where(bm >= 0, w * V[np.maximum(bm, 0)], 0); per_b = np.where(bm >= 0, w * Pp[np.maximum(bm, 0)], 0)
sin_edif = int(((tot == 0) & (V > 0)).sum())
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()
    dentro = (cen[:, 0] > x0) & (cen[:, 0] < x1) & (cen[:, 1] > y0) & (cen[:, 1] < y1)
    o = {'corrida': f.split('/')[-1][:-4], 'sens': {}}
    hb = {}
    for rr in (4.0, 8.0, 16.0):
        rad = np.sqrt(area / np.pi) + rr; hv = np.zeros(len(B))
        for i in np.nonzero(dentro)[0]:
            ii = tree.query_ball_point(cen[i], rad[i]); hv[i] = h[ii].max() if ii else 0.0
        hb[rr] = hv
    # ponderación areal por manzana
    fr = {u: np.zeros(len(M)) for u in (0.02, 0.05, 0.10, 0.30, 0.5, 1.0)}
    for j, g in enumerate(MG):
        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
        for u in fr: fr[u][j] = A[ins][h[ins] > u].sum() / A[ins].sum()
    for u in (0.02, 0.05, 0.10, 0.30):
        for rr in (4.0, 8.0, 16.0):
            e = dentro & (hb[rr] > u)
            o['sens'][f'h{u}_r{int(rr)}'] = {'edificios': int(e.sum()), 'viv_dasim': round(float(viv_b[e].sum())), 'per_dasim': round(float(per_b[e].sum())),
                                             'viv_areal': round(float((fr[u] * V).sum())), 'per_areal': round(float((fr[u] * Pp).sum()))}
    h8 = np.where(dentro, hb[8.0], 0)
    for u in (0.5, 1.0): o[f'viv_areal_h{u}'] = round(float((fr[u] * V).sum())); o[f'edif_h{u}'] = int((h8 > u).sum()); o[f'viv_dasim_h{u}'] = round(float(viv_b[h8 > u].sum()))
    fu, pk = FU(h8), PK(h8)
    o['danio'] = {'fuchs_equivalentes': round(float(fu.sum()), 1), 'pk_equivalentes': round(float(pk.sum()), 1),
                  'fuchs_V>=0.2': int((fu >= 0.2).sum()), 'pk_V>=0.2': int((pk >= 0.2).sum()), 'fuchs_V>=0.5': int((fu >= 0.5).sum()), 'pk_V>=0.5': int((pk >= 0.5).sum()),
                  'viv_dasim_pk_V>=0.2': round(float(viv_b[pk >= 0.2].sum())), 'viv_dasim_fuchs_V>=0.2': round(float(viv_b[fu >= 0.2].sum()))}
    o['hist_h_edificios'] = np.histogram(h8[h8 > 0.05], bins=[0.05, 0.1, 0.3, 0.5, 1, 1.5, 2, 3, 10])[0].tolist()
    return o
out = {'manzanas_con_viviendas_sin_huellas': sin_edif, 'edificios_asignados': int((bm >= 0).sum()), 'res': [calcular(f) for f in sys.argv[1:]]}
json.dump(out, open('res/exposicion2.json', 'w'), indent=1, ensure_ascii=False)
for o in out['res']: print(o['corrida'], o['sens']['h0.05_r8'], o['danio'])
print('manzanas con viviendas y sin huellas:', sin_edif, '· edificios asignados a manzana:', out['edificios_asignados'])