Giacinto Manfron
email: manfron.g@irea.cnr.it

Sylvestre Delmotte

Lorenzo Busetto

Laure Hossard

Luigi Ranghetti

Pietro Alessandro Brivio

Mirco Boschetti

Estimating inter-annual variability in winter wheat sowing dates from satellite time series in Camargue, France

Keywords: Durum wheat, seasonal distribution, MODIS, time series analyses, weather conditions, preceding crop

Crop simulation models are commonly used to forecast the performance of cropping systems under different hypotheses of change. Their use on a regional scale is generally constrained, however, by a lack of information on the spatial and temporal variability of environment-related input variables (e.g., soil) and agricultural practices (e.g., sowing dates) that influence crop yields. Satellite remote sensing data can shed light on such variability by providing timely information on crop dynamics and conditions over large areas. This paper proposes a method for analyzing time series of MODIS satellite data in order to estimate the inter-annual variability of winter wheat sowing dates. A rule-based method was developed to automatically identify a reliable sample of winter wheat field time series, and to infer the corresponding sowing dates.

The method was designed for a case study in the Camargue region (France), where winter wheat is characterized by vernalization, as in other temperate regions. The detection criteria were chosen on the grounds of agronomic expertise and by analyzing high-confidence time-series vegetation index profiles for winter wheat. This automatic method identified the target crop on more than 56% (four-year average) of the cultivated areas, with low commission errors (11%). It also captured the seasonal variability in sowing dates with errors of ± 8 and ± 16 days in 46% and 66% of cases, respectively. Extending the analysis to the years 2002-2012 showed that sowing in the Camargue was usually done on or around November 1 st (±4 days). Comparing inter-annual sowing date variability with the main local agro-climatic drivers showed that the type of preceding crop and the weather conditions during the summer season before the wheat sowing had a prominent role in influencing winter wheat sowing dates.

Introduction

The management of agricultural systems in European regions varies considerably in space and time due to differences in environmental conditions (e.g., pedo-climatic conditions), available technologies (e.g., crop varieties), agricultural policies (e.g. subsidies, environmental regulations), and market prices [START_REF] Bakker | Variability in regional wheat yields as a function of climate, soil and economic variables: Assessing the risk of confounding[END_REF][START_REF] Dury | Models to support cropping plan and crop rotation decisions[END_REF][START_REF] Pettorelli | Using the satellitederived NDVI to assess ecological responses to environmental change[END_REF]. Monitoring and understanding the diversity and dynamics of agricultural systems on a regional scale is crucial to support their evolution towards a more sustainable future [START_REF] Zheng | Remote sensing of crop residue cover using multi-temporal Landsat imagery[END_REF].

Various tools and approaches can be adopted to monitor agricultural systems and support their adaptation [START_REF] Basso | Review of Crop Yield Forecasting Methods and Early Warning Systems[END_REF], including crop simulation models, which are commonly used to understand the current performance of cropping systems, and to predict future trends under different hypotheses of change [START_REF] Soltani | Comparison of three modelling approaches to simulate regional crop yield : a case study of winter wheat yield in western Germany[END_REF]. These hypotheses may relate to: (1) the development and adoption of new crop varieties more resistant to diseases or thermal stresses, or hybrids, for instance [START_REF] Bregaglio | A set of software components for the simulation of plant airborne diseases[END_REF][START_REF] Webber | Simulating canopy temperature for modelling heat stress in cereals[END_REF]; (2) new cropping systems, such as direct seeding, intercropping, or agroforestry [START_REF] Khaledian | Adapting PILOTE model for water and yield management under direct seeding system: The case of corn and durum wheat in a Mediterranean context[END_REF][START_REF] Miao | Modeling water use, transpiration and soil evaporation of spring wheat-maize and spring wheat-sunflower relay intercropping using the dual crop coefficient approach[END_REF][START_REF] Van Der Werf | Yield-SAFE: A parameter-sparse, process-based dynamic model for predicting resource capture, growth, and production in agroforestry systems[END_REF]; and (3) the impact of climate change on cropping systems, and consequent adaptation strategies, such as changing the sowing dates [START_REF] Holzkämper | Projecting climate change impacts on grain maize based on three different crop model approaches[END_REF][START_REF] Nendel | Testing farm management options as climate change adaptation strategies using the MONICA model[END_REF].

When simulating the future of cropping systems, particularly on a regional scale, a key issue concerns the availability of appropriate information on crucial crop cultivation variables to be provided as input in crop models [START_REF] Moulin | Combining agricultural crop models and satellite observations: From field to regional scales[END_REF][START_REF] Therond | Using a cropping system model at regional scale: Low-data approaches for crop management information and model calibration[END_REF][START_REF] Yuping | Monitoring winter wheat growth in North China by combining a crop model and remote sensing data[END_REF]. Retrieving details on crop management is often challenging [START_REF] Clavel | Integrating expert knowledge and quantitative information for mapping cropping systems[END_REF][START_REF] Jiang | Application of Crop Model Data Assimilation With a Particle Filter for Estimating Regional Winter Wheat Yields[END_REF] because they vary considerably within and between regions, and year by year [START_REF] Therond | Using a cropping system model at regional scale: Low-data approaches for crop management information and model calibration[END_REF][START_REF] Vyas | Extracting Regional Pattern of Wheat Sowing Dates Using Multispectral and High Temporal Observations from Indian Geostationary Satellite[END_REF][START_REF] Yuping | Monitoring winter wheat growth in North China by combining a crop model and remote sensing data[END_REF].

Among the possible sources of information, censuses conducted by governmental agencies or regional institutions provide aggregated data that may not be detailed enough for crop modelling purposes [START_REF] Grassini | How good is good enough? Data requirements for reliable crop yield simulations and yield-gap analysis[END_REF]. Data coming from interviews with farmers are more complete, but limited in number and related to a given year, so they are not always sufficiently representative of the spatial and temporal diversity of a region [START_REF] Chen | Derivation and validation of Canada-wide coarseresolution leaf area index maps using high-resolution satellite imagery and ground measurements[END_REF]. The sowing date is an essential information to the accuracy of simulations obtained with crop models [START_REF] Folberth | Impact of input data resolution and extent of harvested areas on crop yield estimates in large-scale agricultural modeling for maize in the USA[END_REF][START_REF] Van Wart | Creating long-term weather data from thin air for crop simulation modeling[END_REF]. For most crops, the sowing date is crucial to explain crop performance variables, such as yield [START_REF] Delmotte | On farm assessment of rice yield variability and productivity gaps between organic and conventional cropping systems under Mediterranean climate[END_REF][START_REF] Kogan | Winter wheat yield forecasting in Ukraine based on Earth observation, meteorological data and biophysical models[END_REF]. Sowing date variability can span several months for some plant species, depending on the interannual climate variability in the region and on local agricultural practices.

