Lab 25 - Reading the Rhythm of Crops

Before We Start

• Prerequisites: Time Series Charting module, Agricultural Monitoring module
• Estimated time: 40 minutes
• You will need: GEE Code Editor access

Let's Get Started

Part 1: What's growing here?

Let's start by loading the 2023 CDL to see what crops are planted across central Iowa. The CDL uses numeric codes for each crop type, and it comes with a built-in color palette that makes the map look like a patchwork quilt.

// ===== PART 1: Load the Cropland Data Layer =====

// Define study area in central Iowa
var studyArea = ee.Geometry.Point([-93.5, 42.0]);

// Load CDL for 2023
var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.eq('system:index', '2023'))
  .first()
  .select('cropland');

// Add CDL to map (uses built-in palette)
Map.addLayer(cdl, {}, 'Cropland Data Layer 2023');
Map.centerObject(studyArea, 12);

// Print crop class names for reference
print('CDL loaded - Key codes: 1=Corn, 5=Soybeans, 24=Winter Wheat');

Understanding the code:

• ee.Filter.eq('system:index', '2023') - selects the 2023 CDL by its index property (a string, not a date)
• .select('cropland') - extracts the crop classification band

Part 2: Dropping pins in the fields

Now let's use the CDL to place sample points in fields of known crop type. The trick is to choose points in the center of large, uniform fields to avoid mixed pixels. Think of it like taking a soil sample: we want a clean reading, not a mix of two fields.

// ===== PART 2: Define crop sample points =====

// Corn field in central Iowa (CDL code 1)
var cornPoint = ee.Geometry.Point([-93.50, 42.00]);

// Soybean field nearby (CDL code 5)
var soyPoint = ee.Geometry.Point([-93.45, 42.05]);

// Winter wheat in central Kansas (CDL code 24)
var wheatPoint = ee.Geometry.Point([-97.50, 38.80]);

// Verify crop types at each point
var cornType = cdl.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: cornPoint,
  scale: 30,
  maxPixels: 1e9
});
var soyType = cdl.reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: soyPoint,
  scale: 30,
  maxPixels: 1e9
});

print('Corn point crop code:', cornType);
print('Soybean point crop code:', soyType);

// Add points to map
Map.addLayer(cornPoint, {color: 'yellow'}, 'Corn Sample');
Map.addLayer(soyPoint, {color: 'green'}, 'Soybean Sample');
Map.addLayer(wheatPoint, {color: 'orange'}, 'Wheat Sample');

Understanding the code:

• ee.Reducer.first() - extracts the pixel value at the exact point location
• Verify the printed values: corn should print 1, soybeans should print 5
• If the values do not match, use the Inspector tool to click on corn and soybean fields and update coordinates

Part 3: Loading the satellite heartbeat monitor

Let's load a full year of MODIS EVI data and apply the scale factor. MODIS MOD13A2 gives us 16-day composites at 1 km resolution with pre-computed EVI values. Think of each image as a check-up on the fields' health.

// ===== PART 3: Prepare MODIS EVI time series =====

// Load MODIS 16-day EVI for 2023
var modisEVI = ee.ImageCollection('MODIS/061/MOD13A2')
  .filterDate('2023-01-01', '2023-12-31')
  .select('EVI');

// Scale EVI values (MODIS stores EVI * 10000)
var eviScaled = modisEVI.map(function(image) {
  return image.multiply(0.0001)
    .copyProperties(image, ['system:time_start']);
});

// Check how many images we have
print('Number of MODIS EVI images:', eviScaled.size());
print('Date range:', eviScaled.first().get('system:time_start'));

Understanding the code:

• .multiply(0.0001) - converts raw MODIS EVI (stored as integers) to the standard 0 to 1 range
• .copyProperties(image, ['system:time_start']) - preserves the timestamp so charts display correct dates on the x-axis
• You should see approximately 23 images (one every 16 days for a full year)

Part 4: The phenology showdown

Here is the moment we have been building toward. Let's chart EVI over time for all three crop types on a single graph. This comparison reveals how each crop's growing season differs in timing, intensity, and duration. It is like watching three different runners in a race.

// ===== PART 4: Create phenology charts =====

// Create a FeatureCollection with labeled crop points
var cropPoints = ee.FeatureCollection([
  ee.Feature(cornPoint, {label: 'Corn'}),
  ee.Feature(soyPoint, {label: 'Soybean'}),
  ee.Feature(wheatPoint, {label: 'Winter Wheat'})
]);

