Google Earth Engine Code

// ============================================================================
// PCA SNOW INDEX - NDSI-LIKE VISUALIZATION (FINAL FIX)
// Resolves: Map.addLayer encoding errors & server/client type mismatches
// ============================================================================

var aoi = table;
Map.centerObject(aoi, 10);
Map.addLayer(aoi.style({color: "red", fillColor: "00000000"}), {}, "Study Area");

// ===== 1. SETUP =====
var start = '2022-11-01';
var end   = '2023-03-30';
var scale = 20;

// ===== 2. SENTINEL-2 =====
var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(aoi)
  .filterDate(start, end)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 50));

function prepS2(img) {
  var refl = img.select(['B2','B3','B4','B8','B11']).divide(10000);
  var ndsi = refl.normalizedDifference(['B3', 'B11']).rename('NDSI');
  var ndvi = refl.normalizedDifference(['B8', 'B4']).rename('NDVI');
  var ndwi = refl.normalizedDifference(['B3', 'B8']).rename('NDWI');
  return refl.addBands([ndsi, ndvi, ndwi]).copyProperties(img, img.propertyNames());
}

var s2_prep = s2.map(prepS2);
var s2_comp = s2_prep.mean().clip(aoi);
Map.addLayer(s2_comp.select(['B4', 'B3', 'B2']), {min:0, max:0.3}, "Sentinel-2 RGB");

// ===== 3. SENTINEL-1 =====
var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(aoi)
  .filterDate(start, end)
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  .select(['VV','VH']);

function toDB_safe(img) {
  var vv = img.select('VV').where(img.select('VV').lte(0), 1e-6);
  var vh = img.select('VH').where(img.select('VH').lte(0), 1e-6);
  return ee.Image([
    vv.log10().multiply(10).rename('VV_dB'),
    vh.log10().multiply(10).rename('VH_dB')
  ]).copyProperties(img, img.propertyNames());
}

var s1_db = s1.map(toDB_safe);
var s1_comp = s1_db.mean().clip(aoi);

// ===== 4. STACK BANDS =====
var stack = s2_comp.select(['B2','B3','B4','B8','B11','NDSI','NDVI','NDWI'])
  .addBands(s1_comp.select(['VV_dB','VH_dB']))
  .clip(aoi);

var nBands = stack.bandNames().length();

// ===== 5. PCA FUNCTION (OFFICIAL GEE PATTERN) =====
function getPrincipalComponents(image, scale, region) {
  var bandNames = image.bandNames();
  
  var means = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9,
    bestEffort: true
  });
  var centered = image.subtract(ee.Image.constant(bandNames.map(function(b) {
    return ee.Number(means.get(b));
  })).rename(bandNames));
  
  var arrayImage = centered.toArray();
  var arrayImage2D = arrayImage.toArray(1);
  
  var covar = arrayImage.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9,
    bestEffort: true
  });
  var covarArray = ee.Array(covar.get('array'));
  var eigens = covarArray.eigen();
  var eigenVectors = eigens.slice(1, 1);
  
  // [bands × components] * [bands × 1] = [components × 1]
  return {
    pcImage: ee.Image(eigenVectors)
      .matrixMultiply(arrayImage2D)
      .arrayProject([0])
      .arrayFlatten([ee.List.sequence(1, nBands).map(function(i) {
        return ee.String('pc').cat(ee.Number(i).format('%d'));
      })])
      .clip(region),
    eigenVectors: eigenVectors
  };
}

// ===== 6. RUN PCA =====
var pcaResult = getPrincipalComponents(stack, scale, aoi);
var pcImage = pcaResult.pcImage;

// ===== 7. FIND BEST SNOW PC =====
var ndsi = s2_comp.select('NDSI');
var pcNames = ee.List.sequence(1, nBands).map(function(i) {
  return ee.String('pc').cat(ee.Number(i).format('%d'));
});

var corrValues = pcNames.map(function(name) {
  var r = pcImage.select([name]).addBands(ndsi).reduceRegion({
    reducer: ee.Reducer.pearsonsCorrelation(),
    geometry: aoi,
    scale: scale,
    maxPixels: 1e9,
    bestEffort: true
  }).get('correlation');
  return ee.Algorithms.If(r, ee.Number(r).abs(), ee.Number(0));
});

var maxCorr = ee.Number(corrValues.reduce(ee.Reducer.max()));
var bestIdx = corrValues.indexOf(maxCorr);
var bestPCName = ee.String(pcNames.get(bestIdx));

