Article

Multi-Temporal Land Cover Classification with
Sequential Recurrent Encoders
Marc Rußwurm 1
1

*

ID

and Marco Körner 1

ID

{marc.russwurm, marco.koerner}@tum.de; Technical University of Munich, Arcisstraße 21, 80333 Munich
Correspondence: marc.russwurm@tum.de; Tel.: +49-172-81-70-121

arXiv:1802.02080v3 [] 20 Feb 2018

submitted to IJGI January 22nd 2018; revised Febrary 20th 2018

Abstract: Earth observation (EO) sensors deliver data at daily or weekly intervals. Most land use and
land cover classification (LULC) approaches, however, are designed for cloud-free and mono-temporal
observations. The increasing temporal capabilities of today’s sensors enable the use of temporal,
along with spectral and spatial features. Domains such as speech recognition or neural machine
translation, work with inherently temporal data and, today, achieve impressive results by using
sequential encoder-decoder structures. Inspired by these sequence-to-sequence models, we adapt an
encoder structure with convolutional recurrent layers in order to approximate a phenological model
for vegetation classes based on a temporal sequence of Sentinel 2 (S2) images. In our experiments,
we visualize internal activations over a sequence of cloudy and non-cloudy images and find several
recurrent cells which reduce the input activity for cloudy observations. Hence, we assume that our
network has learned cloud-filtering schemes solely from input data, which could alleviate the need for
tedious cloud-filtering as a preprocessing step for many EO approaches. Moreover, using unfiltered
temporal series of top-of-atmosphere (TOA) reflectance data, our experiments achieved state-of-the-art
classification accuracies on a large number of crop classes with minimal preprocessing, compared to
other classification approaches.
Keywords: deep learning; multi-temporal classification; land use and land cover classification;
recurrent networks; sequence encoder; crop classification; sequence-to-sequence; sentinel 2

1. Introduction
Land use and land cover classification (LULC) has been a central focus of Earth observation (EO)
since first air- and spaceborne sensors began to provide data. For this purpose, optical sensors
sample the spectral reflectivity of objects on the earth’s surface in a spatial grid at repeated intervals.
Hence, LULC classes can be characterized by spectral, spatial, and temporal features. Today, most
classification tasks focus on spatial and spectral features, while utilizing the temporal domain had long
proven challenging. This is mostly due to limitations on data availability, the cost of data acquisition,
infrastructural challenges regarding data storage and processing, and the complexity of model design
and feature extraction over multiple timeframes.
Some LULC classes, such as urban structures, are mostly invariant to temporal changes and,
hence, are suitable for mono-temporal approaches. Others, predominantly vegetation-related classes,
change their spectral reflectivity based on biochemical processes initiated by phenological events related
to the type of vegetation and to environmental conditions. These vegetation-characteristic phenological
transitions have been utilized for crop yield prediction and, to some extent, for classification [1,2].
However, to circumvent previously mentioned challenges, the dimensionality of spectral bands has
often been compressed by calculating task-specific indices, such as normalized difference vegetation index
(NDVI), normalized difference water index (NDWI), or the enhanced vegetation index (EVI).
Today, most of these temporal data limitations have been alleviated by technological advances.
Reasonable spatial and temporal resolution data of multi-spectral earth observation sensors is available
ISPRS Int. J. Geo-Inf. 2018, , 1–19; doi:

www.mdpi.com/journal/ijgi

ISPRS Int. J. Geo-Inf. 2018, , 1–19

2 of 19

at no cost. Moreover, new services inexpensively provide high temporal and spatial resolution imagery.
The cost of data storage has decreased and data transmission has become sufficiently fast to allow
gathering and processing all available images over a large area and multiple years. Finally, new
advances in machine learning, accompanied by GPU accelerated hardware, have made it possible to
learn complex functional relationships, solely from provided data.
Since now data is available at high resolutions and processing is feasible, the temporal domain
should be exploited for EO approaches. However, this exploitation requires suitable processing
techniques utilizing all available temporal information at reasonable complexity. Other domains, such
as machine translation [3], text summarization [4–6], or speech recognition [7,8], handle sequential
data naturally. These domains have popularized sequence-to-sequence learning, which transforms
a variable-length input sequence to an intermediate representation. This representation is then
decoded to a variable-length output sequence. From this concept, we adopt the sequential encoder
structure and extract characteristic temporal features from a sequence of Sentinel 2 (S2) images using a
straightforward, two-layer network.
Thus, the main contributions of this work are
i) the adaptation of sequence encoders from the field of sequence-to-sequence learning to Earth
observation (EO),
ii) a visualization of internal gate activations on a sequence of satellite observations, and
iii) the application of crop classification over two seasons.
2. Related Work
As we aim to apply our network to vegetation classes, we first introduce common crop
classification approaches, to which we will compare our results in Section 6. Then, we motivate
data-driven learning models and cover the latest work on recurrent network structures in the EO
domain.
Many remote sensing approaches have achieved adequate classification accuracies for
multi-temporal crop data by using multiple preprocessing steps in order to improve feature separability.
Common methods are atmospheric correction [9–13], calculation of vegetation indices [9–13], or
the extraction of sophisticated phenological features [12]. Additionally, some approaches utilize
expert knowledge, for instance, by introducing additional agro-meteorological data [9], by selecting
suitable observation dates for the target crop-classes [13], or by determining rules for classification [10].
Pixel-based [9,12] and object-based [10,11,13] approaches have been proposed. Commonly, decision
trees (DTs) [9,10,13] or random forests (RFs) [11,12] are used as classifiers, of whose rules are sometimes
aided by additional expert knowledge [10].
These traditional approaches generally trade procedural complexity and the use of region-specific
expert knowledge for good classification accuracies in the respective areas of interest (AOIs). However,
these approaches are, in general, difficult to apply to other regions. Furthermore, the processing
structure requires supervision to varying degrees (e.g., product selection, visual image inspection,
parameter tuning), which impedes application at larger scales.
Today, we are experiencing a change in paradigm: away from the design of physically
interpretable, human-understandable models, which require task-specific expert knowledge, towards
data-driven models, which are encoded in internal weight parameters and derived solely from
observations. In that regard, hidden Markov models (HMMs) [14] or conditional random fields (CRFs) [15]
have shown promising classification accuracies with multi-temporal data. However, the underlying
Markov property limits long-term learning capabilities, as Markov-based approaches assume that the
present state only depends on the current input and one previous state.
In recent approaches, deep learning methods have proven particularly beneficial for modeling
physical relationship which are complicated, cannot be generalized, or are not well-understood [16].
Thus, deep learning is potentially well suited to approximate models of phenological changes, which
depend on complex internal biochemical processes of which only the change of surface reflectivity

ISPRS Int. J. Geo-Inf. 2018, , 1–19

3 of 19

Table 1. Update formulas of the convolutional variants of standard recurrent neural networks (RNNs),
long short-term memory (LSTM) cells and gated recurrent units (GRUs). A convolution between matrices a
and b is denoted by a ∗ b, element-wise multiplication by the Hadamard operator a b, and concatenation
on the last dimension is marked by [ akb]. The activation functions sigmoid σ( x ) and tangens
hyperbolicus tanh( x ) are used for non-linear scaling.
Gate

Variant
RNN

LSTM [25]

GRU [26]

h t ← x t , h t −1

h t , c t ← x t , h t −1 , c t −1

h t ← x t , h t −1

Forget/Reset

f t ← σ([ xt k ht−1 ] ∗ W f + 1)

rt ← σ([ xt k ht−1 ] ∗ Wr )

