Gerando Força Meteorológica com Google Earth Engine

HydroPol2D requer dados de forçamento meteorológico distribuídos espacialmente para conduzir simulações hidrológicas. Este script Google Earth Engine (GEE) automatiza a extração e formatação de entradas de precipitação e evapotranspiração (ETP) de vários conjuntos de dados globais. O fluxo de trabalho gera arquivos CSV no formato exato exigido pelo HydroPol2D, com dados extraídos em locais de medição virtuais distribuídos por todo o domínio de modelagem.

O que o fluxo de trabalho faz

O script processa três conjuntos de dados meteorológicos globais e exporta:

Entradas de chuva

• CHIRPS (Precipitação infravermelha do grupo de riscos climáticos com dados da estação)
• ERA5-Land (Reanálise ECMWF)
• NLDAS (Sistema de Assimilação de Dados Terrestres da América do Norte)

Cada conjunto de dados produz um arquivo CSV separado com a intensidade da precipitação (mm/h) em todos os locais do medidor.

Insumos climáticos/ETP

Para cálculos de evapotranspiração, o script extrai:

• Temperatura máxima (Tmax, °C)
• Temperatura mínima (Tmin, °C)
• Temperatura média (Tmed, °C)
• Velocidade do vento em 2m (U2, m/s)
• Umidade relativa (UR, %)
• Radiação solar (G, MJ/m²/dia)

Os insumos climáticos são gerados a partir de:- ERA5-Terreno

• NLDAS

---

Principais recursos

Manuseio Robusto de Unidades

O script inclui detecção e conversão automática de unidades:

• Temperatura: detecta automaticamente Kelvin vs Celsius (evita valores abaixo de zero errados)
• Pressão: detecta automaticamente Pa vs kPa para cálculos precisos de umidade
• Precipitação: converte todas as fontes em mm/h consistentes
• Umidade relativa: calculada corretamente a partir do ponto de orvalho (ERA5) ou umidade específica (NLDAS)

---

Rede de medidores virtuais

O script cria uma grade regular de pluviômetros virtuais em toda a bacia hidrográfica:

• Densidade de grade definida pelo usuário (padrão: 10×10 = 100 medidores)
• Medidores cortados no limite da bacia hidrográfica
• Cada medidor recebeu um ID exclusivo e coordenadas projetadas
• Interpolação espacial tratada por HydroPol2D durante a simulação

---

Garantia de qualidade

Os gráficos de validação integrados mostram:

• Comparação da intensidade da precipitação entre conjuntos de dados (primeiros 30 dias)
• Comparação da velocidade do vento (primeiros 30 dias)
• Série temporal no primeiro local do medidor

Isso ajuda a identificar lacunas de dados, valores discrepantes ou inconsistências nos conjuntos de dados.

Estrutura de saída

Formato CSV de chuva

Cada arquivo de precipitação contém: Linhas de cabeçalho:- Linha de título

• Linha em branco
• Números de índice do medidor
• Nomes de medidores
• Coordenadas leste (m)
• Coordenadas norte (m)
• Nomes de variáveis

Linhas de dados (uma por dia):- Coluna 1: Dias decorridos

• Coluna 2: Data (minutos desde o início)
• Colunas 3+: Intensidade da precipitação (mm/h) em cada medidor

Exemplo de estrutura:

// ======================================================================
// 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
  })
);