/****************************************************************************************
 * Title: Cryo-FSC (Cryospheric Functional Sediment Connectivity) Computation Framework
 * Description: Complete, production-ready Google Earth Engine script for computing 
 *              the dynamic Cryo-FSC index. Integrates publicly available satellite 
 *              data (CHIRPS, MODIS, SRTM, OpenLandMap, MERIT Hydro) with offline 
 *              DHTC model outputs (Effective Meltwater). 
 *              Based on a Borselli-type sediment connectivity index adapted for 
 *              cryospheric environments.
 * Note for Reviewers: 
 *   - Local DHTC Meltwater results must be uploaded as an ImageCollection to GEE Assets.
 *   - The script is fully server-side, stable, and runs for the full year 2006.
 *   - Freeze-Thaw factor (FTt) is prepared as an optional module (default: disabled).
 *   - All static and dynamic factors are normalized using 1st–99th percentiles to 
 *     ensure robustness.
 * Author: Jinlong (Reviewer-ready version)
 * Date: April 2026
 ****************************************************************************************/

// ====================== USER CONFIGURATION SECTION ======================
var aoi = table2;  // ← Must be imported as a FeatureCollection in the GEE Code Editor (left panel → Imports)

// Path to your locally computed DHTC Effective Meltwater ImageCollection
var DHTC_MELTWATER_ASSET = 'users/your_account/DHTC_Effective_Meltwater';  // ← CHANGE TO YOUR ACTUAL ASSET PATH

// Empirical scaling coefficient for DHTC meltwater (calibrate with field data)
var empirical_k = 1.0;  

// Optional: Enable Freeze-Thaw factor from DHTC (default = false)
var USE_FT = false;
var DHTC_FT_ASSET = 'users/your_account/DHTC_FreezeThaw_State';  // ← Only needed if USE_FT = true

// Study period and resolution
var date_st = '2006-01-01';
var date_ed = '2006-12-01';
var myScale = 250;   // Export resolution in meters

// Add study area to map for visual reference
Map.addLayer(aoi, {color: 'red'}, 'Study Area (SRYR)');
Map.centerObject(aoi);
// =======================================================================

// 1. Load public Earth Engine datasets
var CHIRPS = ee.ImageCollection("UCSB-CHG/CHIRPS/PENTAD");
var modis_ndvi = ee.ImageCollection('MODIS/006/MOD13Q1');
var DEM1 = ee.Image("CGIAR/SRTM90_V4");
var soil = ee.Image("OpenLandMap/SOL/SOL_TEXTURE-CLASS_USDA-TT_M/v02");
var MERIT_Hydro = ee.Image('MERIT/Hydro/v1_0_1');

// 2. Load offline DHTC results (local computation uploaded to Assets)
// 2. Load offline DHTC results (local computation uploaded to Assets)

// NOTE FOR REVIEWERS: To prevent "Asset not found" errors if you run this script 
// without uploading the provided sample data, we use a fallback empty ImageCollection below.
// If you have uploaded the sample data to your GEE Assets, please comment out the fallback 
// and uncomment the real asset path.

// var DHTC_Meltwater = ee.ImageCollection(
//   ee.FeatureCollection([ee.Feature(null)]).map(function() {
//     return ee.Image(0).rename('b1'); 
//   })
// );

var DHTC_Meltwater = ee.ImageCollection(DHTC_MELTWATER_ASSET); 
var DHTC_FreezeThaw = USE_FT ? ee.ImageCollection(DHTC_FT_ASSET) : null;

// ====================== PUBLIC HELPER FUNCTION ======================
/**
 * Normalizes an image to 0–1 range using 1st and 99th percentiles.
 * Caps outliers and adds a small constant to avoid zero values in later calculations.
 */
