// ===== 0. Insumos =====
var cuencas = ee.FeatureCollection('projects/chalabeqdaunesco/assets/Cuenca_unesco');
var drenaje = ee.FeatureCollection('projects/chalabeqdaunesco/assets/Corrientes_agua_unesco');
var aoi = cuencas.geometry();
Map.centerObject(cuencas, 9);
print('N° de cuencas:', cuencas.size());
print('N° de tramos de drenaje:', drenaje.size());

// Modelo de elevación Copernicus 30 m
var glo = ee.ImageCollection('COPERNICUS/DEM/GLO30').filterBounds(aoi).select('DEM');
var dem = glo.mosaic().setDefaultProjection(glo.first().projection()).clip(aoi);
var pendiente = ee.Terrain.slope(dem);

// MERIT Hydro (90 m): área de aporte y HAND
var merit = ee.Image('MERIT/Hydro/v1_0_1').clip(aoi);
var upa = merit.select('upa');
var hand = merit.select('hnd');
var tanB = ee.Terrain.slope(merit.select('elv'))
             .multiply(Math.PI / 180).tan().max(0.001);
var twi = upa.multiply(1e6).divide(90).divide(tanB).log().rename('twi');

// Distancia a la red de drenaje (m)
var distDren = drenaje.distance(5000).unmask(5000).clip(aoi);

// Lluvia: Rx1day medio 1991-2020
var chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY').filterBounds(aoi);
var rx1day = ee.ImageCollection(ee.List.sequence(1991, 2020).map(function (y) {
  return chirps.filter(ee.Filter.calendarRange(y, y, 'year')).max();
})).mean().clip(aoi).rename('rx1day');

// NDVI mediano Sentinel-2
var ndvi = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
  .filterBounds(aoi).filterDate('2024-01-01', '2025-12-31')
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
  .median().normalizedDifference(['B8', 'B4']).clip(aoi).rename('ndvi');

// ===== 1. Morfometría por cuenca =====
var cuencasM = dem.reduceRegions({
  collection: cuencas, reducer: ee.Reducer.minMax(), scale: 30
}).map(function (f) {
  var area = f.geometry().area(1);
  var relieve = ee.Number(f.get('max')).subtract(f.get('min'));
  var largo = drenaje.filterBounds(f.geometry())
                .map(function (l) { return l.intersection(f.geometry(), 1); })
                .geometry().length(1);
  return f.set({
    area_km2: area.divide(1e6),
    relieve_m: relieve,
    melton: relieve.divide(area.sqrt()),
    dd_km_km2: largo.divide(1000).divide(area.divide(1e6))
  });
});
print('Morfometría por cuenca', cuencasM);
var melton = ee.Image().float().paint(cuencasM, 'melton').rename('melton');
var dd     = ee.Image().float().paint(cuencasM, 'dd_km_km2').rename('dd');

// ===== 2. Reclasificación a 1-5 =====
function porUmbrales(img, cortes, invertir) {
  var c = ee.Image(1);
  cortes.forEach(function (b) { c = c.add(img.gt(b)); });
  return invertir ? ee.Image(6).subtract(c) : c;
}
function porPercentiles(img, invertir) {
  var p = ee.List(img.reduceRegion({
    reducer: ee.Reducer.percentile([20, 40, 60, 80]),
    geometry: aoi, scale: 90, bestEffort: true, maxPixels: 1e13
  }).values());
  var c = ee.Image(1);
  for (var i = 0; i < 4; i++) c = c.add(img.gt(ee.Number(p.get(i))));
  return invertir ? ee.Image(6).subtract(c) : c;
}

// ===== 3a. Susceptibilidad ALUVIONAL =====
var sAluv = porUmbrales(pendiente, [5, 15, 25, 35]).multiply(0.25)
  .add(porUmbrales(melton, [0.2, 0.3, 0.45, 0.6]).multiply(0.15))
  .add(porPercentiles(rx1day).multiply(0.15))
  .add(porUmbrales(upa.log10(), [-1, 0, 1, 2]).multiply(0.10))
  .add(porPercentiles(dd).multiply(0.10))
  .add(porUmbrales(distDren, [50, 150, 300, 600], true).multiply(0.10))
  .add(porUmbrales(ndvi, [0.1, 0.2, 0.35, 0.5], true).multiply(0.15))
  .rename('susc_aluvional');

// ===== 3b. Susceptibilidad a ANEGAMIENTO =====
var sAneg = porUmbrales(hand, [1, 3, 5, 10], true).multiply(0.30)
  .add(porUmbrales(pendiente, [0.5, 1, 2, 5], true).multiply(0.20))
  .add(porPercentiles(twi).multiply(0.20))
  .add(porUmbrales(distDren, [100, 250, 500, 1000], true).multiply(0.15))
  .add(porPercentiles(rx1day).multiply(0.15))
  .rename('susc_anegamiento');

// ===== 4. Visualización =====
var pal = ['1a9850', '91cf60', 'fee08b', 'fc8d59', 'd73027'];
Map.addLayer(drenaje, {color: '0000ff'}, 'Red de drenaje', false);
Map.addLayer(sAneg, {min: 1, max: 5, palette: pal}, 'Susceptibilidad anegamiento', false);
Map.addLayer(sAluv, {min: 1, max: 5, palette: pal}, 'Susceptibilidad aluvional');
Map.addLayer(ee.Image().paint(cuencas, 0, 2), {palette: ['000000']}, 'Límite de cuencas');

// ===== 5. Exportación =====
Export.image.toDrive({
  image: sAluv.addBands(sAneg).toFloat(),
  description: 'susceptibilidad',
  folder: 'GEE_geologia_ambiental',
  region: aoi, scale: 30, crs: 'EPSG:4326', maxPixels: 1e13
});
Export.table.toDrive({
  collection: cuencasM,
  description: 'morfometria_cuencas',
  folder: 'GEE_geologia_ambiental',
  fileFormat: 'SHP'
});