Sébastien Giordano

Simon Bailly

Loic Landrieu

Nesrine Chehata

Improved Crop Classification with Rotation Knowledge using Sentinel 1 and 2 Time Series

Keywords: classification, time series, temporal regularization, conditional random fields, agriculture, crop rotation, Sentinel, radar

Leveraging the recent availability of accurate, frequent and multimodal (radar and optical) Sentinel 1 & 2 acquisitions, this paper investigates the automation of Land Parcel Identification System (LPIS) crop type classification. Our approach allows for the automatic integration of temporal knowledge i.e. crop rotations using existing parcel-based land cover databases and multimodal Sentinel 1 and 2 time series. The temporal evolution of crop types was modeled with a linear-chain Conditional Random Field, trained with time series of multi-modal (radar and optical) satellite acquisitions and associated LPIS. Our model was tested on two study areas in France (≥1250 km 2 ) which show different crop types, various parcel sizes, and agricultural practices: the Seine et Marne and the Alpes de Haute-Provence classified accordingly to a fine national 25-class nomenclature.

We first trained a Random Forest classifier without temporal structure to achieve 89.0% overall accuracy in Seine et Marne (10 classes) and 73% in Alpes de Haute-Provence (14 classes). We then demonstrated experimentally that taking into account the temporal structure of crop rotation with our model resulted in an increase of 3% to +5% in accuracy.

This increase was especially important (+12%) for classes which were poorly classified without using the temporal structure. A stark positive impact was also demonstrated on permanent crops while it was fairly limited or even detrimental for annual crops.

Introduction

The Sentinel 1&2 satellites provide open and free acquisitions exhibiting unprecedented characteristics which are well-suited to agriculture monitoring: high temporal frequency (5-6 days), the complementary C-Band Sentinel-1 radar images and multispectral Sentinel-2 images with relevant spectral bands to crop mapping, and high spatial resolution (10-20 m).In Europe, several cases of agricultural monitoring using Sentinel images have been proposed [START_REF]Towards future Copernicus service components in support to agriculture[END_REF], such as monitoring crops (crop area estimates, crop map products, crop phenology indicators) , and controlling Common Agricultural Policy (CAP) payments (permanent grasslands, greening measures, . . .). Sentinel Images have also been used for updating and controlling the quality of the Land Parcellar Identification System (LPIS),. a geographical information system on agricultural parcels, at a national scale , updated annually [START_REF] Boryan | Monitoring us agriculture: the us department of agriculture, national agricultural statistics service, cropland data layer program[END_REF].

This paper introduces a tool for automated LPIS crop type classification from multimodal Sentinel time series which incorporate knowledge from existing LPIS editions to improve its accuracy.

Multi-temporal satellite images for crop mapping

Satellite time series are particularly well-suited for identifying different crop types, as they allow for the monitoring of the evolution of the plants phenology. This is particularly crucial in the growing or harvest seasons. Synthetic Aperture Radar (SAR) data are crucial as well, as they mitigate the effect of cloud cover.

Many studies have demonstrated the potential of multi-temporal Sentinel and Landsat-8 images for crop type mapping [START_REF] Palchowdhuri | Classification of multi-temporal spectral indices for crop type mapping: a case study in coalville, uk[END_REF][START_REF] Veloso | Understanding the temporal behavior of crops using sentinel-1 and sentinel-2-like data for agricultural applications[END_REF][START_REF] Vuolo | How much does multi-temporal sentinel-2 data improve crop type classification?[END_REF][START_REF] Belgiu | Sentinel-2 cropland mapping using pixel-based and object-based time-weighted dynamic time warping analysis[END_REF][START_REF] Ottosen | Remote sensing of cropping practice in northern italy using time-series from sentinel-2[END_REF][START_REF] Defourny | Near real-time agriculture monitoring at national scale at parcel resolution: Performance assessment of the sen2-agri automated system in various cropping systems around the world[END_REF] and the contribution of SAR time series for crop monitoring [START_REF] Whelen | Use of time-series l-band uavsar data for the classification of agricultural fields in the san joaquin valley[END_REF][START_REF] Li | Full year crop monitoring and separability assessment with fully-polarimetric l-band uavsar: A case study in the sacramento valley, california[END_REF].

