The Challenge: Monitoring 3.5 Million Hectares

Indonesia is the world's third-largest rice producer. Java Island alone has millions of hectares of paddy fields, and the government needs to know: how much rice is planted, how much is growing, and when will it be harvested?

Sending people to every field is impossible. Optical satellites (like regular cameras) are blocked by clouds — and Java is cloudy most of the year. The solution? SAR radar from the Sentinel-1 satellite.

How Radar "Sees" Rice Growth

The satellite transmits a vertical (V) microwave pulse and records the horizontal (H) echo — called VH backscatter. The strength of this echo (measured in decibels) changes as rice grows:

The 6 Growth Stages

The system classifies every pixel into one of these stages. Think of them as chapters in a rice plant's biography:

Check Your Understanding

Slicing a Year into 31 Snapshots

Imagine a time-lapse camera that photographs every rice field on Java once every 12 days. That's essentially what Sentinel-1 does. The system divides each year into 31 periods of 12 days each:

Looking Backward Through Time

When the system predicts the growth stage for Period 15, it doesn't just look at one snapshot. Like a doctor reviewing your medical history before making a diagnosis, it examines the current period plus the 6 preceding periods — 7 snapshots total.

How the Code Gets the Right Bands

The GeoTIFF file has 31 bands (one per period). The code uses 0-based indexing internally, meaning Band 1 = Index 0, Band 15 = Index 14.

for i in range(n_periods):
    band_idx = current_period - 1 - i

    if band_idx < 0:
        return None

    value = float(s1_data[band_idx, row_idx, col_idx])
    if np.isnan(value) or np.isinf(value):
        return None

    vh_values.append(value)

Source: utils.py lines 199-212

Check Your Understanding

The 29 Features: 4 Categories

Think of feature engineering as giving the model different lenses to examine the same data. Each category reveals something different about the rice growth pattern:

Building the Differences and Ratios

After extracting the 7 raw VH values, the code computes how the signal changes between each pair of consecutive periods:

# Temporal differences
for i in range(len(vh_values)-1):
    diff = vh_values[i] - vh_values[i+1]
    features.append(diff)

# Temporal ratios with safe division
for i in range(len(vh_values)-1):
    denominator = abs(vh_values[i+1])
    if denominator < 1e-10:
        ratio = 0.0
    else:
        ratio = vh_values[i] / denominator
    features.append(ratio)

Source: utils.py lines 220-236

Phenology Detection: Expert Knowledge in Code

The most interesting features come from phenology detection — where domain expertise about rice growth is encoded as rules. The code looks at the VH values and asks: "Does this pattern look like flooding? Like peak biomass? Like a harvest?"

# Flooding: very low VH signal
if vh_smooth[vh_min_idx] < -2000:
    phenology['flooding_detected'] = True

# Early Vegetative: rising after flood
if phenology['flooding_detected'] \
   and np.mean(recent_slopes) > 200:
    if -2000 <= current_vh <= -1700:
        phenology['early_vegetative'] = True

# Post-harvest: sharp decline
if slopes[-1] < -200:
    phenology['post_harvest'] = True

Source: utils.py lines 316-353

Research Finding: Which Features Matter Most?

An experiment tested whether absolute VH values or differential features (differences/ratios) contribute more to accuracy. The results were surprising:

Check Your Understanding

Architecture 1: The Multi-Layer Perceptron

The base model is an MLP — a stack of dense layers that progressively compress the 29 input features into a 6-class prediction. Imagine a funnel: wide at the top (512 neurons), narrow at the bottom (6 output classes).

The Model in Code

model = Sequential([
    Input(shape=(input_shape,)),

    Dense(512, kernel_initializer='he_normal'),
    BatchNormalization(),
    LeakyReLU(negative_slope=0.1),
    Dropout(0.3),

    Dense(num_classes),
    TemperatureLayer(temperature=0.5),
    Activation('softmax')
])

Source: utils.py lines 395-429

Architecture 3: CNN-LSTM — Reading Time Like a Story

