Extra Exercises - Learn Python with BNHR

import ee
import geemap

# Trigger the authentication flow.
ee.Authenticate()

# Initialize the library.
ee.Initialize(project='bnhr-gee') #replace your-gee-project w/ your actual GEE project name

Building Seasonal Composites (Dry vs. Wet Season)¶

Concept: Your GEE Sample Process workflow includes a toggle for “Seasonal Composites,” which generates two distinct map layers: one for the dry season (Jan-Jun) and one for the wet season (Jul-Dec). This is highly useful for agricultural classification because crops look vastly different depending on the season.

# 1. Set up the interactive map
Map2 = geemap.Map()

# 2. Define an Area of Interest based on the survey data coordinates (Mindanao region)
aoi = ee.Geometry.Polygon([
  [[124.2, 7.3], [124.6, 7.3], [124.6, 7.7], [124.2, 7.7], [124.2, 7.3]]
])

# 3. Create the Dry Season Composite (January to June)
dry_season = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
              .filterBounds(aoi)
              .filterDate('2024-01-01', '2024-06-30') # First half of the year
              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
              .median()
              .clip(aoi))

# 4. Create the Wet Season Composite (July to December)
wet_season = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
              .filterBounds(aoi)
              .filterDate('2024-07-01', '2024-12-31') # Second half of the year
              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
              .median()
              .clip(aoi))

# 5. Add both layers to the map
# We use True Color bands (B4, B3, B2) for visualization
Map2.centerObject(aoi, 11)
Map2.addLayer(dry_season.select(['B4', 'B3', 'B2']), {'min': 0, 'max': 0.3}, 'Dry Season (Jan-Jun)')
Map2.addLayer(wet_season.select(['B4', 'B3', 'B2']), {'min': 0, 'max': 0.3}, 'Wet Season (Jul-Dec)')

print("Seasonal composites generated. Use the layer manager (top right of map) to toggle between Dry and Wet seasons.")
Map2

Seasonal composites generated. Use the layer manager (top right of map) to toggle between Dry and Wet seasons.

Exercise 3: Crop NDVI Phenology Smoothing¶

Concept: You monitor “Crop NDVI Phenology” to track the health of plants over time and estimate planting/harvesting dates. Raw satellite data is often noisy, so your JS script uses a Gaussian spatial filter (ee.Kernel.gaussian) to “smooth” the NDVI values.

Here, we will grab a single Tomato farm parcel from your survey data, calculate the NDVI for a specific date, apply the Gaussian smoother, and visually compare the “Raw” vs. “Smoothed” data.

# 1. Initialize Map
Map3 = geemap.Map()

# 2. Extract a single Tomato farm polygon from the PSA Survey data
farm_polygon = ee.Geometry.Polygon([
    [[121.4543512, 14.0900565], [121.4547529, 14.0899463], 
     [121.4552447, 14.0897183], [121.4553856, 14.0899251], 
     [121.4543512, 14.0900565]]
])

# 3. FIXED: Widened dates and used median() to ensure the image is not null
farm_image = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
              .filterBounds(farm_polygon)
              .filterDate('2024-01-01', '2024-06-30')
              .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
              .median()
              .clip(farm_polygon))

# 4. Calculate raw NDVI
ndvi_raw = farm_image.normalizedDifference(['B8', 'B4']).rename('NDVI')

# 5. Apply the Gaussian Smoothing Filter
smoothing_kernel = ee.Kernel.gaussian(radius=2, sigma=1)
ndvi_smoothed = ndvi_raw.convolve(smoothing_kernel).rename('NDVI_Smooth')

# 6. Map the results
Map3.centerObject(farm_polygon, 18) 

# Add Raw NDVI layer
Map3.addLayer(ndvi_raw, {'min': 0, 'max': 0.8, 'palette': ['white', 'green']}, 'Raw NDVI (Noisy)')

# Add Smoothed NDVI layer
Map3.addLayer(ndvi_smoothed, {'min': 0, 'max': 0.8, 'palette': ['white', 'green']}, 'Smoothed NDVI (Cleaned)')