Insert/Update

it ← σ([ xt k ht−1 ] ∗ Wi )
jt ← σ([ xt kht−1 ] ∗ Wj )

ut ← σ([ xt k ht−1 ] ∗ Wu )

Output

ot ← σ([ xt k ht−1 ] ∗ Wo )

h̃t ← [ xt krt

c t ← c t −1
ht ← σ([ xt k ht−1 ] ∗ W )

h t ← ot

f t + it
tanh(ct )

ht−1 ] ∗ Wh̃

jt
ht ← ut

h t −1 + (1 − u t )

tanh(h̃t )

can be observed by EO sensors. A purely data-driven approach might alleviate the need to manually
design a functional model for this complex relationship. However, caution is required, as external
and non class-relevant factors, such as seasonal weather or observation configurations, are potentially
incorporated into the model, which might remain undetected if these factors constantly bias the
dataset.
In remote sensing, convolutional networks have gained increasing popularity for mono-temporal
observation tasks [17–20]. However, for sequential tasks, recurrent network architectures, which
provide an iterative framework to process sequential information, are generally better suited. Recent
approaches utilize recurrent architectures for change detection [21], identification of sea level anomalies
[22], and land cover classification [23]. For multi-seasonal, long-term dependencies, Jia et al. [21]
proposed a new cell architecture, which maintains two separate cell states for single- and multi-seasonal
long-term dependencies. However, the calculation of an additional cell state requires more weights,
which may prolong training and require more training samples.
In previous work, we have experimented with recurrent networks for crop classification [24]
and achieved promising results. Based on this, we propose a network structure using convolutional
recurrent layers and the aforementioned adaptation of a many-to-one classification scheme with
sequence encoders.
3. Methodology
Section 3.1 incrementally introduces the concepts of artificial neural networks (ANNs), feed-forward
networks (FNNs), and recurrent neural networks (RNNs) and illustrates the use of RNNs in
sequence-to-sequence learning. We then describe details of the proposed network structure in
Section 3.3.
3.1. Network Architectures and Sequential Encoders
Artificial neural networks (ANNs) approximate a function ŷ = f ( x; W ) of outputs ŷ (e.g., class
labels) from input data x given a large set of weights W. This approximation is commonly referred to as
the inference phase. These networks are typically composed of multiple cascaded layers with hidden
vectors h as intermediate layer outputs. Analogous to the biological neural cortex, single elements in
these vectors are often referred to as neurons. The quality of the approximation ŷ with respect to ground
truth y is determined by the loss function L(ŷ, y). Based on this function, gradients are back-propagated
through the ANN and adjust network weights W at each training step.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

4 of 19

xt

xt

h t −1

ht
f t it jt ot

h t −1

ht
rt

ut h̃t

a

a

[ akb]
c t −1

+

ct

+

a
a

b
concat

copy

Figure 1. Schematic illustration of long short-term memory (LSTM) and gated recurrent unit (GRU)
cells analog to the cell definitions in Table 1. The cell output ht is calculated via internal gates and based
on the current input xt combined with prior context information ht−1 , ct−1 . LSTM cells are designed to
separately accommodate long-term context in the internal cell state ct−1 , from short-term context ht−1 .
GRU cells combine all context information in a single, but more sophisticated output ht−1 .

Popular feed-forward networks often utilize convolutional or fully-connected layers at which the
input data is propagated through the network once. This is realized by an affine transformation
(fully-connected) h = σ (W x) or a convolution h = σ(W ∗ x) followed by an element-wise, non-linear
transformation σ : R 7→ R.
However, domains like translation [3], text summarization [4–6], or speech recognition [7,8]
formulate input vectors naturally as a sequence of observations x = { x0 , . . . , x T }. In these domains,
individual samples are generally less expressive and the overall model performance is based largely
on contextual information.
Sequential data is commonly processed with recurrent neural network (RNN) layers, in which the
hidden layer output ht is determined at time t by current input xt in combination with the previous
output ht−1 . In theory, the iterative update of ht enables RNNs to simulate arbitrary procedures [27],
since these networks are Turing complete [28]. The standard RNN variant performs the update step ht =
σ (W x̃) by an affine transformation of the concatenated vector x̃ = [ xt kht−1 ] followed by a non-linearity
σ. Consequently, the internal weight matrix is multiplied at each iteration step, which essentially
raises it to a high power [29]. At gradient back-propagation, this iterative matrix multiplication leads
to vanishing and exploding gradients [30,31]. While exploding gradients can be avoided with gradient
clipping, vanishing gradients impede the extraction of long-term feature relationships. This issue has
been addressed by Hochreiter and Schmidhuber [25], who introduced additional gates and an internal
state vector ct in long short-term memory (LSTM) cells to control the gradient propagation through time
and to enable long-term learning, respectively. Analog to standard RNNs, the output gate ot balances
the influence of the previous cell output ht−1 and the current input xt . At LSTMs, the cell output ht is
further augmented by an internal state vector ct , which is designed to contain long-term information.
To avoid aforementioned vanishing gradients, reading and writing to the cell state is controlled by
three additional gates. The forget gate f t decreases previously stored information by element-wise
multiplication ct−1 f t . New information is added by the product of input gate it and modulation gate
jt . Illustrations of the internal calculation can be seen in Figure 1, and the mathematical relations are
shown in Table 1. Besides LSTMs, gated recurrent units (GRUs) [26] have gained increasing popularity,
as these cells achieve similar accuracies to LSTMs with fewer trainable parameters. Instead of separate
vectors for long and short-term memory, GRUs formulate a single, but more sophisticated, output
vector.
To account for the more complicated design, recurrent layers are conventionally referred to as a
collection of cells with a single cell representing the set of elements at one vector-index.
The common output of recurrent layers provides a many-to-many relation by generating an output
vector at each observation ht given previous context ht−1 and ct−1 , as shown in Figure 2(a). However,

ISPRS Int. J. Geo-Inf. 2018, , 1–19

x1

5 of 19

xT

x2

I

live

in

Munich

...
0

intermediate representation c

0
y1

y2

...

yT
Ich

(a) Network structure employed in previous work [24].

lebe

in

München

(b) Illustration of a sequence-to-sequence network [7] as often
used in neural translation tasks.

Figure 2. Illustrations of recurrent network architectures which inspired this work. The network of
previous work [24] shown in Figure 2(a) creates a prediction yt at each observation t based on spectral
input information xt and the previous context ht−1 , ct−1 . Sequence-to-sequence networks, as shown in
Figure 2(b), aggregate sequential information to an intermediate state c T which is a representation of
the entire series.

