# Context-self contrastive pre-training for crop type semantic segmentation

Michail Tarasiou, Riza Alp Güler, and Stefanos Zafeiriou

###### Abstract

In this paper, we propose a fully supervised pre-training scheme based on
contrastive learning particularly tailored to dense classification tasks. The
proposed Context-Self Contrastive Loss (CSCL) learns an embedding space that
makes semantic boundaries pop-up by use of a similarity metric between every
location in a training sample and its local context. For crop type semantic
segmentation from Satellite Image Time Series (SITS) we find performance at
parcel boundaries to be a critical bottleneck and explain how CSCL tackles the
underlying cause of that problem, improving the state-of-the-art performance
in this task. Additionally, using images from the Sentinel-2 (S2) satellite
missions we compile the largest, to our knowledge, SITS dataset densely
annotated by crop type and parcel identities, which we make publicly available
together with the data generation pipeline. Using that data we find CSCL, even
with minimal pre-training, to improve all respective baselines and present a
process for semantic segmentation at super-resolution for obtaining crop
classes at a more granular level. The code and instructions to download the
data can be found in  `ttps://git`ub.com/michaeltrs/DeepSatModels .

###### Index Terms:

IEEE, IEEEtran, journal, LaTeX, paper, template.

Figure 1: Pixel size is coarse enough for the signal at boundaries to mix.
Here the 10m resolution Sentinel-2 grid is overlaid on a high-resolution image
from Google Earth. The magnified regions show grid locations that contain
signal from multiple crop types and objects. The proposed Context-Self
Contrastive Loss directly compares computed embeddings in local
neighbourhoods. As such, interior regions naturally have a small contribution
to the overall loss while at the same time the network learns to disambiguate
boundary pixels by comparing them to nearby locations. Figure 2: In satellite
images interior regions of agricultural fields consist of large homogeneous
regions with little signal variation. (top-left) spectral band B08 for a
location in France, (bottom-left) ground truth crop type maps for the same
location. (right) histograms of intensities over whole parcel regions (blue)
and interior points (orange). Each histogram presents the intensities for a
specific field as indicated by the respective mask (red outline). Interior,
interior+boundary and exterior locations are shown in white, grey and black
colors respectively. The distribution of intensities for other spectral bands
is similar to that of band B08.

##  I Introduction

The large availability of imagery made available by Earth Observation (EO)
satellites has enabled the development of tools for the automated monitoring
of food supplies and the design and control of policies aiming at agricultural
development. In the European Union (EU) agricultural subsidies provided as
part of the Common Agricultural Policy (CAP) comprised approximately 39$\%$ of
last year’s overall budget at €59 billion. Distribution of these subsidies to
over 10 million farms is currently performed by use of crop maps which are
compiled by parcel topographies and crop types declared by EU farmers. New
reforms of the CAP are expected to place particular emphasis on the use of
automated tools for EO. In particular, the Sentinels for Common Agricultural
Policy (Sen4CAP) project aims at making use of the European Space Agency (ESA)
Sentinel satellites for providing stakeholders with validated algorithms,
products and tools relevant for the management of the CAP. Provided examples
for CAP monitoring of agricultural practices 111http://esa-
sen4cap.org/content/agricultural-practices and crop diversification
requirements 222http://esa-sen4cap.org/content/crop-diversification both
impose thresholds on the area of declared arable land dedicated to different
crop types. Their accuracy should benefit greatly from the development of
automated crop identification processes at high granularity which is the goal
of this work.

However, a performance analysis of state-of-the-art large-scale crop-type
identification systems presented in section V-B showed us that the
classification accuracy varies with the location within a parcel and drops
significantly near the parcel boundaries. This is in part expected for all
dense labelling tasks and particularly so for the segmentation of S2 images in
which location can be a few pixels off with respect to the ground truths,
introducing label noise at object boundaries. We identify an additional
difficulty in correctly classifying parcel boundaries to the fact that
satellites dedicated to monitoring land surface variability trade off
resolution with high revisit time. As a result pixel size is coarse enough for
the signal from different crop types or objects to mix. This is demonstrated
in Fig.1, where we overlay the resolution at hand on high resolution images.
The highest resolution band of S2 images covers a 10m$\times$10m region per
pixel. It is worth pointing out that this pixel-size is large enough to
include complete objects such as roads, which adds further complexity to the
problem. Additionally, the temporal patterns of plant growth during a growth
year can offer significant cues for crop type identification. Models employed
for this task have been shown to benefit from spending their computational
budget on capturing the temporal dimension of the input data [1]. However, as
discussed above, the use of SITS poses a limit on spatial resolution, making
medium resolution images captured at a high temporal frequency an
indispensable ingredient for crop type identification. For crop classification
at the parcel level [2] argue that S2 pixel size is coarser than the typical
agricultural textural information and show that using spatial modelling is not
at all critical for good performance. What would be the best strategy though
if we are interested in predicting crop types for every pixel and moreover for
boundary pixels which are the most severely affected by the coarse size? We
argue that class attribution near boundaries can be improved by learning to
detect semantic boundaries by contrasting each pixel with their neighbourhood
pixels. Furthermore, interior locations of agricultural parcels in satellite
images include large homogeneous regions with little structure and small
signal variation as shown in Fig.2. While the Cross-Entropy (CE) criterion
still needs to assign a class to all interior pixels a contrastive criterion
constrained to local neighbourhoods naturally contains a small loss
contribution from these regions as they are encoded by similar embeddings. The
above two observations constitute the motivation behind our proposed method.
Our key contributions are the following:

  1. 1.

Inspired by recent developments on visual self-attention we design a
contrastive learning pre-train scheme, targeting performance at the parcel
boundaries (section III). Applying the proposed approach on a publicly
available crop segmentation dataset, we raise the state-of-the-art mIoU from
$78.75\%$ to $81.03\%$ (section V-A).

  2. 2.

Regional variations of plant growth patterns is a major challenge hindering
generalization. We present the largest, measured by Area of Interest (AOI),
Time Period of Interest (POI) and number of crop types, dataset for crop type
segmentation which we make publicly available together with the data
generation pipeline to facilitate future research (section IV-A). Our dataset
also includes dense annotations for parcel identities and is, to our
knowledge, the first dataset to include labels at super-resolution. We
benchmark this dataset with state-of-the-art models which we consistently
outperform using our pre-training scheme (section V-A).

  3. 3.

We present a simple method for segmenting crops from satellite images at a
higher resolution than the input by leveraging available ground truths and
show performance improvements in segmenting crops at $\times 4$ the input
resolution over training at base resolution (section V-D).

##  II Background

Notation: we follow the notation of [3, 4] and denote scalars as lower-case
letters $x$, vectors as bold lower case letters $\mathbf{x}$, matrices as bold
upper-case letters $\mathbf{X}$, and higher order tensors as bold calligraphic
letters $\mathbfcal{X}$. The number of dimensions of a tensor is the order of
that tensor. An element $(i,j)$ of an order 3 tensor
$\mathbfcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is indicated as a
vector $\mathbf{x}_{ij}$. An element $(i,j)$ of an order-4 tensor
$\mathbfcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}\times I_{4}}$ is
indicated as a matrix $\mathbf{X}_{ij}$. The tensor-dot product ($\bullet$) of
two tensors is defined as a contraction with respect to the last index of the
first one and the first index of the second one. For example, the tensor-dot
product of an order-3 tensor $\mathbfcal{W}\in\mathbb{R}^{I_{1}\times
I_{2}\times d}$ and a vector $\mathbf{x}\in\mathbb{R}^{d}$ is matrix
$\mathbf{Y}=\mathbfcal{W}\bullet\mathbf{x}\in\mathbb{R}^{I_{1}\times I_{2}}$
with elements $y_{ij}=\sum_{k=1}^{d}\mathbf{w}_{ijk}\cdot x_{k}$.

Crop type identification from satellite images is a challenging task that
involves assigning one of $C$ crop categories at a set of desired locations on
a geospatial grid, most commonly at the agricultural parcel or pixel level.
Rather than using single images as inputs crop identification has been shown
to work best by modelling the temporal patterns of plant growth during a
growing season [1]. Remote Sensing approaches to crop type identification have
included multiple preprocessing steps, adding expert knowledge to correct
satellite image capturing and extract hand-crafted features, e.g. various
vegetation indices [5] used as input to machine learning classifiers [6, 7, 8,
9, 10, 11]. More recently Deep neural Networks (DNN) have achieved
increasingly competitive performance outperforming previous approaches by a
large margin [12, 13, 14, 15, 16, 17]. These methods make use of minimal data
preprocessing and are shown to handle challenges such as cloudy images, image
noise and atmospheric corrections automatically during network feature
extraction [13, 18, 19]. Some works involve temporal modelling of single pixel
or parcel level aggregated features [12, 20, 2] while others jointly capture
temporal and spatial patterns [13, 21, 18, 1]. Fewer works have tackled
semantic segmentation of crop types [13, 18]. In particular [18] propose two
types of models: a fully convolutional 3D-UNET [22] and a 2D-UNET feature
extractor coupled with a CLSTM decoder [23] both achieving state-of-the-art
performance. In our work we use a novel contrastive pre-training scheme to
significantly improve these results. On super-resolution with satellite images
[24] present a method for semantic segmentation at the resolution of input
images utilizing only low-resolution ground-truth labels and assuming a known
joint distribution between low and high-resolution ground truths. Our method
utilizes vector representations of parcel geometries, a common output of
agricultural field studies, making no further assumptions. It is also
different in that we predict classes at a higher resolution than inputs. To
our knowledge, there exist no other pre-training methods for improving the
performance of pixel level semantic segmentation of crop types.

