Day 4, Session 1: LSTMs for Earth Observation Time Series

๐ŸŽฏ Learning Objectives

By the end of this session, you will be able to: 1. Understand the importance of time series analysis in Earth Observation 2. Explain why RNNs face challenges with long sequences (vanishing/exploding gradients) 3. Describe LSTM architecture and how gates solve RNN limitations 4. Implement an LSTM model for drought prediction using NDVI time series 5. Apply LSTMs to Philippine EO challenges (Mindanao drought monitoring)

๐Ÿ“‹ Session Overview

  • Duration: 1.5 hours
  • Prerequisites: Basic understanding of neural networks (from Day 3)
  • Application Focus: Drought monitoring in Mindanao agricultural regions
  • Key Dataset: Simulated Sentinel-2 NDVI time series (2019-2024)

๐ŸŒ Philippine Context

The Philippines experiences significant climate variability, with El Niรฑo events causing severe droughts, particularly affecting Mindanaoโ€™s agricultural regions. PAGASA reports that drought events have increased in frequency and intensity, making early warning systems critical for: - Food security in Bukidnon and South Cotabato - Water resource management for irrigation systems - Agricultural planning and crop insurance programs

Time series analysis of satellite-derived vegetation indices enables us to detect early drought signals and predict future conditions.

๐Ÿ“š Module 1: Introduction to Time Series in Earth Observation

What are EO Time Series?

Earth Observation time series are sequences of measurements taken at regular intervals from the same location. Common examples include:

  1. NDVI (Normalized Difference Vegetation Index): Measures vegetation health
    • Range: -1 to +1 (higher values = healthier vegetation)
    • Sentinel-2 provides 5-day revisit time
  2. SAR Backscatter: Radar signal strength indicating surface properties
    • Sentinel-1 provides 6-12 day revisit
    • Sensitive to soil moisture and vegetation structure
  3. Land Surface Temperature: Thermal measurements from satellites
    • Critical for drought and heat stress monitoring

Why Time Series Matter for Philippine EO Applications

  • Phenology Tracking: Monitor rice cropping calendars in Central Luzon
  • Drought Detection: Early warning for Mindanao agricultural zones
  • Land Change Detection: Urban expansion in Metro Manila
  • Disaster Impact Assessment: Pre/post typhoon vegetation analysis

๐Ÿ’ญ Think-Through Discussion

Question: How might seasonal patterns in NDVI differ between irrigated rice fields in Nueva Ecija and rainfed corn farms in Bukidnon? What implications does this have for drought monitoring?

Consider: irrigation availability, crop types, planting calendars, climate zones

๐Ÿ› ๏ธ Setup and Environment Configuration

# Import required libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')

# TensorFlow and Keras imports
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout
from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau
from tensorflow.keras.optimizers import Adam

# Scikit-learn for preprocessing
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

# Set random seeds for reproducibility
np.random.seed(42)
tf.random.set_seed(42)

# Configure matplotlib for better visualizations
plt.rcParams['figure.figsize'] = (12, 6)
plt.rcParams['font.size'] = 10
sns.set_style('whitegrid')

print(f"TensorFlow version: {tf.__version__}")
print(f"GPU Available: {tf.config.list_physical_devices('GPU')}")
print("Setup complete! ๐Ÿš€")

๐Ÿ“Š Data Generation: Simulating Realistic NDVI Time Series

Weโ€™ll create a synthetic but realistic NDVI time series that mimics conditions in Bukidnon, Mindanao: - Normal seasonal patterns (wet/dry seasons) - El Niรฑo drought events (2019, 2023) - Random variations and noise

