Lab 23 - Where Did the Water Go?

Before We Begin

• Prerequisites: Spectral Indices module, Hydrology modules (Introduction and JRC Global Surface Water).
• Estimated time: 40 minutes
• You will need: GEE Code Editor access

Let's Get Started

Part 1: Our first water map

Let's start by loading a cloud-free Sentinel-2 image over Lake Mead and calculating MNDWI. We'll filter for summer imagery when water extent is relatively stable.

MNDWI works like a water detector: it compares Green light (which water reflects) to SWIR light (which water absorbs). The bigger the difference, the more confident we are that a pixel is water.

// ===== Part 1: Load Sentinel-2 and calculate MNDWI =====
var lakeMead = ee.Geometry.Point([-114.75, 36.15]);

// Load a low-cloud Sentinel-2 image from summer 2023
var s2 = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterBounds(lakeMead)
  .filterDate('2023-06-01', '2023-08-31')
  .sort('CLOUDY_PIXEL_PERCENTAGE')
  .first()
  .divide(10000); // scale to reflectance

// Calculate MNDWI = (Green - SWIR) / (Green + SWIR)
var mndwi = s2.normalizedDifference(['B3', 'B11']).rename('MNDWI');

// Display the MNDWI layer
Map.centerObject(lakeMead, 11);
Map.addLayer(s2, {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3}, 'True Color');
Map.addLayer(mndwi, {
  min: -0.5, max: 0.8,
  palette: ['saddlebrown', 'white', 'cyan', 'blue']
}, 'MNDWI');

Understanding the code:

• .sort('CLOUDY_PIXEL_PERCENTAGE').first() grabs the least cloudy image. Nice and clean.
• .divide(10000) converts Sentinel-2 DN values to surface reflectance.
• .normalizedDifference(['B3', 'B11']) computes (B3 - B11) / (B3 + B11).

Part 2: Drawing the line between water and land

Now let's turn our continuous MNDWI into a simple yes/no water map. We do this by picking a threshold: anything above it is water, everything else is not.

Here's the fun part: try using the Inspector tool to click around the water's edge and see what values you get. That will help us pick the right cutoff.

// ===== Part 2: Apply threshold for water mask =====
var waterMask = mndwi.gt(0);

// Use selfMask to show only water pixels (hides 0 values)
Map.addLayer(waterMask.selfMask(), {palette: ['0066FF']}, 'Water Mask (threshold = 0)');

// Try a refined threshold
var waterRefined = mndwi.gt(0.15);
Map.addLayer(waterRefined.selfMask(), {palette: ['003399']}, 'Water Mask (threshold = 0.15)');

Understanding the code:

• .gt(0) creates a binary image: 1 where MNDWI > 0, and 0 elsewhere.
• .selfMask() masks out pixels with value 0, making non-water pixels transparent. This keeps our map clean.
• A higher threshold (0.15) removes some false positives from wet soil or shadows. Think of it like tightening the definition of what counts as "water."

Part 3: How big is the lake right now?

Here's where remote sensing gets powerful. We can go from a pretty map to an actual measurement. Let's use ee.Image.pixelArea() and reduceRegion() to calculate the total water surface area.

// ===== Part 3: Calculate water area in km² =====
var roi = lakeMead.buffer(50000); // 50 km buffer

var waterAreaImage = ee.Image.pixelArea().updateMask(waterMask);

var stats = waterAreaImage.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 10,
  maxPixels: 1e9
});

var areaKm2 = ee.Number(stats.get('area')).divide(1e6);
print('Water surface area (km²):', areaKm2);

// For reference: Lake Mead at full capacity is ~640 km²
// Recent levels have been ~200-300 km²

Understanding the code:

• ee.Image.pixelArea() returns an image where each pixel holds its area in m². Think of it like a grid of tiny measuring tapes.
• .updateMask(waterMask) keeps only water pixels in the area calculation.
• ee.Reducer.sum() adds up all water pixel areas within the ROI.
• Dividing by 1,000,000 converts m² to km².

Part 4: What does 40 years of history tell us?

Our MNDWI gives us a snapshot, a single moment in time. But what if we could see every moment from the last four decades?

That's exactly what the JRC Global Surface Water dataset gives us. It analyzed the entire Landsat archive (1984 to present) and tells us how often each pixel has been water. Let's load it and compare.

// ===== Part 4: Load JRC GSW occurrence =====
var gsw = ee.Image('JRC/GSW1_4/GlobalSurfaceWater');
var occurrence = gsw.select('occurrence');

Map.addLayer(occurrence, {
  min: 0, max: 100,
  palette: ['ffffff', 'ffcccc', 'ff0000', '0000ff']
}, 'JRC Water Occurrence (%)');

// Compare: pixels with > 75% occurrence vs. your MNDWI mask
var jrcPermanent = occurrence.gte(75);
Map.addLayer(jrcPermanent.selfMask(), {palette: ['purple']}, 'JRC Permanent Water (>75%)');

// Overlay your MNDWI mask for comparison
Map.addLayer(waterMask.selfMask(), {palette: ['00FF00']}, 'Your MNDWI Water Mask');

Understanding the code:

• The JRC occurrence layer shows how often a pixel was water from 1984 to present.
• Here's the key insight: areas where occurrence is high but our MNDWI shows no water? That's where water has been lost.
• Toggle layers on and off to compare the two approaches. The gaps tell a powerful story.

Part 5: Reading the story of change

The JRC transition layer is like a biography for each pixel. It tells us how that pixel's relationship with water has changed over the decades. Let's see where water has been permanently lost around Lake Mead.

// ===== Part 5: Water transitions and loss =====
var transition = gsw.select('transition');

Map.addLayer(transition, {
  min: 0, max: 10,
  palette: [
    'ffffff', '0000ff', '22b14c', 'd1102d', '99d9ea',
    'b5e61d', 'e6a1aa', 'ff7f27', 'ffc90e', '7f7f7f', 'c3c3c3'
  ]
}, 'Water Transitions');

// Isolate lost permanent water (class 3)
var lostWater = transition.eq(3);
Map.addLayer(lostWater.selfMask(), {palette: ['FF0000']}, 'Lost Permanent Water');

// Isolate permanent-to-seasonal (class 8)
var permToSeasonal = transition.eq(8);
Map.addLayer(permToSeasonal.selfMask(), {palette: ['FFA500']}, 'Permanent to Seasonal');

Part 6: Putting a number on the loss

Let's quantify the difference between the lake's historical maximum extent and its current state. This is where we turn our visual observations into hard numbers.

// ===== Part 6: Historical vs. current extent =====
var maxExtent = gsw.select('max_extent');
var currentWater = occurrence.gte(75);

// Calculate historical maximum area
var maxArea = ee.Image.pixelArea().updateMask(maxExtent).reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e9
});
var maxKm2 = ee.Number(maxArea.get('area')).divide(1e6);

// Calculate current permanent water area
var currentArea = ee.Image.pixelArea().updateMask(currentWater).reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 30,
  maxPixels: 1e9
});
var currentKm2 = ee.Number(currentArea.get('area')).divide(1e6);

print('Historical maximum extent (km²):', maxKm2);
print('Current permanent water (km²):', currentKm2);
print('Area lost (km²):', maxKm2.subtract(currentKm2));
print('Percentage lost:', maxKm2.subtract(currentKm2).divide(maxKm2).multiply(100));

// Visual comparison
Map.addLayer(maxExtent.selfMask(), {palette: ['AADDFF']}, 'Historical Max Extent');
Map.addLayer(currentWater.selfMask(), {palette: ['0000CC']}, 'Current Permanent Water');