encoding information of the entire sequence in a many-to-one relation is favored in many applications.
Following this idea, sequence-to-sequence learning, illustrated in Figure 2(b), has popularized the use
of the cell state vector c T at the last-processed observation T as a representation of the entire input
sequence. These encoding-decoding networks transform an input sequence of varying length to an
intermediate state representation c of fixed size. Subsequently, the decoder generates a varying length
output sequence from this intermediate representation. Further developments in this domain include
attention schemes. These provide additional intermediate connections between encoder and decoder
layers, which are beneficial for translations of longer sequences [3].
In many sequential applications, the common input form is xt ∈ Rd with a given depth d. The
output vectors ht ∈ Rr are computed by matrix multiplication with internal weights W ∈ R(r+d)×r and
r recurrent cells. However other fields, such as image processing, commonly handle raster data xt ∈
Rh×w×d of specific width w, height h, and spectral depth d. To account for neighborhood relationships and
to circumvent the increasing complexity, convolutional variants of LSTMs [32] and GRUs have been
introduced. These variants convolve the input tensors with weights W ∈ Rk×k×(r+d)×r augmented by
the convolutional kernel size k, which is a hyper-parameter determining the perceptive field.
3.2. Prior Work
Given recurrent networks as popular architectures for sequential data processing, we
experimented with recurrent layers for multi-temporal vegetation classification prior to this work [24].
In conducted experiments, we used a network architecture similar to the illustration in Figure 2(a).
Following the input dimensions of standard recurrent layers, an input sequence x ∈ { x0 , ..., x T } of
observations xt ∈ Rd was introduced to the network. Based on contextual information from previous
observations, a classification for each observation yt was produced. We evaluated the effect of this
information gain by comparing the recurrent network with convolutional neural networks (CNNs) and a
support vector machine (SVM). Standard RNNs and LSTMs, outperformed their non-sequential SVMs
and CNNs counterparts. Further, we observed an increase in accuracy at sequentially later observations,
which were classified with more context information available. Overall, we concluded that recurrent
network architectures are well suited for the extraction of temporal features from multi-temporal EO
imagery, which is consistent with other recent findings [21–23].
However, the experimental setup introduced some limitations regarding applicability in
real-world scenarios. We followed the standard formulation of recurrent networks which process
a d dimensional input vector. This vector included the concatenated bottom-of-atmosphere (BOA)

ISPRS Int. J. Geo-Inf. 2018, , 1–19

6 of 19

reflectances of nine pixels neighboring one point-of-interest. The point-wise classification was sufficient
for quantitative accuracy evaluation, but could not produce areal classification maps. Since a class
prediction was performed on every observation, we introduced additional covered classes for cloudy
pixels at single images. These were derived from the scene classification of the SEN 2 COR atmospheric
correction algorithm, which required additional preprocessing. A single representative classification
for the entire time-series would have required additional post-processing to further aggregate the
predicted labels for each observation. Finally, the mono-directional iterative processing introduced a
bias towards last observations. With more contextual information available, later observation showed
better classification accuracies compared to observations earlier in the sequence.
3.3. This Approach
To address the limitations of previous work, we redesigned and streamlined the network structure
and processing pipeline. Inspired by sequence-to-sequence structures described in Section 3 the
proposed network aggregates the information encoded in the cell state ct within the recurrent cell.
Since one class prediction for the entire temporal series is produced, atmospheric perturbations can be
treated as temporal noise. Hence, explicitly introduced cloud-related labels are not required, which
alleviates the need for prior cloud classification. Without the need for prior scene classification to obtain
these classes, the performance on atmospherically uncorrected top-of-atmosphere (TOA) reflectance
data can be evaluated. We further implemented convolutional recurrent cell variants, as formulated
in Table 1 to process input tensors xt of given height h, width w and depth d. Hence, the proposed
network produces areal prediction maps as shown in the qualitative results Section 5.3. Finally, we
introduce the input sequence in a bidirectional manner to eliminate any bias towards the later elements
in the observation sequence.
Overall, we employ a bidirectional sequential encoder for the task of multi-temporal land cover
classification. As earth observation data is gathered in a periodic manner, many observations of
the same area at consecutive times are available, which may contribute to the classification decision.
Inspired by sequence-to-sequence models, the proposed model encodes this sequence of images to
a fixed-length representation. Compared to previous work, this is an elegant way to condense the
available temporal dimension without further post-processing. A classification map for each class is
derived from this sequence representation. Many optical observations are covered by clouds and prior
cloud classification is often required as additional preprocessing step. As clouds do not contribute to
the classification decision, these observations can be treated as temporal noise and may be potentially
ignored by this encoding scheme. In Section 5.1 we investigate this by visualizing internal activation
states on cloudy and non-cloudy observations.
Figure 3 presents the proposed network structure schematically. The input image sequence
x = { xt , . . . , x T } of observations x ∈ Rh×w×d is passed to gated recurrent layers at each observation
time t. The index T denotes the maximum length of the sequence and d the input feature depth.
In practice, sequence lengths are often shorter than T, as the availability of satellite acquisitions is
variable over larger scales. If less than T observations are present, sequence elements are padded
with a constant value and are subsequently ignored at the iterative encoding steps. To eliminate bias
towards the last observations in the sequence, the data is passed to the encoder in both sequential
(seq) and reversed (rev) order. Network weights are shared between both passes. The initial cell states
seq
h×w×r and output hseq , hrev ∈ Rh×w×r are initialized with zeros. The concatenated final
c0 , crev
0
T ∈R
T
seq
states c T = [c T kcinv
]
are
the
representation
of the entire sequence and are passed to a convolutional
0
layer for classification. A second convolutional classification layer projects the sequence representation
c T to softmax-normalized activation maps ŷ for n classes: c T ∈ Rh×w×2r 7→ ŷ ∈ Rh×w×n . This layer is
composed of a convolution with a kernel size of kclass , followed by batch normalization and a leaky relu
[33] or relu [34] non-linear activation function.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

x0

x1

7 of 19

in sequence

xT

label

...
seq
h0
seq

c0

=0

H (y, ŷ)
seq

cT

=0

crev
0

crev
T =0

classification argmax

hrev
T =0

...

xT

x T −1

reversed

prediction

x0

Figure 3. Schematicial illustration of our proposed bidirectional sequential encoder network. The input
seq
sequence x ∈ { x0 , . . . , x T } of observations xt ∈ Rh×w×d is encoded to a representation c T = [c T kcinv
0 ].
The observations are passed in sequence (seq) and reversed (rev) order to the encoder to eliminate bias
towards recent observations. The concatenated representation of both passes c T is then projected to
softmax-normalized feature maps for each class using a convolutional layer.

At each training step, the cross-entropy loss
H (ŷ, y) = − ∑ ŷi log(yi )

(1)

i

between the predicted activations ŷ and an one-hot representation of the ground truth labels y evaluates
the prediction quality. Consequently, based on this loss function and the Adam optimizer [35],
gradients are back-propagated through the network layers and adjust the model weights.
Tunable hyper-parameters are the number of recurrent cells r and the sizes of the convolutional
kernel krnn and the classification kernel kclass .
4. Dataset
For the evaluation of our approach, we defined a large area of interest (AOI) of 102 km × 42 km
north of Munich, Germany. An overview of the AOI at multiple scales is shown in Figure 4. The
AOI was further subdivided into squared blocks of 3.84 km × 3.84 km (multiples of 240 m and 480 m)
to ensure dataset independence while maintaining similar class distributions. These blocks were
then randomly assigned to partitions for network training, hyper-parameter validation, and model
evaluation in a ratio of 4 : 1 : 1 similar to previous work [24]. The spatial extent of single samples
x are determined by tile-grids of 240 m and 480 m. We bilinearly interpolated the 20 m and 60 m S2
bands to 10 m ground sampling distance (GSD) to harmonize the raster data dimensions. To provide
additional temporal meta information, the year and day-of-year of the individual observations were
added as matrices to the input tensor. Hence, the input feature depth d = 15 is composed of four 10 m
(B4,B3,B2,B8), six 20 m (B5, B6, B7, B8A), and three 60 m (B1,B11, B12) bands combined with year
and day-of-year.
With ground truth labels of two growing seasons 2016 and 2017 available, we gathered 274 (108
in 2016; 166 in 2017) Sentinel 2 products at 98 (46 in 2017; 52 in 2017) observation dates between 3rd
January, 2016 and 15th November, 2017. The obtained time series represents all available S2 products