// Chart EVI for all three crops
var phenologyChart = ui.Chart.image.seriesByRegion(
  eviScaled, cropPoints, ee.Reducer.mean(), 'EVI', 1000,
  'system:time_start', 'label')
  .setChartType('LineChart')
  .setOptions({
    title: 'Crop Phenology Comparison (2023)',
    vAxis: {
      title: 'EVI',
      viewWindow: {min: 0, max: 1}
    },
    hAxis: {title: 'Date'},
    lineWidth: 2,
    pointSize: 4,
    series: {
      0: {color: '#FFD700'}, // Corn - gold
      1: {color: '#228B22'}, // Soybean - green
      2: {color: '#D2691E'}  // Wheat - brown
    }
  });
print(phenologyChart);

// Also chart corn individually for a closer look
var cornChart = ui.Chart.image.series(eviScaled, cornPoint, ee.Reducer.mean(), 1000)
  .setOptions({
    title: 'Corn Phenology - Central Iowa (2023)',
    vAxis: {title: 'EVI', viewWindow: {min: 0, max: 1}},
    hAxis: {title: 'Date'},
    lineWidth: 2, pointSize: 4,
    colors: ['#FFD700']
  });
print(cornChart);

Part 5: Reading the growth curves

Now let's analyze the timing differences between our crops. We will look at when each one greens up, hits its peak, and starts to fade. This information is exactly what crop forecasters use to predict yields and spot problems early.

// ===== PART 5: Analyze phenology timing =====

// Create individual charts to study timing more closely
// Corn phenology
var cornSeries = ui.Chart.image.series(eviScaled, cornPoint, ee.Reducer.mean(), 1000)
  .setOptions({
    title: 'Corn Phenology Detail',
    vAxis: {title: 'EVI'},
    hAxis: {title: 'Date'},
    lineWidth: 2, pointSize: 5,
    colors: ['#FFD700'],
    trendlines: {0: {type: 'polynomial', degree: 5, color: 'gray'}}
  });
print(cornSeries);

// Soybean phenology
var soySeries = ui.Chart.image.series(eviScaled, soyPoint, ee.Reducer.mean(), 1000)
  .setOptions({
    title: 'Soybean Phenology Detail',
    vAxis: {title: 'EVI'},
    hAxis: {title: 'Date'},
    lineWidth: 2, pointSize: 5,
    colors: ['#228B22'],
    trendlines: {0: {type: 'polynomial', degree: 5, color: 'gray'}}
  });
print(soySeries);

// Print interpretation guide
print('--- Phenology Interpretation Guide ---');
print('Corn: Green-up ~May, Peak ~mid-July, Senescence ~Sep');
print('Soybean: Green-up ~June, Peak ~late July, Senescence ~Oct');
print('Winter Wheat: Green-up ~Mar, Peak ~May, Harvest ~June');

What to look for:

• Green-up date: When does EVI first rise above 0.3? This indicates crop emergence.
• Peak date: What date shows the highest EVI value? This is maximum canopy development.
• Senescence: When does EVI drop back below 0.3? This coincides with maturity and harvest.
• Growing season length: Count the days between green-up and senescence.

Part 6: Putting numbers on the curves

Let's extract some hard numbers from our time series. Peak EVI tells us about crop vigor, while growing season length varies by crop type and latitude. These are the metrics that matter for yield forecasting.

// ===== PART 6: Quantitative phenology metrics =====

// Calculate peak (maximum) EVI for each crop point
var peakCorn = eviScaled.max().reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: cornPoint,
  scale: 1000,
  maxPixels: 1e9
});

var peakSoy = eviScaled.max().reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: soyPoint,
  scale: 1000,
  maxPixels: 1e9
});

var peakWheat = eviScaled.max().reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: wheatPoint,
  scale: 1000,
  maxPixels: 1e9
});

print('Peak EVI - Corn:', peakCorn);
print('Peak EVI - Soybean:', peakSoy);
print('Peak EVI - Winter Wheat:', peakWheat);

// Calculate growing season length using a threshold approach
// Count the number of 16-day periods where EVI > 0.3
var growingCorn = eviScaled.map(function(img) {
  return img.gt(0.3).copyProperties(img, ['system:time_start']);
}).sum().reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: cornPoint,
  scale: 1000,
  maxPixels: 1e9
});

var growingSoy = eviScaled.map(function(img) {
  return img.gt(0.3).copyProperties(img, ['system:time_start']);
}).sum().reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: soyPoint,
  scale: 1000,
  maxPixels: 1e9
});

var growingWheat = eviScaled.map(function(img) {
  return img.gt(0.3).copyProperties(img, ['system:time_start']);
}).sum().reduceRegion({
  reducer: ee.Reducer.first(),
  geometry: wheatPoint,
  scale: 1000,
  maxPixels: 1e9
});

// Each count * 16 days = approximate growing season length
print('Growing season (16-day periods above 0.3 EVI):');
print('  Corn:', growingCorn);
print('  Soybean:', growingSoy);
print('  Winter Wheat:', growingWheat);
print('Multiply count by 16 to get approximate days');