Contrastive learning aims at mapping data points to an embedding space in
which a distance metric encodes the similarity between inputs, with similar
samples brought closer together and dissimilar samples moved far apart. It can
be viewed as a subcategory of metric learning in which the loss function
operates directly on the embedding space and is driven by known (supervised
setting) or assumed (self-supervised setting) similarities between samples
[25]. A classification loss, e.g. CE loss, can also be viewed as a metric
style loss as it constructs an embedding space which preserves class
similarity as distance. However, it does so indirectly by transforming the
embedding vector into logits through matrix multiplication with a learnt
weight matrix and maximizing the likelihood that the sample belongs in the
ground truth class. In comparison, a contrastive loss compares the embeddings
produced by different samples and aims at representing sample similarity
directly as distance in embedding space.

The seminal work of [26] for the face verification task used a Siamese Neural
Network [27] architecture to process two input samples with known semantic
similarity and applied a contrastive loss function directly on extracted
features. Similarly, a Triplet Loss [28, 29] uses three samples at a time, an
anchor, a positive and a negative sample and aims at making the anchor-
positive distances smaller and the anchor-negative distances larger. Several
works have employed contrastive learning to enhance the discriminatory
capacity of neural networks in downstream classification tasks. This can be
achieved by either learning a linear classifier directly on the latent space
[30, 31, 32, 33, 34, 35], or by using network weights as a good initialization
point for further end-to-end training [36, 37]. In terms of supervision [30,
31, 32, 33, 35, 36] use instance discrimination coupled with strong data
augmentation as a self-supervision pretext task while [38, 39, 40, 41] use
available annotations to drive sample similarities in a fully supervised
setting. [35, 36] both use self-supervision to improve performance on dense
classification. Our method is a fully supervised contrastive learning scheme
specifically designed for dense classification tasks. Also in contrast to
other approaches we draw both positive and negative pairs from the same sample
using a single model pass and do not require the use of large batch sizes or
data queues.

Self-attention was introduced in [42] for the 1D task of machine translation
and involves learning useful similarities between features at different
positions which are later used to derive weights for averaging them. For
vision tasks [43] were the first to use self-attention as an additional layer
to Convolutional Neural Network (CNN) architectures. For an input map
$\mathbfcal{X}\in\mathbb{R}^{H\times W\times D_{in}}$, where $H,W$ are the
spatial dimensions and $D_{in}$ is the input number of features, queries, keys
and values
$\mathbf{q}_{ij}=\mathbf{W_{q}}\mathbf{x}_{ij}\in\mathbb{R}^{D_{q}}$,
$\mathbf{k}_{ij}=\mathbf{W_{k}}\mathbf{x}_{ij}\in\mathbb{R}^{D_{q}}$,
$\mathbf{v}_{ij}=\mathbf{W_{v}}\mathbf{x}_{ij}\in\mathbb{R}^{D_{out}}$ are
calculated as linear transformations of the input at every spatial location.
Learnt matrices $\mathbf{W_{q}},\mathbf{W_{k}}\in\mathbb{R}^{D_{q}\times
D_{in}}$, $\mathbf{W_{v}}\in\mathbb{R}^{D_{out}\times D_{in}}$ control the
number of features $D_{q}$ (queries, keys) and $D_{out}$ (values)
respectively. The similarity matrix
$\mathbf{S}_{ij}=\mathbfcal{K}_{\mathcal{N}}\bullet\mathbf{q}_{ij}$ contains
the pairwise affinities between location $(i,j)$ and all other locations in
the layer’s receptive field $\mathcal{N}$. These affinities are subsequently
softmax-normalized such that their values sum to one and are used to aggregate
the projected inputs (values).

$\mathbf{y}_{ij}=\sum_{p\in\mathcal{N}}softmax_{p}(\mathbf{q}_{ij}^{T}\mathbf{k}_{p})\mathbf{v}_{p}$
(1)

The self-attention formulation in eq.1 is permutation invariant with respect
to the input, meaning that any change in the spatial ordering of the inputs
does not affect the output. To exploit positional information, which is
critical to capture the spatial structures in vision data, [44] introduced the
use of two dimensional learnt relative positional encodings to remove the
permutation invariance property while remaining translation equivariant. A
learned relative positional encoding term $\mathbfcal{R}$ is incorporated into
the affinities
$\mathbf{S}_{ij}=(\mathbfcal{K}_{\mathcal{N}}+\mathbfcal{R})\bullet\mathbf{q}_{ij}$,
$\mathbf{r}_{ij}\in\mathbb{R}^{D_{q}}$ introducing a prior of where to look at
in the receptive field $\mathcal{N}$ relative to location $(i,j)$.

$\mathbf{y}_{ij}=\sum_{p\in\mathcal{N}}softmax_{p}(\mathbf{q}_{ij}^{T}(\mathbf{k}_{p}+\mathbf{r}_{p}))\mathbf{v}_{p}$
(2)

We note that in addition to relative positional encodings, absolute positional
encodings have also been employed [42] but their performance has been shown to
be inferior in both language [45] and vision tasks [46]. We also note that
$\mathcal{N}$ is abstractly defined as the receptive field of the attention
mechanism in eqs.1, 2. Early works [42, 43] defined $\mathcal{N}$ as the full
extent of the input allowing each location to attend to all other locations.
This could be very desirable for tasks that require modelling of long range
dependencies, however, it makes the calculation of pair-wise affinities
$\mathbfcal{S}$ prohibitive for high resolution images as its complexity is
quadratic with respect to the total number of input locations
$\mathcal{O}(H^{2}W^{2})$. To address this problem [46] restrict the receptive
field of the self-attention module to a local $(h\times w)$ region and reduce
the computational complexity to $\mathcal{O}(HWhw)$ allowing them to use this
module as stand-alone primitive in place of convolution layers to construct
deep architectures. The formulation of the constrained self-attention module
with relative positional encodings takes the following form:

$\mathbf{y}_{ij}=\sum_{p\in{\mathcal{N}_{h\times
w}(ij)}}softmax_{p}(\mathbf{q}_{ij}^{T}(\mathbf{k}_{p}+\mathbf{r}_{p}))\mathbf{v}_{p}$
(3)

Where $\mathcal{N}_{h\times w}(ij)$ denotes the local $h\times w$ region
around location $(i,j)$. In practice instead of applying the attention
function once with $D_{q}$ and $D_{out}$ number of features, all queries, keys
and values are projected into $n$ different spaces with dimensions
$\frac{D_{q}}{n}$ and $\frac{D_{out}}{n}$ and the self-attention module is
applied $n$-times in parallel. The outputs are further concatenated in the
features dimension and are potentially linearly projected to obtain the final
values.

Our method, inspired by the self-attention module and particularly the work of
[46] calculates a local affinity matrix $\mathbfcal{S}$ using the final layer
features of a CNN but rather than using it for feature aggregation we apply
supervision directly on the calculated affinities. Also, because we are
interested only in the total pair-wise affinities we only make use of a single
attention head.

Figure 3: Proposed pre-training scheme using the Context-Self Contrastive Loss
(CSCL). For a CNN feature map $\mathbfcal{Y}$ we define a similarity metric
$\mathbf{S}_{ij}$ between features extracted at all locations $(i,j)$ and
every location in their $w_{d}$, $w_{r}$-dilated neighborhood (top branch).
Similarly, we derive training labels $\mathbf{L}^{s}_{ij}$ from dense
annotations $\mathbfcal{L}$ for fully supervised training (bottom branch).

##  III Context-Self Contrastive pre-training

Applied on an N-dimensional embedding space of a CNN, the affinity matrix of a
self-attention layer naturally encodes semantic similarity between features
extracted at different locations. Using densely annotated data, we derive
class agnostic labels for providing full supervision on the values of the
similarity matrix and train a network on the task of predicting class
agreement between each location paired with all locations in its local
neighbourhood. In its most general form, the proposed pre-training method
includes the following components:

  1. 1.

