Generando forzamiento meteorológico con Google Earth Engine

HydroPol2D requiere datos de forzamiento meteorológico distribuidos espacialmente para impulsar simulaciones hidrológicas. Este script Google Earth Engine (GEE) automatiza la extracción y el formato de entradas de lluvia y evapotranspiración (ETP) de múltiples conjuntos de datos globales. El flujo de trabajo genera archivos CSV en el formato exacto requerido por HydroPol2D, con datos extraídos en ubicaciones de medidores virtuales distribuidos en todo el dominio de modelado.

Qué hace el flujo de trabajo

El script procesa y exporta tres conjuntos de datos meteorológicos globales:

Entradas de lluvia

• CHIRPS (Precipitación infrarroja del grupo de peligros climáticos con datos de la estación)
• ERA5-Land (Reanálisis del ECMWF)
• NLDAS (Sistema de asimilación de datos terrestres de América del Norte)

Cada conjunto de datos produce un archivo CSV independiente con la intensidad de la lluvia (mm/h) en todas las ubicaciones de los medidores.

Entradas climáticas/ETP

Para los cálculos de evapotranspiración, el script extrae:

• Temperatura máxima (Tmáx, °C)
• Temperatura mínima (Tmin, °C)
• Temperatura media (Tmed, °C)
• Velocidad del viento en 2m (U2, m/s)
• Humedad relativa (UR, %)
• Radiación solar (G, MJ/m²/día)

Los insumos climáticos se generan a partir de:- ERA5-Tierra

• NLDAS

---

Características clave

Manejo robusto de la unidad

El script incluye detección y conversión automática de unidades:

• Temperatura: Detecta automáticamente Kelvin frente a Celsius (evita valores bajo cero erróneos)
• Presión: Detecta automáticamente Pa frente a kPa para realizar cálculos precisos de la humedad.
• Precipitación: convierte todas las fuentes a mm/h consistentes
• Humedad relativa: Calculada correctamente a partir del punto de rocío (ERA5) o la humedad específica (NLDAS)

---

Red de ancho virtual

El guión crea una cuadrícula regular de pluviómetros virtuales en toda la cuenca:

• Densidad de cuadrícula definida por el usuario (predeterminada: 10×10 = 100 calibres)
• Medidores recortados al límite de la cuenca
• A cada medidor se le asignó una identificación única y coordenadas proyectadas.
• Interpolación espacial manejada por HydroPol2D durante la simulación

---

Seguro de calidad

Los gráficos de validación incorporados muestran:

• Comparación de la intensidad de las precipitaciones entre conjuntos de datos (primeros 30 días)
• Comparación de la velocidad del viento (primeros 30 días)
• Serie temporal en la ubicación del primer medidor

Estos ayudan a identificar lagunas de datos, valores atípicos o inconsistencias en los conjuntos de datos.

Estructura de salida

Formato CSV de lluvia

Cada archivo de lluvia contiene: Filas de encabezado:- Fila de título

• fila en blanco
• Números de índice de calibre
• Nombres de calibre
• Coordenadas este (m)
• Coordenadas norte (m)
• Nombres de variables

Filas de datos (una por día):- Columna 1: Días transcurridos

• Columna 2: Fecha (minutos desde el inicio)
• Columnas 3+: Intensidad de lluvia (mm/h) en cada pluviómetro

Estructura de ejemplo:

// ======================================================================
// HydroPol2D: Export Rainfall + ETP/Climate-Inputs as "sheet-like" CSVs
// FULLY UPDATED (CHIRPS + ERA5-Land + NLDAS) + ROBUST UNITS FIXES
//
// FIXES:
//  ✅ NLDAS temperature: auto-detect Kelvin vs Celsius (prevents ~ -260°C bug)
//  ✅ UR (%) computed properly for ERA5-Land + NLDAS (not fixed 50, not ~1e-70)
//     - ERA5: RH from T and dewpoint Td
//     - NLDAS: RH from specific humidity q, pressure p, and temperature T
//       * pressure auto-detect Pa vs kPa
//
// QA PLOTS:
//  ✅ Only first 30 days plotted (rain + wind), so charts stay responsive
//
// NOTES:
//  - Rainfall exports remain full period (CSV exports full time range)
//  - Charts are limited to first 30 days ONLY
// ======================================================================