Inglada et al. [START_REF] Inglada | Assessment of an operational system for crop type map production using high temporal and spatial resolution satellite optical imagery[END_REF] assessed state-of-the-art methods for automatic crop mapping with multi-temporal and high resolution optical images. Five different classification approaches using SPOT4 and Landsat-8 images were compared, for 6 annual crops, over 12 different study areas with the best results (OA= 80%) obtained using the Random Forest classifier. In Kussul et al. [START_REF] Kussul | Parcel-based crop classification in Ukraine using Landsat-8 data and Sentinel-1a data[END_REF], Landsat-8 and Sentinel-1 time series were used on a study area in Ukraine. A pixel-based classification combined with a parcel-based regularization (majority voting) was proposed using LPIS ancillary data. An Overall Accuracy (OA) of 89% was reached, but only on a nomenclature comprised of 6 annual crops and large parcels (> 250 ha). [START_REF] Steinhausen | Combining sentinel-1 and sentinel-2 data for improved land use and land cover mapping of monsoon regions[END_REF] combined Sentinel-1 and 2 to improve land cover mapping in cloud-prone regions. [START_REF] Veloso | Understanding the temporal behavior of crops using sentinel-1 and sentinel-2-like data for agricultural applications[END_REF] showed the importance of radar data for crop mapping. More recently, the Sen2-Agri consortium has led experiments at the country level (Czech Republic) using Sentinel time series for crop mapping [START_REF]Sen2-Agri. Czech agriculture national demonstrator -final report[END_REF]. A multi-sensor (Sentinel-1, Sentinel-2) pixel-based supervised classification was performed. The LPIS was used for both learning and validation steps. Monthly cropland maps were produced with an overall accuracy greater than 80%, and each land cover type had a F-score greater than 60%. In [START_REF] Defourny | Near real-time agriculture monitoring at national scale at parcel resolution: Performance assessment of the sen2-agri automated system in various cropping systems around the world[END_REF], three entire countries (Ukraine, Mali and South Africa) and five local cities were mapped using Sen2-agri. Overall accuracy values were higher than 90%, and already as high as 80% mid-season. However, only the five major crops were considered for each site. For the Sen2-agri framework, the nomenclature was generally limited to 5-7 classes and did not fully integrate temporal knowledge from existing data.

Crop rotation integration

Crop rotations knowledge can be used to improve agricultural yields [START_REF] Berzsenyi | Effect of crop rotation and fertilisation on maize and wheat yields and yield stability in a long-term experiment[END_REF] and soil quality [START_REF] Douglas L Karlen | Crop rotation effects on soil quality at three northern corn/soybean belt locations[END_REF]. Crop type prediction can also be improved using prior knowledge on crop rotations per parcel since a crop type is strongly correlated to past crop types. Modeling such temporal structures from Sentinel image time series can lead to significant gains in classification accuracy. To take into account crop rotations in crop mapping, two issues have to be addressed: (i) how to model the crop rotations? (ii) how to integrate crop rotations in a land cover classification process?

Two different approaches can be used to model rotations. The first one uses a priori agronomist expert knowledge. The second one consists in automatically learning crop rotations from the statistical analysis of past practices, as found in the LPIS annual archives. This crop rotation knowledge can then be modeled in a probabilistic framework by estimating the transition probabilities between each crop types from past years to the next. Castellazzi et al. [START_REF] Ms Castellazzi | A systematic representation of crop rotations[END_REF] introduced a mathematical framework modeling such transition probabilities to predict crop rotations at the landscape scale. Few studies have yet focused on the integration of crop rotation information into classification pipelines. Osman et al. [START_REF] Osman | Assessment of a Markov logic model of crop rotations for early crop mapping[END_REF] studied early crop mapping using Markov logic, but not in combination to remote sensing observations.This model proved efficient for early crop type predictions at the beginning of the growing season, when few satellite images are available and crops are hard to distinguish. Other studies proposed to introduce a temporal structure, using Hidden Markov Chains, in a classification pipeline but aimed at modeling phenology instead of crop rotations [START_REF] Aurdal | Use of hidden Markov models and phenology for multitemporal satellite image classification: Applications to mountain vegetation classification[END_REF][START_REF] Beatriz | Hidden Markov Models for crop recognition in remote sensing image sequences[END_REF][START_REF] Siachalou | A hidden Markov models approach for crop classification: Linking crop phenology to time series of multi-sensor remote sensing data[END_REF]. Kenduiywoa et al. [START_REF] Bk Kenduiywoa | Spatial-temporal Conditional Random Fields crop classification from Terrasar-X images[END_REF] modeled phenology information into a Conditional Random Field (CRF), but the classification was performed at different dates through the year. The CRFs were used for classifying land cover classes and crop types on mono-temporal Landsat data [START_REF] Roscher | Kernel discriminative random fields for land cover classification[END_REF]. Hoberg et al. [START_REF] Hoberg | Conditional random fields for multitemporal and multiscale classification of optical satellite imagery[END_REF] used CRF on multi-temporal and multi-scale classification for change detection.

Objectives

This paper focuses on crop type prediction using LPIS and crop rotation knowledge learned from Sentinel-1 & 2 time series. This raises three main applications and methodological questions: 1) what are the respective contribution of optical and radar time series for crop type prediction ? 2) How to combine crop rotation and satellite observations into an unified classification pipeline? 3) What is the contribution of the temporal structure with regard to observation-based classification?

To answer these questions, we propose a series of experiments on the area of interests using a detailed nomenclature with 25 classes with no assumption on parcel sizes. We test the contribution of both optical and radar time series separately and jointly. We then compare the accuracy of classifying these time series each year independently and integrating the temporal structure into a probabilistic model (linear chain Conditional Random Fields) representing the influence of crop rotations.