In this setting, analyzing long-term time series of satellite images is an efficient way to elucidate inter-and intra-annual sowing date variability [START_REF] Bradley | A curve fitting procedure to derive interannual phenologies from time series of noisy satellite NDVI data[END_REF][START_REF] Guyet | Long term analysis of time series of satellite images[END_REF]. Earth Observation satellite data archives can be used to estimate this variable as they provide regular and synoptic information on crop characteristics [START_REF] Justice | An overview of MODIS Land data processing and product status[END_REF][START_REF] Kumar | Remote sensing of crop growth[END_REF][START_REF] Rembold | Using low resolution satellite imagery for yield prediction and yield anomaly detection[END_REF].

In particular, time series of spectral vegetation indexes (VIs), such as the Enhanced Vegetation Index (EVI) [START_REF] Huete | Overview of the radiometric and biophysical performance of the MODIS vegetation indices[END_REF], can be used to characterize vegetation dynamics and retrieve information on key phenological stages by various means, such as thresholding, curve fitting, derivatives' analysis, etc. [START_REF] Curnel | Agrophenology Indicators From Remote Sensing: State of the Art[END_REF]. This approach was first used to monitor natural ecosystems and is often termed land surface phenology (LSP) [START_REF] Ganguly | Remote Sensing of Environment Land surface phenology from MODIS : Characterization of the Collection 5 global land cover dynamics product[END_REF].

Many studies have focused on estimating the timing of crop phenological stages like emergence, flowering or senescence, for example (e.g. [START_REF] Sakamoto | A crop phenology detection method using time-series MODIS data[END_REF][START_REF] Boschetti | Multi-year monitoring of rice crop phenology through time series analysis of MODIS images[END_REF][START_REF] Manfron | Testing automatic procedures to map rice area and detect phenological crop information exploiting time series analysis of remote sensed MODIS data[END_REF][START_REF] Pan | Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data[END_REF]. In agricultural applications, a further and more challenging step is to use the same approach to investigte more specific information about farming practices. For instance, agricultural flooding and the duration of rice paddy flooding were investigated by [START_REF] Sakamoto | Detecting temporal changes in the extent of annual flooding within the Cambodia and the Vietnamese Mekong Delta from MODIS time-series imagery[END_REF] and by [START_REF] Ranghetti | International Journal of Applied Earth Observation and Geoinformation Testing estimation of water surface in Italian rice district from MODIS satellite data[END_REF], while crop management practices such as forage cutting were examined by [START_REF] Halabuk | Towards Detection of Cutting in Hay Meadows by Using of NDVI and EVI Time Series[END_REF].

As claimed by [START_REF] Jin | Mapping Irrigated and Rainfed Wheat Areas Using Multi-Temporal Satellite Data[END_REF] it is generally accepted that this is virtually impossible to detect sowing dates on the basis of remote sensing data alone as there is a period after crop sowing and before crop emergence in which the crop cannot be detected with remote sensing. Sowing dates must therefore be inferred from satellite data on the basis of crop-and region-specific assumptions. Successful applications can be found for rice sown/transplanted in flooded fields because the presence of the water prior to any vegetative growth can be detected on satellite images [START_REF] Sakamoto | A crop phenology detection method using time-series MODIS data[END_REF]Boschetti et al., 2015b).

For other crops, previous studies often focused on identifying the so-called "green-up" date (i.e., coinciding with the earliest reliable evidence of vegetation on satellite images), from which sowing dates could be inferred using various methods. In an application for monitoring wheat in northern India, [START_REF] Lobell | Satellite detection of earlier wheat sowing in India and implications for yield trends[END_REF] assumed that green-up coincided with crop emergence and calculated sowing dates by going back a fixed time (three weeks) from this date. Green-up was established as the point where the curve fitted to the VIs reached 10% of its maximum amplitude for the year in question. A similar approach was taken by [START_REF] Vyas | Extracting Regional Pattern of Wheat Sowing Dates Using Multispectral and High Temporal Observations from Indian Geostationary Satellite[END_REF] in the same geographical area, exploiting daily normalized difference VI (NDVI) produced by the Indian geostationary satellite INSAT 3A CCD. Green-up was identified using a threshold method and the sowing date was set at 7 days prior to green-up. These methods proved efficient for the subtropical areas of northern India, where wheat is sown between October and November, then undergoes a rapid, monotonic increase in biomass up to the flowering phase, and the total growing period lasts from 100 to 170 days.

In temperate climates, winter wheat is sown in autumn and takes from 180 to 300 days to mature [START_REF] Asseng | Wheat. In "Crop yield response to water[END_REF]. After emerging, winter wheat undergoes a tillering stage, then requires a period of cold for up to 90 days (until the end of winter), during which it remains dormant. This is followed by a rapid stem elongation and subsequent plant growth. Hence the "double hump" pattern observed in winter wheat time series by [START_REF] Pan | Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index[END_REF] and by [START_REF] Chu | Estimation of winter wheat phenology under the influence of cumulative temperature and soil salinity in the Yellow River Delta , China , using MODIS time-series data 1161[END_REF], makes it rather difficult to use the previously mentioned solutions based on identifying green-up and then backtracking a fixed number of days to estimate the sowing date. The green-up typically detectable from coarse-resolution satellite data relates to the post-dormancy vegetative phase, and the duration of the dormancy period varies both spatially and temporally, depending on crop variety, for instance, and winter weather conditions. [START_REF] Jin | Mapping Irrigated and Rainfed Wheat Areas Using Multi-Temporal Satellite Data[END_REF] analyzed winter wheat in Shanxi province (China), based on images with a 30 m resolution acquired every two days by the HJ-1 a/b multispectral sensor. They could identify two green-up times, one after sowing (before dormancy), and another more robust one preceding the head development phase.

They were thus able to estimate winter wheat sowing dates by applying the relation proposed by [START_REF] Lobell | Satellite detection of earlier wheat sowing in India and implications for yield trends[END_REF] to the first green-up identified. This was made possible by their use of data with a very high spatial and temporal resolution. Such data are not readily available worldwide and none dating back long enough for use it in medium-to long-term analyses. The feasibility of reliably detecting the weaker greenup preceding dormancy is however not demonstrated, particularly when moderate-resolution satellite data are used.

Employing the same satellite data, [START_REF] Pan | Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data[END_REF] used metrics derived from NDVI signals to map the phenological stages of winter and summer crops in China. These authors estimated only the green-up after dormancy for winter wheat, without retrieving sowing dates.