The most advanced model treats the time series like a sentence. A CNN first scans for local patterns (like detecting a sudden jump or dip), then an LSTM reads the whole sequence to understand the narrative arc.

Check Your Understanding

From CSV to Trained Model

Training starts with a CSV file containing field observations: for each point, the date it was visited, its GPS coordinates, and the growth stage a human observer recorded (1-6).

5-Fold Cross-Validation: Trust but Verify

Instead of training once and hoping for the best, the system uses 5-fold cross-validation. It's like giving the model 5 different exams to ensure it genuinely learned, not just memorized.

The Imbalance Problem

Not all growth stages are equally represented in the training data. Stage 2 (Early Vegetative) has 6.3× more samples than Stage 5 (Late Generative). This is like teaching someone to identify 6 bird species using 500 photos of pigeons but only 80 photos of eagles.

The model tends to "cheat" by defaulting to Stage 2 for ambiguous cases — it's right 37.6% of the time just by guessing! Three strategies combat this:

Check Your Understanding

Processing Millions of Pixels in Chunks

Java Island at 10-meter resolution has hundreds of millions of pixels. Loading all of them at once would overwhelm any computer's memory. The solution is chunked processing — the system processes 20,000 pixels at a time.

Spatial Smoothing: Cleaning Up Noise

Raw pixel-wise predictions can be "noisy" — a single pixel might be classified as Stage 1 (Flooding) while everything around it is Stage 3 (Late Vegetative). That's probably an error. Median filtering fixes this by looking at each pixel's 3×3 neighborhood.

def apply_spatial_smoothing(predictions, kernel_size=3):
    valid_mask = (predictions != -32768)
    smoothed = predictions.copy()

    classes = np.unique(predictions[valid_mask])
    for class_val in classes:
        class_mask = (predictions == class_val)
        smoothed_mask = ndimage.median_filter(
            class_mask, size=kernel_size)
        smoothed[smoothed_mask & valid_mask] = class_val

    smoothed[~valid_mask] = -32768
    return smoothed

Source: utils.py lines 503-515

The 10-50x Speed Optimization

The original code extracts features one pixel at a time using a Python loop. The optimized version uses vectorized NumPy operations to process all pixels simultaneously.

# Loop over each pixel one at a time
for idx in range(n_pixels):
    row, col = valid_indices[idx]
    vh_values = []
    for i in range(7):
        vh_values.append(s1_data[band-i, row, col])
    features[idx] = extract_features(vh_values)
# Time: ~500ms per million pixels

# Process ALL pixels simultaneously
vh_stack = s1_data[band_indices]
vh_values = vh_stack[:, rows, cols].T

vh_diff = vh_values[:, :-1] - vh_values[:, 1:]

denominator = np.abs(vh_values[:, 1:])
vh_ratio = vh_values[:, :-1] / denominator
# Time: ~10-50ms per million pixels

Temporal Filtering: Consistency Across Time

After spatial smoothing, an optional step enforces temporal consistency. If a pixel was classified as Stage 4 in both the previous two periods but suddenly jumped to Stage 1, that's likely an error. Temporal filtering uses confidence-weighted voting across 3 consecutive periods.

Check Your Understanding

The Accuracy Landscape

Overall accuracy sits around 64.6% — but that number hides dramatic per-class variation. Some stages are identified reliably; others are barely detected at all.

Root Cause: The Three-Stage Overlap

Stages 3, 4, and 5 have extremely similar VH values. The model can't distinguish between "plants growing quickly" (Stage 3), "plants at peak height" (Stage 4), and "plants starting to dry" (Stage 5) because the absolute backscatter values are nearly identical.

Optimization Roadmap

Based on the accuracy analysis, here are the highest-impact improvements ranked by effort and expected gain:

Final Assessment

You Made It!

You now understand how satellite radar, time series feature engineering, and neural networks combine to map rice growth stages across 3.5 million hectares of Java Island. From raw microwave echoes to 6-class growth stage maps — and the gaps that still need closing.