This paper is organized as follows: we present the study sites and data in Section 2. In Section 3, an observation-based classification at the parcel scale is presented as well as a temporal-structured framework to integrate crop rotation information. Results are given and discussed in Sections 4 and 5, respectively.

Sites and material

Study sites

Two complementary sites were chosen with French territory. Both sites are research observatories where in-field crop type annotations are made annually. The location and characteristics of each site are provided in Figure 1

Land Parcellar Identification system

In France, the Land Parcellar Identification System is called Registre Parcellaire Graphique (RPG). It is available on the whole territory since 2002. For cultivated areas, the RPG gathers the geometric information (i.e., the parcel delimitation) and the corresponding semantic information such as the owner, the operator, the area and the crop type.

Both geometric and semantic information are updated by the farmers annually. Until 2014, the declarations were made at a block scale that corresponds to contiguous parcels with the same operator. Since 2015, the declarations have been made at the parcel scale, which simplifies machinelearning based approaches to crop prediction. The crop type is specified among more than 300 sub-classes which are organized into 25 classes. 14 and 10 of these classes are present on Site04 and Site77, respectively (cf. Table 1). ). As for Site77, two dominant crops are present: cereals (57.7%), meadows (28.3%) followed by vegetables (5.1%). In this latter case, the data is highly imbalanced, making Site77 classification task more complex.

In this study, to be in tune with the first Sentinel-2 images availability, only the 2016 edition of parcel-based RPG was used for the training and the validation of the supervised classification model (Section 3.1). The 2015 parcels were necessary to train the temporal structured method (Section 3.2). For learning crop rotations, only the geometrically stable blocks of parcels from 2010 to 2014 were used. The number of geometrically stable parcels for both sites is given in Table 1.

Multimodal Sentinel-1 & 2 images

We use both optical and radar Sentinel time-series for crop mapping. Sentinel-2 (S2) provides 10 multispectral bands for earth observation on the VIS-SWIR domain at 10 m and 20 m spatial resolution. Near infra-red (NIR) and red-edge bands allow a fine characterization of crops. Sentinel-1 (S1) is a C-band SAR. The available mode on the studied sites was the Interferometric Wide (IW) mode that presents a dual polarization VV and VH.