// --------------------
// USER INPUTS
// --------------------
var scale_of_image = 1000;
var startyear = 2024;
var endyear   = 2025;
var Folder_Name = 'Pune';
var crsOut = 'EPSG:32643';   // output CRS for x/y in CSV

// --------------------
// ROI
// --------------------
var basins = ee.FeatureCollection("WWF/HydroSHEDS/v1/Basins/hybas_5");
var catchmentId = 7050585120;
var roi = basins.filter(ee.Filter.eq('HYBAS_ID', catchmentId));

// Adding a shapefile (you already have this variable defined in your session)
roi = pune;

var roiGeom = roi.geometry();

Map.addLayer(roi, {color: 'red'}, 'Catchment Area');
Map.centerObject(roi, 9);

// --------------------
// DATASETS
// --------------------

// CHIRPS (daily, precip in mm/day)
var rainfall_chirps = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY');

// ERA5-Land Daily Aggregates (daily)
var era5_land = ee.ImageCollection("ECMWF/ERA5_LAND/DAILY_AGGR")
  .filterDate(ee.Date.fromYMD(startyear, 1, 1), ee.Date.fromYMD(endyear, 12, 31))
  .select([
    'temperature_2m', 'temperature_2m_min', 'temperature_2m_max',
    'dewpoint_temperature_2m',
    'u_component_of_wind_10m', 'v_component_of_wind_10m',
    'total_precipitation_sum'
  ]);

// NLDAS (your band list variant)
var nldas = ee.ImageCollection("NASA/NLDAS/FORA0125_H002")
  .filterDate(ee.Date.fromYMD(startyear, 1, 1), ee.Date.fromYMD(endyear, 12, 31))
  .select([
    'temperature',          // K or °C depending on variant -> we auto-detect
    'specific_humidity',    // kg/kg
    'pressure',             // Pa or kPa depending on variant -> we auto-detect
    'wind_u',               // m/s
    'wind_v',               // m/s
    'total_precipitation'
  ]);


// --------------------
// GAUGE POINT GRID
// --------------------
var num_x = 10;
var num_y = 10;

// bounds coords: [[xmin,ymin], [xmax,ymin], [xmax,ymax], [xmin,ymax], [xmin,ymin]]
var ring = ee.List(roiGeom.bounds().coordinates().get(0));
var p0 = ee.List(ring.get(0));
var p2 = ee.List(ring.get(2));

var xmin = ee.Number(p0.get(0));
var ymin = ee.Number(p0.get(1));
var xmax = ee.Number(p2.get(0));
var ymax = ee.Number(p2.get(1));

var xStep = xmax.subtract(xmin).divide(ee.Number(num_x - 1));
var yStep = ymax.subtract(ymin).divide(ee.Number(num_y - 1));

var points = ee.List.sequence(0, num_x - 1).map(function(i) {
  i = ee.Number(i).int();
  return ee.List.sequence(0, num_y - 1).map(function(j) {
    j = ee.Number(j).int();
    var x = xmin.add(xStep.multiply(i));
    var y = ymin.add(yStep.multiply(j));
    return ee.Feature(ee.Geometry.Point([x, y]));
  });
}).flatten();

var gauge_points = ee.FeatureCollection(points).filterBounds(roiGeom);
Map.addLayer(gauge_points, {color: 'blue'}, 'Gauge Locations');


// --------------------
// Helpers
// --------------------
function nn(v, fallback) {
  return ee.Algorithms.If(ee.Algorithms.IsEqual(v, null), fallback, v);
}

function safeGet(list, i, fallback) {
  list = ee.List(list);
  i = ee.Number(i).int();
  var ok = i.lt(list.size());
  return ee.Algorithms.If(ok, nn(list.get(i), fallback), fallback);
}

function pad3(n) {
  n = ee.Number(n).int();
  return ee.Algorithms.If(
    n.lt(10),  ee.String('00').cat(n.format('%.0f')),
    ee.Algorithms.If(n.lt(100), ee.String('0').cat(n.format('%.0f')),
                     n.format('%.0f'))
  );
}