an encoder function $\mathbfcal{Y}=f(\mathbfcal{X})$ mapping an input tensor
$\mathbfcal{X}$ to an output feature map $\mathbfcal{Y}$

  2. 2.

a self-similarity function $\mathbfcal{S}=g(\mathbfcal{Y})$

  3. 3.

a ground truth generator module $\mathbfcal{L}^{\mathbf{s}}=h(\mathbfcal{L})$
reformatting dense class annotations $\mathbfcal{L}$ into class-agnostic
$\mathbfcal{L}^{\mathbf{s}}$

  4. 4.

a pre-train loss function $L_{pretrain}(\mathbfcal{S},\mathbfcal{L}^{s})$

In this paper, the encoder function is a deep CNN, $\mathbfcal{Y}$ and
$\mathbfcal{L}$ are order-3 tensors while $\mathbfcal{S}$ and
$\mathbfcal{L}^{\mathbf{s}}$ are order-4 tensors.

###  III-A Similarity function

Drawing inspiration from [46] we use a local affinity matrix of embeddings
$\mathbfcal{Y}$ as our similarity metric. Formally, for a feature map
$\mathbfcal{Y}\in\mathbb{R}^{H\times W\times D_{in}}$ we calculate queries and
keys tensors as linear projections of the final layer embeddings.

$\mathbf{q}_{ij}=\mathbf{W_{q}}\mathbf{y}_{ij}\in\mathbb{R}^{D_{q}}$ (4)
$\mathbf{k}_{ij}=\mathbf{W_{k}}\mathbf{y}_{ij}\in\mathbb{R}^{D_{q}}$ (5)

Where $\mathbf{W_{q}},\mathbf{W_{k}}\in\mathbb{R}^{D_{q}\times D_{in}}$ are
learnt matrices and $D_{in},D_{q}$ define the number of features of the input
and keys, queries tensors. In all experiments we use $D_{in}=128$ and
$D_{q}=128$ although similar results are obtained with different choices for
these parameters.

We calculate the similarity tensor $\mathbfcal{S}$ as the tensor-dot product
between unit-sphere normalized queries
$\hat{\mathbf{q}_{ij}}=\frac{\mathbf{q}_{ij}}{\|\mathbf{q}_{ij}\|}$ and
normalized local keys which are first augmented by relative positional
encodings. The inclusion of relative position information takes place by
directly adding $\mathbfcal{R}\in\mathbb{R}^{w_{d}\times w_{d}\times D_{q}}$
to the local keys. $\mathbfcal{R}$ is composed of learnt parameters
$\mathbf{R}^{H}\in\mathbb{R}^{w_{d}\times\frac{D_{q}}{2}}$ and
$\mathbf{R}^{W}\in\mathbb{R}^{w_{d}\times\frac{D_{q}}{2}}$ representing
relative height and width information respectively. We calculate
$\mathbf{r}_{ab}=[\mathbf{r}^{H}_{a},\mathbf{r}^{W}_{b}]$ for each location
$(a,b)$ within the local window. Here $[\cdot,\cdot]$ denotes vector
concatenation.

$\mathbfcal{K^{R}}_{\mathcal{N}_{w_{d},w_{r}}(ij)}=\mathbfcal{K}_{\mathcal{N}_{w_{d},w_{r}}(ij)}+\mathbfcal{R}$
(6)

$\mathbf{S}_{ij}=\hat{\mathbfcal{K^{R}}}_{\mathcal{N}_{w_{d},w_{r}}(ij)}\bullet\hat{\mathbf{q}_{ij}}\in\mathbb{R}^{w_{d}\times
w_{d}}$ (7)

$\mathcal{N}_{w_{d},w_{r}}(ij)$ denotes the $(w_{d}\times w_{d})$,
$w_{r}$-dilated neighbourhood centered around location $(i,j)$. For keys
normalization we normalize each spatial element separately by each own norm.
Normalizing both keys and queries is a key step for the method to work as not
doing so would result in an ever increasing vector norm when training with the
objective function presented in eq.11. The process for calculating local
affinities is presented schematically at the top branch of Fig.3.

Figure 4: CSCL ground truth generator example on a 2D crop type label map. To
generate labels we compare the class at the center location with all locations
in the defined neighbourhood. We use windows with parameters (a, b)
$w_{d}=3,w_{r}=1$, (c) $w_{d}=3,w_{r}=3$, (d) $w_{d}=5,w_{r}=1$, (e)
$w_{d}=3,w_{r}=2$, (f) $w_{d}=5,w_{r}=2$.

###  III-B Loss function

Before defining the pre-training loss function we need to reformat the ground
truth data such that they encode the similarity between each location in a 2D
lattice and all locations in its local neighbourhood. We start from dense one-
hot ground truths $\mathbfcal{L}\in\\{0,1\\}^{H\times W\times
C},\sum_{k}\mathbfcal{L}_{ijk}=1\>\forall\>i,j$ of spatial extent $H\times W$
and $C$ the number of classes. In essence we need pre-training ground truths
$\mathbfcal{L}^{s}\in\\{0,1\\}^{H\times W\times w_{d}\times w_{d}}$ to be a
binary tensor of the same dimensions as $\mathbfcal{S}$ with values

$l^{s}_{ijkm}=\begin{cases}1&,\mathbf{l}_{w_{s}i,w_{s}j}=\mathbf{l}_{w_{s}i-w_{r}(c-k),w_{s}j-w_{r}(c-m)}\\\
0&,otherwise\\\ \end{cases}$ (8)

where $c=\lfloor\frac{w_{d}}{2}\rfloor$ is the index for the center location
of the local window and $(k,m)\in\\{0,...,w_{d}-1\\}^{2}$ define a location
within the spatial extent of the local window. Formally, we obtain suitable
labels that encode class agreement between each location $(i,j)$ and its
$(w_{d}\times w_{d})$, $w_{r}$-dilated neighbourhood by application of the
tensor-dot product.

$\mathbf{L}^{s}_{ij}=\mathbfcal{L}_{\mathcal{N}_{w_{d},w_{r}}(ij)}\bullet\mathbf{l}_{ij}\in\\{0,1\\}^{w_{d}\times
w_{d}}$ (9)

$\mathcal{N}_{w_{d},w_{r}}(ij)$ denotes the $(w_{d}\times w_{d})$,
$w_{r}$-dilated neighbourhood centered around location $(i,j)$. The process is
shown at the bottom branch of Fig.3 while an example of generated labels for
various size and dilation parameters is shown in Fig.4.

For $L_{pretrain}$ we use a contrastive loss function extending [38] to a
dense embedding space and using cosine similarities between elements of
$\mathbfcal{Q}$ and $\mathbfcal{K}$ instead of Euclidean distances. In order
to mask supervision at the zero-padded regions, at the centers of extracted
windows and any other location we might want to exclude from learning, e.g.
background class, we use a mask tensor $\mathbfcal{M}\in\\{0,1\\}^{H\times
W\times w_{d}\times w_{d}}$ of the same dimensions as $\mathbfcal{S}$,
$\mathbfcal{L}^{\mathbf{s}}$.

$\begin{split}L_{pretrain}=-\frac{1}{\sum\mathbfcal{M}}\sum(\mathbfcal{M}\odot(\mathbfcal{L}^{\mathbf{s}}\odot\mathbfcal{S}+(1-\mathbfcal{L}^{\mathbf{s}})\\\
\odot\min(0,m-\mathbfcal{S})))\end{split}$ (10)

Here, summation is over all dimensions of
$\mathbfcal{M},\mathbfcal{L}^{\mathbf{s}},\mathbfcal{S}$ with indices not
shown to avoid cluttering. $m$ is a scalar hyperparameter used to exclude easy
negatives. With $m>0$ the cosine similarities for easily recognised dissimilar
pairs that are smaller than $m$ and have no contribution to the total loss.
Several studies have highlighted the importance of hard negative pairs in
contrastive learning [38, 47, 48]. In our method all negative pairs originate
from nearby locations in the same training sample featuring high overlaps in
their receptive fields. As such there is no possibility for easy negatives and
we can discard the $m$ parameter originally employed to exclude the
contribution of these pairs. For all our experiments we set $m=0$ and weight
the contribution of positive pairs by $\lambda$ leading to the following
equivalent expression.

$L_{pretrain}=-\frac{1}{\sum\mathbfcal{M}}\sum(\mathbfcal{M}\odot((\lambda+1)\mathbfcal{L}^{\mathbf{s}}-1)\odot\mathbfcal{S})$
(11)

###  III-C Interpreting the sliding window parameters

Figure 5: Mean prediction accuracy for positive (left) and negative (right)
pairs w.r.t location in the sliding window ($w_{d}=5$). The prediction
accuracy drops for positive and increases for negative pairs with increasing
distance from the center.