For the year 2016, Sentinel-2 images were automatically downloaded from the Theia platform (https://www.theia-land.fr/) in tiled format, calibrated as Top of Canopy (TOC) reflectance [START_REF] Hagolle | A multi-temporal method for cloud detection, applied to formosat-2, venÂţs, landsat and sentinel-2 images[END_REF] and accompanied with robust cloud mask information [START_REF] Hagolle | A multi-temporal and multi-spectral method to estimate aerosol optical thickness over land, for the atmospheric correction of formosat-2, landsat, venÎijs and sentinel-2 images[END_REF]. The 20m Sentinel-2 images were resampled to 10 m spatial resolution. Radar Sentinel-1 images were downloaded from the Peps platform (https://peps.cnes.fr) in the Ground Range Detected format (GRD) which corresponds to the average of approximately five Single Look Complex ((SLC) acquisitions corrected by the incidence angle and resampled at 10 m spatial resolution.

The total number of images is given in Table 2 and confirms the complementarity between Sentinel-1 and Sentinel-2 images. On Site77, a particularly high cloud coverage led to only 12 Sentinel-2 images. In addition, S2 orbit over Site77 had particularly many acquisition problems in 2016. On the contrary, on this site, ascendant S1 images overlap, leading to more available radar images.

Sentinel images pre-processing

Figure 5 illustrates the optical and radar pre-processing steps to obtain parcel-based features for the classification workflow.

The dual polarization GRD S1 images were first calibrated to σ 0 radar backscattering coefficient. The orthorectification was performed using the SRTM digital terrain model. The speckle is partly removed using a simple 5 × 5 Lee filter [START_REF] Lee | Digital image enhancement and noise filtering by use of local statistics[END_REF]. In addition to VV and VH radar features, an extra radar feature (

σ 0 V H σ 0 V V
) was derived. This ratio is known to be more robust to acquisition system errors or environmental factors such as soil moisture leading to a more stable temporal indicator [START_REF] Veloso | Understanding the temporal behavior of crops using sentinel-1 and sentinel-2-like data for agricultural applications[END_REF]. Average and standard deviation of these three features were then computed for each date and for each parcel. The number of radar features is shown in Table 2 .

Sentinel-2 images were already orthorectified and calibrated in TOC reflectance. On Site77, only 12 Sentinel-2 optical images were obtained in 2016 as shown in Figure 6 with corresponding cloud cover whereas 23 images were available on Site04 (cf Figure 7). The missing data (clouds) were filled using a multi-temporal spline interpolation [START_REF] Inglada | OTB Gapfilling, a temporal gapfilling for image time series library[END_REF]. Average and standard deviation of the 10 spectral bands and the NDVI (Normalized Difference Vegetation Index) per optical image were then computed for each date and for each parcel. The number of optical features is shown in Table 2.

F M A M J J A S O

Methodology

Our method proceeds in two steps: parcel-wise classification and temporal modeling. The first step aims to predict the crop types per parcel using Sentinel time-series for each year independently. The second step integrates the temporal structure into a probabilistic structured model representing the influence of crop rotations. Modeling this temporal dependency may help correct erroneous classification made in the first step, and may also help classify ambiguous parcels by considering crops from previous years.

Parcel-wise multi-source classification

We first compute discriminative parcel-based features from satellite timeseries. For each parcel, we consider all available optical and radar images for one year. We consider the average and standard deviation of each spectral feature over the pixels composing the parcel extent. We then concatenate the observations over the span of a year of acquisitions. A Random Forest classifier provides a parcel-wise prediction under the form of pseudo-probabilities.

For a given parcel i and a given year t, we denote X (t) i ∈ R D the tensor of combined selected features, with D the selected feature size. Since we compute both mean and standard deviation for each channel at each time step, D = 2 × C × S, with C the number of channels, and S the number of acquisitions per year. To counterbalance the over-representation of certain classes in our data sets, we set class weights inversely proportional to the square root of the number of instances in each class. This class weights are used by the random forest classifier to give more importance to rare classes and recover them more easily.

Temporal-structured classification

We now consider the year-by-year temporal structure of each parcel independently. We denote by X i ∈ R T ×D the sequence of observed features X (t) i ∈ R D for parcel i and for all years t = 1, • • • , T available for training. Likewise, we denote Y i ∈ K T the ground truth labels for parcel i for each observed year t = 1, • • • , T and with K the set of all possible labels.

In section 3.2.1, we present the graphical model chosen to capture crop rotation. In section 3.2.2, we explain how the parameters of this model can be learned from previous LPIS editions. In section 3.2.3, we explain how to use our model to compute prediction of the label of a parcel at a given date.

Temporal structure

The aim of this step is to model the yearly crop rotations in order to improve crop type prediction. We model this dependency with a linear chain Conditional Random Field (CRF) of order m, as shown in Figure 8. For a parcel i, we model the posterior distribution P (Z i | X i ) of predicted labels Z i ∈ K T given the observed features X i as:

P (Z i | X i ) = 1 A exp T t=1 O(Z (t) i , X i ) + T t=m+1 I(Z (t-m) i , • • • , Z (t) i , X) , ( 1
)
where A is a normalizing factor, O the observation potentials, and I the interaction potentials, described below.

Observation potential:

The observation potentials model the link between the observed features and the label of each parcel. O(Z (t)

i , X i ) is taken as the logarithm of P RF (Z (t) i , X i ), the pseudo-probability for parcel i at year t to be class Z

(t) i

given by the random forest classifier, described in Section 3.1:

O(Z (t) i , X i ) = log P RF (Z (t) i , X i ).
(

) 2
Interaction potential: This potential models the temporal dependencies between the parcel's labels at a given year t given the labels at the m previous years. We model this potential as the logarithm of the transition probability M (Z

(t-m) i , • • • , Z (t) i ) from the sequence of m previous labels Z (t-m) i , • • • , Z (m) i
to a label Z t i at the current date [START_REF] Liu | Using local transition probability models in Markov random fields for forest change detection[END_REF]. For the sake of simplicity, we choose a temporally homogeneous parametrization, independent of the observed features, and shared by all parcels and years:

I(Z (t-m) i , • • • , Z (t) i , X) = log M (Z (t-m) i , • • • , Z (t) i )) , (3)
with

M ∈ R K m+1 a tensor such that z 1 ,••• ,z m-1 ∈K m-1 M z 1 ,••• ,z m-1 ,zm = 1 for all z m ∈ K, ie.
a stochastic tensor of order m. The tensor M is referred to as the transition tensor.

Z (1) i ... Z (t-2) i Z (t-1) i Z (t) i X (1) i ... X (t-2) i X (t-1) i X (t) i

Learning

The observation potential is obtained by training the random forest classifier. A transition tensor M can be learned from labeled data over past years. Indeed, maximizing the log-likelihood in Equation 1 with respect to M yields the following tensor M , defined for all

z 1 , • • • , z m ∈ K m : Mz 1 ,••• ,zm = N z 1 ,••• ,zm N z 1 ,••• ,z m-1 , (4)
with

N k 1 ,••• ,km the number of sequences k 1 , • • • , k m observed
in the labeled data for all parcels and all years, and

N k 1 ,••• ,k m-1 the number of sequences k 1 , • • • , k m-1 observed in the first T -1 years
, where T is the total number of years available for training. Excluding the last year is necessary to ensure that M is indeed a stochastic tensor.

To account for the large size of this matrix (|K| m+1 ), and to prevent numeric issues, we perform a Laplacian smoothing with α = 1 as described in Manning et al. [20, 11.3.2].

Inference

The aim of this step is to predict the label Z (t) i of a given parcel at year t from the observation of the current year, and knowing its labeling in the m previous years. Once the random forest yielding the observation potential is trained on all available data, and the transition tensor M estimated, the prediction is given by: [START_REF] Ms Castellazzi | A systematic representation of crop rotations[END_REF] and normalizing the results over k ∈ K to obtain a probability.