function colName(k) {
  k = ee.Number(k).int();
  return ee.String('C').cat(k.format('%.0f'));
}

function repeatN(s, N) {
  return ee.List.repeat(s, ee.Number(N).int());
}


// ======================================================================
// ROBUST UNIT HANDLING
// ======================================================================

// Return dict {Tk, Tc} from a temperature image that might be K or °C.
// Heuristic: values > 150 are Kelvin (safe globally for air temp).
function tempToTkTc(Traw) {
  Traw = ee.Image(Traw);
  var isKelvin = Traw.gt(150);
  var Tk = ee.Image(ee.Algorithms.If(isKelvin, Traw, Traw.add(273.15)));
  var Tc = Tk.subtract(273.15);
  return ee.Dictionary({Tk: Tk, Tc: Tc});
}

// Pressure to Pa (auto-detect): if p < 2000 assume kPa, else Pa.
function pressureToPa(p) {
  p = ee.Image(p);
  var isKPa = p.lt(2000);
  return ee.Image(ee.Algorithms.If(isKPa, p.multiply(1000), p));
}


// ======================================================================
// RH helpers
// ======================================================================

// Saturation vapor pressure (Pa) from temperature in °C (Tetens/Magnus)
function esPa_from_Tc(Tc) {
  Tc = ee.Image(Tc);
  return ee.Image(610.78).multiply(
    (ee.Image(17.2694).multiply(Tc).divide(Tc.add(237.3))).exp()
  );
}

// ERA5 RH (%): RH = 100 * es(Td) / es(T)
function RH_from_T_Td(Tk, Tdk) {
  var Tc  = ee.Image(Tk).subtract(273.15);
  var Tdc = ee.Image(Tdk).subtract(273.15);
  var esT  = esPa_from_Tc(Tc);
  var esTd = esPa_from_Tc(Tdc);
  return esTd.divide(esT).multiply(100).clamp(0, 100).rename('UR');
}

// NLDAS RH (%): compute actual vapor pressure from q + p
// e = (q*p) / (0.622 + 0.378 q)  (Pa)   q=kg/kg, p=Pa
function RH_from_q_p_T(q, p, Tk_or_Tc) {
  q = ee.Image(q);
  var pPa = pressureToPa(p);

  // handle T in either unit via tempToTkTc
  var tmp = tempToTkTc(Tk_or_Tc);
  var Tc = ee.Image(ee.Dictionary(tmp).get('Tc'));

  var e  = q.multiply(pPa).divide(ee.Image(0.622).add(ee.Image(0.378).multiply(q))); // Pa
  var es = esPa_from_Tc(Tc); // Pa
  return e.divide(es).multiply(100).clamp(0, 100).rename('UR');
}


// --------------------
// Attach gauge IDs + projected coords
// --------------------
var gpList  = gauge_points.toList(gauge_points.size());
var idxList = ee.List.sequence(0, gauge_points.size().subtract(1));

var gauge_points_id = ee.FeatureCollection(
  idxList.map(function(i) {
    i = ee.Number(i).int();
    var f = ee.Feature(gpList.get(i));

    var gid   = ee.String('A').cat(ee.String(pad3(i.add(101))));
    var gname = ee.String('Gauge ').cat(i.add(1).format('%.0f'));

    var pt = f.geometry().transform(crsOut, 1);
    var xy = pt.coordinates();

    return f.set({
      gauge_id: gid,
      gauge_name: gname,
      x: nn(xy.get(0), ''),
      y: nn(xy.get(1), '')
    });
  })
);

var N = ee.Number(gauge_points_id.size()).int();
print('N gauges:', N);
print('Example gauge:', gauge_points_id.first());

// Client-side count ONLY for building selectors (NOT inside map)
var Nclient = gauge_points_id.size().getInfo();

// These lists are used multiple times (build once)
var gaugeNames = ee.List(gauge_points_id.aggregate_array('gauge_name'));
var gids       = ee.List(gauge_points_id.aggregate_array('gauge_id'));
var xs         = ee.List(gauge_points_id.aggregate_array('x'));
var ys         = ee.List(gauge_points_id.aggregate_array('y'));


