Crop Yield Modeling

Crop Yield Modeling & Forecasting

Overview

Develop predictive models for crop yield using weather data, soil conditions, satellite imagery, historical yields, and agronomic practices. Covers machine learning approaches (random forest, XGBoost, LSTM), feature engineering, model validation, and integration with farm management systems.

When to Use

• "Predict crop yields using weather/soason data"
• "Build ML model for yield forecasting"
• "Analyze satellite imagery for crop health"
• "Calibrate yield models with field data"
• "Integrate yield predictions into farm planning"

Essential Input Features

Weather Variables

Soil Variables

Management Variables

ML Model Implementation

import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, r2_score

def build_yield_model(farm_data_path, weather_data_path):
    """
    Build ML model for yield prediction
    
    Args:
        farm_data_path: historical farm data (yield, management)
        weather_data_path: weather data for growing season
    
    Returns:
        Trained model and performance metrics
    """
    # Load and merge data
    farm_df = pd.read_csv(farm_data_path)
    weather_df = pd.read_csv(weather_data_path)
    
    # Feature engineering for weather
    weather_df = weather_df.groupby('date').agg({
        'temp_max': 'max',
        'temp_min': 'min',
        'precipitation': 'sum',
        'gdd': 'sum',  # Growing degree days
        'vpd': 'mean'  # Vapor pressure deficit
    }).reset_index()
    
    weather_df['avg_temp'] = (weather_df['temp_max'] + weather_df['temp_min']) / 2
    
    # Key growing season periods (before harvest)
    critical_periods = {
        'emergence_to_vegetative': ('2023-05-01', '2023-06-15'),
        'vegetative_to_flowering': ('2023-06-16', '2023-07-30'),
        'flowering_to_grain_fill': ('2023-07-31', '2023-09-15')
    }
    
    # Calculate accumulated weather for each period
    features = {}
    for period_name, (start, end) in critical_periods.items():
        period_weather = weather_df[
            (weather_df['date'] >= start) & 
            (weather_df['date'] <= end)
        ]
        
        features[f'{period_name}_gdd'] = period_weather['gdd'].sum()
        features[f'{period_name}_precip'] = period_weather['precipitation'].sum()
        features[f'{period_name}_avg_temp'] = period_weather['avg_temp'].mean()
    
    # Merge with farm management data
    model_data = farm_df.merge(
        pd.DataFrame([features]), 
        left_index=True, right_index=True
    )
    
    # Prepare features for ML
    feature_cols = [
        'planting_density', 'fertilizer_n_rate', 'hybrid_year',
        'emergence_to_vegetative_gdd', 'flowering_to_grain_fill_precip',
        'avg_soil_moisture', 'organic_matter_pct'
    ]
    
    X = model_data[feature_cols]
    y = model_data['actual_yield']
    
    # Train/test split
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, random_state=42
    )
    
    # Try multiple models
    models = {
        'Random Forest': RandomForestRegressor(n_estimators=100, random_state=42),
        'Gradient Boosting': GradientBoostingRegressor(n_estimators=100, random_state=42)
    }
    
    best_model = None
    best_mae = float('inf')
    
    for name, model in models.items():
        model.fit(X_train, y_train)
        predictions = model.predict(X_test)
        mae = mean_absolute_error(y_test, predictions)
        
        if mae < best_mae:
            best_mae = mae
            best_model = model
    
    return {
        'model': best_model,
        'feature_importance': dict(zip(feature_cols, best_model.feature_importances_)),
        'mae': round(best_mae, 2),
        'r2_score': round(best_model.score(X_test, y_test), 3)
    }

Model Validation Strategy

Cross-Validation by Year

def year_stratified_cv(model_data, feature_cols, n_years=10):
    """
    Use leave-year-out cross-validation to test model generalization
    """
    years = sorted(model_data['year'].unique())
    results = []
    
    for test_year in years[-5:]:  # Last 5 years as test
        train_years = [y for y in years if y < test_year]
        
        train_data = model_data[model_data['year'].isin(train_years)]
        test_data = model_data[model_data['year'] == test_year]
        
        model = RandomForestRegressor(n_estimators=100, random_state=42)
        model.fit(train_data[feature_cols], train_data['actual_yield'])
        
        predictions = model.predict(test_data[feature_cols])
        mae = mean_absolute_error(test_data['actual_yield'], predictions)
        
        results.append({
            'test_year': test_year,
            'training_years': train_years,
            'mae': round(mae, 2),
            'r2': round(r2_score(test_data['actual_yield'], predictions), 3)
        })
    
    return results