p(Z (t) i = k | Z (t-m,••• ,t-1) i , X i ) ∝ P RF (Z (t) i , X i ) × M [Z (t-m) i , • • • , Z (t-1) i ],

Experimental setup

The random forest classifier is composed of 100 decision trees. The meta-parameters of the forest, such as the maximum number of attributes considered at each node, are chosen by k-fold cross-validation with k = 4.

For the temporal structure, spatio-temporal homogeneity hypothesis allows us to estimate the transition tensor M . For each study site, we use the geometrically stable parcel blocks over a period of 5 years (2010-2014). To decrease the number of parameters, only first order transitions were used (transition from one year to the next).

The data is randomly split equally into two distinct training and testing sets. The model is trained and validated on the training set while the quality of the model is estimated on the testing set. The overall accuracy (OA) is used to assess the general performance of the model. The F-score combines User accuracy (UA or precision) and Producer accuracy (PA or recall) and allows estimating the per-class quality. The F-score for a class C is defined as follows:

F-score(C) = 2 * U A c * P A c U A c + P A c (6)
To sum up this information, the F-scores are averaged, with and without weighting by each class cardinality. The weighted F-scores reduce biases due to imbalanced data. The results are averaged over 10 runs.

Results

The results are presented on two distinct agricultural sites (≥ 1250 km 2 ), showing different crop types with highly imbalanced data, various parcel sizes, different agricultural practices and acquired by different Sentinel image distributions. Prediction accuracies are presented using different feature combinations and both with and without temporal modeling. The impact of parcel size on the classification accuracy is also studied.

Transition matrix assessment

Figure 9 shows the estimated transitions between crop types as Hinton diagrams, for both sites. First order transitions are normalized by the number of parcels of the initial class (year n-1) which ensures to be non-sensitive to imbalanced data. Besides, in case of missing crops, a smoothing is processed to avoid zeros in the transition matrix (cf. Section 3.2.2). This smoothing may lead to some biases on minor classes but do not impact the results. On Site04, the most probable transitions are to and from permanent crops, such as olive groves, vineyards, orchards, permanent meadows and fruit trees reaching 98.34%, 93.87%, 92.72%, 91.89% and 84.23%, respectively. From figure 9, we can observe that the standard rotation patterns of annual crops are generally not applied in this area. For instance, the rape seed, proteins and sun flowers have probabilities of 76.53%, 66.78% and 64.25%, respectively to be transformed to other cereals the following year.

Site04 Site77

On Site77, more transitions are observed for the annual crops. Agricultural rules for annual crop rotations seem to be better followed in this area. The rape seed and proteins have probabilities of 97.09% and 94.85%, respectively, to be transformed to other cereals the following year. Indeed, rapeseed → winter wheat (in other cereals) → barley is a well-known 3-year rotation for farmers of this area. Permanent crops such as meadows and fruit trees have a probability of being carried over the next year of 94.45% and 81.39%, respectively.

Optical Vs. Radar Sentinel Time Series

Overall accuracy and F-scores, using different configurations of optical and radar data, are displayed in Tables 3 and4 for Site04 and Site77, respectively.

As reported in Site04 (cf. Table 3), optical data lead to better results than radar data (+9% for OA and +10% for weighted F-score). This may be explained by a low cloud cover in this area (Figure 7) and a finer native resolution of optical imagery that is more suited to small parcel sizes.Table 3 confirms that optical and radar combination led to the best results when not modeling the temporal structure.

Contrarily to the previous site, on Site77 (cf. Table 4), radar attributes improved the results of optical ones by 7% achieving an overall accuracy of 89%. This can be explained by a combination of frequent acquisition problems and a high cloud cover in 2016, leading to many missing optical Sentinel-2 data (cf. Figure 6). In addition, the parcels on Site77 are larger and thus more compatible with radar Sentinel-1 image spatial resolution. Consequently, using radar imagery solely led to similar results when combining optical and radar attributes.

Weighted F-scores, using combined radar and optical images, reached 88% and 71% on Site77 and Site04, respectively. The crop type mappings and prediction errors are illustrated on Figures 10,11 and 12,[START_REF] Douglas L Karlen | Crop rotation effects on soil quality at three northern corn/soybean belt locations[END_REF] for Site04 and Site77, respectively for test parcels. 5 and Table 6 display the F-score, the user and producer accuracy measures per class for both approaches with and without temporal modeling using combined radar and optical images on Site04 and Site77, respectively.

