Multi-Temporal Land Cover Classification with Sequential Recurrent Encoders

In plain words

A satellite sees the same field again and again over a season. Crops change how they reflect light as they grow, so the sequence of pictures tells you more about what is planted than any single picture does. Most land-cover methods throw that sequence away and work from one picture, or they first need clouds removed and the images corrected for the atmosphere.

This paper feeds the whole sequence into a recurrent network, the kind that reads one item after another and keeps a running memory. The network is convolutional, so it works on small image patches rather than single pixels. It reads the pictures forwards and backwards, so that neither the start nor the end of the season dominates. At the end it has one summary of the whole season for each spot on the map, and a small final layer turns that summary into a crop label for every pixel.

The authors also looked inside the network and found a few memory cells that react to clouds. The network seems to have learned to ignore cloudy pictures on its own, which would remove a preprocessing step. That observation rests on visual inspection of a few cells, not on a measured test.

TL;DR

A bidirectional convolutional recurrent encoder (ConvGRU as the main model, ConvLSTM also trained) reads a Sentinel-2 time series, and a convolutional softmax layer classifies every pixel into one of 17 crop classes §3.1 §3.3 fig. 3. Training uses a random subset of 30 observations per step, sampled from up to 52 available dates §5.2. Inputs are uncorrected top-of-atmosphere reflectance with clouds still present §3.3 §7. Author-reported: overall accuracy 89.7% (2016) and 89.5% (2017), overall kappa 0.870 for both years, and class recall between 41.5% (peas, 2017) and 96.8% (maize, 2017) tab. 2. Training used both seasons, so these numbers do not measure generalisation to an unseen year. The released code differs from the paper in several defaults and does not compute kappa (see Paper vs code).

Problem

Land-cover classification from optical satellites has mostly used a single date, or a fixed set of index features such as NDVI, NDWI and EVI, which keep only a few spectral bands §1 §6. Vegetation changes its spectral reflectance with phenological events, so a time series carries class information that single dates lack §1. Working with the full time series was hard because of data availability, storage, processing cost and the design of features across many dates §1. Cloud-free composites and atmospheric correction were usually required as preprocessing §2 §7. Earlier recurrent approaches by the same authors needed a per-pixel setup with nine neighbouring pixels packed into one vector, gave only point predictions, and needed extra cloud classes §3.2.

Method

The model is a sequence encoder in the style of sequence-to-sequence learning: an input sequence of variable length is compressed into one fixed-size state, which here is classified §3.1. The recurrent cells are convolutional, so the state is a map with spatial size h × w and r channels, not a vector §3.1 tab. 1.

1. Sample construction. Each sample is a tile of 240 m (24 × 24 pixels at 10 m) or 480 m §4. Per observation, the network input has 15 channels: four 10 m bands (B4, B3, B2, B8), six 20 m bands (B5, B6, B7, B8A) and three 60 m bands (B1, B11, B12), plus two constant maps for year and day-of-year §4. The 20 m and 60 m bands are bilinearly resampled to 10 m §4. Only observations with cloud cover below 80% are included; no cloud mask is applied §4.
2. Padding. A tile may have fewer than T observations. Missing slots are filled with a constant and ignored by the encoder §3.3. The paper does not state the constant.
3. Temporal subsampling (training only). At each training step, 30 observations are drawn at random from all available dates, while their temporal order is kept §5.2.
4. Forward and reverse encoding. The same recurrent cell is run over the observations in sequence order and again in reverse order, with shared weights and zero initial states §3.3 fig. 3. Each pass produces a final cell state. The forward final state is the state after the last observation; the reverse final state is the state after the first observation §3.3.
5. Concatenation. The two final states are concatenated along the channel axis, giving 2r channels §3.3 fig. 3.
6. Per-pixel classification. A convolution with kernel k_class = 3 maps the 2r-channel state to n = 17 class channels, followed by batch normalisation and a ReLU or leaky ReLU, and then a softmax over the classes §3.3 §5.2.
7. Loss. The cross-entropy between the softmax output and the one-hot label is computed for each training pixel §3.3 eq. 1.

