Lab 26 - Watching Cities Grow from Space

Before We Dive In

• Prerequisites: Spectral Indices module, Urban Analysis module
• Estimated time: 45 minutes
• You will need: GEE Code Editor access

Lab Instructions

Part 1: Can we see concrete from space?

Let's find out! We'll load Landsat 8 imagery, apply the scale factors, and calculate NDBI to highlight built-up surfaces across the Houston metro area.

NDBI works a lot like NDVI, but instead of highlighting vegetation, it highlights impervious surfaces. Concrete, asphalt, and rooftops reflect more shortwave infrared light than near-infrared, and NDBI exploits that difference.

// ===== PART 1: Calculate NDBI =====

// Define the Houston metro study area
var houstonCenter = ee.Geometry.Point([-95.37, 29.76]);
var roi = ee.Geometry.Rectangle([-95.8, 29.5, -95.0, 30.0]);

// Load Landsat 8 SR, summer 2023
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterDate('2023-06-01', '2023-08-31')
  .filterBounds(roi)
  .filter(ee.Filter.lt('CLOUD_COVER', 20))
  .median();

// Apply Landsat Collection 2 scale factors
var scaled = l8.multiply(0.0000275).add(-0.2);

// Calculate NDBI = (SWIR1 - NIR) / (SWIR1 + NIR)
// Landsat 8: SWIR1 = SR_B6, NIR = SR_B5
var ndbi = scaled.normalizedDifference(['SR_B6', 'SR_B5']).rename('NDBI');

// Visualize NDBI
Map.addLayer(ndbi.clip(roi), {
  min: -0.5, max: 0.5,
  palette: ['blue', 'white', 'red']
}, 'NDBI - Houston');
Map.centerObject(houstonCenter, 10);
print('NDBI range: -1 (vegetation) to +1 (built-up)');

Let's break that down:

• .normalizedDifference(['SR_B6', 'SR_B5']) computes (B6 - B5) / (B6 + B5), which is the NDBI formula. Same math we used for NDVI, just different bands!
• .multiply(0.0000275).add(-0.2) is the required scale/offset for Landsat Collection 2 surface reflectance
• Red pixels have high NDBI (built-up); blue pixels have low NDBI (vegetation). The city should pop right out.

Part 2: Drawing the line between city and countryside

Now let's turn our continuous NDBI image into a simple yes-or-no question: is this pixel built-up or not? We do this by applying a threshold. It's like setting a cutoff grade on an exam.

// ===== PART 2: Threshold NDBI to create urban mask =====

// Pixels with NDBI > 0 are classified as built-up
var urbanMask = ndbi.gt(0).rename('urban');

// Self-mask to show only built-up pixels
Map.addLayer(urbanMask.selfMask().clip(roi), {
  palette: ['red']
}, 'Urban Mask (NDBI > 0)');

// Also try a stricter threshold
var urbanStrict = ndbi.gt(0.1).rename('urban_strict');
Map.addLayer(urbanStrict.selfMask().clip(roi), {
  palette: ['darkred']
}, 'Urban Mask (NDBI > 0.1)');

print('Compare the two thresholds by toggling layers on/off');

Let's break that down:

• .gt(0) creates a binary image where 1 = NDBI above threshold, 0 = below. Simple and elegant.
• .selfMask() makes 0 pixels transparent so only built-up areas are visible on the map
• A threshold of 0 captures most built-up areas but may include some bare soil; 0.1 is stricter but may miss low-density residential areas. Toggle both layers on and off to see the difference!

Part 3: What does a city look like at night?

Here's where it gets really interesting. Let's look at Houston from a completely different perspective: nighttime lights from the VIIRS satellite.

Think of it this way: NDBI tells us where the concrete is, but nighttime lights tell us where the people are. Lights reveal economic activity and population density in a way that spectral indices simply can't.

// ===== PART 3: VIIRS Nighttime Lights =====

// Load VIIRS monthly composite for June 2023
var ntl2023 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG')
  .filterDate('2023-06-01', '2023-06-30')
  .first()
  .select('avg_rad');

// Visualize with a night-sky color palette
Map.addLayer(ntl2023.clip(roi), {
  min: 0, max: 80,
  palette: ['black', 'navy', 'blue', 'purple', 'magenta', 'orange', 'yellow', 'white']
}, 'Nighttime Lights - 2023');

// Print statistics
var ntlStats = ntl2023.reduceRegion({
  reducer: ee.Reducer.minMax(),
  geometry: roi,
  scale: 500,
  maxPixels: 1e9
});
print('Nighttime light range (nW/cm2/sr):', ntlStats);

Part 4: Where did the city grow in the last decade?

Now let's use nighttime lights as a time machine. We'll load data from 2014, subtract it from 2023, and see exactly where Houston has been expanding. Areas that got brighter over the decade are areas that urbanized.