We note two ways in which the values of the sliding window parameters $w_{d}$,
$w_{r}$ affect training. First, they introduce a trade-off between the
difficulty of fitting positive and negative pairs. Small $w_{d}$, $w_{r}$
values result in comparing locations that are closer together making it easier
to recognise positive pairs and harder for negative pairs, the other way being
the case for large window size and dilation values. This can be seen in Fig.5
where we plot the mean prediction accuracy w.r.t. the location on the sliding
window. We assume a positive prediction if the calculated affinity is greater
than a specified threshold. We note how the accuracy is shown to be lower for
positive pairs the farther away from the center of the sliding window and
larger for negative pairs.

Secondly, they affect the statistics of the generated labels. Small $w_{d}$,
$w_{r}$ values results in comparing locations that are more likely to belong
to the same object resulting in a higher proportion of positive pairs. Class
imbalance in class-agnostic ground truths - in favour of positive pairs - has
been identified as a key issue. This effect is more pronounced for high
resolution outputs with fewer distinct objects, in which case we have found it
beneficial to increase $w_{d}$, $w_{r}$ and reduce $\lambda$. We note the
option of applying eqs.(7, 9) at a stride $w_{s}$ to reduce the amount of
computation in high resolution outputs. This reduces the attention complexity
by a factor of $w_{s}^{2}$ to $\mathcal{O}(\frac{HWhw}{w_{s}^{2}})$ compared
to the non-strided case. Finally, for some special case of parameters
$w_{r}>1$, $w_{s}>1$ some locations $(i,j)$ are never used for calculating
$\mathbfcal{S}$, facilitating a more efficient implementation. This is further
elaborated in Appendix A.

country |  $A_{AOI}$ (ha) | $N_{parc}$ | $N_{years}$ | $N_{bands}$ | $N_{T}$ | $CC_{max}$ | $N_{cl}$ | dense annt. | ids annt. | SR labels
---|---|---|---|---|---|---|---|---|---|---
Romania [49] | 64k | 2.5k | 1 | 4 | 13 | $10\%$ | 4 |  |  |
France [2] | 1.21M | 200k | 1 | 10 | 24 | - | 20 |  |  |
France [19] $\dagger$ | 2.7M | 768k | 1 | 13 | 51-102 | $80\%$ | 9 |  |  |
Austria [50] $\ddagger$ | - | - | 1 | 9 | 26-37 | $80\%$ | 16 |  |  |
Germany [13] | 428k | 137k | 2 | 14 | 2-51 | $80\%$ | 17 | ✓ |  |
France [51] | 400k | 124k | 1 | 10 | 33-61 | - | 19 | ✓ | ✓ |
France (T31TFM-1618) | 1.21M | 200k | 3 | 14 | 14-33 | $70\%$ | 166 | ✓ | ✓ | ✓
Table I: Large scale datasets of SITS for crop type recognition. $A_{AOI}$:
AOI size in hectares, $N_{parc}$: number of parcels in AOI, $N_{years}$:
number of growing seasons in dataset, $N_{bands}$: number of satellite image
bands, $N_{T}$: number of observations for each growing season, $CC_{max}$:
maximum cloud coverage, $N_{cl}$: number of crop types, dense annt.: dense
annotations included, ids annt.: parcel ids included, SR labels: labels at
super-resolution included. Input data are: timeseries of mean values over the
extent of an agricultural field $\dagger$, sampled per pixel timeseries
$\ddagger$. Figure 6: Pixel counts for all 166 crop types available in the
T31TFM-1618 dataset. Red column corresponds to background class. Black columns
indicate the 20 classes used in experiments. Figure 7: Day of year (doy) of
satellite observations included in the T31TFM-1618 dataset. Figure 8: Example
of data included in the T31TFM-1618 dataset. (left) Satellite image inputs. We
retain four 10m bands, six 20m bands and three 60m bands for variable number
of time steps T depending on location. (right) Ground truth data. From left to
right we show crop types, parcel ids and doubly assigned pixels. The same
ground truth data are plotted at 10m resolution (top row) and 2.5m $\times 4$
super-resolution (bottom row).

##  IV Experiments

###  IV-A T31TFM-1618 dataset for crop type segmentation

For generating our densely annotated crop segmentation dataset we are using
ground truth crop type data available in the records of the French Land Parcel
Identification System (RPG) 333https://www.data.gouv.fr/en/datasets/registre-
parcellaire-graphique-rpg-contours-des-parcelles-et-ilots-culturaux-et-leur-
groupe-de-cultures-majoritaire. In total, the RGP includes approximately $10M$
parcels in the form of parcel geometries. As a source of satellite images we
use S2 products available through the Copernicus Open Access Hub (COAH)
444https://scihub.copernicus.eu/. The new dataset which we will refer to as
T31TFM-1618 covers the full extent of the T31TFM S2 tile in France for the
years 2016-18. This is the same AOI as [2], however, we extend the time period
of interest to three years and provide dense annotations for crop types and
parcel identities. Our AOI also has some overlap with the AOI of [51] which is
concurrent work to ours. In Table I we compare our dataset with other large
scale crop type identification datasets. The T31TFM-1618 dataset is the
largest, to our knowledge, publicly available dataset for crop type semantic
segmentation. Its AOI covers the full extent of a densely cultivated S2 tile
and by including three years it spans the longest time period among similar
datasets. Additionally, it includes dense annotations for parcel identities
enabling the training of instance segmentation models and is the first dataset
to include semantic labels at a higher resolution than the input images. In
Fig.6 we show pixel counts per label type for all labels. Days in the form of
”day-of-the-year” (doy) for the observations included in the dataset are shown
in Fig.7 separate for each year.

The automated process for generating the dataset includes the following steps:
defining an AOI and a POI, downloading all relevant satellite images,
extracting and rasterizing all relevant ground truth data, grouping and
sorting satellite images to create timeseries objects and matching ground
truths with timeseries by location. As a final step, we create boolean masks
for regions that should be masked out during training and evaluation. These
include non-agricultural regions, crop types we do not wish to include as well
as pixels contained to more than one parcels due to geocoding errors. Samples
that contain no classes other than the background class are removed from the
training and evaluation sets.

We split inputs into 480m square regions corresponding to $48\times 48$ pixels
for the highest resolution satellite band. During the ground truth
rasterization step we need to make a choice for the size of the grid. We
proceed with the rasterization of ground truth polygons using two different
grid sizes. A grid size of 10m matches the maximum resolution of S2 images. In
this manner, ground truth samples consist of $48\times 48$ pixels and each
pixel from the highest resolution S2 bands corresponds to one pixel from the
raster ground truth data. Additionally, we use a 2.5m grid size to extract the
$\times 4$-super-resolution ground truths which leads to a $4\times 4$ label
map per pixel for the highest resolution S2 bands. Super-resolution patches
consist of $192\times 192$ pixels. Example ground truth data rasterized using
10m and 2.5m grids are presented in Fig.8.

To facilitate further research, we make the T31TFM-1618 dataset publicly
available together with evaluation codes.

###  IV-B Network architectures

Figure 9: UNET3D ($a\rightarrow b$), UNET3Df ($a\rightarrow c\rightarrow d$)
and UNET3Df-x4 ($a\rightarrow c\rightarrow e$) models. All models share the
same backbone. Dashed lines indicate how different combinations of components
lead to different architectures. Box size represents spatial-temporal
dimensions, feature dimensions are shown by color coding. Indicators $\times
2$, $\times 4$ show super-resolution feature maps which are not drawn in
scale.

All models explored here learn a mapping from image timeseries to a dense
label map , i.e. $f:\mathbb{R}^{T\times H\times W\times
D}\rightarrow\mathbb{R}^{H\times W\times C}$ where $D$ equals the number of
bands in a multispectral satellite image plus a spatially-constant doy feature
and $C$ is the number of crop types.

In the experiments that follow we use the backbone architectures from [13,
18]. We have modified the UNET3D architecture [22] to make it more suitable
for our contrastive pre-training scheme. Towards this goal we only modified
network components after the last residual connection which leads to a feature
map of size $(T\times H\times W\times 128)$. In the subsequent layers of
UNET3D the number of channels is further reduced to $32$ and $C$ while all
spatial-temporal dimensions are not affected. Finally, the classifier layer
collapses the temporal dimension resulting to an output logits map of size
$(H\times W\times C)$. The issue we encounter here is that the final 2D
embedding space contains as few as $C$ channels which can be small to encode
meaningful similarities between locations in eq.11. We thus aim to construct a
richer final 2D embedding space. In UNET3Df we maintain the number of features
from the last unaffected feature map to $128$ for the following two layers and
then collapse the temporal dimension. This way the final feature map
$y\in\mathbb{R}^{H\times W\times 128}$ contains $128$ features ($128>C$) which
is typical of the number of features used in deep contrastive learning [52,
33, 34]. To obtain dense logits we transform $\mathbfcal{Y}$ into a feature
map with $C$ number of features by use of a linear classifier layer. A
schematic comparison between the UNET3D and the modified UNET3Df architectures
is presented in Fig.9.

