# Cuenca aportante a los puntos de salida H1 (abanico, aguas arriba de Av. La Plaza) y H2 (cruce del cauce con Fco. Bulnes Correa).
# FABDEM 1" → relleno de depresiones → direcciones D8 → acumulación → delimitación. Área, relieve (R = z_max − z_salida),
# pendiente media y longitud del cauce principal. Salida: cuenca.json + cuenca_H1.geojson / cuenca_H2.geojson
import json, numpy as np, rasterio
if not hasattr(np, "in1d"): np.in1d = np.isin   # pysheds con numpy 2
from rasterio.windows import from_bounds
from pysheds.grid import Grid
from pyproj import Geod
SRC = '/home/linuxoficina/Proyectos/dem/fabdem/S34W071_FABDEM_V1-2.tif'
W, S, E, N = -70.53, -33.445, -70.41, -33.385
with rasterio.open(SRC) as d:
    win = from_bounds(W, S, E, N, d.transform); prof = d.profile.copy(); a = d.read(1, window=win); tr = d.window_transform(win)
prof.update(width=a.shape[1], height=a.shape[0], transform=tr, nodata=-9999, dtype='float32')
with rasterio.open('dem/fabdem_cuenca.tif', 'w', **prof) as o: o.write(a.astype('float32'), 1)
g = Grid.from_raster('dem/fabdem_cuenca.tif'); dem = g.read_raster('dem/fabdem_cuenca.tif')
fd = g.resolve_flats(g.fill_depressions(g.fill_pits(dem))); fdir = g.flowdir(fd); acc = g.accumulation(fdir)
geod = Geod(ellps='WGS84'); dlat = abs(tr.e); dlon = tr.a
cell_m2 = abs(geod.geometry_area_perimeter(__import__('shapely.geometry', fromlist=['box']).box(-70.47, -33.41, -70.47 + dlon, -33.41 + dlat))[0])
out = {}
for nombre, (lat, lon) in {'H1': (-33.40030, -70.49960), 'H2': (-33.40590, -70.51630)}.items():
    x, y = g.snap_to_mask(acc > 300, (lon, lat))
    c = g.catchment(x=x, y=y, fdir=fdir, xytype='coordinate'); m = np.asarray(c, bool)
    area_km2 = m.sum() * cell_m2 / 1e6; z = np.asarray(dem)[m]
    r, cidx = ~tr * (x, y); zs = float(np.asarray(dem)[int(r if False else int((y - tr.f) / tr.e)), int((x - tr.c) / tr.a)])
    shapes = list(g.polygonize(c.astype(np.uint8) if hasattr(c, 'astype') else c))
    poly = [s for s, v in shapes if v == 1]
    json.dump({'type': 'FeatureCollection', 'features': [{'type': 'Feature', 'properties': {'nombre': nombre}, 'geometry': p} for p in poly]}, open(f'cuenca_{nombre}.geojson', 'w'))
    out[nombre] = {'salida_lat_lon_ajustada': [round(y, 5), round(x, 5)], 'area_km2': round(area_km2, 3), 'z_max_m': round(float(z.max()), 1),
                   'z_salida_m': round(zs, 1), 'relieve_R_m': round(float(z.max()) - zs, 1), 'z_media_m': round(float(z.mean()), 1),
                   'celdas': int(m.sum()), 'acumulacion_en_salida': int(np.asarray(acc)[int((y - tr.f) / tr.e), int((x - tr.c) / tr.a)])}
json.dump(out, open('cuenca.json', 'w'), indent=1); print(json.dumps(out, indent=1))