On Site77, from Table 6, one can see that high F-scores are obtained for annual crops (Corn (94%), Barley (90%), Other cereals (95%), Rape seed (96%), . . . ). However, forage crops and fruit trees were not identified by the classifier and were instead classified as meadows (cf. Table 8). Forage crops class is harder to classify using only satellite observations since it is more 8 shows that the meadows class is often confused with other classes, particularly fruit trees and other cereals. The quality of the prediction is lower for Site04 , with a weighted F-scores varying between 61% and 71%. This is due to the presence of more classes [START_REF] Bk Kenduiywoa | Spatial-temporal Conditional Random Fields crop classification from Terrasar-X images[END_REF], small size of parcels, and highly imbalanced classes where meadows and other cereals represented almost 50% of the area as shown in Figure 3. Similarly to Site77, the confusion matrix for combined optical and radar data only large parcels (area >3 Ha), overall accuracy improves by 15%, 5.9% and 8.7% for radar, optical, and combined optical/ radar attributes, respectively (Figure 14). Indeed, due to the limited spatial resolution of Sentinel-1 images, radar attributes are less robust on small parcel sizes. When considering parcels larger than >3 Ha, radar images achieve similar results as optical images, with an overall accuracy of 79.3%. Finally, combining optical and radar data lead to better results in all cases, and especially for large parcels.

As for Site77, similarly to Site04, when considering only parcels larger than >3 Ha, the overall accuracy is greatly improved by radar attributes (+ 8.1% reaching 97.1%), which confirms that radar images are less robust for smaller parcels.

Optical Vs. Radar Sentinel Time Series

Depending on the image distribution and the cloud cover, optical data may not lead to good results for crop type mapping. In our case, temporal missing data interpolation was used, which led to some uncertainties and decreased the classification accuracy. On the other hand, radar data are less robust for small parcels (≤ 1.5 Ha). However, this issue can be reduced by refining the preprocessing framework of radar data.

These results confirm that combining optical and radar data ensures higher crop type prediction accuracy, and lead to more robust prediction, independently from the study site.

Impact of parcel size

Weighted F-scores on both sites are highly dependent on parcel sizes as detailed hereby. For both sites, when keeping only large parcels (area >3 Ha), overall accuracies are improved by radar attributes (+15% and + 8.1% for Site04 and Site77, respectively). Indeed, due to the limited spatial resolution of Sentinel-1 images, radar attributes are less robust on small parcel sizes.