Cell update rules (convolutions with kernel k_rnn = 3; the notation is the paper's, with conv() denoting a convolution with the cell's own weights) tab. 1:

# One ConvGRU step  [tab. 1]
def convgru_step(x, h_prev):
    # Reset and update gates read the input and the previous state  [tab. 1]
    stacked = concat(x, h_prev)
    reset = sigmoid(conv(stacked, W_r))
    update = sigmoid(conv(stacked, W_u))
    # Candidate uses the reset-gated previous state  [tab. 1]
    candidate = conv(concat(x, reset * h_prev), W_h)
    # Blend previous state and candidate; update gate weights h_prev  [tab. 1]
    return update * h_prev + (1 - update) * tanh(candidate)

# One ConvLSTM step  [tab. 1]
def convlstm_step(x, h_prev, c_prev):
    stacked = concat(x, h_prev)
    # Forget gate has +1 added before the sigmoid  [tab. 1]
    forget = sigmoid(conv(stacked, W_f) + 1)
    input_gate = sigmoid(conv(stacked, W_i))
    # Table 1 writes sigmoid here; §5.1 says the range is [-1, 1]  [§5.1]
    modulation = sigmoid(conv(stacked, W_j))
    output_gate = sigmoid(conv(stacked, W_o))
    c = c_prev * forget + input_gate * modulation  # cell state  [tab. 1]
    h = output_gate * tanh(c)  # [tab. 1]
    return h, c

Full encoder and classifier, from the description above §3.3 fig. 3:

# Bidirectional sequential encoder with per-pixel classifier  [§3.3]
# Input:  images (T, h, w, 15), n valid observations (n <= T)
# Output: class logits (h, w, 17)

def forward(images, n):
    # Forward pass over the valid observations  [§3.3]
    state_fwd = zeros(h, w, r)
    for t in range(n):
        state_fwd = cell_step(images[t], state_fwd)  # [tab. 1]
    # Reverse pass; the same cell weights are reused  [§3.3]
    state_rev = zeros(h, w, r)
    for t in reversed(range(n)):
        state_rev = cell_step(images[t], state_rev)
    # Concatenate the two final states  [§3.3] [fig. 3]
    representation = concat(state_fwd, state_rev)  # (h, w, 2r)
    # Conv k=3, then batch norm and leaky ReLU  [§3.3] [§5.2]
    scores = leaky_relu(batch_norm(conv(representation, 3, 17)))
    return scores

# Training loss for one step  [§3.3] [eq. 1]
def loss(scores, one_hot_labels, labelled):
    # Softmax over the class channel
    probs = softmax(scores, axis=class)
    # Per-pixel cross-entropy  [eq. 1]
    per_pixel = -sum(one_hot_labels * log(probs), axis=class)
    # Average over labelled pixels; the paper's text does not state the mask
    return mean(per_pixel[labelled])

Cross-entropy, as the paper writes it eq. 1:

H(\hat{\mathbf{y}}, \mathbf{y}) = -\sum_i y_i \log \hat{y}_i eq. 1

Tunable hyperparameters are the number of recurrent cells r, the recurrent kernel size k_rnn, and the classification kernel size k_class §3.3.

Hyperparameters

Results

Author-reported. Setup: Sentinel-2 top-of-atmosphere time series over one study area north of Munich, 17 crop classes, a convolutional GRU encoder with 256 cells trained for 60 epochs on both 2016 and 2017 §4 §5.2 tab. 2. Table 2 reports precision, recall, F-measure, conditional kappa and pixel counts per class and year. Precision is the user's accuracy and recall the producer's accuracy tab. 2. Values are copied as printed; the source's bold highlighting is omitted.

tab. 2. Caption: pixel-wise accuracies after 60 epochs, GRU encoder, on data of both growth seasons.

Overall accuracy is 89.7% for 2016 and 89.5% for 2017. Overall kappa is 0.870 for both years tab. 2.

Headline results, all author-reported:

• Weighted precision, recall and F-measure lie between 89.3% and 89.9% across both years tab. 2 §5.2.
• Class recall ranges from 41.5% (peas, 2017) to 96.8% (maize, 2017). The text calls these "class-wise accuracies" §5.2 tab. 2.
• Conditional kappa ranges from 0.414 (peas, 2017) to 0.957 (rapeseed, 2017). These extremes come from the 2017 column; the 2016 minimum is 0.466 (rye) tab. 2 §5.2.
• Frequent classes tend to score higher than rare ones, with exceptions such as peas in 2016 §5.2.
• Confusions that appear in only one season, such as soybeans with potato in 2017, are attributed to seasonal factors such as the timing of phenology and local weather §5.2.
• Qualitative examples A–D are good classifications; E and F show misclassifications, and in F the activations are blurred and low in amplitude fig. 8 §5.3.
• Internal activations: of 256 cells, three are shown. Cell 47 has input and modulation gates that change on cloudy observations; the authors infer cloud-related cells from this visual evidence fig. 6 §5.1.
• Comparison with other crop-classification studies is tabulated tab. 3, but the authors state that numerical comparison is difficult because the data, areas and class counts differ §6.

Paper vs code

Official code: TUM-LMF/MTLCC, MIT licence (2018, Marc Rußwurm), at commit f6e5aed5c59af084a91428fdd285a17fcf6344f4. The paper links it in its supplementary note §7; the arXiv version history (v4) records the updated link. The repository owner is the authors' lab, and the first author is a listed author of the paper.

• Recurrent cell count. The paper uses r = 256 for the final evaluation §5.2. The code's default is 24 seqencmodel.py, and the repository's examples use 128 readme.md. The configuration that produced Table 2 is not in the repository.
• Epochs and batch size. The paper uses 60 epochs with batch size 32 (GRU) or 28 (LSTM) §5.2. The training script defaults to batch size 16 train.py and 1 epoch train.py. The example command uses 30 epochs and batch size 4 readme.md readme.md. These are demo settings, not the paper's.
• Optimizer. The paper does not state one. The code uses Adam with learning rate 0.01, β1 0.9, β2 0.999 and ε 0.9 seqencmodel.py seqencmodel.py. The value ε = 0.9 is far from the common default, and the paper does not confirm it was used for the reported runs.
• Masking of unlabelled pixels. The active text of the paper does not mention this. The code removes class-0 (unknown) pixels from the loss seqencmodel.py seqencmodel.py and from the metrics seqencmodel.py evaluate.py.
• Classification output. The paper describes a convolution with batch norm and ReLU or leaky ReLU, then softmax §3.3. The code applies the batch norm and activation before the softmax, so the logits are non-negative for ReLU or slightly negative for leaky ReLU seqencmodel.py seqencmodel.py. The code default is leaky ReLU seqencmodel.py.
• Convolution details. The code uses reflection padding and no bias term in its convolutions seqencmodel.py seqencmodel.py. The paper does not mention either. Batch norm decay is 0.9 seqencmodel.py.
• Cell implementation. The ConvGRU and ConvLSTM cells come from the git submodule utils/convrnn (MarcCoru/tensorflow-convlstm-cell) .gitmodules. That submodule is not part of this clone, so the cell equations of Table 1 could not be checked against the code. The code calls the cells with batch normalisation on by default and no recurrent dropout seqencmodel.py seqencmodel.py, neither of which the paper mentions.
• Shared weights and state selection. The same cell object is passed to both directions seqencmodel.py, which matches the paper's shared weights §3.3. For LSTM, the code uses the cell state c as the representation seqencmodel.py, as the paper does §3.3. For GRU, the representation is the single GRU state.
• Sequence length. The code counts an observation as present when the first pixel of band B4 is positive seqencmodel.py. This assumes that missing observations are stored at the end and have non-positive values. The paper does not state the padding value.
• Normalisation. Reflectance is multiplied by 1e-4, day-of-year is divided by 365 and year minus 2016 is used Dataset.py. The paper does not state any scaling.
• Augmentation. The code implements random up-down and left-right flips Dataset.py, but the dataset object is created with augmentation off by default Dataset.py, and the training script does not enable it train.py. The paper does not mention augmentation.
• Validation during training. The training script passes the same temporal subsampling to every partition train.py Dataset.py. Validation accuracy during training is therefore measured on random 30-observation subsets. The evaluation script uses all observations evaluate.py. The paper does not say which is used for validation.
• Metrics. The released scripts compute precision, recall and F1 with scikit-learn, averaged with class weights evaluate.py. They do not compute the conditional or overall kappa reported in Table 2. Those values cannot be reproduced from this code.
• Evaluation partition and fold. Evaluation reads the eval partition at fold 0 evaluate.py. Training defaults to fold 0 train.py. The paper's active text does not name the fold.
• Hardware and framework. The paper uses a P100 GPU §5.2. The code targets TensorFlow 1.7 and Python 2.7 and was developed in TensorFlow 1.4 readme.md readme.md.