Integration with Farm Management

class YieldForecast:
    def __init__(self, model, current_season_data):
        self.model = model
        self.data = current_season_data
        
    def get_forecast(self, confidence_level=0.8):
        """
        Generate yield forecast with confidence intervals
        """
        point_estimate = self.model.predict([self.data])[0]
        
        # Bootstrap confidence interval
        predictions = []
        for _ in range(1000):
            # Resample training data with replacement
            bootstrap_sample = np.random.choice(
                self.training_yields, 
                size=len(self.training_yields), 
                replace=True
            )
            predictions.append(np.random.normal(
                point_estimate, 
                np.std(bootstrap_predictions)
            ))
        
        ci_lower = np.percentile(predictions, (1-confidence_level)*50)
        ci_upper = np.percentile(predictions, (1+confidence_level)*50)
        
        return {
            'point_estimate': round(point_estimate, 2),
            'confidence_interval': [round(ci_lower, 2), round(ci_upper, 2)],
            'risk_level': self.classify_risk(point_estimate)
        }

def classify_risk(self, predicted_yield):
    """
    Classify yield risk level
    """
    historical_avg = self.historical_yields.mean()
    deviation = (predicted_yield - historical_avg) / historical_avg
    
    if deviation > 0.15:
        return 'HIGH_YIELD'
    elif deviation > 0.05:
        return 'MODERATE_HIGH'
    elif deviation > -0.05:
        return 'NORMAL'
    elif deviation > -0.15:
        return 'MODERATE_LOW'
    else:
        return 'LOW_YIELD'

Satellite & Remote Sensing Integration

NDVI-Based Growth Staging

def ndvi_growth_staging(ndvi_timeseries, planting_date, harvest_date):
    """
    Use NDVI to identify crop growth stages and correlate with yield impact
    """
    stages = {
        'emergence': {'ndvi_range': (0.1, 0.3), 'days_after_planting': '10-30'},
        'vegetative': {'ndvi_range': (0.3, 0.6), 'days_after_planting': '30-70'},
        'reproductive': {'ndvi_range': (0.6, 0.8), 'days_after_planting': '70-110'},
        'grain_fill': {'ndvi_range': (0.4, 0.7), 'days_after_planting': '110-140'},
        'senescence': {'ndvi_range': (0.1, 0.4), 'days_after_planting': '140+'}
    }
    
    # Identify peak NDVI timing
    peak_ndvi_date = ndvi_timeseries.loc[
        ndvi_timeseries['ndvi'].idxmax(), 'date'
    ]
    days_to_peak = (peak_ndvi_date - planting_date).days
    
    return {
        'peak_ndvi_date': peak_ndvi_date,
        'days_to_peak': days_to_peak,
        'growth_stage_at_peak': find_stage(days_to_peak, stages),
        'yield_impact_score': calculate_yield_impact(ndvi_timeseries, stages)
    }

Common Pitfalls

1. Weather data not aligned to growth stages — accumulating all season's weather instead of critical periods
2. Not accounting for year-to-year variability — climate change shifts optimal timing
3. Using too few training years — models need 8-10+ years of data
4. Ignoring interaction effects — fertilizer effectiveness depends on weather
5. Satellite cloud cover gaps — missing data during critical growth stages
6. Not validating with independent farms — models overfit to local conditions
7. Assuming linear relationships — crop responses are often non-linear
8. Not updating models — conditions change annually
9. Ignoring spatial variability — field zones need zone-specific models
10. Over-relying on satellite data alone — combine with ground truth

Verification Checklist

• Weather features calculated per growth stage, not entire season
• Model trained on ≥8 years of historical data
• Cross-validation by year (not random split)
• Feature importance analysis completed
• Confidence intervals calculated via bootstrapping
• Satellite NDVI aligned with growth staging
• Independent validation on separate farms
• Model updates scheduled for each harvest season
• Spatial variability addressed with zone-based modeling
• Risk classification aligned with farm management decisions