def generate_mindanao_ndvi_timeseries(start_date='2019-01-01', end_date='2024-12-31', 
                                      location='Bukidnon'):
    """
    Generate synthetic NDVI time series mimicking Mindanao agricultural patterns.
    Includes seasonal variations, drought events, and realistic noise.
    
    Parameters:
    -----------
    start_date : str
        Start date of the time series
    end_date : str
        End date of the time series
    location : str
        Location name for reference
    
    Returns:
    --------
    pd.DataFrame
        DataFrame with date, NDVI, precipitation, and drought index
    """
    
    # Create date range (10-day composites, similar to Sentinel-2)
    dates = pd.date_range(start=start_date, end=end_date, freq='10D')
    n_samples = len(dates)
    
    # Initialize arrays
    ndvi = np.zeros(n_samples)
    precipitation = np.zeros(n_samples)
    
    # Base NDVI for healthy vegetation in Mindanao
    base_ndvi = 0.75
    
    for i, date in enumerate(dates):
        # Seasonal component (wet season: Nov-Apr, dry season: May-Oct)
        month = date.month
        if month in [11, 12, 1, 2, 3, 4]:  # Wet season
            seasonal_factor = 0.85 + 0.1 * np.sin(2 * np.pi * month / 12)
            precip_base = 250 + 50 * np.random.randn()  # mm/month
        else:  # Dry season
            seasonal_factor = 0.7 + 0.05 * np.sin(2 * np.pi * month / 12)
            precip_base = 100 + 30 * np.random.randn()  # mm/month
        
        # El Niรฑo drought events (2019 Q2-Q3, 2023 Q1-Q2)
        drought_factor = 1.0
        if (date.year == 2019 and 4 <= month <= 9) or \
           (date.year == 2023 and 2 <= month <= 7):
            drought_factor = 0.6 + 0.2 * np.random.random()
            precip_base *= 0.4  # Reduced precipitation during drought
        
        # Calculate NDVI with noise
        ndvi[i] = base_ndvi * seasonal_factor * drought_factor
        ndvi[i] += np.random.normal(0, 0.03)  # Add noise
        ndvi[i] = np.clip(ndvi[i], 0.1, 0.95)  # Realistic bounds
        
        # Calculate precipitation
        precipitation[i] = max(0, precip_base)
    
    # Smooth the time series (moving average)
    window = 3
    ndvi_smooth = pd.Series(ndvi).rolling(window=window, center=True).mean()
    ndvi_smooth = ndvi_smooth.fillna(method='bfill').fillna(method='ffill')
    
    # Calculate drought index (simplified: based on NDVI deviation from normal)
    ndvi_mean = ndvi_smooth.rolling(window=36, center=True).mean().fillna(method='bfill').fillna(method='ffill')
    drought_index = (ndvi_smooth - ndvi_mean) / ndvi_mean.std()
    
    # Create DataFrame
    df = pd.DataFrame({
        'date': dates,
        'ndvi': ndvi_smooth.values,
        'precipitation': precipitation,
        'drought_index': drought_index.values,
        'location': location
    })
    
    return df

# Generate data for Bukidnon, Mindanao
df_mindanao = generate_mindanao_ndvi_timeseries()
print(f"Generated {len(df_mindanao)} time points of NDVI data")
print(f"Date range: {df_mindanao['date'].min()} to {df_mindanao['date'].max()}")
print(f"\nFirst 5 rows:")
print(df_mindanao.head())
print(f"\nBasic statistics:")
print(df_mindanao[['ndvi', 'precipitation', 'drought_index']].describe())

๐Ÿ“ˆ Visualizing the Time Series Data

Letโ€™s visualize our NDVI time series to understand the patterns, including the drought events.

# Create comprehensive visualization
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)

# Plot 1: NDVI Time Series
ax1 = axes[0]
ax1.plot(df_mindanao['date'], df_mindanao['ndvi'], 'g-', linewidth=1.5, label='NDVI')
ax1.axhline(y=df_mindanao['ndvi'].mean(), color='gray', linestyle='--', alpha=0.5, label='Mean NDVI')

# Highlight drought periods
drought_periods = [
    ('2019-04-01', '2019-09-30', '2019 El Niรฑo'),
    ('2023-02-01', '2023-07-31', '2023 El Niรฑo')
]
for start, end, label in drought_periods:
    ax1.axvspan(pd.to_datetime(start), pd.to_datetime(end), alpha=0.2, color='orange', label=label)