Pitfalls

• Modulation gate activation. Table 1 writes a sigmoid for the LSTM modulation gate, but §5.1 describes it as ranging over [-1, 1] tab. 1 §5.1. A tanh is the consistent reading. The paper does not clarify.
• GRU update convention. The GRU blend in Table 1 gives weight u to the previous state and (1 − u) to the candidate tab. 1. Some GRU formulations swap these, so check against the cell you use.
• Forget bias. The forget gate of the LSTM has +1 inside the sigmoid tab. 1. Omitting it changes early training.
• Padding and subsampling. Missing observations are padded with a constant that is not stated. Subsampling picks random slots among all stored slots Dataset.py, so the number of real observations in a training sample varies. Trailing padding is needed for the length count to work seqencmodel.py.
• Unlabelled pixels. The loss and the metrics ignore class 0 in the code seqencmodel.py. The paper's text does not state this; without it, unlabelled areas would contribute to the loss.
• Input scale. Missing from the paper: the reflectance scaling of 1e-4 and the day-of-year and year scalings Dataset.py. Without them, inputs will have a very different range.
• Batch norm in training versus evaluation. The code switches batch norm to its stored statistics at evaluation time evaluate.py. Evaluate with the same mode, or the results will change.
• Same-year training. Both seasons are used for training, and Table 2 reports results on held-out blocks of the same two seasons §4 §5.2. The paper's own note on temporal overfitting, where a network trained on one year fails on the other, is in commented-out text. The results therefore do not show cross-year transfer.
• Cloud claim. The statement that the network "learned to internally filter clouds" is supported by visual inspection of a few cells on one tile, not by a quantitative cloud test §5.1.
• Spatial offset. Some Sentinel-2 images are offset by about one pixel. The authors did not correct this and treated it as noise §4.
• Compute cost. Training ran for about 51–58 h on one P100 GPU §5.2.
• Source inconsistencies.
  • Section 4 gives the observation count as "98 (46 in 2017; 52 in 2017)"; §5.2 says 46 in 2016 and 52 in 2017. The 2016 figure is 46.
  • Section 5.2 says the GRU run took 58 h and the LSTM 51 h, then says the GRU was seven hours faster because of its larger batch size.
  • Text refers to Figure 2a and 2b through an undefined label fig2 §3.1 §3.2. Use the subfigure captions instead.
• Unidentified prior work. The authors' earlier crop-classification study §3.2 is cited only by key in the source, and it is not identified on this card.

Lineage

• Builds on: sequence-to-sequence encoders, attention and convolutional recurrent cells (ConvLSTM, ConvGRU) cited in §3.1. The change here is using the encoder's final state as a whole-series representation for per-pixel classification, and running the sequence in both directions §3.3.
• Earlier version: the authors' earlier recurrent crop-classification study, which used point-wise input and needed cloud classes §3.2.
• Compared against: the approach cited as Hao15 (random forest on MODIS phenological features from NDVI and NDWI) and the one cited as Siachalou2015 (hidden Markov model on Landsat and RapidEye) §6 tab. 3.