ISPRS Int. J. Geo-Inf. 2018, , 1–19

8 of 19

tile (240 m)

training

validation

evaluation

40 km

Berlin

margin (480 m)
102 km

3840 m blocks (14 240 m tiles, 480 m margin)

Munich

Figure 4. Area of interest (AOI) north of Munich containing 430 kha and 137 k field parcels. The AOI is
further tiled at multiple scales into datasets for training, validation and evaluation and footprints of
individual samples.

labeled with cloud coverage less than 80%. In some S2 images, we noticed a spatial offset in the scale
of one pixel. However, we did not perform additional georeferencing and treated the spatial offset
as data-inherent observation noise. Overall, we relied on the geometrical and spectral reference as
provided by C OPERNICUS ground segment.
Ground truth information was provided by the Bavarian Ministry of Food, Agriculture and Forestry
(StMELF) in form of geometry and semantic labels of 137 k field parcels. The crop-type is reported
by farmers to the ministry as mandated by the European crop subsidy program. We selected and
aggregated 17 crop-classes from approximately 200 distinct field labels, occurred at least 400 times in
the AOI. With modern agriculture, centered around few predominant crops, the distribution of classes
is not uniform, as can be observed from Figure 5(a). This non-uniform class distribution is generally
not optimal for the classification evaluation as it skews the overall accuracy metric towards classes of
high frequency. Hence, we additionally calculated kappa metrics [36] for the quantitative evaluation
in Section 5.2 to compensate for unbalanced distributions.
5. Results

field parcels

In this section, we first visualize internal state activations in Section 5.1 to gain a visual
understanding of the sequential encoding process. Further findings on internal cloud masking are
presented before the classification results on crop classes are quantitatively and qualitatively evaluated
in Sections 5.2 and 5.3.
30,000
10,000
3,000
1,000
300

2016

2017

e at w y e d y p le at ye et lt us s as s
aiz e o rle ato ee rle ho ca o r r be sperag beanpe bean
m wheadr bapot pesr ba triti
a
y
ga asp
m te
rame
so
su
n
m
wi
su
crop classes
(a) Non-uniform distribution of field classes in the AOI

S2B
S2A
01
2016

04

07

10

02

05

08

12

2017

(b) Acquired Sentinel 2 (S2) observations of the twin
satellites S2A and S2B

Figure 5. Information of the area of interest containing location, division schemes, class distributions
and dates of acquired satellite imagery.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

9 of 19

5.1. Internal Network Activations
In Section 3.1, we gave an overview of the functionality of recurrent layers and discussed
the property of LSTM state vectors ct ∈ Rh×w×r to encode sequential information over a series of
observations. The cell state is updated by internal gates it , jt , f t ∈ Rh×w×r , which in turn are calculated
based on previous cell output ht−1 and cell state ct−1 (see Table 1). To assess prior assumptions
regarding cloud filtering and to visually assess the encoding process, we visualized internal LSTM
cell tensors for a sequence of images and show representative activations of three cells in Figure 6.
The LSTM network, from which these activations are extracted, was trained on 24 px × 24 px tiles
with r = 256 recurrent cells and krnn = kclass = 3px. Additionally, we inferred the network with tiles
of height h and width w of 48 px. Experiments with the input size of 24 px show similar results and
are included in the complementary material to this work. At the first row, a 4σ band-normalized
RGB image represents the input satellite image xt ∈ Rh=48 × Rw=48 × Rd=15 at each time-frame t. The
next rows show the activations of input gate iti , modulation gate jti , forget gate f ti , and cell state cit at
three selected recurrent cells, which are denoted by the raised index i ∈ {3, 22, 47}. After iteratively
processing the sequence, the final cell state c T =36 is used to produce activations for each class, as
described in Section 3.3.
In the encoding process, the detail of structures at the cell state tensor increased gradually. This
may be interpreted as additional information written to the cell state. It further appeared that the
structures visible at the cell states resembled shapes, which were present in cloud-free RGB images
(3)

(22)

(e.g., ct=15 or ct=28 ). Some cells (e.g., cell 3 or cell 22) changed their activations gradually over the
span of multiple observations, while others (e.g., 48) changed more frequently. Since forget gate f
activations are element-wise multiplied with the previous cell state ct−1 and range between zero and
one. Low values in this gate numerically reduce the cell state which can be potentially interpreted
as a change of decision. The input i and modulation gate j control the degree of new information
written to the cell state. While the input gate is scaled between zero and one, the modulation gate
j ∈ [−1, 1] determines the sign of change. In general, we found the activity of a majority of cells (e.g.,
cell 3 or cell 22) difficult to associate with distinct events in the current input. However, we assumed
that classification-relevant features were expressed as a combination of cell activations similar to other
neural network approaches. Nevertheless, we could identify a proportionally small number of cells, in
which the shape of clouds visible in the image was projected on the internal state activations. One of
these was cell i = 47. At cloudy observations, the input gate approached zero either over the entire tile
(e.g., t = {10, 18, 19, 36}) or at patches of cloudy pixels (e.g., t = {11, 13, 31, 33}). At some observation
(47)

times (e.g., t = {13, 31, 32}) the modulation gate jt additionally changed the sign.
In a similar fashion, Karpathy [37] evaluated cell activations for the task of text processing. He
could associate a small number of cells to a set of distinct tasks, such as monitoring the lengths of a
sentence or maintaining a state-flag for text inside and outside of brackets.
Summarizing this experiment, the majority of cells showed increasingly detailed structures when
new information was provided in the input sequence. It is likely that the grammar of crop-characteristic
phenological changes was encoded in the network weights and we suspect that a certain amount of
these cells were sensitive to distinct events relevant for crop identification. However, these events
may be encoded in multiple cells and were difficult to visually interpret. A small set of cells could be
visually associated with individual cloud covers and may be used for internal cloud masking. Based
on these findings, we are confident that our network has learned to internally filter clouds without
explicitly introduced cloud-related labels.

t2

t3

t10

t4

t11

t12

t13

t14

t15

t16

t17

t18

t19

t20

t28

x

...

...

f (3)

...

...

i (3)

...

...

j (3)

...

...

c (3)

...

...

f (22)

...

...

i(22)

...

...

j(22)

...

...

c(22)

...

...

f (47)

...

...

i(47)

...

...

j(47)

...

...

c(47)

...

...

t29

t30

t31

t32

t33

t34

t35

t36

1

0
1

0
1

ISPRS Int. J. Geo-Inf. 2018, , 1–19

t1

-1
1

-1
1

0
1

0
1

-1
1

-1
1

0
1

0
1

-1
1

-1

10 of 19

Figure 6. Internal LSTM cell activations of input gate i(i) , forget gate f (i) , modulation gate j(i) and cell state c(i) at three (of r = 256) selected cells i ∈ {3, 22, 47} given
the current input xt over the sequence of observations t = {1, .., 36}. The detail of features at the cell states increased gradually, which indicated the aggregation of
information over the sequence. While most cells likely contribute to the classification decision, only some cells are visually interpretable with regard to the current
input xt . One visually interpretable cell i = 47 has learned to identify cloud, as input and modulation gates show different activation patterns on cloudy and
non-cloudy observations.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

11 of 19