ax1.set_ylabel('NDVI', fontsize=11)
ax1.set_title('NDVI Time Series for Bukidnon, Mindanao (2019-2024)', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_ylim([0.2, 0.9])

# Plot 2: Precipitation
ax2 = axes[1]
ax2.bar(df_mindanao['date'], df_mindanao['precipitation'], color='blue', alpha=0.6, width=8)
ax2.set_ylabel('Precipitation (mm)', fontsize=11)
ax2.set_title('Precipitation Patterns', fontsize=11)
ax2.grid(True, alpha=0.3)

# Plot 3: Drought Index
ax3 = axes[2]
colors = ['red' if x < -1 else 'orange' if x < 0 else 'green' for x in df_mindanao['drought_index']]
ax3.scatter(df_mindanao['date'], df_mindanao['drought_index'], c=colors, alpha=0.6, s=10)
ax3.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
ax3.axhline(y=-1, color='red', linestyle='--', alpha=0.5, label='Drought threshold')
ax3.set_ylabel('Drought Index', fontsize=11)
ax3.set_xlabel('Date', fontsize=11)
ax3.set_title('Drought Index (Negative values indicate drought stress)', fontsize=11)
ax3.legend(loc='upper right')
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print drought statistics
drought_days = (df_mindanao['drought_index'] < -1).sum()
print(f"\n๐Ÿ“Š Drought Statistics:")
print(f"Total days with severe drought (index < -1): {drought_days}")
print(f"Percentage of time in drought: {drought_days/len(df_mindanao)*100:.1f}%")

๐Ÿง  Module 2: Understanding RNNs and Their Limitations

Recurrent Neural Networks (RNNs) Basics

RNNs are designed to work with sequential data by maintaining a hidden state that acts as memory:

\[h_t = \tanh(W_{hh} \cdot h_{t-1} + W_{xh} \cdot x_t + b_h)\] \[y_t = W_{hy} \cdot h_t + b_y\]

Where: - \(h_t\): hidden state at time \(t\) - \(x_t\): input at time \(t\) - \(W_{hh}, W_{xh}, W_{hy}\): weight matrices - \(b_h, b_y\): bias terms

The Vanishing Gradient Problem

During backpropagation through time, gradients are multiplied repeatedly:

Example: If gradient = 0.5 at each step: - After 10 steps: \(0.5^{10} โ‰ˆ 0.001\) - After 50 steps: \(0.5^{50} โ‰ˆ 8.9 \times 10^{-16}\) (essentially zero!)

This means the network cannot learn long-term dependencies.

The Exploding Gradient Problem

Conversely, if gradient = 1.5 at each step: - After 10 steps: \(1.5^{10} โ‰ˆ 58\) - After 50 steps: \(1.5^{50} โ‰ˆ 6.4 \times 10^{8}\) (numerical overflow!)

This causes unstable training and NaN values.

๐ŸŽฏ Mini-Challenge 1

Task: Calculate how many time steps it takes for a gradient of 0.9 to shrink below 0.01. What does this mean for analyzing a year of monthly NDVI data?

# TODO: Mini-Challenge 1 - Calculate gradient vanishing
# Hint: Use a loop or logarithm to find when 0.9^n < 0.01

gradient_factor = 0.9
threshold = 0.01

# Your code here:
# steps = ???

# Uncomment to test:
# print(f"Gradient shrinks below {threshold} after {steps} steps")
# print(f"For monthly data, this means we can only learn patterns from the last {steps} months")

๐Ÿ—๏ธ Module 3: LSTM Architecture - The Solution

The LSTM Cell: A Smart Memory System

LSTMs solve the gradient problems through a sophisticated gating mechanism. Think of an LSTM as a smart student taking notes during a lecture:

  1. Forget Gate: Decides what old information to discard
  2. Input Gate: Determines what new information to store
  3. Output Gate: Controls what information to pass forward
  4. Cell State: The โ€œnotebookโ€ carrying information through time

Mathematical Formulation

The LSTM operations at time step \(t\):

Forget Gate: \(f_t = \sigma(W_f \cdot [h_{t-1}, x_t] + b_f)\)

Input Gate: \(i_t = \sigma(W_i \cdot [h_{t-1}, x_t] + b_i)\)

Candidate Values: \(\tilde{C}_t = \tanh(W_C \cdot [h_{t-1}, x_t] + b_C)\)

Cell State Update: \(C_t = f_t * C_{t-1} + i_t * \tilde{C}_t\)

Output Gate: \(o_t = \sigma(W_o \cdot [h_{t-1}, x_t] + b_o)\)

Hidden State: \(h_t = o_t * \tanh(C_t)\)

Where \(\sigma\) is the sigmoid function (outputs 0-1) and \(*\) is element-wise multiplication.

Why LSTMs Work

The cell state \(C_t\) acts as a โ€œconveyor beltโ€ that can carry information unchanged across many time steps. The gates (using sigmoid activation) produce values between 0 and 1, acting as โ€œvalvesโ€ that control information flow without causing gradient explosion or vanishing.

๐Ÿ’ญ Think-Through Discussion

Question: In drought monitoring, what kind of information might the โ€œforget gateโ€ discard and what might the โ€œinput gateโ€ preserve? Think about seasonal patterns vs. anomalies.

# Visualize LSTM Architecture Conceptually
def visualize_lstm_concept():
    """
    Create a conceptual visualization of LSTM information flow.
    """
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
    
    # Left plot: RNN vs LSTM gradient flow
    time_steps = np.arange(0, 50)
    rnn_gradient = 0.9 ** time_steps
    lstm_gradient = 0.95 ** time_steps  # LSTMs maintain gradients better
    
    ax1.plot(time_steps, rnn_gradient, 'r-', label='Vanilla RNN', linewidth=2)
    ax1.plot(time_steps, lstm_gradient, 'b-', label='LSTM', linewidth=2)
    ax1.axhline(y=0.01, color='gray', linestyle='--', alpha=0.5, label='Effective threshold')
    ax1.set_xlabel('Time Steps', fontsize=11)
    ax1.set_ylabel('Gradient Magnitude', fontsize=11)
    ax1.set_title('Gradient Flow: RNN vs LSTM', fontsize=12, fontweight='bold')
    ax1.set_yscale('log')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Right plot: LSTM gates behavior simulation
    time = np.linspace(0, 24, 100)  # 24 months
    
    # Simulate gate activations during drought event
    normal_period = (time < 6) | (time > 18)
    drought_period = (time >= 6) & (time <= 18)
    
    forget_gate = np.where(normal_period, 0.8, 0.3)  # Forget more during normal times
    input_gate = np.where(drought_period, 0.9, 0.4)   # Store more during drought
    output_gate = np.where(drought_period, 0.95, 0.6) # Output more during drought
    
    # Add some smooth transitions
    from scipy.ndimage import gaussian_filter1d
    forget_gate = gaussian_filter1d(forget_gate, sigma=2)
    input_gate = gaussian_filter1d(input_gate, sigma=2)
    output_gate = gaussian_filter1d(output_gate, sigma=2)
    
    ax2.plot(time, forget_gate, 'r-', label='Forget Gate', linewidth=2)
    ax2.plot(time, input_gate, 'g-', label='Input Gate', linewidth=2)
    ax2.plot(time, output_gate, 'b-', label='Output Gate', linewidth=2)
    ax2.axvspan(6, 18, alpha=0.2, color='orange', label='Drought Period')
    ax2.set_xlabel('Time (months)', fontsize=11)
    ax2.set_ylabel('Gate Activation (0-1)', fontsize=11)
    ax2.set_title('LSTM Gates During Drought Monitoring', fontsize=12, fontweight='bold')
    ax2.legend(loc='right')
    ax2.grid(True, alpha=0.3)
    ax2.set_ylim([0, 1])
    
    plt.tight_layout()
    plt.show()

visualize_lstm_concept()

print("๐Ÿ” Key Insights:")
print("1. LSTMs maintain gradient flow much better than vanilla RNNs")
print("2. During drought events, the input and output gates open more to capture and propagate anomaly information")
print("3. The forget gate reduces during drought to preserve important drought signals")

๐Ÿ”ง Data Preparation for LSTM

LSTMs require data in a specific format: sequences of fixed length. Weโ€™ll create sliding windows from our time series.

Sliding Window Approach

For drought prediction, weโ€™ll use: - Input: 12 months of historical NDVI values - Output: Next monthโ€™s NDVI value

Example:

Window 1: Months 1-12 โ†’ Predict Month 13
Window 2: Months 2-13 โ†’ Predict Month 14
Window 3: Months 3-14 โ†’ Predict Month 15
...
def create_sequences(data, seq_length=12, prediction_horizon=1):
    """
    Create sequences for LSTM training.
    
    Parameters:
    -----------
    data : np.array
        Time series data
    seq_length : int
        Number of time steps to use as input
    prediction_horizon : int
        Number of time steps ahead to predict
    
    Returns:
    --------
    X, y : np.arrays
        Input sequences and targets
    """
    X, y = [], []
    
    for i in range(len(data) - seq_length - prediction_horizon + 1):
        # Input sequence
        X.append(data[i:i + seq_length])
        # Target value
        y.append(data[i + seq_length + prediction_horizon - 1])
    
    return np.array(X), np.array(y)

# Prepare NDVI data
ndvi_values = df_mindanao['ndvi'].values.reshape(-1, 1)

# Normalize data (important for neural networks)
scaler = MinMaxScaler(feature_range=(0, 1))
ndvi_scaled = scaler.fit_transform(ndvi_values)

# Create sequences
SEQ_LENGTH = 12  # Use 12 time steps (120 days) to predict next time step
X, y = create_sequences(ndvi_scaled, seq_length=SEQ_LENGTH)

print(f"๐Ÿ“ฆ Data Shape:")
print(f"Input sequences (X): {X.shape}")
print(f"Target values (y): {y.shape}")
print(f"\nExample:")
print(f"First input sequence (scaled): {X[0].flatten()[:5]}... (showing first 5 values)")
print(f"Corresponding target: {y[0]}")

# Split into training and validation sets
split_index = int(0.8 * len(X))
X_train, X_val = X[:split_index], X[split_index:]
y_train, y_val = y[:split_index], y[split_index:]

print(f"\n๐Ÿ“Š Dataset Split:")
print(f"Training samples: {len(X_train)}")
print(f"Validation samples: {len(X_val)}")

๐ŸŽฏ Mini-Challenge 2: Data Exploration

Before we build the model, letโ€™s understand our sequences better.

# TODO: Mini-Challenge 2 - Analyze the sequences
# Task: Find and visualize a sequence that leads to a drought (low NDVI) prediction

# Hint: Look for sequences where y < 0.3 (scaled value)
# Your code here:

# drought_indices = np.where(y_train < 0.3)[0]
# if len(drought_indices) > 0:
#     drought_idx = drought_indices[0]
#     plt.figure(figsize=(10, 4))
#     plt.plot(range(SEQ_LENGTH), X_train[drought_idx], 'b-o', label='Input sequence')
#     plt.axvline(x=SEQ_LENGTH-0.5, color='gray', linestyle='--')
#     plt.plot(SEQ_LENGTH, y_train[drought_idx], 'ro', markersize=10, label='Predicted drought')
#     # Add your visualization code here

๐Ÿ›๏ธ Building the LSTM Model

Now letโ€™s build our LSTM model for drought prediction. Weโ€™ll use a architecture suitable for time series forecasting.

def build_lstm_model(seq_length, n_features=1, lstm_units=[64, 32], 
                    dropout_rate=0.2, learning_rate=0.001):
    """
    Build an LSTM model for time series prediction.
    
    Parameters:
    -----------
    seq_length : int
        Length of input sequences
    n_features : int
        Number of features per time step
    lstm_units : list
        Number of units in each LSTM layer
    dropout_rate : float
        Dropout rate for regularization
    learning_rate : float
        Learning rate for Adam optimizer
    
    Returns:
    --------
    model : keras.Model
        Compiled LSTM model
    """
    
    model = Sequential([
        # First LSTM layer with return sequences
        LSTM(lstm_units[0], 
             activation='tanh',
             return_sequences=True,  # Return full sequence for next LSTM
             input_shape=(seq_length, n_features),
             name='lstm_1'),
        
        # Dropout for regularization
        Dropout(dropout_rate, name='dropout_1'),
        
        # Second LSTM layer
        LSTM(lstm_units[1], 
             activation='tanh',
             return_sequences=False,  # Only return last output
             name='lstm_2'),
        
        # Dropout
        Dropout(dropout_rate, name='dropout_2'),
        
        # Dense layer for final prediction
        Dense(16, activation='relu', name='dense_1'),
        
        # Output layer
        Dense(1, activation='linear', name='output')
    ])
    
    # Compile model
    optimizer = Adam(learning_rate=learning_rate)
    model.compile(
        optimizer=optimizer,
        loss='mse',
        metrics=['mae']
    )
    
    return model

# Build the model
model = build_lstm_model(
    seq_length=SEQ_LENGTH,
    n_features=1,
    lstm_units=[64, 32],
    dropout_rate=0.2,
    learning_rate=0.001
)

# Display model architecture
model.summary()

# Visualize model architecture
tf.keras.utils.plot_model(
    model, 
    to_file='lstm_architecture.png',
    show_shapes=True,
    show_layer_names=True,
    dpi=100
)

print("\nโœ… Model built successfully!")
print(f"Total parameters: {model.count_params():,}")

๐Ÿš€ Training the LSTM Model

Weโ€™ll train our model with early stopping to prevent overfitting and learning rate reduction for better convergence.

Training Best Practices

  1. Early Stopping: Stop training when validation loss stops improving
  2. Learning Rate Reduction: Reduce LR when loss plateaus
  3. Batch Size: Balance between stability (large) and generalization (small)
  4. Monitoring: Track both training and validation metrics
# Define callbacks
early_stopping = EarlyStopping(
    monitor='val_loss',
    patience=15,
    restore_best_weights=True,
    verbose=1
)

reduce_lr = ReduceLROnPlateau(
    monitor='val_loss',
    factor=0.5,
    patience=5,
    min_lr=1e-6,
    verbose=1
)

# TODO: Complete the training code
# Train the model
EPOCHS = 100
BATCH_SIZE = 32

print("๐Ÿ‹๏ธ Training LSTM model...")
print(f"Epochs: {EPOCHS}")
print(f"Batch size: {BATCH_SIZE}")
print(f"Training samples: {len(X_train)}")
print(f"Validation samples: {len(X_val)}\n")

# TODO: Uncomment and complete the training
# history = model.fit(
#     X_train, y_train,
#     epochs=???,
#     batch_size=???,
#     validation_data=(???, ???),
#     callbacks=[???, ???],
#     verbose=1
# )

# print("\nโœ… Training completed!")

๐Ÿ“Š Training Visualization and Analysis

def plot_training_history(history):
    """
    Visualize training and validation metrics.
    """
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
    
    # Plot loss
    ax1.plot(history.history['loss'], label='Training Loss', linewidth=2)
    ax1.plot(history.history['val_loss'], label='Validation Loss', linewidth=2)
    ax1.set_xlabel('Epoch')
    ax1.set_ylabel('Loss (MSE)')
    ax1.set_title('Model Loss During Training')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Plot MAE
    ax2.plot(history.history['mae'], label='Training MAE', linewidth=2)
    ax2.plot(history.history['val_mae'], label='Validation MAE', linewidth=2)
    ax2.set_xlabel('Epoch')
    ax2.set_ylabel('MAE')
    ax2.set_title('Mean Absolute Error During Training')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Print final metrics
    final_train_loss = history.history['loss'][-1]
    final_val_loss = history.history['val_loss'][-1]
    final_train_mae = history.history['mae'][-1]
    final_val_mae = history.history['val_mae'][-1]
    
    print("๐Ÿ“ˆ Final Training Metrics:")
    print(f"Training Loss: {final_train_loss:.6f}")
    print(f"Validation Loss: {final_val_loss:.6f}")
    print(f"Training MAE: {final_train_mae:.6f}")
    print(f"Validation MAE: {final_val_mae:.6f}")
    
    if final_val_loss > final_train_loss * 1.5:
        print("\nโš ๏ธ Warning: Model might be overfitting!")
    else:
        print("\nโœ… Model generalization looks good!")

# TODO: Uncomment after training
# plot_training_history(history)

๐Ÿ”ฎ Model Evaluation and Predictions

Letโ€™s evaluate our modelโ€™s performance on the validation set and make predictions.

# TODO: Complete the prediction code
# Make predictions on validation set
# y_pred_scaled = model.predict(X_val)

# Inverse transform to get actual NDVI values
# y_pred = scaler.inverse_transform(y_pred_scaled)
# y_val_actual = scaler.inverse_transform(y_val)

# Calculate metrics
# mse = mean_squared_error(y_val_actual, y_pred)
# mae = mean_absolute_error(y_val_actual, y_pred)
# r2 = r2_score(y_val_actual, y_pred)

# print("๐Ÿ“Š Validation Set Performance:")
# print(f"Mean Squared Error: {mse:.6f}")
# print(f"Mean Absolute Error: {mae:.6f}")
# print(f"Rยฒ Score: {r2:.4f}")
# print(f"\nIn NDVI terms:")
# print(f"Average prediction error: ยฑ{mae:.3f} NDVI units")

๐Ÿ“ˆ Visualizing Predictions

def visualize_predictions(y_true, y_pred, dates=None, n_points=100):
    """
    Visualize actual vs predicted NDVI values.
    """
    fig, axes = plt.subplots(2, 1, figsize=(14, 8))
    
    # Plot 1: Time series comparison
    ax1 = axes[0]
    x_axis = range(len(y_true[:n_points]))
    
    ax1.plot(x_axis, y_true[:n_points], 'g-', label='Actual NDVI', linewidth=2, alpha=0.7)
    ax1.plot(x_axis, y_pred[:n_points], 'b--', label='Predicted NDVI', linewidth=2, alpha=0.7)
    
    # Highlight areas with large errors
    errors = np.abs(y_true[:n_points].flatten() - y_pred[:n_points].flatten())
    large_errors = errors > 0.1
    if np.any(large_errors):
        ax1.scatter(np.where(large_errors)[0], y_true[:n_points][large_errors], 
                   color='red', s=30, alpha=0.5, label='Large errors (>0.1)')
    
    ax1.set_xlabel('Time Steps')
    ax1.set_ylabel('NDVI')
    ax1.set_title('LSTM Predictions vs Actual NDVI Values', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Plot 2: Scatter plot
    ax2 = axes[1]
    ax2.scatter(y_true, y_pred, alpha=0.5, s=10)
    
    # Add perfect prediction line
    min_val = min(y_true.min(), y_pred.min())
    max_val = max(y_true.max(), y_pred.max())
    ax2.plot([min_val, max_val], [min_val, max_val], 'r--', alpha=0.5, label='Perfect prediction')
    
    ax2.set_xlabel('Actual NDVI')
    ax2.set_ylabel('Predicted NDVI')
    ax2.set_title('Prediction Accuracy Scatter Plot', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

# TODO: Uncomment after making predictions
# visualize_predictions(y_val_actual, y_pred)

๐ŸŒพ Module 4: LSTM Applications in Philippine Earth Observation

Real-World Applications

  1. Drought Forecasting (Mindanao)
    • Predict drought 1-3 months ahead
    • Enable early warning for farmers
    • Support irrigation planning
  2. Crop Yield Prediction
    • Combine NDVI with weather data
    • Forecast rice/corn yields
    • Support food security planning
  3. Phenology Analysis
    • Track cropping calendars
    • Detect planting/harvest dates
    • Monitor seasonal shifts due to climate change
  4. Land Cover Change Detection
    • Identify deforestation patterns
    • Monitor urban expansion
    • Track agricultural conversion

Practical Deployment Considerations

For operational use in Philippine agencies:

  1. Data Requirements
    • Minimum 2-3 years of historical data
    • Regular updates (weekly/bi-weekly)
    • Cloud-free observations critical
  2. Model Updates
    • Retrain quarterly with new data
    • Validate against ground truth
    • Account for seasonal variations
  3. Integration with Existing Systems
    • DOST-ASTI DATOS platform
    • PAGASA seasonal forecasts
    • PhilSA Space+ Dashboard

๐Ÿ’ญ Think-Through Discussion

Question: How would you modify this LSTM approach to integrate multiple data sources (e.g., Sentinel-1 SAR, Sentinel-2 optical, weather data) for improved drought prediction? What challenges might arise?

๐ŸŽฏ Mini-Challenge 3: Drought Alert System

Create a simple drought alert function based on LSTM predictions.

# TODO: Mini-Challenge 3 - Implement drought alert system
def drought_alert_system(predicted_ndvi, historical_mean=0.7, threshold_mild=0.6, 
                         threshold_severe=0.4):
    """
    Generate drought alerts based on predicted NDVI values.
    
    Parameters:
    -----------
    predicted_ndvi : float
        Predicted NDVI value
    historical_mean : float
        Historical mean NDVI for the location
    threshold_mild : float
        Threshold for mild drought alert
    threshold_severe : float
        Threshold for severe drought alert
    
    Returns:
    --------
    alert_level : str
        Alert level and message
    """
    
    # TODO: Implement the alert logic
    # Your code here:
    # if predicted_ndvi < threshold_severe:
    #     return "๐Ÿ”ด SEVERE DROUGHT ALERT: Immediate action required"
    # elif predicted_ndvi < threshold_mild:
    #     return "๐ŸŸก MILD DROUGHT WARNING: Monitor closely"
    # else:
    #     return "๐ŸŸข NORMAL CONDITIONS: No drought detected"
    
    pass

# Test the alert system
# test_values = [0.3, 0.5, 0.7, 0.8]
# for val in test_values:
#     alert = drought_alert_system(val)
#     print(f"NDVI: {val:.2f} โ†’ {alert}")

๐Ÿ”ฌ Advanced: Multi-Step Ahead Prediction

For operational drought forecasting, we often need to predict multiple time steps ahead.

def multi_step_prediction(model, initial_sequence, n_steps=3):
    """
    Predict multiple time steps into the future.
    
    Parameters:
    -----------
    model : keras.Model
        Trained LSTM model
    initial_sequence : np.array
        Initial sequence to start predictions from
    n_steps : int
        Number of steps to predict ahead
    
    Returns:
    --------
    predictions : list
        List of predictions for each time step
    """
    
    current_sequence = initial_sequence.copy()
    predictions = []
    
    for _ in range(n_steps):
        # Predict next time step
        next_pred = model.predict(current_sequence.reshape(1, -1, 1), verbose=0)
        predictions.append(next_pred[0, 0])
        
        # Update sequence: remove first element, add prediction
        current_sequence = np.append(current_sequence[1:], next_pred)
    
    return predictions

# TODO: Test multi-step prediction
# Example usage (uncomment after training):
# test_sequence = X_val[0]
# multi_predictions = multi_step_prediction(model, test_sequence, n_steps=3)
# 
# print("Multi-step ahead predictions:")
# for i, pred in enumerate(multi_predictions, 1):
#     pred_actual = scaler.inverse_transform([[pred]])[0, 0]
#     print(f"Step {i}: {pred_actual:.3f}")

๐Ÿ“‹ Key Takeaways

What Weโ€™ve Learned

  1. Time Series in EO: Critical for monitoring environmental changes and predicting future conditions

  2. RNN Limitations: Vanilla RNNs suffer from vanishing/exploding gradients, limiting their ability to learn long-term dependencies

  3. LSTM Architecture: Gates (forget, input, output) and cell state enable learning of both short and long-term patterns

  4. Implementation:

    • Data preparation with sliding windows
    • Normalization is crucial
    • Early stopping prevents overfitting

  5. Philippine Applications: Drought monitoring in Mindanao is a critical use case with immediate practical value

Best Practices

โœ… DO: - Normalize your data before training - Use sufficient historical data (2+ years) - Validate predictions against ground truth - Consider ensemble approaches for operational systems - Account for data gaps and cloud cover

โŒ DONโ€™T: - Ignore seasonal patterns in your data - Use too short sequences (< 6 time steps) - Deploy without thorough validation - Forget to retrain with new data periodically

Next Steps

  1. Experiment: Try different sequence lengths and LSTM architectures
  2. Enhance: Add weather data and other indices (EVI, NDWI)
  3. Scale: Apply to your area of interest using Google Earth Engine
  4. Integrate: Connect with Philippine EO platforms (DATOS, Space+)

๐Ÿš€ Your Turn!

Complete the TODOs in this notebook and experiment with: - Different LSTM architectures (more layers, different units) - Alternative optimizers (RMSprop, SGD) - Bidirectional LSTMs - Attention mechanisms - Multi-variate inputs (NDVI + temperature + precipitation)

๐Ÿ“š References and Further Reading

Scientific Papers

  1. Hochreiter, S., & Schmidhuber, J. (1997). โ€œLong Short-Term Memory.โ€ Neural Computation.
  2. RuรŸwurm, M., & Kรถrner, M. (2020). โ€œSelf-attention for raw optical Satellite Time Series Classification.โ€ ISPRS.
  3. Interdonato, R., et al. (2019). โ€œDuPLO: A DUal view Point deep Learning architecture for time series classificatiOn.โ€ ISPRS.

Philippine EO Resources

Tutorials and Documentation

End of Session 1: LSTMs for Earth Observation Time Series

Proceed to Session 2: Foundation Models and Transfer Learning ๐Ÿš€