print('🎯 Best snow-correlated PC:', bestPCName);
print('📊 Max |correlation| with NDSI:', maxCorr.getInfo());

// ===== 8. CREATE & FLIP PCA SNOW INDEX =====
var pcaSnowIndex = pcImage.select([bestPCName]).rename('PCA_Snow_Index');

// Get raw correlation as JS number for safe client-side branching
var rawCorrVal = ndsi.addBands(pcImage.select([bestPCName])).reduceRegion({
  reducer: ee.Reducer.pearsonsCorrelation(),
  geometry: aoi,
  scale: scale,
  maxPixels: 1e9,
  bestEffort: true
}).get('correlation').getInfo();

// FIX: Use client-side if/else to avoid ee.Algorithms.If encoding issues
if (rawCorrVal < 0) {
  pcaSnowIndex = pcaSnowIndex.multiply(-1).rename('PCA_Snow_Index');
}

// ===== 9. AUTO-VISUALIZATION (CLIENT-SIDE NUMBERS) =====
var stats = pcaSnowIndex.reduceRegion({
  reducer: ee.Reducer.percentile([1, 99]),
  geometry: aoi,
  scale: scale,
  maxPixels: 1e9,
  bestEffort: true
});

// CRITICAL: Convert ee.Number to JS numbers for Map.addLayer
var p1 = ee.Number(stats.get('PCA_Snow_Index_p1')).getInfo() || -1;
var p99 = ee.Number(stats.get('PCA_Snow_Index_p99')).getInfo() || 1;

var snowVis = {
  min: p1,
  max: p99,
  palette: ['brown', 'orange', 'yellow', 'white', 'cyan', 'blue']
};

// ===== 10. DISPLAY =====
Map.addLayer(pcaSnowIndex, snowVis, '🌨️ PCA Snow Index (NDSI-like)');

Map.addLayer(ndsi, {min: -0.2, max: 1.0, palette: ['brown', 'white', 'blue']}, '🔍 Traditional NDSI');

// Binary mask
var meanStd = pcaSnowIndex.reduceRegion({
  reducer: ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true),
  geometry: aoi,
  scale: scale,
  maxPixels: 1e9,
  bestEffort: true
});
var threshold = ee.Number(meanStd.get('PCA_Snow_Index_mean'))
  .add(ee.Number(meanStd.get('PCA_Snow_Index_stdDev')).multiply(0.5));
var snowMask = pcaSnowIndex.gt(threshold);
Map.addLayer(snowMask.updateMask(snowMask), {palette: 'blue'}, '❄️ Snow Cover');

// ===== 11. WET/DRY CLASSIFICATION =====
var wetIndicators = ee.Image.cat([
  s2_comp.select('NDWI').gt(0.05),
  s1_comp.select('VV_dB').lt(-10),
  s1_comp.select('VH_dB').divide(s1_comp.select('VV_dB').add(0.001)).gt(0.4)
]).reduce(ee.Reducer.sum());

var snowClass = ee.Image(1)
  .where(snowMask.and(wetIndicators.gte(2)), 2)
  .where(snowMask.not(), 0)
  .rename('Snow_Class');

Map.addLayer(snowClass, {min:0, max:2, palette:['gray','white','cyan']}, '🌧️ Wet/Dry Snow');

// ===== 12. EXPORT =====
Export.image.toDrive({
  image: pcaSnowIndex,
  description: 'PCA_Snow_Index_' + start.replace(/-/g, ''),
  folder: "Hohe-Tauren Snow Study",
  region: aoi,
  scale: 10,
  crs: "EPSG:4326",
  maxPixels: 1e13,
  fileFormat: 'GeoTIFF'
});

Export.image.toDrive({
  image: snowClass,
  description: 'Snow_Class_' + start.replace(/-/g, ''),
  folder: "Hohe-Tauren Snow Study",
  region: aoi,
  scale: 10,
  crs: "EPSG:4326",
  maxPixels: 1e13,
  fileFormat: 'GeoTIFF'
});

// ===== 13. LEGEND =====
var legend = ui.Panel({style: {position: 'bottom-left', padding: '8px 15px'}});
legend.add(ui.Label('🌨️ PCA Snow Index', {fontWeight: 'bold'}));
legend.add(ui.Thumbnail({
  image: ee.Image.pixelLonLat().select(0),
  params: {
    bbox: [0,0,1,0.1], dimensions: '200x30', format: 'png', 
    min: p1, max: p99, palette: snowVis.palette
  },
  style: {stretch: 'horizontal'}
}));
Map.add(legend);