5.2. Quanititative Classificaton Evaluation
For the quantitative evaluation of our approach, we trained networks with bidirectional
convolutional LSTM and GRU cells with r ∈ {128, 256} recurrent cells. Kernel sizes of krnn = kclass = 3
were used for the evaluation since previous tests with larger kernel sizes showed similar accuracies.
For these initial experiments, we predominantly tested network variants with r = 128 recurrent cells,
as these networks could be trained within a reasonable time-frame. We decided to use networks with
r = 256 recurrent cells for the final accuracy evaluation, as we found that these variants achieved
slightly better results in prior tests. The convolutional GRU and LSTM networks were trained on a
P100 GPU for 60 epochs (3.51m 24 px × 24 px tiles seen) and took 58 h and 51 h, respectively. However,
reasonable accuracies were achieved within the first twelve hours and further training increased the
accuracies on validation data only marginally. At each training step, a subset of 30 observations was
randomly sampled from all available 46 (2016) and 52 (2017) observations to randomize the sequence
while the sequential order is maintained. For all our tests, the performance of LSTM and GRU networks
was similar. The fewer weights of GRU cells, however, allowed using a slightly larger batch size of
32 samples compared to 28 samples of the LSTM variant. This lead to a seven-hour faster training
compared to the LSTM variant.
For these reasons, we decided to report evaluation results of the GRU network at Table 2. Precision
and recall are common accuracy measures which normalize the sum of correctly predicted samples
with the total number of predicted and reference samples of a given class, respectively. These measures
are equivalent to user’s and producer’s accuracies and inverse to errors of commission and omission,
which are popular metrics in the remote sensing community. We further calculated the f -measure as
the harmonic average of precision and recall and the overall accuracy, as the sum of correctly classified
samples normalized by the total number of samples. These metrics weight each sample equally. This
introduces a bias towards classes frequent in the dataset, such as maize or wheat. To compensate
for the non-uniform class distribution, we additionally report the conditional [38] and overall kappa
[36] coefficients, which are normalized by the probability of a hypothetical correct classification by
chance. The kappa coefficient κ is a measure of agreement and typically ranges between κ = 0 for
no- and κ = 1 for complete agreement. McHugh [39] provides an interpretative table in which values
0.4 ≤ κ < 0.6 are considered ’weak’, values 0.6 ≤ κ < 0.8 ’moderate’, 0.8 ≤ κ ≤ 0.9 considered ’strong’
and values beyond 0.9 ’almost perfect’.
The provided table of accuracies shows precision, recall, f -measure and the conditional kappa
coefficient for each class over the two evaluated seasons. Further overall accuracy and overall kappa
coefficients indicate the quality for the classification and report good accuracies. The pixel-averaged
achieved precision, recall, and f-score accuracies were consistent and ranged between 89.3 % and 89.9 %.
The kappa coefficients of 0.870 and overall accuracies of 89.7 % and 89.5 % show similar consistency.
While these classification measures reported good performances, the class-wise accuracies varied
largely between 41.5 % (peas) and 96.8 % (maize). For better visibility, we emphasized the best and
worst metrics by boldface. The conditional kappa scores are similarly variable and range between
0.414 (peas) and 0.957 (rapeseed ).
Frequent classes (e.g., maize, meadow) have been in general more confidently classified than less
frequent classes (e.g., peas, summer oat, winter spelt, winter triticale). Nonetheless this relation
has exceptions. The least frequent class, peas, performed relatively well on data of 2016 and
other less frequent classes, such as asparagus or hop, showed good performances despite their
underrepresentation in the dataset.
To investigate the causes of varying accuracies, we calculated confusion matrices for both seasons
as shown in Figure 7. These error matrices are two-dimensional histograms of classification samples
aggregated by class prediction and ground truth reference. To account for the non-uniform class
distribution, the absolute number of samples for each row-column pair is normalized. We decided to
normalize the confusion matrices by row to obtain recall (producer’s) accuracies, due to their direct
relation to available ground truth labels. The diagonal elements of the matrices represent correctly

ISPRS Int. J. Geo-Inf. 2018, , 1–19

12 of 19

Table 2. Pixel-wise accuracies of the trained convolutional GRU sequential encoder network after
training over 60 epochs on data of both growth seasons. The conditional kappa metrics [38] for each
class and the overall kappa [36] measure are given for both growth seasons.
Class

Year
2016

2017

Precision
(Users Acc.)

Recall
(Prod. Acc.)

f -Meas.

Kappa

#Pixels

Precision
(Users Acc.)

Recall
(Prod. Acc.)

f -Meas.

Kappa

#Pixels

sugar beet
oat
meadow
rapeseed
hop
spelt
triticale
beans
peas
potato
soybeans
asparagus
wheat
winter barley
rye
summer barley
maize

94.6
86.1
90.8
95.4
96.4
55.1
69.4
92.4
93.2
90.9
97.7
89.2
87.7
95.2
85.6
87.5
91.6

77.6
67.8
85.7
90.0
87.5
81.1
55.7
87.1
70.7
88.2
79.6
78.8
93.1
87.3
47.0
83.4
96.3

85.3
75.8
88.2
92.6
91.7
65.6
61.8
89.6
80.4
89.5
87.7
83.7
90.3
91.0
60.7
85.4
93.9

.772
.675
.845
.896
.873
.807
.549
.869
.706
.876
.795
.787
.902
.861
.466
.830
.944

59k
36k
233k
125k
51k
38k
65k
27k
9k
126k
21k
20k
806k
258k
43k
73k
919k

89.2
63.8
88.1
96.2
92.5
75.3
62.4
92.8
60.9
95.2
75.9
81.6
90.1
92.5
76.7
77.9
92.3

78.5
62.8
85.0
95.9
74.7
46.7
57.2
63.2
41.5
73.8
79.9
77.5
95.0
92.2
61.9
88.5
96.8

83.5
63.3
86.5
96.1
82.7
57.6
59.7
75.2
49.3
83.1
77.8
79.5
92.5
92.4
68.5
82.9
94.5

.779
.623
.837
.957
.743
.463
.563
.630
.414
.728
.798
.773
.930
.915
.616
.880
.953

94k
38k
242k
114k
53k
31k
64k
28k
6k
140k
26k
19k
783k
255k
30k
91k
876k

weight. avg

89.9

89.7

89.5

89.5

89.5

89.3

Overall Accuracy

Overall Kappa

Overall Accuracy

Overall Kappa

89.7

.870

89.5

.870

classified samples with values equivalent to Table 2. Structures outside the diagonal indicate systematic
confusions between classes and may give insight into the reasoning behind varying classification
accuracies.
Some crops likely share common spectral or phenological characteristics. Hence, we expected
some symmetric confusion between classes, which would be expressed as diagonal symmetric
confusions consistent in both years. Examples of this were triticale and rye, or oat and summer
barley. However, these relations were not frequent in the dataset, which indicates that the network
had sufficient capacity to separate the classes by provided features. In some cases, one class may
share characteristics with another class. This class may be further distinguished by additional unique
features, which would be expressed by asymmetric confusions between these two classes in both
seasons. Relations of this type were more dominantly visible in the matrices and included confusions
between barley and triticale, triticale and spelt, or wheat confused with triticale and spelt. These types
of confusion were consistent over both seasons and may be explained by a spectral or phenological
similarity between individual crop-types.
More dominantly, many confusions were not consistent over the two growing seasons. For
instance, confusions occurring only in the 2017 season were soybeans with potato, or peas with
meadow and potato. Since the cultivated crops are identical in these years and the class distributions
were consistent, seasonally variable factors were likely responsible for these relations. As reported in
Table 2, peas have been classified well in 2016, but poorly in 2017, due to aforementioned confusions
with meadow and potato. These results indicate that external and not crop-type related factors had a
negative influence on classification accuracies, which appeared unique to one season. One of these
might be the variable onset of phenological events, which are indirectly observed by the change of
reflectances by the sensors. These events are influenced by local weather and sun exposure which may
vary over large regional scales or multiple years.