Based on this literature review, we decided to develop a specific algorithm to estimate winter wheat sowing dates in temperate climates from moderate-resolution satellite data. The study was conducted in the Camargue region (France). Durum wheat sowing conditions vary greatly in this region, influenced by environmental conditions, agricultural practices and subsidies, and sowing may take place in appropriate conditions over a time frame of about ten weeks (Mouret J.C., pers. comm., Dec. 5, 2015). No data are currently available on the local intra-and inter-annual variability of winter wheat sowing dates, but this variability needs to be analyzed, and the factors influencing the timing of this operation need to be clarified in order to improve the crop models used to characterize the present conditions and simulate future scenarios (including e.g., climate change). In this context, the aims of the present study are twofold: (1) developing an approach based on remote sensing capable of identifying a robust sample of areas where winter wheat is grown from MODIS time series with a 250 m spatial resolution, and estimating the corresponding sowing dates; (2) analyzing the inter-annual variability of winter wheat sowing dates in the Camargue over the years 2002 to 2012 in relation to meteorological and anthropic drivers.

Materials and methods

The study involved three main phases to develop, apply and validate the proposed method (Figure 1). MODIS data were first preprocessed and representative temporal crop signatures were extracted by exploiting reference information (sections 2.2.1 and 2.2.2). A rule-based algorithm for detecting a robust sample of areas where winter wheat is grown, and identifying the corresponding sowing dates was developed, applied to the years 2002 to 2012 (section 2.2.3) and validated (Section 2.3). The inter-annual variability of the sowing dates thus estimated was then analyzed in relation to meteorological data and main farming practices (section 2.4).

Study area