// --------------------
// DAILY CALENDAR
// --------------------
var startDate = ee.Date.fromYMD(startyear, 1, 1);
var endDate   = ee.Date.fromYMD(endyear, 12, 31).advance(1, 'day'); // exclusive
var nDays = endDate.difference(startDate, 'day');
var dayOffsets = ee.List.sequence(0, nDays.subtract(1));


// ======================================================================
// BUILD DAILY RAINFALL IMAGECOLLECTIONS (mm/h) FOR EACH SOURCE
// ======================================================================

// 1) CHIRPS: precip is mm/day -> mm/h (divide by 24)
var dailyRain_CHIRPS = ee.ImageCollection.fromImages(
  dayOffsets.map(function(d) {
    d = ee.Number(d).int();
    var date = startDate.advance(d, 'day');

    var rain = rainfall_chirps
      .filterDate(date, date.advance(1, 'day'))
      .sum()
      .divide(24)
      .clip(roiGeom)
      .rename('rain');

    return rain.set('system:time_start', date.millis());
  })
);

// 2) ERA5-Land: total_precipitation_sum is meters/day -> mm/h
var era5DummyRainMasked = ee.Image.constant(0).rename('rain').updateMask(ee.Image(0));

var dailyRain_ERA5 = ee.ImageCollection.fromImages(
  dayOffsets.map(function(d) {
    d = ee.Number(d).int();
    var date = startDate.advance(d, 'day');

    var era0 = era5_land.filterDate(date, date.advance(1, 'day')).first();

    var rain = ee.Image(
      ee.Algorithms.If(
        ee.Algorithms.IsEqual(era0, null),
        era5DummyRainMasked,
        ee.Image(era0).select('total_precipitation_sum')
          .multiply(1000)
          .divide(24)
          .rename('rain')
      )
    ).clip(roiGeom);

    return rain.set('system:time_start', date.millis());
  })
);

// 3) NLDAS: total_precipitation -> assume mm/day -> mm/h
var nldasDummyRainMasked = ee.Image.constant(0).rename('rain').updateMask(ee.Image(0));

var dailyRain_NLDAS = ee.ImageCollection.fromImages(
  dayOffsets.map(function(d) {
    d = ee.Number(d).int();
    var date = startDate.advance(d, 'day');

    var im0 = nldas.filterDate(date, date.advance(1, 'day')).first();

    var rain = ee.Image(
      ee.Algorithms.If(
        ee.Algorithms.IsEqual(im0, null),
        nldasDummyRainMasked,
        ee.Image(im0).select('total_precipitation')
          .divide(24)
          .rename('rain')
      )
    ).clip(roiGeom);

    return rain.set('system:time_start', date.millis());
  })
);


// ======================================================================
// HEADER ROW BUILDER (list only)
// ======================================================================
function makeHeaderRowList(c0, c1, valuesList) {
  valuesList = ee.List(valuesList);
  var base = ee.Dictionary({C0: c0, C1: c1});

  var out = ee.Dictionary(
    ee.List.sequence(0, N.subtract(1)).iterate(function(i, acc) {
      i = ee.Number(i).int();
      acc = ee.Dictionary(acc);
      var v = safeGet(valuesList, i, '');
      return acc.set(colName(i.add(2)), v);
    }, base)
  );

  return ee.Feature(null, out);
}