13 of 19

2016

2017

predicted class

predicted class

sugar beet
oat
meadow
rapeseed
hop
spelt
triticale
beans
peas
potato
soybeans
asparagus
wheat
w. barley
rye
s. barley
maize

recall
100 %
80 %
60 %
40 %

sugar beet
oat
meadow
rapeseed
hop
spelt
triticale
beans
peas
potato
soybeans
asparagus
wheat
w. barley
rye
s. barley
maize

20 %

sugar beet
oat
meadow
rapeseed
hop
spelt
triticale
beans
peas
potato
soybeans
asparagus
wheat
w. barley
rye
s. barley
maize

ground truth class

ISPRS Int. J. Geo-Inf. 2018, , 1–19

0%

Figure 7. Confusion matrix of the trained convolutional GRU network on data of the seasons 2016
and 2017. While the confusion of some classes was consistent over both seasons (e.g., winter triticale
to winter wheat ), other classes are classified at different accuracies at consecutive years (e.g., winter
barley to winter spelt ).

5.3. Qualitative Classification Evaluation
For the qualitative evaluation, we used the same network structure as in the previous section. We
inferred the network with 48 px tiles from the evaluation dataset of 2017 for better visibility. At Figure 8
a series of good (A-D) and bad (E, F) classification examples are shown. The first column represents
the input sequence x as band-normalized RGB images from one selected cloud-free observation
x RGB,t . Further columns show the available ground truth labels y, predictions ŷ, and the cross-entropy
loss H (y, ŷ). Additionally, four selected softmax-normalized class activations are displayed in the
last columns. These activations can be interpreted as classification confidences for each class. The
prediction map contains the index of the most activated class at each pixel, which may be interpreted
as the class of highest confidence. The cross-entropy loss is the measure the agreement between the
one-hot representation of the ground truth labels and the activations per class. It is used as the objective
function for network training indicates disagreement between ground truth and prediction even if the
final class prediction is correct. This relation can be observed at fields of several examples, such as
peas in example A, spelt in example B, and oat in example C. However, most classifications at these
examples were accurate, which is expressed by well-defined activation maps.
Often, classifiers use low-pass filters in the spatial dimensions to compensate for high-frequent
noise. These filters typically limit the ability to classify small objects. To evaluate to which degree
the network has learned to apply low-pass filtering we show a tile with a series of narrow fields in
example D. Two thin wheat and maize fields have been classified correctly. However, some errors
occurred on the southern end of an adjacent potato field, as indicated by the loss map. It appears that
the network was able to resolve high-frequent spatial changes and did not apply smoothing the class
activations, as in example F.
Two misclassified fields are shown in example E. The upper wheat field has been confidently
misclassified to summer barley. Underneath, the classification at a second rye field was uncertain
between rye, wheat, and triticale. While triticale, as the least activated class, was not present in the
prediction map, the mixture of rye and wheat is visible in the class predictions.
Example F shows a mostly misclassified tile. Only a few patches of meadow and winter barley
were correctly predicted. The activations of these classes were, compared to previous examples,
generally more blurred and of lower amplitude. Similar to example D, the most activated classes are
also the most frequent in the dataset. In fact, the entire region around the displayed tile seemed

ISPRS Int. J. Geo-Inf. 2018, , 1–19

x RGB,t

labels y

14 of 19

pred. ŷ

loss H (y, ŷ)

activation

activation

activation

activation

maize

meadow

peas

rape

A

1

0

B

spelt

wheat

s. barley

maize

C

meadow

wheat

oat

maize

1

0
1

0

D

meadow

wheat

potato

1

maize

0

E

rye

wheat

triticale

1

s. barley

0

F

wheat

meadow

maize

1

w.barley

0

aspar.

bean

hop

maize

mead.

peas

pot.

rape

soyb.

beet

s.barl

oat

w.barl

rye

spelt

tritic

wheat

Figure 8. Qualitative results of the convolutional GRU sequential encoder. Examples A-D show good
classification results. At example E the network misclassified one maize parcel with high confidence,
which is indicated by incorrect, but well defined activations. At a second field the class activations
reveal a confusion between wheat, meadow and maize. At example F most pixels are misclassified.
However, the class activations show uncertainty in the classification decision.

to be classified poorly. This region was located on the north-west border of the AOI. Further
examination showed that for this region fewer satellite images were available. The lack of temporal
information likely explains the poor classification accuracies. However, this example illustrates that
the class activations give an indication on classification confidence independently from ground truth
information.
6. Discussion
In this section, we compare our approach with other multi-temporal classifications. Unfortunately,
to our best knowledge, no multi-temporal benchmark dataset is available to compare remote sensing
approaches on equal footing. Nevertheless, we provide some perspective of the study domain by
gathering multi-temporal crop classification approaches in Table 3 and categorizing these by their
applied methodology and achieved overall accuracy. However, the heterogeneity of data sources, the
varying extents of their evaluated areas and the number of classes used in these studies impedes a

ISPRS Int. J. Geo-Inf. 2018, , 1–19

15 of 19

Table 3. Overview over recent approaches for crop classification.
Approach

Details
Sensor

Preprocessing

Features

Classifier

Accuracy

# Classes

this work

S2

none

TOA reflect.

ConvRNN

90

17

Rußwurm and Körner [24, 2017]

S2

atm. cor. (SEN 2 COR)

BOA reflect.

RNN

74

18

Siachalou et al. [14, 2015]

LS, RE

geometric
correction,
image registration

TOA reflect.

HMM

90

6

Hao et al. [12, 2015]

MODIS

image reprojection,
atm. cor. [41]

statistical
phen. features

RF

89

6

Conrad et al. [11, 2014]

SPOT,
RE, QB

segmentation,
atm. cor. [41]

vegetation
indices

OBIA+RF

86

9

Foerster et al. [9, 2012]

LS

phen.
normalization,
atm. cor. [41]

NDVI
statistics

DT

73

11

Peña-Barragán et al. [13, 2011]

ASTER

segmentation,
atm. cor. [42]

vegetation
indices

OBIA+DT

79

13

Conrad et al. [10, 2010]

SPOT
ASTER

segmentation,
atm. cor. [41]

vegetation
indices

OBIA+DT

80

6

numerical comparison of achieved accuracies. Despite this, we hope that this table will provide an
overview of the state-of-the-art in multi-temporal crop identification.
Earth observation (EO) data is acquired in periodic intervals at high spatial resolutions. From
an information theoretical perspective, utilizing additional data should lead to better classification
performance. However, the large quantity of data requires methods which are able to process
this information and are robust with regard to observation noise. Optimally, these approaches are
scalable with minimal supervision so that data of multiple years can be included over large regions.
Existing approaches in multi-temporal EO tasks often use multiple separate processing steps, such
as preprocessing, feature extraction, and classification, as summarized by Ünsalan and Boyer [40].
Generally, these steps require manual supervision or the selection of additional parameters based on
region-specific expert knowledge, a process which impedes applicability on large scales. The cost of
data acquisition is an additional barrier, as multiple and potentially expensive satellite images are
required. Commercial satellites, such as RapidEye (RE), Satellite Pour l’Observation de la Terre (SPOT), or
QuickBird (QB), provide images at excellent spatial resolution. However, predominantly inexpensive
sensors, such as Landsat (LS), Sentinel 2 (S2), Moderate-resolution Imaging Spectroradiometer (MODIS),
or Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), can be applied at large
scales, since the decreasing information gain of additional observations must justify image acquisition
costs. Many approaches use spectral indices, such as normalized difference vegetation index (NDVI),
normalized difference water index (NDWI), or enhanced vegetation index (EVI), to extract statistical features
from vegetation related signals and are invariant to atmospheric perturbations. Commonly decision
trees (DTs) or random forests (RFs) are used for classification. The exclusive use of spectral indices
simplifies the task of feature extraction. However, these indices utilize only a small number of available
spectral bands (predominantly blue, red and near-infrared). Thus, methods that utilize all reflectance
measurements, either at top-of-atmosphere (TOA), or atmospherically corrected to bottom-of-atmosphere
(BOA), are favorable, since all potential spectral information can be extracted.
In general, a direct numerical comparison of classification accuracies is difficult, since these are
dependent on the number of evaluated samples, the extent of evaluated area, and the number of
classified categories. Nonetheless, we compare our method with the approaches of Siachalou et al.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