Our super-resolution network builds on the UNET3Df architecture without the
classifier layer. The final 2D embedding space $\mathbfcal{Y}$ is further
processed by two 2D upsampling blocks each consisting of a deconvolution layer
[53] with kernel size $3\times 3$ and stride $2$, a Batch Normalization (BN)
layer [54] and Leaky ReLu activation [55]. We denote
$\mathbfcal{Y}^{\prime}\in\mathbb{R}^{4H\times 4W\times 128}$ the feature map
after the two upsampling layers. For CSCL pre-training we now use
$\mathbfcal{Y}^{\prime}$ in the same manner that we used $\mathbfcal{Y}$ for
the base resolution case. For semantic segmentation training we apply a linear
classifier in every spatial location leading to a $C$-dimensional map of
logits. A schematic representation of the super-resolution component and its
relationship with the UNET3Df architecture is presented in Fig.9.

###  IV-C Implementation details

Evaluation protocol.  Our main goal is to test the benefit of the proposed
pre-training method as an alternative to random initialization for crop type
semantic segmentation. All experiments dedicated to this purpose follow the
same two step protocol:

  1. 1.

we obtain the baseline performance for a model of choice and dataset by
training the network to convergence starting from random initialization of
model parameters

  2. 2.

we pre-train the same model using the methodology presented in section III
starting from random initialization of model parameters. The pre-trained model
is subsequently used as the initialization point for supervised end-to-end
training following the same training schedule and hyperparameters as step 1.

Presented evaluation metrics are pixel-level overall accuracy (averaged over
pixels) and mIoU, macro F1 scores (averaged over classes) in the evaluation
sets for Germany and France. We run all experiments five times at base
resolution and three times at super-resolution and report the mean value and
the $95\%$ confidence interval around the mean for each metric.

Datasets. We split the T31TFM-1618 dataset into training and evaluation sets
consisting of $120,438$ and $18,749$ samples respectively with no overlapping
locations. We do not apply any label grouping but rather select classes
directly from the $349$ RPG codes such that there are at least $20,000$ such
parcels in the train set and $2,000$ in the evaluation set for each year
resulting in $C=20$ classes in total. Pixel counts for these classes are
presented in the black columns of Fig.6. From publicly available crop
segmentation datasets we make use of the dataset presented in [13] which
covers an AOI of $102km\times 42km$ north of Munich, Germany containing $C=17$
classes. We train all models using the provided train split and report results
for the evaluation split for a fair comparison with literature.

| Germany | France
---|---|---
model | Acc. | mIoU | F1 | Acc. | mIoU | F1
BiCGRU [13] | 0.897 | - | 0.831 | $88.4\pm 0.1$ | $57.3\pm 0.1$ | $68.7\pm 0.1$
UNET3D [18] | $92.5\pm 0.1$ | $76.1\pm 0.2$ | $85.0\pm 0.2$ | $88.2\pm 0.1$ | $57.1\pm 0.2$ | $68.6\pm 0.2$
UNET2D-CLSTM [18] | $92.9\pm 0.2$ | $78.8\pm 0.2$ | $87.1\pm 0.2$ | $89.0\pm 0.2$ | $58.3\pm 0.2$ | $68.9\pm 0.2$
UNET2D-CLSTM-CSCL | $\mathbf{94.2\pm 0.1}$ | $\mathbf{81.0\pm 0.1}$ | $\mathbf{88.7\pm 0.2}$ | $\mathbf{89.9\pm 0.1}$ | $\mathbf{60.0\pm 0.1}$ | $\mathbf{70.9\pm 0.2}$
UNET3Df | $92.5\pm 0.2$ | $76.8\pm 0.2$ | $85.5\pm 0.1$ | $88.5\pm 0.1$ | $57.2\pm 0.2$ | $68.8\pm 0.3$
UNET3Df-CSCL | $\mathbf{93.7\pm 0.2}$ | $\mathbf{79.4\pm 0.2}$ | $\mathbf{87.6\pm 0.3}$ | $\mathbf{89.3\pm 0.1}$ | $\mathbf{59.8\pm 0.1}$ | $\mathbf{70.4\pm 0.2}$
Table II: Model comparison for Germany and France. Presented metrics are mean
and $95\%$ confidence intervals among five runs. Models from [18] were
reimplemented using provided code as reported results in the original study
used a custom evaluation set. A significant improvement over the baselines is
achieved through CSCL pre-training using no additional data.  Figure 10:
Confusion matrices for models trained in Germany. From left to right: (top)
UNET3Df, UNET3Df-CSCL, (bottom) UNET2D-CLSTM, UNET2D-CLSTM-CSCL. For diagonal
and all other terms $\geq 0.05$ we explicitly show numerical values. Figure
11: Confusion matrices for models trained in France (T31TFM-1618). From left
to right: (top) UNET3Df, UNET3Df-CSCL, (bottom) UNET2D-CLSTM, UNET2D-CLSTM-
CSCL. For diagonal and all other terms $\geq 0.05$ we explicitly show
numerical values. Figure 12: Qualitative comparison of models in Germany
(left) and France (right). Top to bottom rows: UNET3Df, UNET3Df-CSCL, ground
truth crop types, UNET2D-CLSTM, UNET2D-CLSTM-CSCL. White ”x” indicates false
prediction at a particular location.

Pre-training and training details.  For all experiments, including pre-
training and fine-tuning, at base and super-resolution, we use the Adam
optimizer [56] with initial learning rate $0.0001$ and decay rate $0.975$
applied every second epoch, beta parameters $(0.9,0.999)$ and batch size $32$.
All models are trained for $50$ epochs for France and $300$ epochs for
Germany, which is enough for the optimization process to converge. For Germany
we reset the learn rate back to its initial value at epoch $150$, and continue
with the same decay rate as before, which leads to improved performance. This
strategy makes no difference, and is not used, when training with T31TFM-1618.

In terms of data augmentation we take advantage of the bilateral symmetry in
EO imagery and perform random horizontal and vertical flipping with
probability 0.5 in each dimension. We do not perform random cropping but
rather use the original patch dimensions per dataset ($24\times 24$ for
Germany and $48\times 48$ for France).

For pre-training we use a sliding window with $w_{d}$=3, $w_{r}$=1, $w_{s}$=1
and $\lambda$=0.125. Following the completion of pre-training we discard all
components apart from the encoder which is used as an initialization point for
fine-tuning. We add a randomly initialized linear classifier mapping encoded
features to logits $\mathbfcal{O}\in\mathbb{R}^{H\times W\times C}$, where $C$
is the number of classes, and proceed with end-to-end training of the network
using the Masked CE loss function.

##  V Results

###  V-A Comparison with state-of-the-art

We compare our method’s performance over random initialization in training
recently proposed models using the T31TFM-1618 and German datasets. Table II
and Figs.10, 11 respectively present performance metrics and confusion
matrices while Fig.12 visually compares predictions. Both quantitative and
qualitative comparisons indicate a clear advantage of CSCL over the randomly
initialized baselines, achieving new state-of-the-art results. In Fig.14 we
present train and evaluation loss curves for Germany. Both training and
evaluation losses are lower for the pre-trained models compared to their
random initialization counterparts. We note that performance is significantly
lower for our experiments using the T31TFM-1618 compared to the German data.
This is mainly attributed to a few bad performing classes as shown in Fig.11.
As a sanity check we run further tests with the label groupings from [2] who
reported overall accuracy $0.942$ and mIoU $0.509$ reaching similar
performance ($0.963$ and $0.543$) and although not directly comparable it
indicates the added difficulty of training in this AOI.

Figure 13: Model performance at parcel interior vs boundary locations. We
demonstrate that the boost in segmentation performance mostly occurs near the
parcel boundaries. Figure 14: Training and evaluation loss curves for Germany.

###  V-B Performance at parcel boundaries

For our analysis we define a semantic boundary as any location for which not
all ground truths in a $3\times 3$ local neighbourhood share the same value.
By plotting model performance in terms of interior and boundary locations in
Fig.13 we notice a clear performance drop at boundaries also evidenced by
qualitative analysis of results in Fig.12. The same conclusion can be reached
by inspecting the confusion matrices for interior and boundary points in
Figs.10, 11. We also notice smaller but consistent improvements at interior
locations for all models, however, using CSCL there is a significant
improvement at boundaries over the baselines showcasing the desired effect of
the method.