// ======================================================================
// RAINFALL SHEET BUILDER (MATCH TEMPLATE) - PARAMETERIZED
// ======================================================================
function buildRainfallSheet(dailyRainIC) {
  var rainRows = dailyRainIC.map(function(img) {
    img = ee.Image(img);

    var date = ee.Date(img.get('system:time_start'));
    var dateMin = date.difference(startDate, 'minute');

    // HydroPol2D convention you used
    var elapsedDays = ee.Number(dateMin).divide(1400);

    var fc = img.reduceRegions({
      collection: gauge_points_id,
      reducer: ee.Reducer.mean(),
      scale: scale_of_image
    });

    var vals = ee.List(fc.aggregate_array('mean'));

    var row = ee.Dictionary({ C0: elapsedDays, C1: dateMin });

    var row2 = ee.Dictionary(
      ee.List.sequence(0, N.subtract(1)).iterate(function(i, acc) {
        i = ee.Number(i).int();
        acc = ee.Dictionary(acc);
        var v = safeGet(vals, i, 0);
        return acc.set(colName(i.add(2)), v);
      }, row)
    );

    return ee.Feature(null, row2);
  });

  var header0 = ee.Feature(null, {
    C0: 'HydroPol2D - Spatial Rainfall Parameters: We assume the Durations Match with the Running Control Begin',
    C1: ''
  });

  var header1 = ee.Feature(null, {C0: '', C1: ''});

  var header2 = makeHeaderRowList('', 'Index', ee.List.sequence(1, N));
  var header3 = makeHeaderRowList('', 'Rain Gauge', gaugeNames);
  var header4 = makeHeaderRowList('', 'Easting [m]', xs);
  var header5 = makeHeaderRowList('', 'Northing [m]', ys);
  var header6 = makeHeaderRowList('Elapsed Days', 'Date (min)', repeatN('Rainfall Intensity (mm/h)', N));

  return ee.FeatureCollection([header0, header1, header2, header3, header4, header5, header6])
    .merge(rainRows);
}

// selectors: C0, C1, C2..C(N+1)
var selectorsRain = ['C0', 'C1'];
for (var i = 0; i < Nclient; i++) selectorsRain.push('C' + (i + 2));


// ======================================================================
// 1) RAINFALL SHEETS EXPORTS (FULL PERIOD)
// ======================================================================
Export.table.toDrive({
  collection: buildRainfallSheet(dailyRain_CHIRPS),
  description: 'Rainfall_HP2D_CHIRPS',
  folder: Folder_Name,
  fileFormat: 'CSV',
  selectors: selectorsRain
});

Export.table.toDrive({
  collection: buildRainfallSheet(dailyRain_ERA5),
  description: 'Rainfall_HP2D_ERA5',
  folder: Folder_Name,
  fileFormat: 'CSV',
  selectors: selectorsRain
});

Export.table.toDrive({
  collection: buildRainfallSheet(dailyRain_NLDAS),
  description: 'Rainfall_HP2D_NLDAS',
  folder: Folder_Name,
  fileFormat: 'CSV',
  selectors: selectorsRain
});


// ======================================================================
// 2) ETP / CLIMATE INPUTS SHEET (MATCH TEMPLATE)
// ======================================================================

function buildEtpHeaderFeatures() {
  var titleETP = ee.Feature(null, {C0: 'HydroPol2D - Spatial ETP Parameters', C1: ''});

  var metaRow = ee.Dictionary({C0:'', C1:''});
  metaRow = ee.Dictionary(
    ee.List.sequence(0, N.subtract(1)).iterate(function(i, acc){
      i = ee.Number(i).int();
      acc = ee.Dictionary(acc);
      var baseCol = i.multiply(6).add(2);

      acc = acc.set(colName(baseCol.add(0)), 'Index');
      acc = acc.set(colName(baseCol.add(1)), safeGet(gids, i, ''));
      acc = acc.set(colName(baseCol.add(2)), 'x easting (m)');
      acc = acc.set(colName(baseCol.add(3)), safeGet(xs, i, ''));
      acc = acc.set(colName(baseCol.add(4)), 'y northing (m)');
      acc = acc.set(colName(baseCol.add(5)), safeGet(ys, i, ''));
      return acc;
    }, metaRow)
  );
  var metaETP = ee.Feature(null, metaRow);

  var varRow = ee.Dictionary({C0:'', C1:'Date '});
  varRow = ee.Dictionary(
    ee.List.sequence(0, N.subtract(1)).iterate(function(i, acc){
      i = ee.Number(i).int();
      acc = ee.Dictionary(acc);
      var baseCol = i.multiply(6).add(2);

      acc = acc.set(colName(baseCol.add(0)), 'Tmax (oC)');
      acc = acc.set(colName(baseCol.add(1)), 'Tmin (oC)');
      acc = acc.set(colName(baseCol.add(2)), 'Tmed (oC)');
      acc = acc.set(colName(baseCol.add(3)), 'U2 [m/s]');
      acc = acc.set(colName(baseCol.add(4)), 'UR (%)');
      acc = acc.set(colName(baseCol.add(5)), 'G (MJ/(m2.dia))');
      return acc;
    }, varRow)
  );
  var varsETP = ee.Feature(null, varRow);

  return ee.FeatureCollection([titleETP, metaETP, varsETP]);
}