In order to make radar data more robust to parcel of limited sizes, some improvements could be undertaken on radar data pre-processing. To this end, we used speckle filtering [START_REF] Lee | Digital image enhancement and noise filtering by use of local statistics[END_REF] on a restricted local neighborhood (5×5). This is suitable for large parcels as the radar scattering coefficients are averaged afterwards at the parcel level. However, when the parcel area is too small with respect to the Sentinel-1 spatial resolution, this method is no longer suitable. Adaptive radar speckle filtering to small objects should be investigated [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolutionpreserving (pol)(in) SAR denoising[END_REF].

For optical attributes, the relation between parcel size and accuracy is less pronounced. Many factors may impact the overall accuracy such as cloud cover, data imbalance, and the parcel size. Indeed, some classes are more represented in small parcels (0.5-1.5 ha) and are well identified such as Rape seed, Protein, meadows and Fiber plants. Removing these small parcels may decrease the overall accuracy. Combining optical and radar data lead to better results in all cases, and especially for large parcels.Finally, this sensitivity study confirms the robustness of combined radar and optical data to the parcel size.

Impact of temporal structure

The modeling of temporal structure, i.e crop rotation modeling, improved the global prediction accuracies on both sites. On Site04 (cf. Table 5), temporally-structured classification improved the overall accuracy and the weighted F-score by 5% and 1% respectively. When considering only radar data, unstructured accuracies were very low with an overall accuracy of 64% and a weighted F-score of 61%. The temporal modeling approach improved corresponding OA and weighted F-score by 12% and 9%, respectively, confirming the contribution of temporal structure even if the accuracy of the parcel-wise prediction was low. As for Site77, the structured approach slightly improved the weighted F-scores by 3%, 4% and 3% for radar, optical, and combined optical/radar attributes, respectively. The contribution of temporal structure is lower than that in Site04, as the initial parcel-wise accuracies were already high (weighted F-scores >0.88).

As seen in Table 5, temporal structure modeling significantly improved per-class accuracies of permanent crops (fruit trees +36%, vineyards +30%, olive groves +25%, aromatic groves +19% for Site04 ), which reached F-scores higher than 93%. The meadows class F-score was improved by 4%. These results were expected since the permanent crops have the highest transition probability as shown in Section 4.1. However, F-scores of annual crops classes decreased when using temporal structure; only slightly so for corn (-11%), other cereals (-7%), and sunflowers(-8%) but rape seed(-21%), barley (-22%) , other oilseeds (-16%), protein(-49%) and forage crops(-46%) are significantly more often misclassified (cf. Table 7). This may be explained, in our opinion, by two facts: First, the crop rotations are modeled by a temporal regularization between the observation-based term (classification without temporal structure) and the crop transition probabilities. The prediction is a trade-off between both data and regularization terms. If the data-term is high (good prediction with observations), adding crop rotation information does not impact the results significantly ( as for Corn, rape seed and sun flower classes). Secondly, this may be due to the fact that the first order transitions of annual crops are less stable and highly variable with agricultural practices and operators in this area. Hence temporal modeling does not add useful information, and may even wrongly correct an initially correct parcel-wise prediction.

As for Site77, similarly to Site04, the best improvements occur on permanent crops such as meadows and fruit trees (cf. Table 5). Moreover, the temporal structure improved the prediction of some annual crops such as other cereals, Rape seeds and proteins since they have a high first order transition probability to other cereals. The prediction of forage crops is also highly improved using crop rotations information.

On both sites, including rotation knowledge improved the overall accuracy of crop classification. The proposed model is a trade-off between observation-based classification and temporal regularization using learned rotation kwnoledge. If the precision of observation-based classification is already high and the transition patterns inconclusive or poorly followed, integrating rotation knowledge may decrease the accuracy. However, ambiguous observation-based prediction can be improved by modeling the temporal structure, especially if the temporal aspect is very influential, as for permanent crops or crops alternating with other cereals. The detrimental effect of temporal modeling on some annual crops were can also be explained by the limited availability LPIS (only one edition-2016-was available at the time), which might not be sufficient to model crop rotation schemes occurring over 2 or 3 years. Indeed, in this paper, only first order crop rotation were modeled. However, our approach could be straightforwardly extended to rotations over multiple years, provided more data is available. Further tests should be processed with a newer LPIS edition and over larger areas in order to assess the effect of modeling the temporal structure.

Conclusion and perspectives

This study focused on improving the automatic prediction of crop types using Sentinel 1&2 time series and learned rotation knowledge. This study demonstrated the efficiency of multi-temporal and multi-modal Sentinel (optical and radar) images for crop type classification using a fine nomenclature ( > 10 classes) and without filtering small parcels. The joint use of optical and radar features ensured more stable and accurate results. However, results varied highly depending on sites depending on cloud cover, crop types and parcel size.

We modeled the temporal structure (i.e. crop knowledge) with Conditional random fields and automatically learning the probability of crop rotations from previous LPIS editions. This rotation knowledge markedly improved the prediction of crop types. However, while a positive impact is demonstrated on permanent crops using first order crop transitions, this impact is fairly limited or even detrimental for some annual crops. Higher transition orders should be investigated to confirm the interest of temporal structure for annual crops and larger areas with more representative classes should be used. Finally, thanks to the growing volume of available LPIS data, and the free availability of numerous Sentinel images, deep learning approaches for both parcel-wise feature extraction and temporal modeling should be investigated.

Figure 2 :

 2 Figure 2: Normalized histogram of parcel areas for Site04 and Site77.

Figure 3 :Figure 4 :

 34 Figure 3: Site04 : 2016 RPG parcellar superimposed to a very high resolution Digital Terrain Model.

Figure 5 :

 5 Figure 5: Sentinel-1 & 2 pre-processing steps

Figure 7 :

 7 Figure 7: S2 optical images over the year 2016 and corresponding cloud cover on Site77.

Figure 8 :

 8 Figure 8: Graph structure of the temporal dependency at order 2.

Figure 9 :

 9 Figure 9: Representation of the transition matrices with a Hinton diagram.

Figure 10 :

 10 Figure 10: Site04 : Crop type predictions

Figure 11 :

 11 Figure 11: Site04 : Prediction errors

Figure 12 :Figure 13 :

 1213 Figure 12: Site77 : Crop type predictions

Figure 14 :

 14 Figure 14: Site04 : Impact of parcel size on the overall accuracy of the classification without temporal modeling. x-axis: only parcels whose surface area exceeds the threshold (in Ha) are considered.

Figure 15 :

 15 Figure 15: Site77 : Impact of parcel size on the overall accuracy of the classification without temporal modeling; x-axis: only parcels whose surface area exceeds the threshold (in Ha) are considered.

Table 1 ,

 1 and respectively. The site name refers to the national number of the corresponding administrative department. Site04 is located in South Eastern France, in the Alpes de Haute-Provence region, in the Durance river Valley. It is a representative of Mediterranean cultivated areas. It covers 1050 km2 and is characterized by a highly variable topography, a very fragmented landscape and a high diversity of crop types. Site77 is located near Paris, in the Seine

	-5°0′	0°0′	5°0′	10°0′
	50°0′			50°0′
		Site77
	-5°0′
	45°0′			45°0′
			Site04
				10°0′
		0°0′	5°0′
	Figure 1: Localization of Site04 and Site77
	et Marne region.It covers 233 km 2 . Contrarily to Site04, it is characterized
	by a flat relief, with large parcels and a majority of cereal crops. Figure 2
	shows the distribution of parcel sizes on both sites. One can see that Site04
	is much more fragmented with very small parcels while Site 77 has larger
	parcels reaching 20 ha.

Table 1 :

 1 Comparison of both study sites in terms of areas and crop types.

Class nb parcels -Site04 nb parcels -Site77 Corn 147 350

	on Site04 and Site77 respectively, with the corresponding crop types. For
	Site04, dominant crops are: cereals (23.8%), meadows (30.7%), aromatic crops
	(12.2%), forage crops (10.2%) and olive groves (8.6%
	Barley	517	158
	Other cereals	2176	889
	Rape seed	154	85
	Sun flower	293
	Other oilseeds	116
	Protein(peas)	87	76
	Fiber plants		76
	Forage crops	1215	46
	Meadows	3652	725
	Fruit trees	298	30
	vineyards	249
	Olive groves	1029
	Aromatic crops	1452
	Vegetables	520	131
	Total nb classes	14	10
	Total nb stable parcels (2015-16)	9230	1902
	Site area (km 2 )	1050	233
	Figures 3 and 4 show the 2016 RPG edition i.e., the ground truth data

Table 2 :

 2 Characteristics of the parcel-based features for both sites. Features are detailed in the text.

	Site

Nb of dates Optical features Radar features Total

	04	Optical: 23	22 per image	6 per image	Optical: 460
		Radar: 28	σ, µ of	(σ, µ of	Radar: 168
			(10 bands + NDVI) 3 radar features)	total: 628
	77	Optical: 12	22 per image	6 per image	Optical: 240
		Radar: 85	σ, µ of	(σ, µ of	Radar: 509
			(10 bands + NDVI) 3 radar features)	total: 749

Table 3 :

 3 Site04 -Global prediction accuracies, using optical and radar imagery.

		No temporal modeling
	Config	OA F-score Weighted F-score
	Radar	0.64	0.59	0.61
	Optical	0.73	0.67	0.71
	Radar+Optical 0.73	0.68	0.71
		Temporal modeling
	Radar	0.76	0.60	0.7
	Optical	0.78	0.63	0.72
	Radar+Optical 0.78	0.64	0.72

Table 4 :

 4 Site77 -Global prediction accuracy, using optical and radar imagery.

		No temporal modeling
	Config	OA F-score Weighted F-score
	Radar	0.89	0.73	0.88
	Optical	0.82	0.62	0.81
	Radar+Optical 0.89	0.74	0.88
		Temporal modeling
	Radar	0.92	0.78	0.91
	Optical	0.87	0.67	0.85
	Radar+Optical 0.92	0.76	0.91
	4.3. Impact of temporal structure
	Table

Table 5 :

 5 Site04 -Effect of temporal modeling on classification metrics, using combined radar and optical attributes.

		No temporal modeling	Temporal modeling
	Class	F-score User. Acc Proc. Acc F-score User. Acc Proc. Acc
	Corn	0.89	0.95	0.83	0.78	0.84	0.73
	Barley	0.40	0.85	0.26	0.18	0.67	0.11
	Other cereals	0.85	0.83	0.86	0.78	0.73	0.84
	Rape seed	0.92	1	0.86	0.71	1	0.55
	Sun flower	0.79	0.78	0.79	0.71	0.75	0.68
	Other oilseeds	0.57	0.74	0.46	0.41	1	0.26
	Protein	0.49	1	0.32	0	0.07	0
	Forage crops	0.47	0.73	0.34	0.1	0.88	0.05
	Meadows	0.76	0.67	0.87	0.80	0.67	0.98
	Fruit trees	0.61	0.86	0.47	0.97	0.96	0.97
	vineyards	0.69	0.74	0.65	0.99	0.99	0.98
	Olive groves	0.74	0.82	0.67	0.99	0.99	0.99
	Aromatic crops	0.74	0.65	0.86	0.93	0.94	0.92
	Vegetables	0.61	0.72	0.53	0.67	0.9	0.54

Table 6 :

 6 Site77 -Effect of temporal modeling on accuracy metrics, using combined radar and optical attributes.

		No temporal modeling	Temporal modeling
	Class	F-score User. Acc Proc. Acc F-score User. Acc Proc. Acc
	Corn	0.94	0.93	0.95	0.88	0.83	0.93
	Barley	0.90	0.94	0.86	0.82	0.78	0.85
	Other cereals	0.95	0.96	0.94	0.95	0.94	0.97
	Rape seed	0.96	0.97	0.94	0.97	0.98	0.95
	Protein	0.95	0.93	0.97	0.95	0.97	0.94
	Fiber plants	0.97	1	0.95	0	0.1	0
	Forage crops	0	0.1	0	0.70	0.78	0.65
	Meadows	0.87	0.81	0.93	0.95	0.94	0.97
	Fruit trees	0.01	0.1	0	0.94	1	0.89
	Vegetables	0.89	0.91	0.88	0.45	0.97	0.30

Acknowledgments

The authors want to thank the ASP (French payment agency) and the French agricultural ministry for the rich discussions and their guidance for the nomenclature choice and crop rotation schemes.

23

Structured Optical and radar

Discussion

In the following, we will compare the contribution of optical and Radar Sentinel time series and the impact of parcel size. We then detail the impact of modeling the temporal structure and finally explain the impact of the site characteristics on the results.

Classification Without Temporal Modeling

First, we will compare the results of both sites using parcel-based crop type prediction based on image observations only (without temporal modeling).