model | pre-train | $\mathbfcal{R}=0$ | $\mathbf{W_{k}}=\mathbf{W_{q}}$ | $\mathbf{W_{k}}=\mathbf{W_{q}}=\mathbf{I}$ | Acc. | mIoU | F1
---|---|---|---|---|---|---|---
UNET3Df |  |  |  |  | $88.5\pm 0.1$ | $57.2\pm 0.2$ | $68.8\pm 0.3$
UNET3Df-scaled |  |  |  |  | 88.40 | 57.70 | 69.20
UNET3Df-equal |  |  |  |  | 88.60 | 57.93 | 69.60
UNET3Df-CSCL | ✓ | ✓ |  |  | 89.25 | 59.48 | 70.34
UNET3Df-CSCL | ✓ | ✓ | ✓ |  | 89.30 | 59.22 | 70.23
UNET3Df-CSCL | ✓ | ✓ |  | ✓ | 89.17 | 59.06 | 69.95
UNET3Df-CSCL | ✓ |  |  |  | $\mathbf{89.3\pm 0.1}$ | $\mathbf{59.8\pm 0.1}$ | $\mathbf{70.4\pm 0.2}$
UNET3Df-CSCL | ✓ |  | ✓ |  | 89.26 | 59.40 | 70.27
UNET3Df-CSCL | ✓ |  |  | ✓ | 89.27 | 59.13 | 70.04
Table III: Results of ablation study in France. (top) Alternative
methodologies to pre-training for reducing statistical bias towards interior
locations. ”UNET3Df-scaled” is initialized randomly and trained such that the
loss contribution from interior pixels is scaled down by the ratio (# boundary
pixels)/(# interior pixels). ”UNET3Df-equal” is initialized randomly and
trained on an equal number of interior and boundary pixels by randomly masking
out a variable number of interior pixel locations for every training sample.
(center) Ablation on CSCL components without positional encodings, and
(bottom) ablation on CSCL components with positional encodings.

###  V-C Ablation study

We present an ablation study on CSCL using the T31TFM-1618 dataset. Here we
only use the UNET3Df model mainly because it is faster to train compared to
UNET2D-CLSTM.

First, given the large number of interior compared to boundary pixels, we
present a comparison with simpler methods aiming to reduce the statistical
bias of the model towards interior pixels. Towards this goal, we discard pre-
training and train directly by scaling - separately for each sample - the loss
component from interior pixels by the ratio of $\gamma=\frac{\textrm{num.
boundary pixels}}{\textrm{num. interior pixels}}$ (UNET3Df-scaled). In a
similar spirit we also train on a balanced number of interior and boundary
pixels (UNET3Df-equal). For every training sample we simply mask some interior
(or boundary) locations at random such that an equal number of respective
locations contribute in the loss function. Results, presented in Table III
suggest that scaling down the loss component of interior pixels does improve
performance compared to the baseline case ($57.9\%$ vs $57.3\%$ mIoU) but this
effect is much less significant compared to the improvements obtained by CSCL
($+0.6\%$ vs $+2.5\%$ mIoU).

Next, we perform and ablation on the various CSCL components presented in
section III. In particular, we examine the downstream segmentation performance
with or without the inclusion of positional encodings in eq. 7, assuming that
learnt parameters $\mathbf{W_{q}}=\mathbf{W_{k}}$, i.e. same keys and queries,
and finally assuming $\mathbf{W_{q}}=\mathbf{W_{k}}=\mathbf{I}$, i.e. no
feature projection prior to calculating the pairwise affinities. Results are
presented in Table III. Overall, we observe significant improvements for all
variants of the pre-training method compared to random initialization but
lower performance compared to the general form of CSCL which includes
positional encodings and different projection matrices.

$\lambda$ | 0.01 | 0.05 | 0.125 | 0.25 | 0.50 | 1.0
---|---|---|---|---|---|---
mIoU | 0.551 | 0.556 | 0.567 | 0.545 | 0.541 | 0.534
Table IV: Model mIoU from varying the weight $\lambda$ of loss contribution from positive pairs. Figure 15: Model mIoU w.r.t. pre-training performance. For each presented metric we fit degree 1 (grey) and degree 2 (grey dashed) polynomials and report $R^{2}$ values. Black dashed line $mIoU=0.532$ shows the no pre-training baseline. $w_{r}$ $w_{d}$ | 3 | 5 | 7
---|---|---|---
1 |  0.567 (3) | 0.553 (5) | 0.550 (7)
2 | 0.552 (5) | 0.541 (9) | 0.536 (13)
Table V: Model mIoU for different CSCL window parameters. Corresponding
$w_{d}^{eff}$ in parenthesis.

Finally, we perform three ablations on CSCL hyperparameters using our dataset
in France for 2018. We find that the method is robust to varying parameters as
we improve upon the no pre-training baseline (mIoU=$0.532$) even when
initializing with a model pre-trained for only 1 epoch (mIoU=$0.541$) which is
the least amount of pre-training tested. Overall, our ablation study
highlights the significance of negative pairs in achieving good performance.

In Table IV we compare results varying the weight of loss from positive pairs
$\lambda$. We find smaller weights generally leading to a better performance
than larger ones. However, best performance is achieved by $\lambda=0.125$
which is close to the ratio of negative to positive pairs in our training set
($0.111$).

The relationship between pre-training performance and final segmentation
performance is examined in Fig.15. This is done by sampling initializations
from various pre-taining steps. For each presented metric we fit degree 1 and
2 polynomials, plot the regression line and report the coefficient of
determination $R^{2}$ indicating the proportion of variance in mIoU that can
be predicted by each pre-train metric following these simple regression
models. As expected, segmentation performance is generally found to improve
with more pre-train steps. We further demonstrate that lower loss and more
heavily contrasted negative pairs during pre-training lead to a better final
performance and are better predictors of mIoU than pre-training steps. There
appears to be no clear relationship between the mean cosine similarity for
positive pairs and final performance.

Results for varying CSCL window parameters $w_{d}$, $w_{r}$ are presented in
Table V. We find that windows covering smaller effective receptive fields
($w_{d}^{eff}=(w_{d}-1)*w_{r}+1$), thus containing harder negative pairs on
average, perform better overall.

###  V-D Semantic segmentation at super-resolution

In Table VI we present results for the baseline and CSCL pre-trained models
using the $\times 4$ ground truths. There, the ”train res.” column indicates
which data were used both for pre-training and semantic segmentation training,
e.g. the same dataset ($\times 1$ or $\times 4$) is always used for both
tasks. To show the benefit over training with $\times 4$ ground truths and no
pre-training we upscale the output of the best pre-trained model from Table II
($\times 1$) via nearest neighbour interpolation. We note that $\times 4$ with
pre-training outperforms the upscaled model (mIoU 0.597 vs 0.586) and that
pre-training has a more significant effect that training at $\times 4$
resolution alone with no pre-training (mIoU 0.586 vs 0.581). Fig.16 shows a
qualitative comparison of predictions using the models from Table VI. We note
the clustering of false predictions in blocks for the non pre-trained upscaled
output and the performance improvement at boundaries for the CSCL vs the
randomly initialized model.

model | train res. | Acc. | mIoU | F1
---|---|---|---|---
UNET3Df | $\times 4$ | $66.2\pm 0.5$ | $58.1\pm 0.2$ | $69.1\pm 0.2$
UNET3Df-CSCL | $\times 1$ | $66.9\pm 0.2$ | $58.6\pm 0.2$ | $69.8\pm 0.2$
UNET3Df-CSCL | $\times 4$ | $67.7\pm 0.4$ | $59.7\pm 0.3$ | $70.9\pm 0.2$
Table VI: Super-resolution model performance. The output of models trained at
base ($\times 1$) resolution is upscaled to match $\times 4$ ground truths via
nearest neighbour interpolation. Presented metrics are mean and $95\%$
confidence intervals among 3 runs.  Figure 16: Qualitative comparison of
super-resolution models from Table VI. From left to right: UNet3Df-$\times 4$,
UNet3Df-CSCL-$\times 1$, UNet3Df-CSCL-$\times 4$. White ”x” indicates false
prediction at a particular location.

##  VI Conclusion

In this study, we tackled the problem of crop type semantic segmentation from
satellite images. Motivated by an identified performance drop at object
boundaries we proposed CSCL, a pre-training scheme targeting the ability of a
CNN to model the relationship between each input location and its local
neighborhood. Using our method we improved the state-of-the-art performance in
crop type semantic segmentation. Finally, we create a very large dataset of
SITS for crop type semantic segmentation which also enables training of
semantic segmentation models at super-resolution. We make this data publicly
available to facilitate further research into this domain.

##  Appendix A CSCL implementation for reduced space complexity

⬇

1# f: encoder network

2# X: input to encoder network

3# w_d: CSCL window size

4# w_r: CSCL windowdilation

5# w_s: CSCL window stride

6# W_q: queries weights (Dy x Dq)

7# W_k: keys weights (Dy x Dq)

8# R: relative positional encodings (w_d x w_d x Dq)

9# \------------------------------------------------

10# tensordot: tensor contraction with respect to the

11# last index of the first tensor and the first index

12# of the second tensor

13# batch_tensordot: tensordot product applied

14# separately to each element in a batch

15

16Y = f(X)

17# (H x W x Dy)

18

19if w_s >= w_r:

20 assert stride % dilation == 0

21 w_s0 = w_r

22 w_s1 = w_s // w_r

23 w_d1 = w_d

24 w_r = 1

25elif w_s < w_r:

26 assert w_r % w_s == 0

27 w_s0 = w_s

28 w_s1 = 1

29 w_d1 = (w_d - 1) * w_r + 1

30 w_r = int(w_r / w_s)

31

32Q = tensordot(Y[::w_s, ::w_s], W_q)

33# (H/w_s x W/w_s x Dq)

34Q = Q.reshape(-1, Dq)

35# (H*W/w_s**2 x Dq)

36

37K = tensordot(Y[::w_s0, ::w_s0], W_k)

38# (H/w_s0 x W/w_s0 x Dq)

39K = K.unfold(1, size=w_d1, step=w_s1)\

40 .unfold(2, size=w_d1, step=w_s1)

41# (H/w_s x W/w_s x w_d1 x w_d1 x Dq)

42K = K[:, :, :, ::w_r, ::w_r] + R

43# (H/w_s x W/w_s x w_d x w_d x Dq)

44K = K.reshape(-1, w_d, w_d, Dq)

45# (H*W/w_s**2 x w_d x w_d x Dq)

46

47S = batch_tensordot(norm(K, dim=3), norm(Q, dim=1))

48# (H*W/w_s**2 x w_d x w_d)

49S = S.reshape(H/w_s x W/w_s x w_d x w_d)

50# (H/w_s x W/w_s x w_d x w_d)

Listing 1: Pseudocode. CSCL special case for reduced space complexity

The lack of optimised code for running the operations used by CSCL poses a
major impediment for the implementation of the method on large size output
tensors as the memory required to keep all operations in parallel can be
prohibitively large. We note that when either of $w_{r}$, $w_{s}$ is an
integer multiple of the other and $w_{r}>1$, $w_{s}>1$ there are some
locations in $\mathbfcal{Y}$ which are never used for calculating
$\mathbfcal{S}$, facilitating a more efficient implementation. In essence it
is possible to exclude these locations by striding the input tensor by
$w_{s}^{0}$ before the application of CSCL similarity calculation. This
implementation reduces the space complexity by a factor of
$\max(\frac{w_{r}}{w_{s}},\frac{w_{s}}{w_{r}})^{N}$, where $N$ is the number
of spatial dimensions, with respect to the naive algorithm. The process is
described in python pseudocode in algorithm 1 for a 2D embedding space
$\mathbfcal{Y}\in\mathbb{R}^{H\times W\times D}$. Fig.17 shows a visual
example how this works for the two distinct subcases. The same process can be
applied for generating ground truths $\mathbfcal{L}^{s}$.

Figure 17: Visual representation of the special case for CSCL when
$w_{r}>w_{s}$ (left) and $w_{s}>w_{r}$ (right). Grid points containing numbers
are the ones eventually used during the CSCL computation, empty grid points
never participate in the loss function. The number in each grid point denotes
the $(i,j)$ indices of that point in the spatial extent of the feature map,
e.g. 11 refers to row 1, column 1 of the feature map. Grid points in bold are
used by the top-left CSCL window. Note that following the special case
implementation the same indices are eventually used before (naive
implementation) and after the initial stride.

## References

  * [1] V. S. F. Garnot, L. Landrieu, S. Giordano, and N. Chehata, “Time-space tradeoff in deep learning models for crop classification on satellite multi-spectral image time series,” in _IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium_ , 2019, pp. 6247–6250.
  * [2] V. S. F. Garnot, L. Landrieu, S. Giordano, and N. Chehata, “Satellite image time series classification with pixel-set encoders and temporal self-attention,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020.
  * [3] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM Rev._ , vol. 51, no. 3, p. 455–500, Aug. 2009. [Online]. Available: https://doi.org/10.1137/07070111X
  * [4] F. Babiloni, I. Marras, F. Kokkinos, J. Deng, G. Chrysos, and S. Zafeiriou, “Poly-nl: Linear complexity non-local layers with 3rd order polynomials,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , October 2021, pp. 10 518–10 528.
  * [5] X. Jinru and B. Su, “Significant remote sensing vegetation indices: A review of developments and applications,” _Journal of Sensors_ , vol. 2017, pp. 1–17, 01 2017.
  * [6] C. Conrad, S. Fritsch, J. Zeidler, G. Rücker, and S. Dech, “Per-field irrigated crop classification in arid central asia using spot and aster data,” _Remote Sensing_ , vol. 2, no. 4, pp. 1035–1056, 2010. [Online]. Available: https://www.mdpi.com/2072-4292/2/4/1035
  * [7] J. M. Peña-Barragán, M. K. Ngugi, R. E. Plant, and J. Six, “Object-based crop identification using multiple vegetation indices, textural features and crop phenology,” _Remote Sensing of Environment_ , vol. 115, no. 6, pp. 1301–1316, 2011. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0034425711000290
  * [8] P. Hao, Y. Zhan, L. Wang, Z. Niu, and M. Shakir, “Feature selection of time series modis data for early crop classification using random forest: A case study in kansas, usa,” _Remote Sensing_ , vol. 7, no. 5, pp. 5347–5369, 2015\. [Online]. Available: https://www.mdpi.com/2072-4292/7/5/5347
  * [9] C. Conrad, S. Dech, O. Dubovyk, S. Fritsch, D. Klein, F. Löw, G. Schorcht, and J. Zeidler, “Derivation of temporal windows for accurate crop discrimination in heterogeneous croplands of uzbekistan using multitemporal rapideye images,” _Computers and Electronics in Agriculture_ , vol. 103, pp. 63–74, 2014. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168169914000386
  * [10] S. Siachalou, G. Mallinis, and M. Tsakiri-Strati, “A hidden markov models approach for crop classification: Linking crop phenology to time series of multi-sensor remote sensing data,” _Remote Sensing_ , vol. 7, no. 4, pp. 3633–3650, 2015. [Online]. Available: https://www.mdpi.com/2072-4292/7/4/3633
  * [11] C. Pelletier, S. Valero, J. Inglada, N. Champion, and G. Dedieu, “Assessing the robustness of random forests to map land cover with high resolution satellite image time series over large areas,” _Remote Sensing of Environment_ , vol. 187, pp. 156–168, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0034425716303820
  * [12] M. Rußwurm and M. Körner, “Temporal vegetation modelling using long short-term memory networks for crop identification from medium-resolution multi-spectral satellite images,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , 2017, pp. 1496–1504.
  * [13] M. Rußwurm and M. Körner, “Multi-temporal land cover classification with sequential recurrent encoders,” _ISPRS International Journal of Geo-Information_ , vol. 7, no. 4, p. 129, Mar 2018. [Online]. Available: http://dx.doi.org/10.3390/ijgi7040129
  * [14] N. Audebert, B. Le Saux, and S. Lefèvre, “Semantic segmentation of earth observation data using multimodal and multi-scale deep networks,” in _Computer Vision – ACCV 2016_ , S.-H. Lai, V. Lepetit, K. Nishino, and Y. Sato, Eds. Cham: Springer International Publishing, 2017, pp. 180–196.
  * [15] A. Kamilaris and F. X. Prenafeta-Boldú, “Deep learning in agriculture: A survey,” _Computers and Electronics in Agriculture_ , vol. 147, pp. 70–90, 2018. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168169917308803
  * [16] D. Ienco, R. Gaetano, C. Dupaquier, and P. Maurel, “Land cover classification via multitemporal spatial data by deep recurrent neural networks,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 10, pp. 1685–1689, 2017.
  * [17] N. Kussul, M. Lavreniuk, S. Skakun, and A. Shelestov, “Deep learning classification of land cover and crop types using remote sensing data,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 5, pp. 778–782, 2017.
  * [18] R. Rustowicz, R. Cheong, L. Wang, S. Ermon, M. Burke, and D. B. Lobell, “Semantic segmentation of crop type in africa: A novel dataset and analysis of deep learning methods,” in _CVPR Workshops_ , 2019.
  * [19] M. Rußwurm, C. Pelletier, M. Zollner, S. Lefèvre, and M. Körner, “Breizhcrops: A time series dataset for crop type mapping,” _ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences_ , vol. XLIII-B2-2020, pp. 1545–1551, 08 2020.
  * [20] C. Pelletier, G. I. Webb, and F. Petitjean, “Temporal convolutional neural network for the classification of satellite image time series,” _Remote Sensing_ , vol. 11, no. 5, 2019. [Online]. Available: https://www.mdpi.com/2072-4292/11/5/523
  * [21] R. Interdonato, D. Ienco, R. Gaetano, and K. Ose, “Duplo: A dual view point deep learning architecture for time series classification,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 149, pp. 91 – 104, 2019. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0924271619300115
  * [22] Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger, “3d u-net: Learning dense volumetric segmentation from sparse annotation,” _CoRR_ , vol. abs/1606.06650, 2016. [Online]. Available: http://arxiv.org/abs/1606.06650
  * [23] X. SHI, Z. Chen, H. Wang, D.-Y. Yeung, W.-k. Wong, and W.-c. WOO, “Convolutional lstm network: A machine learning approach for precipitation nowcasting,” in _Advances in Neural Information Processing Systems_ , C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, Eds., vol. 28. Curran Associates, Inc., 2015, pp. 802–810. [Online]. Available: https://proceedings.neurips.cc/paper/2015/file/07563a3fe3bbe7e3ba84431ad9d055af-Paper.pdf
  * [24] K. Malkin, C. Robinson, L. Hou, and N. Jojic, “Label super-resolution networks,” in _International Conference on Learning Representations_ , 2019\. [Online]. Available: https://openreview.net/forum?id=rkxwShA9Ym
  * [25] K. Musgrave, S. Belongie, and S.-N. Lim, “A metric learning reality check,” in _Computer Vision – ECCV 2020_ , A. Vedaldi, H. Bischof, T. Brox, and J.-M. Frahm, Eds. Cham: Springer International Publishing, 2020, pp. 681–699.
  * [26] S. Chopra, R. Hadsell, and Y. LeCun, “Learning a similarity metric discriminatively, with application to face verification,” in _Proc. of Computer Vision and Pattern Recognition Conference_. IEEE Press, 2005.
  * [27] J. Bromley, I. Guyon, Y. LeCun, E. Säckinger, and R. Shah, “Signature verification using a ”siamese” time delay neural network,” in _Proceedings of the 6th International Conference on Neural Information Processing Systems_ , ser. NIPS’93. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1993, p. 737–744.
  * [28] K. Q. Weinberger, J. Blitzer, and L. Saul, “Distance metric learning for large margin nearest neighbor classification,” in _Advances in Neural Information Processing Systems_ , Y. Weiss, B. Schölkopf, and J. Platt, Eds., vol. 18. MIT Press, 2006. [Online]. Available: https://proceedings.neurips.cc/paper/2005/file/a7f592cef8b130a6967a90617db5681b-Paper.pdf
  * [29] F. Schroff, D. Kalenichenko, and J. Philbin, “Facenet: A unified embedding for face recognition and clustering,” _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , Jun 2015. [Online]. Available: http://dx.doi.org/10.1109/CVPR.2015.7298682
  * [30] Z. Wu, Y. Xiong, S. X. Yu, and D. Lin, “Unsupervised feature learning via non-parametric instance discrimination,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 3733–3742.
  * [31] A. van den Oord, Y. Li, and O. Vinyals, “Representation learning with contrastive predictive coding,” _CoRR_ , vol. abs/1807.03748, 2018. [Online]. Available: http://arxiv.org/abs/1807.03748
  * [32] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” 2020.
  * [33] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020\.
  * [34] X. Chen, H. Fan, R. Girshick, and K. He, “Improved baselines with momentum contrastive learning,” 2020.
  * [35] P. O. O. Pinheiro, A. Almahairi, R. Benmalek, F. Golemo, and A. C. Courville, “Unsupervised learning of dense visual representations,” in _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, Eds., vol. 33. Curran Associates, Inc., 2020, pp. 4489–4500. [Online]. Available: https://proceedings.neurips.cc/paper/2020/file/3000311ca56a1cb93397bc676c0b7fff-Paper.pdf
  * [36] X. Wang, R. Zhang, C. Shen, T. Kong, and L. Li, “Dense contrastive learning for self-supervised visual pre-training,” 2020.
  * [37] J. Son Chung, A. Senior, O. Vinyals, and A. Zisserman, “Lip reading sentences in the wild,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , July 2017.
  * [38] R. Hadsell, S. Chopra, and Y. LeCun, “Dimensionality reduction by learning an invariant mapping,” in _2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’06)_ , vol. 2, 2006, pp. 1735–1742.
  * [39] R. Salakhutdinov and G. Hinton, “Learning a nonlinear embedding by preserving class neighbourhood structure,” in _Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, M. Meila and X. Shen, Eds., vol. 2. San Juan, Puerto Rico: PMLR, 21–24 Mar 2007, pp. 412–419. [Online]. Available: http://proceedings.mlr.press/v2/salakhutdinov07a.html
  * [40] P. Khosla, P. Teterwak, C. Wang, A. Sarna, Y. Tian, P. Isola, A. Maschinot, C. Liu, and D. Krishnan, “Supervised contrastive learning,” 2020.
  * [41] K. Kamnitsas, D. C. de Castro, L. L. Folgoc, I. Walker, R. Tanno, D. Rueckert, B. Glocker, A. Criminisi, and A. V. Nori, “Semi-supervised learning via compact latent space clustering,” in _Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018_ , ser. Proceedings of Machine Learning Research, J. G. Dy and A. Krause, Eds., vol. 80. PMLR, 2018, pp. 2464–2473. [Online]. Available: http://proceedings.mlr.press/v80/kamnitsas18a.html
  * [42] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in Neural Information Processing Systems 30_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds. Curran Associates, Inc., 2017, pp. 5998–6008. [Online]. Available: http://papers.nips.cc/paper/7181-attention-is-all-you-need.pdf
  * [43] X. Wang, R. Girshick, A. Gupta, and K. He, “Non-local neural networks,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 7794–7803.
  * [44] I. Bello, B. Zoph, Q. Le, A. Vaswani, and J. Shlens, “Attention augmented convolutional networks,” in _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , 2019, pp. 3285–3294.
  * [45] P. Shaw, J. Uszkoreit, and A. Vaswani, “Self-attention with relative position representations,” _CoRR_ , vol. abs/1803.02155, 2018. [Online]. Available: http://arxiv.org/abs/1803.02155
  * [46] P. Ramachandran, N. Parmar, A. Vaswani, I. Bello, A. Levskaya, and J. Shlens, “Stand-alone self-attention in vision models,” in _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds., vol. 32. Curran Associates, Inc., 2019, pp. 68–80. [Online]. Available: https://proceedings.neurips.cc/paper/2019/file/3416a75f4cea9109507cacd8e2f2aefc-Paper.pdf
  * [47] Y. Kalantidis, M. B. Sariyildiz, N. Pion, P. Weinzaepfel, and D. Larlus, “Hard negative mixing for contrastive learning,” in _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, Eds., vol. 33. Curran Associates, Inc., 2020, pp. 21 798–21 809. [Online]. Available: https://proceedings.neurips.cc/paper/2020/file/f7cade80b7cc92b991cf4d2806d6bd78-Paper.pdf
  * [48] J. D. Robinson, C.-Y. Chuang, S. Sra, and S. Jegelka, “Contrastive learning with hard negative samples,” in _International Conference on Learning Representations_ , 2021. [Online]. Available: https://openreview.net/forum?id=CR1XOQ0UTh-
  * [49] O. Csillik and M. Belgiu, “Cropland mapping from sentinel-2 time series data using object-based image analysis,” 05 2017.
  * [50] G. Weikmann, C. Paris, and L. Bruzzone, “Timesen2crop: A million labeled samples dataset of sentinel 2 image time series for crop-type classification,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 4699–4708, 2021.
  * [51] V. S. F. Garnot and L. Landrieu, “Panoptic segmentation of satellite image time series with convolutional temporal attention networks,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , October 2021, pp. 4872–4881.
  * [52] T. Wang and P. Isola, “Understanding contrastive representation learning through alignment and uniformity on the hypersphere,” in _Proceedings of the 37th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, H. D. III and A. Singh, Eds., vol. 119. PMLR, 13–18 Jul 2020, pp. 9929–9939. [Online]. Available: https://proceedings.mlr.press/v119/wang20k.html
  * [53] E. Shelhamer, J. Long, and T. Darrell, “Fully convolutional networks for semantic segmentation,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 39, no. 4, pp. 640–651, 2017.
  * [54] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” ser. Proceedings of Machine Learning Research, F. Bach and D. Blei, Eds., vol. 37. Lille, France: PMLR, 07–09 Jul 2015, pp. 448–456. [Online]. Available: http://proceedings.mlr.press/v37/ioffe15.html
  * [55] A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” in _in ICML Workshop on Deep Learning for Audio, Speech and Language Processing_ , 2013.
  * [56] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014, cite arxiv:1412.6980Comment: Published as a conference paper at the 3rd International Conference for Learning Representations, San Diego, 2015. [Online]. Available: http://arxiv.org/abs/1412.6980