var selectorsETP = ['C0', 'C1'];
for (var j = 0; j < Nclient * 6; j++) selectorsETP.push('C' + (j + 2));


// --------------------
// ERA5 ETP/CLIMATE ROWS (daily)
// --------------------
var eraDummyMasked = ee.Image.constant([0, 0, 0, 0, 0, 0])
  .rename([
    'temperature_2m', 'temperature_2m_min', 'temperature_2m_max',
    'dewpoint_temperature_2m',
    'u_component_of_wind_10m', 'v_component_of_wind_10m'
  ])
  .updateMask(ee.Image(0));

var etpRows_ERA5 = dailyRain_CHIRPS.map(function(dummyImg) {
  var date = ee.Date(ee.Image(dummyImg).get('system:time_start'));
  var era0 = era5_land.filterDate(date, date.advance(1, 'day')).first();

  var era = ee.Image(
    ee.Algorithms.If(ee.Algorithms.IsEqual(era0, null), eraDummyMasked, era0)
  );

  var tmed = era.select('temperature_2m').subtract(273.15).rename('Tmed');
  var tmax = era.select('temperature_2m_max').subtract(273.15).rename('Tmax');
  var tmin = era.select('temperature_2m_min').subtract(273.15).rename('Tmin');

  var U2  = era.select('u_component_of_wind_10m').pow(2)
    .add(era.select('v_component_of_wind_10m').pow(2))
    .sqrt()
    .rename('U2');

  var UR = RH_from_T_Td(era.select('temperature_2m'), era.select('dewpoint_temperature_2m'));
  var G  = ee.Image.constant(0).rename('G');

  var stack = ee.Image.cat([tmax, tmin, tmed, U2, UR, G]);

  var fc = stack.reduceRegions({
    collection: gauge_points_id,
    reducer: ee.Reducer.mean(),
    scale: scale_of_image
  });

  var tmaxL = ee.List(fc.aggregate_array('Tmax'));
  var tminL = ee.List(fc.aggregate_array('Tmin'));
  var tmedL = ee.List(fc.aggregate_array('Tmed'));
  var u2L   = ee.List(fc.aggregate_array('U2'));
  var urL   = ee.List(fc.aggregate_array('UR'));
  var gL    = ee.List(fc.aggregate_array('G'));

  var row = ee.Dictionary({ C0: '', C1: date.format('yyyy-MM-dd HH:mm:ss') });

  var row2 = ee.Dictionary(
    ee.List.sequence(0, N.subtract(1)).iterate(function(i, acc){
      i = ee.Number(i).int();
      acc = ee.Dictionary(acc);
      var baseCol = i.multiply(6).add(2);

      acc = acc.set(colName(baseCol.add(0)), safeGet(tmaxL, i, ''));
      acc = acc.set(colName(baseCol.add(1)), safeGet(tminL, i, ''));
      acc = acc.set(colName(baseCol.add(2)), safeGet(tmedL, i, ''));
      acc = acc.set(colName(baseCol.add(3)), safeGet(u2L,   i, ''));
      acc = acc.set(colName(baseCol.add(4)), safeGet(urL,   i, ''));
      acc = acc.set(colName(baseCol.add(5)), safeGet(gL,    i, ''));
      return acc;
    }, row)
  );

  return ee.Feature(null, row2);
});

Export.table.toDrive({
  collection: buildEtpHeaderFeatures().merge(etpRows_ERA5),
  description: 'ETP_Inputs_HP2D_ERA5',
  folder: Folder_Name,
  fileFormat: 'CSV',
  selectors: selectorsETP
});