The Camargue is a region in south-eastern France, at the mouth of the Rhône River. The Natural Regional Park of Camargue (PNRC: Parc Naturel Régional de Camargue in French) is a protected area covering about 108,000 ha and including almost all the land between the two branches of the Rhône (Grande Camargue), plus some areas further east and west (Figure 2). Based on the PNRC land use maps, which also provide in vector format information on individual crops on a scale of 1:5000 (http://www.pnrpaca.org/, last accessed September 2016), the areas being cultivated in 2006 and 2011 accounted for 33% and 29% of the park's surface, respectively (Figure 2 While rice is sown within a narrow time frame (usually from the end of April to the middle of May), the winter wheat sowing dates span a longer period, from late September to early December. The timing depends largely on climatic and anthropic factors: it rains more frequently in autumn, and this can delay sowing; and winter wheat is usually less remunerable than rice in terms of market price and/or common agricultural policy (CAP) subsidies, so farmers tend to give priority to their rice farming schedule (e.g., postponing their wheat sowing if they need to delay harvesting their rice crop) (Mouret J.C., pers. comm., Dec. 5, 2015).

Data processing

2.2.1 MODIS data pre-processing MOD13Q1 (MODIS/Terra) and MYD13Q1 (MODIS/Aqua) v.005 Level-3 composite products with a nominal 250 m spatial resolution (the pixel size at the nadir is 231.65 m, corresponding to an area of 5.37 ha) and sinusoidal projection [START_REF] Didan | MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V005[END_REF] were used to identify winter wheat growing areas and to estimate when they were sown. These products are based on a maximum value compositing method that considers a time span of 16 days and retrieves the best observation obtained during the compositing period [START_REF] Solano | MODIS vegetation index user's guide (MOD13 series)[END_REF]. They are provided with an eight-day nominal time shift between TERRA and AQUA, and combining Science (EROS) Center (https://lpdaac.usgs.gov/), using the MODIStsp package [START_REF] Busetto | Computers & Geosciences MODIStsp : An R package for automatic preprocessing of MODIS Land Products time series[END_REF] of the R software (R Core Team, 2016). The proposed approach is based on the analysis of MODIS EVI time series (Eq. 1).

(Eq. 1)

The EVI was chosen for its well-known suitability for monitoring vegetation dynamics (e.g., [START_REF] Dubovyk | Monitoring vegetation dynamics with medium resolution MODIS-EVI time series at sub-regional scale in southern Africa[END_REF][START_REF] Sakamoto | Detecting temporal changes in the extent of annual flooding within the Cambodia and the Vietnamese Mekong Delta from MODIS time-series imagery[END_REF][START_REF] Sakamoto | A crop phenology detection method using time-series MODIS data[END_REF][START_REF] Son | A phenology-based classification of timeseries MODIS data for rice crop monitoring in Mekong Delta[END_REF]. Information on EVI were extracted from the MOD/MYD 13Q1 product and stacked to create a multi-temporal array. Information was also obtained on EVI quality flags, reflectance in the blue band, the DOY of acquisition, and pixel quality. A weighted Savitzky-Golay smoothing filter (SAV-GOL - [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF] was then applied to the EVI time series to remove residual noise from the data [START_REF] Pettorelli | Using the satellitederived NDVI to assess ecological responses to environmental change[END_REF]. The SAV-GOL filter has been widely used for smoothing time series of MODIS-derived VIs (Boschetti et al., 2015a(Boschetti et al., , 2015b;;[START_REF] Chen | A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[END_REF][START_REF] Geng | Comparison of Eight Techniques for Reconstructing Multi-Satellite Sensor Time-Series NDVI Data Sets in the Heihe River Basin, China[END_REF][START_REF] Manfron | Testing automatic procedures to map rice area and detect phenological crop information exploiting time series analysis of remote sensed MODIS data[END_REF][START_REF] Ren | Regional yield estimation for winter wheat with MODIS-NDVI data in Shandong, China[END_REF][START_REF] White | Intercomparison, interpretation, and assessment of spring phenology in North America estimated from remote sensing for 1982-2006[END_REF].

EVI temporal crop signature extraction

The available PNRC land use maps were used to identify winter wheat fields in the MODIS 250 m pixels covering the study area. Geographic Information System (GIS) analyses were run to extract polygons corresponding to agricultural areas, and then calculate the fractions covered by specific crops in each MODIS pixel for the years 2006 and 2011. MODIS pixels with at least 80% of their surface covered by a particular crop (rice, winter wheat, etc.) were considered homogeneous. These were used to calculate statistical descriptors (median and standard deviation) for each date in the EVI time series, and thus obtain representative temporal signatures for the various crops. Our aim was to identify the more important descriptors for the purpose of distinguishing winter wheat from the other crops, and estimating when it was sown.

A preliminary analysis (see Supplementary material S1) showed that the small green-up identified by [START_REF] Jin | Mapping Irrigated and Rainfed Wheat Areas Using Multi-Temporal Satellite Data[END_REF] using high-resolution data, which coincided with the emergence-tillering period (a few weeks after sowing), could not be pointed reliably on the MODIS time series. The green-up automatically identified on the strength of a 10% increase in the EVI signal [START_REF] Lobell | Satellite detection of earlier wheat sowing in India and implications for yield trends[END_REF][START_REF] Vyas | Extracting Regional Pattern of Wheat Sowing Dates Using Multispectral and High Temporal Observations from Indian Geostationary Satellite[END_REF][START_REF] Dong | Estimating winter wheat biomass by assimilating leaf area index derived from fusion of Landsat-8 and MODIS data[END_REF] was mainly related to the faster growth of the vegetation after dormancy, in spring (see Supplementary material S1). It is therefore not easy to compute the sowing date by going back a given number of days from the green-up date because the duration of the dormancy period depends largely on winter temperatures (vernalization requires a certain number of chilling days). We thus developed and tested a rule-based approach for estimating sowing dates that can cope with: i) the peculiarity of the winter wheat signal in temperate areas (long interval between sowing and green-up); and ii) the specific characteristics of the case study area (fragmented agricultural environment) in relation to moderateresolution satellite data.

Algorithm development and implementation

On the basis of expert agronomic knowledge, a set of rules for automatically detecting winter wheat growing areas and estimating sowing dates was obtained by analyzing the crop's temporal signatures (Figure 3 and Supplementary material S2). The main features of the phenological cycle of winter wheat in the Camargue can be summarized as follows:

 the seeding period is in Autumn (late September to early December), followed by crop emergence after 6-10 days;

 the plant reaches the three-leaf stage after a vegetative period of 2-3 weeks, then the tillering period begins, with a long ensuing period of dormancy in winter (December-February, ending when the need for chilling has been met);

 the tillering phase is completed at the end of the winter, and followed by a green-up phase (February-March), during which the plant rapidly increases in biomass and leaf area;

 the crop heading period, when plants reach their maximum Leaf Area Index, and the subsequent reproductive phase take place between April and May;

 grain filling and ripening, leading to the plant's drying and senescence, generally occur in late June.

The detection of some of these distinctive features in EVI time series (Figure 3) enables winter wheat to be distinguished from other crops [START_REF] Huang | Monitoring Winter Wheat Phenology Using Time Series of Remote Sensing Data[END_REF][START_REF] Lu | Detecting winter wheat phenology with SPOT-VEGETATION data in the North China Plain[END_REF][START_REF] Pan | Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index[END_REF] (see Supplementary material S2). The most significant features for this purpose are crop heading in late spring, and the duration of the crop cycle. In fact, heading generally occurs in May for winter wheat, and in July (nearly two months later) for summer crops. Winter wheat also has a much longer growing cycle -lasting up to 8 months, from autumn to late spring -than summer crops (though it is shorter than for perennial crops like alfalfa). Baresoil conditions due to ploughing in late September to early October were also judged to be an important feature to consider in estimating the sowing date. The above-mentioned features were formalized in rules and implemented in an algorithm to automatically generate, from MODIS time series, raster maps of winter wheat fields and sowing dates.

It is widely accepted in the literature that specific local minima (MIN) and maxima (MAX) points in a time series respectively indicate the start of active crop growth [START_REF] Pan | Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data[END_REF][START_REF] Sakamoto | A crop phenology detection method using time-series MODIS data[END_REF] and of the crop heading period [START_REF] Atzberger | Mapping the spatial distribution of winter crops at sub-pixel level using AVHRR NDVI time series and neural nets[END_REF][START_REF] Curnel | Agrophenology Indicators From Remote Sensing: State of the Art[END_REF]. A pixel was thus recognized as winter wheat if its time series featured a MAX point (crop heading) in the second quarter of the year associated with a MIN point preceding the MAX within a time window compliant with the length of a winter wheat crop cycle.

In particular, a MAX point was assumed to correspond to crop heading if the following requirements were simultaneously met:

1. crop calendar: the MAX occurs in the second quarter of the year (early April to late June);

2. rapid growth: the MAX is preceded by a rapid EVI growth typical of the vegetative phase (February -March), identifiable from the presence of at least three positive derivatives in a temporal window of 40 days (5 composites) before the MAX;

3. senescence: the MAX is followed by a relatively fast EVI decrease typical of the senescence phase, and identifiable from the presence of at least three negative derivatives in a temporal window of 40 days

(5 composites) after the MAX, associated with a drop of at least 33% in the EVI vis-à-vis the MAX value;

4. high biomass: the EVI of the MAX is greater than 0.42 (see Supplementary material S4).

A MIN point was assumed to correspond to the wheat sowing period [START_REF] Pan | Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data[END_REF] because it is characterized by field preparation, i.e. bare soil conditions (ploughing) indicated by a low EVI, and crop growth after sowing corresponding to a moderate increase in EVI. The identification of a suitable MIN was based on the following conditions:

1. crop calendar: the MIN point occurs in a temporal window at least three months before the crop heading date identified;

2. bare soil conditions: the MIN has an EVI lower than 0.3, indicative of bare soil conditions [START_REF] Vyas | Extracting Regional Pattern of Wheat Sowing Dates Using Multispectral and High Temporal Observations from Indian Geostationary Satellite[END_REF];

3. subsequent crop growth: the MIN is followed by the presence of at least three positive derivatives in a temporal window of 40 days (5 composites).

If a suitable MIN is identified, the pixel analyzed is then assumed to correspond to winter wheat.

Additional rules were used to estimate the sowing date near the MIN identified as outlined above. The procedure involved finding the most likely sowing date within a time window ranging between September 30 th (the date that agronomists consider the earliest possible sowing date in the study area) and 140 days before the estimated crop heading (MAX), 140 days being considered by the experts the shortest possible interval between sowing and heading for winter wheat (Mouret J.C., pers. comm., Dec. 5, 2015).

All the MIN points falling within this period are checked according to a "flatness" criterion. Previous experiments (data not shown) showed that "false" sowing dates were identified in line with random minima points generated as a result of residual noise in the smoothed EVI time series obtained for periods containing several consecutive low EVI values. The flatness criterion is formalized by checking the 4 EVI values centered on the local MIN points within the sowing dates search window. A MIN point surrounded by EVI values of similar magnitude (±5%) is identified as noise and discarded. Once these points have been detected and discarded, the earliest MIN point retained is considered the final estimate of the sowing date.

The above-described rule-based algorithm (implemented in IDL-Interactive Data Language v8.3) was applied to the MODIS time series for the years 2002-2012 to generate annual maps of wheat growing areas and associated sowing dates.

Validation

The

Landsat validation of winter wheat growing areas

The MODIS-derived winter wheat growing areas were validated against independent reference data obtained from a visual interpretation of Landsat satellite images for the years 2003, 2006, 2009 and 2012.

The visual interpretation of high-resolution satellite images is a widely accepted practice in remote sensing, especially when in situ data suitable for building reference datasets are lacking. The method is particularly useful when applied to large areas (from regional to continental), where collecting data on the ground becomes too expensive [START_REF] Stroppiana | Integration of Optical and SAR Data for Burned Area Mapping in[END_REF][START_REF] Chen | Derivation and validation of Canada-wide coarseresolution leaf area index maps using high-resolution satellite imagery and ground measurements[END_REF][START_REF] Sun | Winter wheat mapping using temporal signatures of MODIS vegetation index data[END_REF][START_REF] Wulder | Remote sensing of forest environments: concepts and case studies[END_REF][START_REF] Zhao | Towards a common validation sample set for global land-cover mapping[END_REF].

A set of Landsat-7 ETM+ images (spatial resolution of 30 m) was therefore used to build an independent reference dataset for validation purposes. Judging from the land cover information provided by the PNRC, winter wheat is the only winter crop grown in the study area, so Landsat images were selected in the most suitable periods for distinguishing it from summer crops. Two Landsat images were selected for the years 2003, 2006, 2009 and 2012, one obtained around mid-April (when winter wheat is nearing its seasonal heading and summer crops are only just starting to grow), and the other in the second half of July (when winter wheat is harvested and summer crops are in the heading phase) (Table 1). Then 300 MODIS pixels were randomly chosen in the study area and labelled as "winter wheat" or "other crops" based on the visual interpretation of the Landsat images. Where fragmented conditions were detected by Landsat at MODIS pixel level, the pixel was skipped. The random selection of MODIS pixels was then repeated recursively, for the four years considered for validation purposes (2003, 2006, 2009, and 2012), if fewer than 50 MODIS pixels were labelled as winter wheat. The final database of 1200 labelled MODIS pixels was used as the reference data set for assessing winter wheat detection performance using standard validation accuracy metrics calculated from error matrices [START_REF] Congalton | A review of assessing the accuracy of classifications of remotely sensed data[END_REF][START_REF] Brivio | Principi e metodi di telerilevamento[END_REF]:

overall accuracy (OA), class omission and commission errors (OE, CE), and Cohen's kappa coefficient (K) [START_REF] Cohen | A coefficient of agreement for nominal scales[END_REF].

Ground reference validation of sowing dates

The

Analysis of inter-annual variability in sowing dates

The results obtained with the algorithm for the years 2002-2012 were analyzed to: (1) describe the interannual variability in sowing dates; and (2) link this variability with the main drivers influencing farmers' sowing practices, i.e. weather conditions (rainfall, temperature) and preceding crops.

For the first point, two analyses were performed to investigate the diversity/similarity in the distribution of the winter wheat sowing dates obtained for the various years. First, the non-parametric Kolmogorov-Smirnov test [START_REF] Young | Proof without prejudice: use of the Kolmogorov-Smirnov test for the analysis of histograms from flow systems and other sources[END_REF] was applied to each pair of years considered in order to highlight the years exhibiting significantly different distributions. Then a more qualitative analysis was conducted to compare the relative earliness and skewness of the cumulative annual sowing dates (see Supplementary material S3).

Regarding the second point, our analyses were based on daily measurements obtained by the weather station at Fourques from 2002 to 2012. The inter-annual variability in the estimated sowing dates was analyzed first in relation to rainfall patterns, then in relation to temperatures recorded during the preceding crop's growth season. For this latter analysis, the preceding crop was identified as wheat (adopting the algorithm described here), or rice, which is the only summer crop to precede wheat in the study area. The rice fields were identified by applying the PhenoRice algorithm on the MODIS time series (Boschetti et al., 2015a;[START_REF] Manfron | Testing automatic procedures to map rice area and detect phenological crop information exploiting time series analysis of remote sensed MODIS data[END_REF].

Finally, a significant relationship was sought between variations in the weather conditions during the ricegrowing season and variations in the winter wheat sowing dates. This was done by characterizing the temperatures recorded during rice growth by cumulating the mean daily temperatures (from May, 1 st to September, 30 th ) for each year, then fitting linear regression models between the cumulative temperatures and the different deciles of the distributions of the MODIS-estimated sowing dates for the year concerned.

The significance of the explanatory variable (the cumulative mean temperature) was assessed with Student's statistical tests. Regression analyses were run separately for each decile, and for areas characterized by different preceding crops, i.e. winter wheat preceded by winter wheat (Wheat->Wheat), or winter wheat preceded by rice (Rice->Wheat).

Results

Accuracy of winter wheat field identification and sowing dates

The outcome of our validation process concerning the detection of areas where winter wheat was grown is 76 (2009). Concerning the winter wheat class, the commission error ranged between 3% and 19% (mean 11%), and the omission error between 28% and 56% (mean 44%); this means that, on average, the method can identify 56% of winter wheat fields, with quite low misclassification rates. Figure 4 shows a comparison between the MODIS-derived sowing dates and information provided by farmers: 66% of the estimates fell within ± 16 days of the true sowing dates, and 46% within ± 8 days (Figure 4A and4B). The errors were centered on zero (Figure 4B and4C). Figure 4C shows the distribution of the residual errors as a function of the true sowing date. Predictions were more likely to be too early for late sowing dates (e.g. from November 25 th to December 19 th ), and too late for early sowings (from October 16 th to 24 th ). Overall, the method estimated the sowing dates with a mean absolute error of 16.5 days, and a root mean square error of 22.05 days.

shown in Table 2,

Analysis of inter-annual variability

The analysis conducted on the years 2002 to 2012 showed that, on average, winter wheat was sown on November 1 st , with 50% of the estimated sowing dates falling within a 25-day interval between October 16 th and November 9 th (see Supplementary material S2). Year-to-year comparisons drawn with the Kolmogorov-Smirnov test revealed significantly different distributions of winter wheat estimated sowing date (at p = 0.05) for almost all years, the exceptions (distributions were not significantly different) being the pairs 2002-2009, 2002-2011[START_REF] Brivio | Principi e metodi di telerilevamento[END_REF]-2009, and 2009-2012 (Figure 5 (Figure 5, and Table S3.2 in Supplementary material S3). Qualitative analyses using logit curves showed that sowing was done earlier in 2003, 2005 and 2008, and later in 2007 and 2012, and also that the sowing period was more concentrated in the middle of the possible time window (late September to early December) in 2004, 2008 and 2011 (See Table S3. .g. 2004, 2006, and 2012), most of the rain fell before the main sowing periods, while in others (e.g. 2007, 2008 and 2011), much of the rainfall occurred after sowing, and in 2003 the rainfall was evenly distributed before and after the estimated sowing period. There was no evidence of any correlation between the estimated sowing dates and the inter-annual variability in rainfall events.

Relationship between inter-annual sowing date variability and external factors

Influence of rainfall distribution

Influence of preceding crop and temperature

Figure 6 shows the average distribution of the estimated sowing dates grouped by the type of crop preceding the winter wheat. Clearly, the Wheat->Wheat areas tended to be sown earlier than the Rice-> Wheat areas (Figure 6). Sowing between September 30 th and October 16 th was done less frequently (-16% the sum of sowings in the period) in Rice->Wheat fields than in Wheat->Wheat fields, whereas late sowing (e.g. November 9 th to 25 th ) was more common in Rice->Wheat areas (+16% the sum of sowings in the period).

A significant relationship emerged between the first deciles of the sowing date distributions in Rice-> Wheat fields and the temperature during the preceding rice growing season. In particular, the earliest 40% of the wheat sowing dates in Rice->Wheat fields showed a significant negative correlation with the cumulative temperatures during the rice growing season (p<0.01; R 2 ranging between 0.6 and 0.9 depending on the decile -Table 3). In other words, wheat was sown earlier in Rice->Wheat fields in warmer years, when the rice was harvested earlier. The effect of temperature on wheat sowing dates in Rice-> Wheat fields was no longer significant, however, beyond the 5 th decile (p>0.05) of the distribution.

In Wheat->Wheat areas, the effect of summer temperatures on winter wheat sowing dates was never significant (p>0.05; Table 3). In this case, farmers were in a position to schedule the sowing of winter wheat in the most suitable period.

Discussion

The method's performance

The method developed for this study was able to identify winter wheat fields with overall omission and commission errors of 44% and 11%, respectively. This means that, although not all the winter wheat fields were identified, the areas detected by our method have a very high likelihood of being accurate. Since the aim of the method is to automatically identify a reliable and representative sample of winter wheat fields, the automatic detection of 56% of the target area with a low probability of commission errors is more than adequate for this purpose.

Our validation of the estimated sowing dates showed that the method detected the sowing date variability with a margin of error of 8 days in 46% of cases, and 16 days in 66% of cases. Bearing in mind that sowing dates in the Camargue can span a period of more than ten weeks, our method provided useful information for the purpose of analyzing inter-annual variability in crop management practices over the course of a decade.

Sowing date variability and consequences for crop modeling

Our analysis identified a marked inter-annual variability in winter wheat sowing dates, with two main sowing seasons, one in early October (as in 2003, 2005 and 2008), the other in late October (as in 2007 and 2012). Rainfall patterns did not appear to influence the sowing dates in the Camargue, though some authors have reported such a relationship in other Mediterranean areas (e.g., [START_REF] Bassu | Optimising sowing date of durum wheat in a variable Mediterranean environment[END_REF].

The winter wheat sowing dates did instead correlate with the preceding crop (rice vs. wheat) and the summer temperatures. With climate change, weather-related constraints may play a fundamental role in changing agricultural practices, strongly influencing the growth cycles of crops such as rice [START_REF] Bregaglio | Identifying trends and asociated uncertainties in potential rice production under climate change in Mediterranean areas[END_REF]. Other sources of variability may affect winter wheat sowing dates including type of soil, differences in soil moisture and field accessibility, as well as organizational constraints at farm level (e.g. [START_REF] Mailly | Un modèle de décision d'assolement en riziculture conventionnelle et biologique pour prédire les usages des sols sous différents scénarios: cas de la Camargue (Sud de la France)[END_REF][START_REF] Schaller | Combining farmers' decision rules and landscape stochastic regularities for landscape modelling[END_REF], grain prices and CAP subsidies.

Our method's ability to provide information on winter wheat sowing dates (hitherto unavailable for the study area considered here) and on their inter-annual variability marks a step forward in the characterization of spatially explicit and yearly variable crop calendars. Such information may be crucial for yield forecasting, for instance, because durum wheat sowing dates strongly influence the final grain yield [START_REF] Bassu | Optimising sowing date of durum wheat in a variable Mediterranean environment[END_REF][START_REF] Ferrise | Sowing date and nitrogen fertilisation effects on dry matter and nitrogen dynamics for durum wheat: An experimental and simulation study[END_REF]. It can have positive implications for the performance of crop models too, because sowing date parameters strongly influence the simulation of leaf appearance dynamics [START_REF] Bassu | Optimising sowing date of durum wheat in a variable Mediterranean environment[END_REF][START_REF] Ferrise | Sowing date and nitrogen fertilisation effects on dry matter and nitrogen dynamics for durum wheat: An experimental and simulation study[END_REF] and nitrogen uptake [START_REF] Ehdaie | Sowing date and nitrogen rate effects on dry matter and nitrogen partitioning in bread and durum wheat[END_REF].

Method's exportability, limits and future prospects

The proposed method relies mainly on expert-derived rules and thresholds. While, in the specific context considered here, the choices made were the only feasible options to ensure the robustness of our data, they might limit the exportability of our method to other areas. The temporal signature of winter wheat on MODIS images may differ, for example, in the presence of different crop species (e.g. Triticum durum vs

Triticum aestivum), with different dormancy and vegetative periods, or due to snow cover in winter. It is worth noting, however, that the method could easily be replicated for other areas because: (i) MODIS data are available worldwide; and (ii) local experts can adapt the agronomic-based rules to their own settings. A source of uncertainty in the results obtained with our model stems from the rather coarse resolution of MODIS data in relation to the dimensions of the fields (the MODIS pixel area of 5.37 ha is almost twice the size of a typical wheat field in the Camargue, which is 2.73-2.83 ha), and to the fragmentation of the landscape (the fields have irregular shapes, and wheat fields are often interspersed with fields with other crops). This leads to the well-known "low resolution bias" and "mixed-pixel" effects, which introduce noise in time series of moderate-resolution satellite images [START_REF] Boschetti | Analysis of the conflict between omission and commission in low spatial resolution dichotomic thematic products: The Pareto Boundary[END_REF][START_REF] Mingwei | Crop discrimination in Northern China with double cropping systems using Fourier analysis of time-series MODIS data[END_REF][START_REF] Pan | Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data[END_REF][START_REF] Sakamoto | A crop phenology detection method using time-series MODIS data[END_REF], making it more difficult to identify reliable time series for the crops of interest. Another source of uncertainty derives from the effect of the bi-directional reflectance distribution function (BRDF), which can be pronounced for wide-angle sensors like MODIS, further adding to the noise in time-series data [START_REF] Hansen | A review of large area monitoring of land cover change using Landsat data[END_REF][START_REF] Huete | Overview of the radiometric and biophysical performance of the MODIS vegetation indices[END_REF][START_REF] Sakamoto | A Two-Step Filtering approach for detecting maize and soybean phenology with time-series MODIS data[END_REF].

For the purpose of analyzing long-term past dynamics of crop practices, the spatial and temporal resolution of archival data could be enhanced by combining heterogeneous satellite data (e.g. from MODIS and Landsat) using "data fusion" techniques to simulate time series with a higher spatio-temporal resolution [START_REF] Bisquert | A Simple Fusion Method for Image Time Series Based on the Estimation of Image Temporal Validity[END_REF][START_REF] Dong | Estimating winter wheat biomass by assimilating leaf area index derived from fusion of Landsat-8 and MODIS data[END_REF]. An interesting opportunity for present and future analyses could also come from examining improved multi-sensor time series obtained from datasets with a higher spatial resolution, by making a combined use of Landsat 8 and Sentinel-2 data, for instance.

Conclusions

In this study, expert agronomic knowledge was combined with an analysis of the temporal signatures of crops on MODIS satellite images to design a method for estimating the inter-annual variability of winter wheat sowing dates. To our knowledge, no studies have yet focused on deriving winter wheat sowing dates from the analysis of moderate-resolution satellite data.

MODIS spatial and temporal resolutions proved adequate for depicting inter-annual winter wheat sowing date variability, providing information previously unavailable for the study area considered, thus representing a step forward vis-à-vis currently available static crop calendars. Our findings revealed a marked variability in sowing date distributions from one year to the other, with differences in the timing and range of the sowing period. They also showed that sowing dates in fields where wheat was sown after rice correlated strongly with summer temperatures (i.e., winter wheat sowing could be delayed by late rice harvests). Identifying the seasonal distribution of sowing dates can also be useful for assessing regionalscale winter wheat yield variability by means of crop modelling. While our proposed method is suitable for depicting variability over large areas, satellite time series with higher spatial and temporal resolutions would be needed for more detailed analyses on crop management practices. A better spatial resolution would reduce the noise in time series due to mixed pixel contamination, while improving the temporal resolution, would reveal further characteristics of winter wheat time signatures that could be associated with the crop's phenological features.

authors are very grateful to Francesco Nutini, Daniela Stroppiana, Paolo Villa and Simone Bregaglio for their support in the preparation of this work.

Figure 1 -

 1 Figure 1 -Flow chart representing the key steps in the method described and the datasets involved in each phase.

  -A and Figure 2-B). The two main crops are: rice (Oryza sativa L.), covering 36% (11,684 ha) and 56% (16,037 ha) of the area under cultivation in 2006 and 2011, respectively; and winter wheat (Triticum durum L.), occupying 22% (7,350 ha), and 10% (2,927 ha) of the cultivated area in 2006 and 2011, respectively. The average size of the winter wheat fields was 2.73 ha (min = 0.21 ha; max = 30.26 ha; median = 2.26 ha) in 2006, and 2.83 ha (min = 0.38 ha; max = 23.69 ha; median = 2.31 ha) in 2011.

Figure 2 -

 2 Figure 2 -Camargue, France (Europe): cultivated areas according to the Parc Naturel Régional de Camargue (PNRC) land use maps. Main crops in 2006 (A), and 2011 (B): rice (green), winter wheat (blue), and other crops (orange). Bodies of water, nonagricultural areas, and areas of no interest are colored black, white and grey, respectively.

  them enables VI time series to be constructed with a nominal eight-day temporal frequency. In particular, a time series of 484 MODIS MOD/MYD 13Q1 images covering the period 2002-2013, from DOY 177 in 2002 (end of June) to DOY 361 in 2013 (end of December) was downloaded for pre-processing from the NASA Land Processes Distributed Active Archive Center (LP DAAC) at the USGS Earth Resources Observation and

Figure 3 -

 3 Figure 3 -Typical EVI temporal signature for winter wheat in the Camargue. The graph shows median values (red line) and the interquartile range (red shadow) for 479 winter wheat time series in 2006 and 2011, considering MODIS pixels with a fractional cover exceeding 80%.

  results obtained with the algorithm (years 2002-2012) were validated by comparing: (1) the winter wheat raster maps obtained for 2003, 2006, 2009 and 2012 with Landsat data; and (2) the sowing dates estimated in 2010, 2011, 2012 with ground reference data.

  reliability of the sowing dates estimated from the MODIS data was assessed by comparison with field information on agronomic practices. The available field survey data included field-scale winter wheat sowing dates for three different farms, which accounted for a total of 1,334 ha, corresponding to 4.6% of the area under cultivation in 2011 according to the PNRC. The information concerned 83 winter wheat fields (261 ha) in 2011, 94 (343 ha) in 2012, and 111 (363 ha) in 2013. The information associated with the different fields was first aggregated at the spatial resolution of MODIS by means of a GIS spatial averaging procedure. The averaged sowing dates were then associated with the time of the nearest MODIS 8-day temporal composite, and a comparison was drawn between the averaged field data and the MODIS-estimated sowing dates for the years 2010, 2011 and 2012. A total of 146 pixels were involved in this comparison: 52 for 2010, 30 for 2011 and 64 for 2012.

  in the form of error matrices for the years 2003, 2006, 2009 and 2012, and for all four years together. The results were good, with annual OA values ranging between 84% (2006) and 93% (2009), and Cohen's K coefficient in the range of 0.54 (2006) to 0.

Figure 4 -

 4 Figure 4 -Comparison between true winter wheat sowing dates and MODIS-based estimates for three farms in the years 2011, 2012 and 2013 (n=146). (A) Scatter plot between observed and estimated dates, with 8-day and 16-day confidence intervals in dashed red lines. The grey scale indicates the number of pixels in each binning square spanning 8 days. (B) Distribution of residual errors. (C) Distribution of residual errors as a function of true sowing date.

  Figure S3.1).

Figure 5 -

 5 Figure 5 -Frequency histograms of estimated winter wheat sowing dates in the Camargue from 2002 to 2012. The lower-right histogram (in dark grey) shows the mean frequencies for the years considered. Black vertical lines superimposed on each histogram indicate the median (dashed line) and the first and third quartiles (continuous lines) of the distributions. Rainfall (mm) recorded at the Fourques weather station is also shown for each year, and the mean for the 11-year period.

Figure 5

 5 Figure 5 also presents the estimated winter wheat sowing dates in relation to the distribution of rainfall events in the Camargue for the years 2002-2012. In some seasons(e.g. 2004, 2006, and 2012), most of the

Figure 6 -

 6 Figure 6 -Frequency histograms of winter wheat sowing dates for Wheat->Wheat areas (light grey histograms, n = 1771) and Rice->Wheat areas (dark grey histograms, n = 1286). Data refer to the years 2002-2012.

Table 1 -List of Landsat-7 ETM+ images used to build the reference dataset for validating the MODIS-derived winter wheat growing area maps. SLC refers to the Scan Line Corrector device of Landsat 7, which failed on May 31, 2003.

 1

	ACQUISITION DATE	DOY	SLC status CLOUD COVER %
	April-22-2003	112	On	2
	August-12-2003	224	Off	1
	April-14-2006	104	Off	0
	July-19-2006	200	Off	15
	April-22-2009	112	Off	2
	August-28-2009	240	Off	2
	April-14-2012	105	Off	18
	July-20-2012	201	Off	23

Table 2 -Accuracy matrix obtained by comparing the reference dataset (based on a visual interpretation of Landsat images) with the winter wheat maps based on MODIS EVI time series analyses. CE: commission error, OE: omission error, OA: overall accuracy, K: coefficient of agreement.

 2

				REFERENCE DATASET		ACCURACY
				winter wheat	other	CE	OE	OA	K
		2003	winter wheat other	44 23	10 223	19% 9%	34% 4%	89%	0.66
	MODIS ESTIMATION	2006 2009 2012	winter wheat other winter wheat other winter wheat other	40 45 39 15 28 35	2 213 5 241 1 236	5% 17% 11% 6% 3% 13%	53% 1% 28% 2% 56% 0%	84% 93% 88%	0.54 0.76 0.55
		All	winter wheat other	151 118	18 913	11% 11%	44% 2%	89%	0.62

Table 3 -Summary statistics of the regression models fitting Rice->Wheat and Wheat ->Wheat sowing dates with the cumulative mean temperatures during the rice-growing season. R 2 : coefficient of determination. D1-D5: 1 st to 5 th decile of the winter wheat sowing date distribution.

 3

	Preceding crop Regression model	D1	D2	D3	D4	D5
	Rice	Slope	-0.068	-0.057	-0.074	-0.054	-0.046
	Wheat	Slope	-0.005	0.007	-0.004	-0.015	0.026
	Rice	p-value	<0.001	0.001	<0.001	0.005	0.076
	Wheat	p-value	0.759	0.681	0.851	0.541	0.172
	Rice	R 2	0.779	0.722	0.91	0.606	0.259
	Wheat	R 2	0	0	0	0	0.122

Acknowledgements

This research was conducted within the Scenarice project (n°1201-008) funded by the Agropolis and Cariplo foundations. It was supported by the French National Research Agency under the program "Investments for the future" ANR-10-LABEX-0001-01. Field data, land cover data and agronomical expertise were kindly shared by Vincent Couderc, Jean-Claude Mouret and Jean-Marc Barbier from INRA UMR Innovation. The

Supplementary material

S1: Preliminary analysis of EVI temporal signature of winter wheat

Figure S1.1 shows the smoothed EVI temporal signature for a sample of pure MODIS wheat pixels (38 pixels with wheat fc > 75%) for which farmers provided sowing dates in 2011, 2012 or 2013. Data are given as box-plots using the MODIS eight-day temporal resolution. EVI time series were aligned with the sowing dates to produce day after sowing (DAS) time series profiles.

Although the average VI signal starts to increase slowly after 20-40 DAS, the green-up signal coinciding with the initial plant leaf development and tillering phases prior to dormancy cannot be identified. A significant VI increase is only evident after dormancy, corresponding to the stem elongation phase. From this moment, up to senescence, the VI profile shows a continuous monotonic increase that plateau at the end of the flowering period. Given this behavior, calculating the sowing dates by backtracking a given number of days from the observable green-up does not seem feasible for the case study considered here. Figure S1.2 shows the frequency histograms for automatic green-up detection on the time series for 479 pure wheat pixels (i.e., where the presence of the crop in the pixels was >80%, judging from PNRC 2006 and 2011 land cover maps). According to [START_REF] Lobell | Satellite detection of earlier wheat sowing in India and implications for yield trends[END_REF], the green-up is defined as the moment when the curve fitted to the VIs reaches 10% of its maximum amplitude for a given year. This preliminary analysis shows that most of the green-ups (>60%) apparent on time series analysis coincided with the postdormancy stem elongation phase, while there was a low prevalence of green-ups in the pre-dormancy tillering phase. In these conditions, it would clearly be inappropriate to adopt a set number of days to establish the sowing date.

Supplementary Material

Click here to download Supplementary Material: Supplementary-Material-R2_JAG-D-16-00323.docx

S2: EVI temporal signature of crops in the Camargue

All the main crops in the Camargue showed different profiles in representative time series (Figure S2.1), and the rules applied in the algorithm prevented winter wheat from being confused with other crops. The temporal signature of winter wheat (panel A) shows an initial sowing phase (from October onwards) characterized by low EVI values (around 0.2) that slowly decrease up until the end of vernalization (in January), followed by a rapid increase in the signal marking rapid growth (green-up) up until crop heading

S3: Statistics for estimating sowing dates and results of year-by-year comparisons

A qualitative analysis was run to compare yearly sowing date distributions paying particular attention to: i) the (relatively) early sowing dates; and ii) the (relative) skewness to a given period of the sowing date distribution. This agronomic expert-based analysis was based on logistic curves fitted to the cumulative yearly distributions of estimated winter wheat sowing dates, expressed as frequencies.

Table S3.1 shows the statistics for the sowing date estimated in each year from 2002 to 2012, together with the overall average for the period analyzed.Table S3.Table S3.2 shows the results of the Kolmogorov-Smirnov coupled test. Only a few years were not significantly different from the others.

Table S3.3 shows the results of the logistic curve fitting procedure on the cumulative frequencies of the yearly distributions of the estimated winter wheat sowing dates. The logistic curves feature three numerical parameters, namely: the asymptote (here forced to 1 because frequencies were used); the x value (date) at the inflection point of the curve (here representing how early sowing took place); and the scale parameter on the x-axis, representing the slope (i.e., if sowing was mainly concentrated on certain consecutive dates).

Figure S3.1 shows the fitted logistic curves. The graphs highlight the years in which the sowing season was early or late, and also the range of the sowing dates.

S4: Threshold on EVI local MAX point for identifying winter wheat heading

During the mapping of winter wheat fields, an EVI threshold of 0.42 was chosen to identify EVI MAX points corresponding to a strong presence of biomass. This value was obtained by means of a recursive partitioning technique implemented in R ("rpart" package) to select the best-performing threshold for separating a total of 1000 samples. This partitioning is based on minimal error fitting of the calibration dataset and identifies the EVI value that best discriminates between pixels of winter wheat and other crops.

Samples were extracted considering the higher EVI value in the second quarter of the year (April 1 st to June 30 th ) in highly-reliable time series. In particular, 500 samples were randomly chosen for winter wheat and another 500 by means of a stratification procedure according to how well they represented other crops in the study area.

Figure S4.1 shows the box plot of the two sets of samples involved in this analysis. Partitioning led to the value of 0.42 being identified as the most appropriate for distinguishing between the samples. Of the 460 samples above said cutoff (46% of the total), 62% belonged to the "winter wheat" group, and 38% to the "other crops" group; and of the other 540 samples below the cutoff (54% of the total), 60% were in the "other crops" group, and 40% in the "winter wheat" group.