/**********************************************
 Wet vs Dry Snow classification (rule-based)
 Inputs: Sentinel-2 (NDSI) + Sentinel-1 (VV,VH in dB)
 Outputs: 0 = No snow, 1 = Dry snow, 2 = Wet snow
**********************************************/

// ---------- USER PARAMETERS ----------
var aoi = table;                       // replace with your imported table or geometry
Map.centerObject(aoi, 9);
Map.addLayer(aoi.style({color:"red", fillColor: "00000000"}), {}, "AOI");

// Target analysis period (when you want to classify snow)
var targetStart = '2017-05-01';
var targetEnd   = '2017-05-30';

// Baseline period used to detect abrupt backscatter drops (e.g., preceding month)
var baselineStart = '2017-04-01';
var baselineEnd   = '2017-04-30';

// Processing scale (meters)
var scale = 10;

// Thresholds (defaults — **tune for your area**)
var ndsiSnowThresh = 0.4;      // NDSI threshold to call snow
var vvAbsoluteThresh = -11;    // If VV (dB) < this -> likely wet (absolute rule)
var deltaVVThresh = -2.0;      // If targetVV - baselineVV < this (dB) -> likely wet (relative rule)
var minSnowAreaPixels = 10;    // minimum connected pixel area to keep (optional cleaning)

// --------------------------------------
// 1) Prepare Sentinel-2 median composite + NDSI
// --------------------------------------
var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterBounds(aoi)
  .filterDate(targetStart, targetEnd)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 80));

function prepS2(img) {
  var scl = img.select('SCL');
  // mask cloud & cloud-shadow classes from SCL (3=cloud shadow, 8-11 clouds/snow etc.)
  var cloudMask = scl.eq(3) //.or(scl.eq(8)).or(scl.eq(9)).or(scl.eq(10)).or(scl.eq(11));
  img = img.updateMask(cloudMask.not());
  var refl = img.select(['B2','B3','B4','B8','B11']).divide(10000);
  var ndsi = refl.normalizedDifference(['B3','B11']).rename('NDSI');
  var ndvi = refl.normalizedDifference(['B8','B4']).rename('NDVI');
  return refl.addBands(ndsi).addBands(ndvi).copyProperties(img, img.propertyNames());
}

var s2_prepped = s2.map(prepS2);
var s2_target = s2_prepped.median().clip(aoi);
Map.addLayer(s2_target.select(['B4','B3','B2']), {min:0, max:1.3}, 'S2 RGB (target)');

// NDSI image
var ndsiImg = s2_target.select('NDSI');
Map.addLayer(ndsiImg, {min:0, max:1, palette:['blue','white']}, 'NDSI (target)');

// Snow mask (binary)
var snowMask = ndsiImg.gt(ndsiSnowThresh).selfMask(); // masked image for display
Map.addLayer(snowMask, {palette:['cyan']}, 'Snow mask (NDSI>' + ndsiSnowThresh + ')');

// --------------------------------------
// 2) Prepare Sentinel-1 composites (target & baseline) in dB
// --------------------------------------
var s1col = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(aoi)
  .filter(ee.Filter.eq('instrumentMode','IW'))
  .filter(ee.Filter.eq('orbitProperties_pass','ASCENDING')) // optional: include both or choose pass
  .select(['VV','VH']); 

function toDB_safe(img) {
  var vv = img.select('VV').where(img.select('VV').lte(0), 1e-6);
  var vh = img.select('VH').where(img.select('VH').lte(0), 1e-6);
  return ee.Image([vv.log10().multiply(10).rename('VV'),
                   vh.log10().multiply(10).rename('VH')])
            .copyProperties(img, img.propertyNames());
}

// target composite
var s1_target = s1col.filterDate(targetStart, targetEnd).map(toDB_safe).median().clip(aoi);
Map.addLayer(s1_target, {bands:['VV'], min:-25, max:0}, 'S1 target VV (dB)');

// baseline composite
var s1_baseline = s1col.filterDate(baselineStart, baselineEnd).map(toDB_safe).median().clip(aoi);
Map.addLayer(s1_baseline, {bands:['VV'], min:-25, max:0}, 'S1 baseline VV (dB)');

// Compute delta VV (target - baseline); negative means drop in backscatter
var deltaVV = s1_target.select('VV').subtract(s1_baseline.select('VV')).rename('dVV');
Map.addLayer(deltaVV, {min:-6, max:6, palette:['red','white','blue']}, 'ΔVV (target - baseline)');