16 of 19

[14] and Hao et al. [12] in detail since their achieved classification accuracies are on a similar level to
ours. Hao et al. [12] used a RFs classifier on phenological features, which were extracted from NDVI
and NDWI time series of MODIS data. Their results demonstrate that good classification accuracies
with hand-crafted feature extraction and classification methods can be achieved if data of sufficient
temporal resolution is available. However, the large spatial resolution (500 m) of the MODIS sensor
limits the applicability of this approach to areas of large homogeneous regions. On a smaller scale,
Siachalou et al. [14] report good levels of accuracy on small fields. For this, they used a hidden Markov
models (HMMs) with a temporal series of four LS images combined with one single RapidEye (RE)
image for field border delineation. Methodically, HMMs and conditional random fields (CRFs) [15] are
closer to our approach since the phenological model is approximated with an internal chain of hidden
states. However, these methods might not be applicable for long temporal series, since Markov-based
approaches assume that only one previous state contains classification-relevant information.
Overall, this comparison shows that our proposed network can achieve state-of-the-art
classification accuracy with a comparatively large number of classes. Furthermore, the S2 data
of non-atmospherically corrected TOA values can be acquired easily and does not require further
preprocessing. Compared to previous work, we were able to process larger tiles by using convolutional
recurrent cells with only a single recurrent encoding layer. Moreover, we neither required atmospheric
correction nor additional cloud classes, since one classification decision is derived from the entire
sequence of observations.
7. Conclusion
In this work, we proposed an automated end-to-end approach for multi-temporal classification,
which achieved state-of-the-art accuracies in crop classification tasks with a large number of crop
classes. Furthermore, the reported accuracies were achieved without radiometric and geometric
preprocessing. The trained and inferred data was atmospherically uncorrected and contained clouds.
In traditional approaches, multi-temporal cloud detection algorithms [43] utilize the sudden positive
change in reflectivity of cloudy pixels and achieve better results than other traditional mono-temporal
remote sensing classifiers, as pointed out by a recent comparison [44]. Results of this work indicate
that cloud masking can be learned jointly together with classification. By visualizing internal gate
activations in our network in Section 5.1, we found evidence that some recurrent cells were sensitive to
cloud coverage. These cells may be used by the network to internally mask cloudy pixels similar to an
external cloud filtering algorithms.
In Sections 5.2 and 5.3, we further evaluated the classification results quantitatively and
qualitatively. Based on several findings, we derived that the network has approximated a
discriminative crop-specific phenological model based on a raw series of TOA S2 observations. Further
inspection revealed that some crops were inconsistently classified in both growing seasons. This may
be caused by seasonally-variable environmental conditions, which may have been implicitly integrated
into the encoded phenological model. We employed our network for the task crop classification
since vegetative classes are well characterized by their inherently temporal phenology. However, the
network architecture is methodically not limited to vegetation modeling and may be employed for
further tasks, which may benefit from the extraction of temporal features. We hope that our results
encourage the research community to utilize the temporal domain for their applications. In this regard,
we plan to publish the T ENSOR F LOW source code of our network along with the evaluations and
experiments from this work.
Supplementary Materials: The code of the network and further material was made publicly available at
https://github.com/TUM-LMF/MultitemporalLandCoverClassification.
Acknowledgments: We would like to thank the Bavarian Ministry of Food, Agriculture and Forestry (StMELF) for
providing ground truth data in excellent semantic and geometric quality. Furthermore, we thank the Leibnitz
Supercomputing Centre (LRZ) for providing access to computational resources, such as DGX-1 and P100 servers
and N VIDIA for providing one T ITAN X GPU.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

17 of 19

ArXiv version history
v1 submitted version to IJGI
v2 added notice of submission date to history
v3 revised manuscript version on review comments. Added section Section 3.2 Prior Work and
improved English language and style.

1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.

12.
13.

14.

15.

16.
17.
18.

19.

20.

Odenweller, J.B.; Johnson, K.I. Crop identification using Landsat temporal-spectral profiles. Remote Sensing
of Environment 1984, 14, 39–54.
Reed, B.C.; Brown, J.F.; VanderZee, D.; Loveland, T.R.; Merchant, J.W.; Ohlen, D.O. Measuring Phenological
Variability from Satellite Imagery. Journal of Vegetation Science 1994, 5, 703–714.
Bahdanau, D.; Cho, K.; Bengio, Y. Neural Machine Translation by Jointly Learning to Align and Translate.
Arxiv 2014, [1409.0473].
Rush, A.; Chopra, S.; Weston, J. A Neural Attention Model for Sentence Summarization. Arxiv 2017.
Shen, S.; Liu, Z.; Sun, M. Neural Headline Generation with Minimum Risk Training. Arxiv 2016,
[arXiv:1604.01904v1].
Nallapati, R.; Zhou, B.; dos Santos, C.N.; Gulcehre, C.; Xiang, B. Abstractive Text Summarization Using
Sequence-to-Sequence RNNs and Beyond. Arxiv 2016, [1602.06023].
Sutskever, I.; Vinyals, O.; Le, Q.V. Sequence to Sequence Learning with Neural Networks. Arxiv 2014,
[1409.3215].
Chorowski, J.; Bahdanau, D.; Serdyuk, D.; Cho, K.; Bengio, Y. Attention-based models for speech
recognition. NIPS, 2015, pp. 577–585, [1506.07503].
Foerster, S.; Kaden, K.; Foerster, M.; Itzerott, S. Crop type mapping using spectral-temporal profiles and
phenological information. Computers and Electronics in Agriculture 2012, 89, 30–40.
Conrad, C.; Fritsch, S.; Zeidler, J.; Rücker, G.; Dech, S. Per-Field Irrigated Crop Classification in Arid
Central Asia Using SPOT and ASTER Data. Remote Sensing 2010, 2, 1035–1056.
Conrad, C.; Dech, S.; Dubovyk, O.; Fritsch, S.; Klein, D.; Löw, F.; Schorcht, G.; Zeidler, J. Derivation
of temporal windows for accurate crop discrimination in heterogeneous croplands of Uzbekistan using
multitemporal RapidEye images. Computers and Electronics in Agriculture 2014, 103, 63–74.
Hao, P.; Zhan, Y.; Wang, L.; Niu, Z.; Shakir, M. Feature Selection of Time Series MODIS Data for Early Crop
Classification Using Random Forest: A Case Study in Kansas, USA. Remote Sensing 2015, 7, 5347–5369.
Peña-Barragán, J.M.; Ngugi, M.K.; Plant, R.E.; Six, J. Object-based crop identification using multiple
vegetation indices, textural features and crop phenology.
Remote Sensing of Environment 2011,
115, 1301–1316.
Siachalou, S.; Mallinis, G.; Tsakiri-Strati, M. A hidden markov models approach for crop classification:
Linking crop phenology to time series of multi-sensor remote sensing data. Remote Sensing 2015,
7, 3633–3650.
Hoberg, T.; Rottensteiner, F.; Feitosa, R.Q.; Heipke, C. Conditional random fields for multitemporal and
multiscale classification of optical satellite imagery. IEEE Transactions on Geoscience and Remote Sensing
2015, 53, 659–673.
Zhu, X.X.; Tuia, D.; Mou, L.; Xia, G.S.; Zhang, L.; Xu, F.; Fraundorfer, F. Deep learning in remote sensing: a
review. Arxiv 2017, [1710.03959].
Hu, F.; Xia, G.S.; Hu, J.; Zhang, L. Transferring Deep Convolutional Neural Networks for the Scene
Classification of High-Resolution Remote Sensing Imagery. Remote Sensing 2015, 7, 14680–14707.
Scott, G.J.; England, M.R.; Starms, W.A.; Marcum, R.A.; Davis, C.H. Training Deep Convolutional Neural
Networks for Land-Cover Classification of High-Resolution Imagery. IEEE Geoscience and Remote Sensing
Letters (GRSL) 2017, 14, 549–553.
Makantasis, K.; Karantzalos, K.; Doulamis, A.; Doulamis, N. Deep Supervised Learning for Hyperspectral
Data Classification through Convolutional Neural Networks. Proceedings GARSS 2015. 2015 IEEE
International Geoscience and Remote Sensing Symposium. IEEE, 2015, pp. 4959–4962.
Castelluccio, M.; Poggi, G.; Sansone, C.; Verdoliva, L. Land Use Classification in Remote Sensing Images
by Convolutional Neural Networks. arXiv preprint arXiv:1508.00092 2015, pp. 1–11, [1508.00092].