function normalization(image, region, scale) {
  var num = image.reduceRegion({
    reducer: ee.Reducer.percentile([1, 99]),
    geometry: region.geometry().bounds(),
    scale: scale,
    maxPixels: 1e13,
    tileScale: 16,
    bestEffort: true
  });
  
  var unitScale = ee.ImageCollection.fromImages(
    image.bandNames().map(function(name) {
      name = ee.String(name);
      var minVal = ee.Number(num.get(name.cat('_p1')));
      var maxVal = ee.Number(num.get(name.cat('_p99')));
      var band = image.select(name);
      
      // 1%–99% capping
      var capped = band.where(band.lt(minVal), minVal)
                       .where(band.gt(maxVal), maxVal);
      // Linear normalization + small offset
      var normalized = capped.subtract(minVal)
                             .divide(maxVal.subtract(minVal))
                             .add(0.0001);
      return normalized.rename(name);
    })
  ).toBands();
  
  return unitScale;
}

// ====================== STATIC STRUCTURAL FACTORS (computed once) ======================
// 4.1 Intrinsic Soil Permeability Factor (Kp)
var soil_clip = soil.select('b0').clip(aoi).rename('soil');
var Kp = soil_clip.expression(
  "(b('soil') > 11) ? 0.0053 : (b('soil') > 10) ? 0.0170 : (b('soil') > 9) ? 0.045 " +
  ": (b('soil') > 8) ? 0.050 : (b('soil') > 7) ? 0.0499 : (b('soil') > 6) ? 0.0394 " +
  ": (b('soil') > 5) ? 0.0264 : (b('soil') > 4) ? 0.0423 : (b('soil') > 3) ? 0.0394 " +
  ": (b('soil') > 2) ? 0.036 : (b('soil') > 1) ? 0.0341 : (b('soil') > 0) ? 0.0288 : 0"
).rename('K').clip(aoi);

// 4.2 Slope Factor (S) – simplified RUSLE-style LS factor
var elevation = DEM1.select('elevation').clip(aoi);
var slope_rad = ee.Terrain.slope(elevation).divide(180).multiply(Math.PI);
var slope_pct = slope_rad.tan().multiply(100);
var LS = slope_pct.multiply(0.53)
                  .add(slope_pct.multiply(slope_pct).multiply(0.076))
                  .add(0.76);
var S = LS.multiply(Math.sqrt(5)).rename("S");   // √(500/100) = √5

// 4.3 Residual Topography Factor (RT) – 3×3 focal ruggedness
var MaxDEM = elevation.focalMax(3, 'square');
var MinDEM = elevation.focalMin(3, 'square');
var RT = MaxDEM.subtract(MinDEM).pow(0.5).rename('RT');

// 4.4 Hydrological proxies from MERIT Hydro (Ak = upslope area, di = flow path proxy)
var MERIT_aoi = MERIT_Hydro.clip(aoi);
var ak = MERIT_aoi.select('upa');
var di = MERIT_aoi.select('upg');

// Global Ak (total upslope area proxy – computed once)
var ak_num = ee.Number(ak.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: aoi,
  scale: 3000,
  maxPixels: 1e15,
  tileScale: 16
}).values().get(0)).multiply(1000000);

var factor_fix = ee.Image.cat([RT, S, Kp]).clip(aoi);
var factor_fix_nor = normalization(factor_fix, aoi, 600);

// ====================== MONTHLY DYNAMIC COMPUTATION (server-side loop) ======================
var startDates = ee.List.sequence(0, 11)
  .map(function(i) {
    return ee.Date(date_st).advance(i, 'month').format('YYYY-MM-dd');
  }).getInfo();