// --------------------
// NLDAS ETP/CLIMATE ROWS (daily)  [ROBUST FIX]
// --------------------
var nldasDummyEtpMasked = ee.Image.constant([0, 0, 0, 0, 0, 0])
  .rename(['Tmax', 'Tmin', 'Tmed', 'U2', 'UR', 'G'])
  .updateMask(ee.Image(0));

var etpRows_NLDAS = dailyRain_CHIRPS.map(function(dummyImg) {
  var date = ee.Date(ee.Image(dummyImg).get('system:time_start'));
  var im0 = nldas.filterDate(date, date.advance(1, 'day')).first();

  var daily = ee.Image(
    ee.Algorithms.If(
      ee.Algorithms.IsEqual(im0, null),
      nldasDummyEtpMasked,
      (function(){
        var im = ee.Image(im0);

        // robust temperature handling
        var tmp = tempToTkTc(im.select('temperature'));
        var Tc = ee.Image(ee.Dictionary(tmp).get('Tc'));
        // NLDAS daily: using same Tc for Tmax/Tmin/Tmed (as before)
        var Tmax = Tc.rename('Tmax');
        var Tmin = Tc.rename('Tmin');
        var Tmed = Tc.rename('Tmed');

        var U2 = im.select('wind_u').pow(2)
          .add(im.select('wind_v').pow(2))
          .sqrt()
          .rename('U2');

        var UR = RH_from_q_p_T(
          im.select('specific_humidity'),
          im.select('pressure'),
          im.select('temperature') // can be K or °C, RH function handles it
        ).rename('UR');

        var G = ee.Image.constant(0).rename('G');

        return ee.Image.cat([Tmax, Tmin, Tmed, U2, UR, G]);
      })()
    )
  ).clip(roiGeom);

  var fc = daily.reduceRegions({
    collection: gauge_points_id,
    reducer: ee.Reducer.mean(),
    scale: scale_of_image
  });

  var tmaxL = ee.List(fc.aggregate_array('Tmax'));
  var tminL = ee.List(fc.aggregate_array('Tmin'));
  var tmedL = ee.List(fc.aggregate_array('Tmed'));
  var u2L   = ee.List(fc.aggregate_array('U2'));
  var urL   = ee.List(fc.aggregate_array('UR'));
  var gL    = ee.List(fc.aggregate_array('G'));

  var row = ee.Dictionary({ C0: '', C1: date.format('yyyy-MM-dd HH:mm:ss') });

  var row2 = ee.Dictionary(
    ee.List.sequence(0, N.subtract(1)).iterate(function(i, acc){
      i = ee.Number(i).int();
      acc = ee.Dictionary(acc);
      var baseCol = i.multiply(6).add(2);

      acc = acc.set(colName(baseCol.add(0)), safeGet(tmaxL, i, ''));
      acc = acc.set(colName(baseCol.add(1)), safeGet(tminL, i, ''));
      acc = acc.set(colName(baseCol.add(2)), safeGet(tmedL, i, ''));
      acc = acc.set(colName(baseCol.add(3)), safeGet(u2L,   i, ''));
      acc = acc.set(colName(baseCol.add(4)), safeGet(urL,   i, ''));
      acc = acc.set(colName(baseCol.add(5)), safeGet(gL,    i, ''));
      return acc;
    }, row)
  );

  return ee.Feature(null, row2);
});

Export.table.toDrive({
  collection: buildEtpHeaderFeatures().merge(etpRows_NLDAS),
  description: 'ETP_Inputs_HP2D_NLDAS',
  folder: Folder_Name,
  fileFormat: 'CSV',
  selectors: selectorsETP
});


// ======================================================================
// QA PLOTS (FIRST 30 DAYS ONLY)
// ======================================================================

// Pick first gauge
var oneGauge = ee.Feature(gauge_points_id.first());
var onePt = oneGauge.geometry();
Map.addLayer(onePt, {color: 'yellow'}, 'Gauge (first)');
print('Selected gauge:', oneGauge);

// Limit IC to first 30 days for plotting ONLY
var plotStart = startDate;
var plotEnd   = startDate.advance(30, 'day'); // first 30 days
function clipTo30Days(ic) {
  return ic.filterDate(plotStart, plotEnd);
}