ISPRS Int. J. Geo-Inf. 2018, , 1–19

21.

22.
23.
24.

25.
26.

27.
28.
29.

30.

31.
32.

33.
34.
35.
36.
37.
38.
39.
40.
41.
42.

43.

44.

18 of 19

Jia, X.; Khandelwal, A.; Nayak, G.; Gerber, J.; Carlson, K.; West, P.; Kumar, V. Incremental Dual-memory
LSTM in Land Cover Prediction. 23rd ACM SIGKDD International Conference on Knowledge Discovery
and Data Mining, 2017, pp. 867–876.
Braakmann-Folgmann, A.; Roscher, R.; Wenzel, S.; Uebbing, B.; Kusche, J. Sea Level Anomaly Prediction
using Recurrent Neural Networks. Arxiv 2017, [1710.07099].
Sharma, A.; Liu, X.; Yang, X. Land Cover Classification from Multi-temporal, Multi-spectral Remotely
Sensed Imagery using Patch-Based Recurrent Neural Networks. Arxiv 2017, [1708.00813].
Rußwurm, M.; Körner, M. Temporal Vegetation Modelling using Long Short-Term Memory Networks for
Crop Identification from Medium-Resolution Multi-Spectral Satellite Images. IEEE/ISPRS Workshop on
Large Scale Computer Vision for Remote Sensing Imagery (EarthVision), 2017, Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition (CVPR) Workshops.
Hochreiter, S.; Schmidhuber, J. Long Short-Term Memory. Neural Computation 1997, 9, 1735–1780.
Cho, K.; van Merrienboer, B.; Gulcehre, C.; Bahdanau, D.; Bougares, F.; Schwenk, H.; Bengio, Y. Learning
Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation. arXiv preprint
arXiv:1406.1078 2014, [1406.1078].
Graves, A.; Wayne, G.; Danihelka, I. Neural Turing Machines. Arxiv 2014, [1410.5401].
Siegelmann, H.; Sontag, E. On the Computational Power of Neural Nets. Journal of Computer and System
Sciences 1995, 50, 132–150.
Rafal, J.; Wojciech, Z.; Ilya, S. An Empirical Exploration of Recurrent Network Architectures. Proceedings
of the 32nd International Conference on International Conference on Machine Learning, 2015, Vol. 7, pp.
2342–2350, [1512.03385].
Hochreiter, S.; Bengio, Y.; Frasconi, P.; Schmidhuber, J. Gradient flow in recurrent nets: the difficulty
of learning long-term dependencies. A Field Guide to Dynamical Recurrent Networks 2001, pp. 237–243,
[arXiv:1011.1669v3].
Yoshua, B.; Patrice, S.; Paolo, F. Learning long-term dependencies with gradient descent is difficult. IEEE
transactions on neural networks 1994, 5, 157–166.
Shi, X.; Chen, Z.; Wang, H.; Yeung, D.Y.; Wong, W.K.; Woo, W.C. Convolutional LSTM Network: A Machine
Learning Approach for Precipitation Nowcasting. Advances in Neural Information Processing Systems 28
2015.
Maas, A.L.; Hannun, A.Y.; Ng, A.Y. Rectifier Nonlinearities Improve Neural Network Acoustic Models.
Proceedings of the 30 th International Conference on Machine Learning 2013, 28, 6.
Hahnloser R.; Sarpeshkar R.; Mahowald M A.; Douglas R. J..; Seung H. S.. Digital selection and analogue
amplification coexist in a cortex-inspired silicon circuit. Nature 2000, 405, 947–951.
Kingma, D.P.; Ba, J.L. Adam: A Method for Stochastic Optimization. ArXiv Preprint (1412.6980 ) 2014,
[arXiv:1412.6980].
Cohen, J. A coefficient of agreeement for nominal scales. Educational and Psychological Measurement 1960,
20, 37–46.
Karpathy, A. The unreasonable effectiveness of recurrent neural networks, 2015.
Fung, T.; Ledrew, E. The Determination of Optimal Threshold Levels for Change Detection Using Various
Accuracy Indices. Photogrammetric Engineering & Remote Sensing 1988, 54, 1449–1454.
McHugh, M.L. Interrater reliability: the kappa statistic. Biochemia medica 2012, 22, 276–82.
Ünsalan, C.; Boyer, K.L. Review on Land Use Classification. In Multispectral Satellite Image Understanding:
From Land Classification to Building and Road Detection; Springer, 2011; pp. 49–64.
Richter, R. A spatially adaptive fast atmospheric correction algorithm. International Journal of Remote
Sensing 1996, 17, 1201–1214.
Matthew, M.W.; Adler-Golden, S.M.; Berk, A.; Richtsmeier, S.C.; Levine, R.Y.; Bernstein, L.S.; Acharya, P.K.;
Anderson, G.P.; Felde, G.W.; Hoke, M.P. Status of Atmospheric Correction using a MODTRAN4-Based
Algorithm. SPIE proceeding. Algorithms for multispectral, hyperspectral, and ultra-spectral imagery VI,
2000, pp. 199–207.
Hagolle, O.; Huc, M.; Villa Pascual, D.; Dedieu, G. A multi-temporal method for cloud detection, applied
to FORMOSAT-2, VENuS, LANDSAT and SENTINEL-2 images. Remote Sensing of Environment 2010,
114, 1747–1755.
Zupanc, A.S. Improving Cloud Detection with Machine Learning, 2017.

ISPRS Int. J. Geo-Inf. 2018, , 1–19

19 of 19

c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).