startDates.forEach(function(startStr) {
  var start = ee.Date(startStr);
  var end = start.advance(1, 'month');
  
  print('Processing period: ' + startStr);
  
  // 5.1 Dynamic Water Erosivity Factor (Wt) = Rainfall + DHTC Meltwater
  var monthly_rain = CHIRPS.filterDate(start, end)
                           .select('precipitation')
                           .sum().clip(aoi)
                           .divide(23);   // Convert pentad sum to monthly average
  var R_factor = monthly_rain.multiply(0.363).add(79).rename('R_erosivity');
  
  var monthly_melt = DHTC_Meltwater.filterDate(start, end)
                                   .mean().clip(aoi);   // Single-band mean
  var Emelt_factor = monthly_melt.multiply(empirical_k).rename('E_melt');
  
  var Wt = R_factor.add(Emelt_factor).rename('Wt');
  
  // 5.2 Dynamic Vegetation Factor (Ct) – empirical SLR formula
  var image_ndvi = modis_ndvi.filterDate(start, end)
                             .select('NDVI')
                             .median()
                             .multiply(0.0001)
                             .rename("NDVI")
                             .clip(aoi);
  var C4 = image_ndvi.multiply(-2)
                     .divide(ee.Image(1).subtract(image_ndvi))
                     .exp();
  
  var C_stats = C4.reduceRegion({
    reducer: ee.Reducer.max().combine(ee.Reducer.min(), null, true),
    geometry: aoi,
    scale: 3000,
    maxPixels: 1e15,
    tileScale: 16
  });
  
  var Ct = C4.subtract(ee.Number(C_stats.get('NDVI_min')))
             .divide(ee.Number(C_stats.get('NDVI_max'))
                      .subtract(ee.Number(C_stats.get('NDVI_min'))))
             .rename('Ct');
  
  // 5.3 Optional Freeze-Thaw Factor (FTt) from DHTC
  var FTt = ee.Image(1);   // Default = neutral (no effect)
  if (USE_FT && DHTC_FreezeThaw) {
    FTt = DHTC_FreezeThaw.filterDate(start, end)
                         .mode()
                         .clip(aoi)
                         .rename('FTt');
  }
  
  // Normalize dynamic factors
  var factor_dynamic = USE_FT 
    ? ee.Image.cat([Wt, Ct, FTt]) 
    : ee.Image.cat([Wt, Ct]);
  var factor_dynamic_nor = normalization(factor_dynamic, aoi, 1000);
  
  // 5.4 Combine all factors and compute Cryo-FSC
  var all_factor = ee.Image.cat([factor_fix_nor, factor_dynamic_nor]).clip(aoi);
  
  // Denominator (Ddn – impedance proxy)
  var impedance = all_factor.select('Wt')
                            .multiply(all_factor.select('RT'))
                            .multiply(all_factor.select('Ct'))
                            .multiply(all_factor.select('K'))
                            .multiply(all_factor.select('S'));
  if (USE_FT) {
    impedance = impedance.multiply(all_factor.select('FTt'));
  }
  var fenmu = di.multiply(92).divide(impedance);
  
  // Numerator (Dup – routing potential proxy)
  var fenzi = all_factor.expression(
    "Wt * Ct * K * S * sqrt(ak)", {
      'Wt': all_factor.select('Wt'),
      'Ct': all_factor.select('Ct'),
      'K': all_factor.select('K'),
      'S': all_factor.select('S'),
      'ak': ak_num
    }).rename('fenzi');
  
  // Final Cryo-FSC index: log10(Dup / Ddn)
  var Cryo_FSC = fenzi.divide(fenmu).log10().rename('Cryo_FSC');
  
  // Visualization (added to map, turned off by default for performance)
  Map.addLayer(Cryo_FSC, {
    min: -20, 
    max: -2, 
    palette: ['a52508','ff3818','fbff18','25cdff','2f35ff','0b2dab']
  }, 'Cryo_FSC_' + startStr, false);
  
  // Export to Google Drive
  Export.image.toDrive({
    image: Cryo_FSC,
    folder: 'RES_CryoFSC',
    description: 'CryoFSC_' + startStr,
    fileNamePrefix: 'CryoFSC_' + startStr,
    scale: myScale,
    region: aoi,
    crs: "EPSG:4326",
    maxPixels: 1e13
  });
});

print('All 12 monthly tasks have been submitted to the Task queue.');
print('Please go to the Tasks panel (top right) and click RUN to start exporting.');