// Compute VH-VV difference and ratio (could help)
var vh = s1_target.select('VH');
var vv = s1_target.select('VV');
var vhMinusVv = vh.subtract(vv).rename('VH_minus_VV');
var vhDivVv = vh.expression('vh - vv', {vh: vh, vv: vv}).rename('VH_minus_VV2'); // same but as expression

// --------------------------------------
// 3) Classify wet vs dry within snow mask
//    Rule: wet if (VV < absolute threshold) OR (ΔVV < delta threshold)
//    Otherwise dry (within snowMask)
// --------------------------------------
var isSnow = ndsiImg.gt(ndsiSnowThresh); // boolean image (0/1)
var vvImg = s1_target.select('VV');
var dVVimg = deltaVV.select('dVV');

var wetRuleAbs = vvImg.lt(vvAbsoluteThresh);            // absolute low VV
var wetRuleDelta = dVVimg.lt(deltaVVThresh);           // sudden negative change
var wetRule = wetRuleAbs.or(wetRuleDelta);

// Apply only within snow
var wetSnow = isSnow.and(wetRule);
var drySnow = isSnow.and(wetRule.not());

// Build classified map: 0=no snow, 1=dry, 2=wet
var classImage = ee.Image(0).updateMask(ee.Image(1).clip(aoi)); // start with zeros inside AOI
classImage = classImage.where(drySnow, 1);
classImage = classImage.where(wetSnow, 2);

// Clean small speckles (optional): remove small objects using connectedPixelCount
var classMasked = classImage.updateMask(classImage.gte(0)); // keep only AOI
var connected = classMasked.connectedPixelCount(256, true); // window radius big enough
var classClean = classMasked.updateMask( connected.gte(minSnowAreaPixels) ).unmask(classMasked);

// Visualization palette
var palette = ['ffffff', 'ffe67f', 'e61a6a']; // no-snow white, dry yellow, wet red
Map.addLayer(classClean, {min:0, max:2, palette: palette}, 'Snow classification (0=none,1=dry,2=wet)');

// Legend helper (simple)
var legend = ui.Panel({style: {position: 'bottom-left'}});
legend.add(ui.Label('Legend: 0=no snow, 1=dry, 2=wet'));
Map.add(legend);

// --------------------------------------
// 4) Diagnostics & tuning aids
//    - Show areas where absolute and delta rules fired
//    - Export small thumbnails for inspection
// --------------------------------------
Map.addLayer(wetRuleAbs.updateMask(wetRuleAbs), {palette:['magenta']}, 'Wet rule (VV < ' + vvAbsoluteThresh + ' dB)');
Map.addLayer(wetRuleDelta.updateMask(wetRuleDelta), {palette:['orange']}, 'Wet rule (ΔVV < ' + deltaVVThresh + ' dB)');

// Print basic stats
print('Area stats (pixel counts) of classes in AOI (approx):',
  classClean.reduceRegion({
    reducer: ee.Reducer.frequencyHistogram(),
    geometry: aoi,
    scale: scale,
    maxPixels: 1e9
  }));

// --- Optional: export classification as GeoTIFF to Drive ---
Export.image.toDrive({
  image: classClean.toUint8(),
  description: 'snow_wet_dry_class',
  folder: 'Hohe-Tauren Snow Study',
  fileNamePrefix: 'snow_wet_dry_'+targetStart+'_'+targetEnd,
  region: aoi.geometry().bounds(),
  scale: 10,
  crs: 'EPSG:4326',
  maxPixels: 1e13
});

// --------------------------------------
// 5) Notes for tuning & validation
// --------------------------------------
/*
- Thresholds are region-dependent. Suggested calibration:
  * Collect a few reference points (field measurements or very high-res imagery) with known wet/dry labels.
  * Print S1 values at those points: print(s1_target.reduceRegions(...)) and adjust vvAbsoluteThresh / deltaVVThresh.

- If you have multiple days of S1, deltaVV can be computed against a multi-year dry-season baseline for robustness.

- Additional heuristics to improve classification:
  * Require NDSI > threshold AND NDVI below vegetation threshold to avoid bright vegetation confusion.
  * Use VH-VV or VH/VV ratio: wet snow may change polarization ratio.
  * Include topographic filters (slope/aspect) because steep shadows may bias S2 and S1.

- Evaluation:
  * If you have ground truth, compute confusion matrix and accuracy metrics in GEE (reduceRegions with confusion matrices).
  * Otherwise visually inspect classification together with S2 RGB, NDSI and S1 composites.

- If patchiness is an issue: apply morphological opening/closing or increase minSnowAreaPixels.
*/