// Helper: IC -> point TS FC
function icToPointTS(ic, bandName, datasetName) {
  return ee.FeatureCollection(
    ic.map(function(img) {
      img = ee.Image(img);
      var date = ee.Date(img.get('system:time_start'));

      var val = img.select(bandName).reduceRegion({
        reducer: ee.Reducer.mean(),
        geometry: onePt,
        scale: scale_of_image,
        bestEffort: true,
        maxPixels: 1e13
      }).get(bandName);

      return ee.Feature(null, {
        time: date.millis(),
        date: date.format('YYYY-MM-dd'),
        value: val,
        dataset: datasetName
      });
    })
  );
}

// Rainfall charts (30 days)
var tsRain_CHIRPS = icToPointTS(clipTo30Days(dailyRain_CHIRPS), 'rain', 'CHIRPS');
var tsRain_ERA5   = icToPointTS(clipTo30Days(dailyRain_ERA5),   'rain', 'ERA5-Land');
var tsRain_NLDAS  = icToPointTS(clipTo30Days(dailyRain_NLDAS),  'rain', 'NLDAS');

var tsRainAll = tsRain_CHIRPS.merge(tsRain_ERA5).merge(tsRain_NLDAS);

print(
  ui.Chart.feature.groups({
    features: tsRainAll,
    xProperty: 'time',
    yProperty: 'value',
    seriesProperty: 'dataset'
  })
  .setChartType('LineChart')
  .setOptions({
    title: 'Rainfall intensity at first gauge (mm/h) - first 30 days',
    hAxis: { title: 'Date' },
    vAxis: { title: 'Rain (mm/h)' },
    lineWidth: 1,
    pointSize: 0,
    interpolateNulls: false
  })
);

// Wind charts (30 days)
var eraU2DummyMasked = ee.Image.constant(0).rename('U2').updateMask(ee.Image(0));
var nldasU2DummyMasked = ee.Image.constant(0).rename('U2').updateMask(ee.Image(0));

var dailyWind_ERA5 = ee.ImageCollection.fromImages(
  dayOffsets.map(function(d) {
    d = ee.Number(d).int();
    var date = startDate.advance(d, 'day');

    var era0 = era5_land.filterDate(date, date.advance(1, 'day')).first();

    var U2 = ee.Image(
      ee.Algorithms.If(
        ee.Algorithms.IsEqual(era0, null),
        eraU2DummyMasked,
        ee.Image(era0)
          .select(['u_component_of_wind_10m', 'v_component_of_wind_10m'])
          .pow(2)
          .reduce('sum')
          .sqrt()
          .rename('U2')
      )
    ).clip(roiGeom);

    return U2.set('system:time_start', date.millis());
  })
);

var dailyWind_NLDAS = ee.ImageCollection.fromImages(
  dayOffsets.map(function(d) {
    d = ee.Number(d).int();
    var date = startDate.advance(d, 'day');

    var im0 = nldas.filterDate(date, date.advance(1, 'day')).first();

    var U2 = ee.Image(
      ee.Algorithms.If(
        ee.Algorithms.IsEqual(im0, null),
        nldasU2DummyMasked,
        ee.Image(im0).select('wind_u').pow(2)
          .add(ee.Image(im0).select('wind_v').pow(2))
          .sqrt()
          .rename('U2')
      )
    ).clip(roiGeom);

    return U2.set('system:time_start', date.millis());
  })
);

var tsWind_ERA5  = icToPointTS(clipTo30Days(dailyWind_ERA5),  'U2', 'ERA5-Land');
var tsWind_NLDAS = icToPointTS(clipTo30Days(dailyWind_NLDAS), 'U2', 'NLDAS');

var tsWindAll = tsWind_ERA5.merge(tsWind_NLDAS);

print(
  ui.Chart.feature.groups({
    features: tsWindAll,
    xProperty: 'time',
    yProperty: 'value',
    seriesProperty: 'dataset'
  })
  .setChartType('LineChart')
  .setOptions({
    title: 'Wind speed at first gauge (m/s) - first 30 days',
    hAxis: { title: 'Date' },
    vAxis: { title: 'Wind (m/s)' },
    lineWidth: 1,
    pointSize: 0,
    interpolateNulls: false
  })
);