// ===== PART 4: Nighttime Light Change Detection =====

// Load 2014 nighttime lights
var ntl2014 = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG')
  .filterDate('2014-06-01', '2014-06-30')
  .first()
  .select('avg_rad');

// Calculate change (positive = brighter = growth)
var ntlChange = ntl2023.subtract(ntl2014).rename('NTL_Change');

// Visualize change
Map.addLayer(ntlChange.clip(roi), {
  min: -10, max: 40,
  palette: ['blue', 'lightblue', 'white', 'orange', 'red']
}, 'NTL Change (2014-2023)');

// Add side-by-side comparison layers
Map.addLayer(ntl2014.clip(roi), {
  min: 0, max: 80,
  palette: ['black', 'yellow', 'white']
}, 'NTL 2014');

Map.addLayer(ntl2023.clip(roi), {
  min: 0, max: 80,
  palette: ['black', 'yellow', 'white']
}, 'NTL 2023');

print('Red = areas that became brighter (urban growth)');
print('Blue = areas that became dimmer (rare in growing cities)');

Part 5: How much did the city actually grow?

Maps are great, but sometimes we need a number. Let's quantify the urban footprint by calculating the total area of built-up pixels for both years. This is where ee.Image.pixelArea() becomes our best friend.

// ===== PART 5: Urban Area Statistics =====

// Define "urban" as nighttime light radiance > 5 nW/cm2/sr
var urban2014 = ntl2014.gt(5).rename('urban_2014');
var urban2023 = ntl2023.gt(5).rename('urban_2023');

// Calculate area of urban pixels using pixelArea()
var urbanArea2014 = urban2014.multiply(ee.Image.pixelArea())
  .reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: roi,
    scale: 500,
    maxPixels: 1e9
  });

var urbanArea2023 = urban2023.multiply(ee.Image.pixelArea())
  .reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: roi,
    scale: 500,
    maxPixels: 1e9
  });

// Print results (values are in square meters)
print('=== Urban Area Statistics ===');
print('Urban area 2014 (sq m):', urbanArea2014);
print('Urban area 2023 (sq m):', urbanArea2023);

// Calculate total ROI area for context
var totalArea = ee.Image.pixelArea().reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: roi,
  scale: 500,
  maxPixels: 1e9
});
print('Total ROI area (sq m):', totalArea);

// Visualize both urban extents
Map.addLayer(urban2014.selfMask().clip(roi), {palette: ['yellow']}, 'Urban Extent 2014');
Map.addLayer(urban2023.selfMask().clip(roi), {palette: ['red']}, 'Urban Extent 2023');

Let's break that down:

• ee.Image.pixelArea() creates an image where each pixel's value equals its ground area in square meters. Think of it like stamping a tiny measuring tape onto every pixel.
• Multiplying the binary urban mask by pixel area gives the area of only the urban pixels (non-urban pixels become zero)
• ee.Reducer.sum() adds up all those pixel areas to get total urban area
• Divide by 1,000,000 to convert square meters to square kilometers for more interpretable numbers

Part 6: Putting it all together with a land use map

Here's the cool part. We can create a simple but powerful land use classification by combining three spectral indices we already know: NDVI, NDBI, and NDWI. Each one is like a specialist detector, and together, they separate Houston's landscape into water, vegetation, and built-up categories.

// ===== PART 6: Combined Index Classification =====

// Calculate NDVI and NDWI from the same Landsat image
var ndvi = scaled.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
var ndwi = scaled.normalizedDifference(['SR_B3', 'SR_B5']).rename('NDWI');

// Classification rules:
// Water: NDWI > 0
// Vegetation: NDVI > 0.3 (and not water)
// Built-up: NDBI > 0 (and not water or vegetation)
// Other: everything else
var water = ndwi.gt(0);
var vegetation = ndvi.gt(0.3).and(water.not());
var builtUp = ndbi.gt(0).and(water.not()).and(vegetation.not());

// Create classified image
var classification = ee.Image(0)
  .where(water, 1)       // 1 = Water (blue)
  .where(vegetation, 2)  // 2 = Vegetation (green)
  .where(builtUp, 3)     // 3 = Built-up (red)
  .selfMask()
  .clip(roi);

Map.addLayer(classification, {
  min: 1, max: 3,
  palette: ['#2196F3', '#4CAF50', '#F44336']
}, 'Land Use Classification');

// Calculate area for each class
var classAreas = ee.Image.pixelArea().addBands(classification)
  .reduceRegion({
    reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'class'
    }),
    geometry: roi,
    scale: 30,
    maxPixels: 1e9
  });
print('Area by class (sq m):', classAreas);
print('Class legend: 1=Water, 2=Vegetation, 3=Built-up');