# Outline the farm parcel
Map3.addLayer(ee.FeatureCollection(farm_polygon).style(color='red', fillColor='00000000'), {}, 'Tomato Farm Boundary')

print("Phenology Smoothing Complete. Toggle between the Raw and Smoothed layers to see how the Gaussian kernel cleans the data.")
Map3

Phenology Smoothing Complete. Toggle between the Raw and Smoothed layers to see how the Gaussian kernel cleans the data.

Processing and Mapping¶

This script integrates Sentinel-2, Landsat 8/9, MODIS, and DEM data, calculates indices (NDVI, NDWI, NDBI, EVI), and displays the result on an interactive geemap widget.

import geemap

# 1. Authenticate and Initialize (Run this once per session)
# ee.Authenticate() # Uncomment if running for the very first time
ee.Initialize(project='bnhr-gee') # Replace with your GCP project ID

# 2. Initialize interactive map (Replaces standard JS Map)
Map = geemap.Map()

# Define a sample AOI (e.g., a province boundary or manual polygon)
aoi = ee.Geometry.Polygon([
  [[124.615, 8.445], [124.610, 8.390], [124.645, 8.390], [124.650, 8.435], [124.615, 8.445]]
])

# 3. Processing Function
def getCombinedImages(aoi_region, maxCloud=20):
    # Sentinel-2 with calculated indices
    s2 = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
          .filterBounds(aoi_region)
          .filterDate('2024-01-01', '2024-12-31')
          .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', maxCloud))
          .median().clip(aoi_region)
          .select(['B2','B3','B4','B8','B11']))
    
    # Landsat 8/9
    l8 = (ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
          .merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2'))
          .filterBounds(aoi_region)
          .filterDate('2024-01-01', '2024-12-31')
          .filter(ee.Filter.lt('CLOUD_COVER', maxCloud))
          .median().clip(aoi_region)
          .multiply(0.0000275).subtract(0.2)
          .select(['SR_B2','SR_B3','SR_B4','SR_B5'])
          .rename(['B2','B3','B4','B8']))
    
    # Combine S2 and L8
    combined = s2.addBands(l8)
    
    # Calculate Indices
    ndvi = combined.normalizedDifference(['B8','B4']).rename('NDVI')
    ndwi = s2.normalizedDifference(['B3','B8']).rename('NDWI')
    ndbi = s2.normalizedDifference(['B11','B8']).rename('NDBI')
    
    # EVI calculation requires an expression mapping
    evi = s2.expression(
        '2.5*((NIR-RED)/(NIR+6*RED-7.5*BLUE+1))',
        {'NIR': s2.select('B8'), 'RED': s2.select('B4'), 'BLUE': s2.select('B2')}
    ).rename('EVI')
    
    combined = combined.addBands([ndvi, ndwi, ndbi, evi])
    
    # Add MODIS
    modis = (ee.ImageCollection('MODIS/006/MOD13Q1')
             .filterBounds(aoi_region)
             .filterDate('2024-01-01', '2024-12-31')
             .select(['NDVI','EVI'])
             .map(lambda img: img.multiply(0.0001))
             .median().clip(aoi_region))
    
    combined = combined.addBands(modis)
    
    # Add DEM & Slope
    dem = ee.Image('USGS/SRTMGL1_003').clip(aoi_region)
    slope = ee.Terrain.slope(dem).rename('slope')
    combined = combined.addBands([dem.rename('DEM'), slope])
    
    return combined

# 4. Execute and add to map
processed_image = getCombinedImages(aoi, maxCloud=20)

Map.centerObject(aoi, 11)
Map.addLayer(processed_image.select(['B4', 'B3', 'B2']), {'min': 0, 'max': 0.3}, 'True Color (S2)')
Map.addLayer(processed_image.select('NDVI'), {'min': 0, 'max': 1, 'palette': ['white', 'green']}, 'NDVI Layer')
Map