# Country-scale crop-specific phenology from disaggregated PROBA-V

Henry Rivas, Nicolas Delbart, Fabienne Maignan, Emmanuelle Vaudour, Catherine Ottlé

► To cite this version:

Henry Rivas, Nicolas Delbart, Fabienne Maignan, Emmanuelle Vaudour, Catherine Ottlé. Country-scale crop-specific phenology from disaggregated PROBA-V. Remote Sensing, 2024, Special Issue Remote Sensing of Land Surface Phenology II, 16 (23), pp.4521. 10.3390/rs16234521. hal-04814690

HAL Id: hal-04814690

https://hal.science/hal-04814690v1

Submitted on 2 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Country-Scale Crop-Specific Phenology from Disaggregated PROBA-V

Henry Rivas <fnref n="1" />,<fnref n="2" />,* , Nicolas Delbart <fnref n="1" /> , Fabienne Maignan <fnref n="3" /> , Emmanuelle Vaudour <fnref n="4" /> and Catherine Ottlé <fnref n="3" />

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Rivas, H.; Delbart, N.; Maignan, F.; Vaudour, E.; Ottlé, C. Country-Scale Crop-Specific Phenology from Disaggregated PROBA-V. *Remote Sens.* **2024**, *16*, 4521.
https://doi.org/10.3390/rs16234521

Academic Editor: Jianxi Huang

Received: 29 October 2024

Revised: 27 November 2024

Accepted: 30 November 2024

Published: 2 December 2024

<figure><img src="image_4.png" /></figure>

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

**Abstract:** Large-scale crop phenology monitoring is essential for agro-ecosystem policy. Remote sensing helps track crop development but requires high-temporal and spatial resolutions. While datasets with both attributes are now available, their large-scale applications require significant resources. Medium-resolution data offer daily observations but lack detail for smaller plots. This study generated crop-specific phenomaps for mainland France (2016–2020) using PROBA-V data. A spatial disaggregation method reconstructed NDVI time series for individual crops within mixed pixels. Then, phenometrics were extracted from disaggregated PROBA-V and Sentinel-2 separately and compared to observed phenological stages. Results showed that PROBA-V-based phenomaps closely matched observations at regional level, with moderate accuracy at municipal level. PROBA-V demonstrated a higher detection rate than Sentinel-2, especially in cloudy periods, and successfully generated phenomaps before Sentinel-2B’s launch. The study highlights PROBA-V’s potential for operational crop monitoring, i.e., wheat heading and oilseed rape flowering, with performance comparable to Sentinel-2. PROBA-V outputs complement Sentinel-2: phenometrics cannot be generated at plot level but are efficiently produced at regional or national scales to study phenological gradients more easily than with Sentinel-2 and with similar accuracy. This approach could be extended to MODIS or SPOT-VGT, to generate historical phenological data, providing that a crop map is available.

**Keywords:** crops phenometrics; mixels; spatial disaggregation; winter wheat heading; oilseed rape flowering; phenomaps

# 1. Introduction

Phenology is the study of the timing of recurrent seasonal events in the life of plants and animals, including its observation and description, the analysis of biotic and abiotic causes, and its role in the interplay between species [1]. In botanics, it concerns the timing of key events such as flowering, leaf bud-break or senescence, the role of explanatory factors such as air temperature [2], soil moisture, photo-period, etc. Phenology metrics are essential parameters for assessing variations in the functioning of ecosystems in the context of climate change [3–10]. In particular, the advance of spring events are considered among the clearest effects of climate change [11,12]. These changes themselves have consequences on carbon [9,13,14] and water exchanges balance [15], biodiversity [16,17], and most eco-logical aspects as detailed by [18]. For agriculture, phenological information plays a role in (i) scheduling agricultural activities such as irrigation, fertilization, agrochemical applica-tion and harvesting, (ii) choosing suitable species and varieties based on local pedoclimatic

<footnote n="1">¹ Laboratoire Interdisciplinaire des Energies de Demain, Université Paris Cité, 75013 Paris, France; nicolas.delbart@univ-paris-diderot.fr</footnote>

<footnote n="2">² CESBIO, Université de Toulouse, CNES/CNRS/INRAE/IRD/UT3-Paul Sabatier, 31401 Toulouse, France</footnote>

<footnote n="3">³ Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, 91190 Gif-sur-Yvette, France; fabienne.maignan@lsce.ipsl.fr (F.M.); catherine.ottle@lsce.ipsl.fr (C.O.)</footnote>

<footnote n="4">⁴ INRAE, AgroParisTech, UMR EcoSys, Université Paris-Saclay, 91120 Palaiseau, France; emmanuelle.vaudour@inrae.fr</footnote>

<footnote>* Correspondence: hrivasullon@gmail.com</footnote>
conditions, and (iii) planning crop rotations [19]. It also contributes to the prediction of
yield by indicating the stages to be evaluated as well as the optimal moment to measure
the parameters (e.g., biomass, number of pods, weight of 100 seeds, etc.) [20], and therefore
supports food security [21,22]. For crops, the main stages are emergence, leaf development,
tillering, bolting, inflorescence/heading, flowering, fruit development, maturation, and
senescence [23].

Hence, it is essential to acquire phenology observations for both natural and agricultural vegetation. These are obtained through networks of naturalist field observations [24–28], by analyzing daily radiometric measurements from field sensors [29–33] or by remote sensing [5,34–39]. In remote sensing, the concept of observing and interpreting phenological information is different from observing specific phenological stages in the field. Land Surface Phenology (LSP) refers to seasonal variations in surface reflectance that correspond to the overall development of vegetation, without being linked to any specific phenological phase [40,41]. These seasonal variations are therefore the result of the accumulation of specific phenological events that are not explicitly detectable by the sensors. Phenometrics usually detected by satellites are the start (SOS), peak (POS), end (EOS), and duration of the growing season, which can be partially linked to biological events [32].

Over the past four decades, daily observations from medium-spatial resolution (MSR) sensors like AVHRR [42], SPOT-VGT [43], MODIS [44], or PROBA-V [45], have significantly contributed in assessing changes in natural vegetation [46]. These observations have allowed to quantify the evolution of the ecosystem's seasonality [7,35,47,48], and provided valuable insights into the consequences of these changes on biogeochemical processes [9,49,50]. In the literature, most remote sensing-based phenological studies have focused on monitoring entire vegetated surfaces [51–53] or large thematically homogeneous areas [54–56]. This focus is mainly driven by the limitations imposed by the spatial resolution of available data and the undesirable effects of mixed pixels (mixels) [57].

The remote sensing-based monitoring of crop phenology, however, presents additional specific challenges. These arise from the limitations in current sensor capabilities, including the needs for high-frequency observations during critical phenological transitions and for adequate spatial resolution to differentiate individual agricultural plots. Hence, satellite-based monitoring of phenology is more complex for crops than for forests or natural grasslands, due to the fragmented nature of agricultural landscapes resulting in mixels, and the diversity of agricultural practices leading to complex and multi-modal signals. Despite this, MSR images (250 m–4 km) allowed to characterize certain agricultural systems by their phenology [58,59].

Today, the recent availability of high-spatio-temporal resolution data (HSR), provided by Sentinel-2 (10 m and 5 days of revisit) and Sentinel-1 (10 m and 6 days of revisit), or harmonized Landsat-8/Sentinel-2 [60], allows the monitoring of phenology at a finer scale [61–64]. High-spatio-temporal resolution time series were used to obtain LSP indicators for various crops in many regions of the world. For example, in Russia, LSP indicators agreed with the timing of tillering, rising, or maturity of crops [65]. In the eastern United States, these data were able to extract crop emergence date as well as cover crop termination date [66,67]. High-spatio-temporal resolution phenology monitoring has allowed to improve crop types mapping at a 30 m resolution by a synergistic use of Sentinel-2 and Landsat 7 and 8 [68], to map rice-cultivated surfaces at 10 m [69], or soya-cultivated surfaces with MODIS for large homogeneous fields or Sentinel-2 for smaller parcels [70]. Similarly, it was shown that considering phenology is necessary when estimating yield from remote sensing time series for sunflower [71], manioc [72], or rice [73]. LSP indicators also allow mapping crop successions [74] that may impact water quality [75]. Additionally, commercial PlanetScope (3 m and 1 day revisit) constellation data have already demonstrated their high potential in phenology studies [76], mainly in areas with short-cycle vegetation and heterogeneous land cover [77]. Similarly, Venus (5 m and 2–3 days revisit) data have also been able to accurately detect crop-specific phenological metrics and have exhibited near-real-time detection capabilities [66,67]. These studies showed increased
performances compared to Sentinel-2 data, largely due to the higher temporal revisit (<3 days). This very high spatio-temporal resolution phenological information has helped to understand how vegetation responds to local weather conditions, whether at the individual level with PlanetScope (3 m) data or at the landscape level with Sentinel-2 (10 m) and MODIS (250 m) data [77]. Refs. [66,67] have demonstrated that phenological information can be extracted in near-real-time from Venus data (5 m) and likely from Sentinel-2 (10 m) and harmonized Landsat/Sentinel-2 (30 m) data, which could support agricultural scheduling or public policy monitoring.

Phenological indicators can thus be extracted at high- and even very high spatio-temporal resolution from optical remote sensing time series. Radar data have been used for phenology monitoring because of their insensitivity to atmospheric conditions and their current cost-free availability. The synergistic use of optical Sentinel-2 data and radar Sentinel-1 data has further increased the possibility of mapping phenological events [78–80]. Altogether, the above-mentioned optical high-spatio-temporal data and Sentinel-1 data offer an unprecedented opportunity in the study of current and future phenology, both for natural and agricultural vegetation.

However, these datasets have only become usable for phenological purposes in the last decade due to improved revisit capabilities. For instance, Sentinel-2 data have been exploitable for phenological monitoring only since the launch of Sentinel-2B in March 2017, while Sentinel-1 data's usefulness improved only after the launch of Sentinel-1B in April 2016. Harmonized Landsat/Sentinel-2 data, which show promise, have been available since 2015. Venus, designed for scientific purposes, offers data from August 2017 for specific sites only, while PlanetScope, among sensors with extensive global observation capabilities, provides data from June 2016, and so forth. With those data, a long-term phenomaps archive cannot be built to assess the climate change's impact on agro-ecosystems. Moreover, their use for large-scale applications requires significant computational resources, which could be an operational limitation, especially when it relies on very high-spatio-temporal resolution data such as PlanetScope data, which indeed could pose financial challenges for large-scale applications in support of public policy initiatives, this sensor being commercial. Thus, such data cannot easily support large-scale phenological mapping that would facilitate the implementation of new adaptation strategies at the agro-ecosystem level. Finally, recent studies using high-spatio-temporal resolution data (3–10 m and 1–5 days revisit, i.e., PlanetScope, Venus, and Sentinel-2) also highlighted that despite acquiring a dense dataset, the number of cloud-free observations might not be enough to ensure accurate detection of phenological stages, as a function of location and temporal period in the growing season [66,67,76,77].

In this context, MSR data have great potential for monitoring crop-specific stages because of their large-scale operational applicability and their long period availability (about 40 years). MSR also offers a daily temporal resolution that increases the probability of cloud-free observations and is suitable for studying plant dynamics during the main development periods. However, in fragmented agricultural landscapes, the problem of mixels must be correctly addressed. In this study, we explore the potential of MSR optical data to provide crop phenometric maps to study the phenological gradients at a national scale. We evaluate this potential comparatively to the performances obtained with Sentinel-2 data. For this purpose, the MSR data need to be disaggregated to each specific crops. The approach consists in isolating the contribution of each crop type in the mixed signal to extract its average phenometrics within a group of mixels, allowing to study the phenology across scales from the municipality to the country, and allowing the comparison of the phenology of the various crops. The disaggregation does not aim to provide the phenometrics for each stand individually, but the average phenology of a crop within an area corresponding to several pixels. As a feasibility test, we generated country-scale crop-specific phenomaps over mainland France (except Corsica) during 2016–2020, from medium-spatial resolution PROBA-V (300 m) data. Additionally, this feasibility test aims to determine whether this PROBA-V's ability remains prior to the launch of Sentinel-2B, as
this would extend in the past the possibility to derive crop phenology. Nevertheless, the description of phenology at plot level is outside the scope of this study due to the nature of PROBA-V's disaggregated data, which synthesizes information from all plots within a neighborhood. Despite synthesizing spatial variations within neighborhoods, disaggregated data could effectively describe large-scale spatial patterns within municipalities or regions, and allow analyzing spatial gradients in phenology within the country limits.

This work is based on the spatial disaggregation method previously developed [81], which provides NDVI time series for individual crops within mixels, provided that the surface fraction of each crop is known. This approach is in continuity with other disaggregation techniques in the literature [82–87]. Following this pre-processing, we applied usual LSP extraction techniques to resulting crop-specific time series. These techniques involved gap-filling methods [38, 88, 89], followed by the extraction of the date when the adjusted curve reaches a threshold or maximizes a curvature [90]. These LSP extraction techniques were also applied identically to Sentinel-2 data. Subsequently, phenometrics were compared to observed phenological stage dates provided by: 1/the Céré'Obs dataset, which gives a regional synthesis of the phenological progression, arising from the analysis by agricultural advisers in *chambres d'agriculture*, and 2/TEMPO observations, coming from visual inspection on the ground of the timing of phenological stages following the *Biologische Bundesanstalt, Bundessortenamt und CChemische Industrie* (BBCH) protocol [23]. Given their significant socio-economic importance and the high performance of the spatial disaggregation algorithm for these major crops compared to other minority ones like maize or spring barley [81], our study was focused on three crops in France: winter wheat, winter barley, and oilseed rape. These crops represent about 50% of the agricultural area excluding grasslands. Here, we added the winter barley plots to the winter wheat ones due to their similar growth patterns and indistinguishable radiometric responses observed in this work. To enhance clarity, this paper organizes the presented results into only two thematic categories: winter wheat and oilseed rape.

The advantages and limitations of phenometrics extracted from the disaggregated PROBA-V data, and MSR data in general, are discussed in comparison with the results obtained from Sentinel-2 data. The results also involve an analysis of the succession of a phenological event among two crops across the study years and the common cultivation area.

## 2. Materials and Methods

We estimated crop-specific phenometrics from both disaggregated PROBA-V (300 m) and Sentinel-2 (10 m) NDVI time series for the period 2016–2020 (Figure 1). From the disaggregated PROBA-V data, phenometrics were derived at the national level through a pixel-by-pixel strategy, while a parcel-by-parcel strategy was employed for the high-resolution data at seven Sentinel-2 tiles (100 km × 100 km) level, here called inter-comparison sites. Then, a comparison was made between the phenometrics derived from each sensor within these inter-comparison sites. Finally, the performance of each sensor was assessed by comparing the estimated phenometrics with ground data from the French phenological observation networks, TEMPO and Céré'Obs. In the following sections, we will explain the different methodological steps carried out.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** General flowchart of this study. Crop-specific phenological mapping procedure using (A) PROBA-V (300 m) and (B) Sentinel-2 (10 m) data. PROBA-V-based phenometrics were extracted at pixel level, while those from Sentinel-2 were extracted at plot level. Both are inter-compared and compared with TEMPO data at the municipal level. PROBA-V-based phenometrics are then compared with Céré'Obs data at the regional level. For phenometrics extraction, thresholds were calibrated in Block A using disaggregated PROBA-V NDVI time series and TEMPO data, and were then applied identically to both disaggregated PROBA-V and Sentinel-2 NDVI time series.</figcaption></figure>

## 2.1. *Study Area*

Our study area covers about 25% (75,000 km²) of the agricultural area declared in mainland France's (except Corsica) Land Parcel Identification System (LPIS, [91]) in 2016–2020 (Figure 2). According to the Köppen–Geiger classification [92], the climate is mainly oceanic, with warm summers throughout mainland France and a Mediterranean climate in its south. The main extensive crops are wheat, barley, maize, sunflower, oilseed rape, sugar beets, grapes, and potatoes (https://agreste.agriculture.gouv.fr, accessed on 1 October 2024). Crop production typically follows a conventional approach, beginning with deep ploughing in November–December, followed by chiseling in March, and subsequent seedbed preparation for spring cereals. Winter crops are sown in October–November, and spring crops in March–April. Annual rainfall is around 800–1000 mm, with a contrast between the western (>1000 mm) and the southeastern regions (600–800 mm). The average annual temperature is about 11–13 °C, with an average value of 20–25 °C in summer and of 5–10 °C in winter (https://meteofrance.com, accessed on 1 October 2024).
<figure><img src="image_6.png" /></figure>

**Figure 2.** Study area location. The validation of the phenomaps was carried out at regional and municipal scales. (A) The regional comparison between the phenometrics (from PROBA-V) and the observed (from Céré'Obs) median phenological dates was made across administrative regions, highlighted in both yellow and orange (n = 14), while the comparison between phenometrics and observed phenological progression (i.e., the percentage of area reaching a given phenological stage, as a function of date) focused on regions highlighted in orange only (n = 3). Regional ground data were not available for the three regions highlighted in gray. (B,C) The comparison at municipal level between the phenometrics (from PROBA-V and Sentinel-2) and the observed (from TEMPO) median phenological stage dates was made across points highlighted in blue (winter wheat) and red (oilseed rape) within inter-comparison sites. Inter-comparison sites were delimited by seven Sentinel-2 tiles: 30UWU, 31UCQ, 31UDP, 31TCN, 31UGQ, 31TFN, and 31TCJ. Additionally, PROBA-V-based phenomaps were also compared to TEMPO data available outside these inter-comparison sites. Finally, the green mask shows the winter wheat and oilseed rape areas declared in 2019 according to LPIS.

Here, we evaluated our results at the regional and municipal administrative scales. Regional evaluation was conducted at the national level and municipal evaluation was conducted only within inter-comparison sites. These sites were chosen to compare the performances of PROBA-V (300 m) and Sentinel-2 (10 m) since high-resolution data analysis demands substantial computing power, posing challenges for nationwide comparisons. Inter-comparison sites are representative of four eco-climatic regions [93] that present important phenological contrasts and host the largest agricultural production basins in mainland France.

## 2.2. Satellite Data

All available scenes from PROBA-V and Sentinel-2 sensors between September 2015 and October 2020 were used. Images were preprocessed by masking non-valid pixels using their quality flags (clouds, shadows, undefined...) for each acquisition date.

### 2.2.1. Disaggregated PROBA-V (300 m)

PROBA-V images were acquired at a 300 m spatial resolution with a daily revisit [94]. Annual NDVI time series, which included 420 scenes, were obtained from the Vito Remote Sensing data center (https://proba-v.vgt.vito.be, accessed on 15 July 2024). Because a 300 m spatial resolution pixel may include several crop types and even other land uses, we generated disaggregated NDVI time series using the spatial disaggregation (SD) algorithm
developed by [81] to derive the NDVI time series for several crop types. SD was based on the knowledge of crop fractions within a mixel. This linear mixing approach assumes that the signal value of a mixel is the weighted average of the contributions from each crop fraction within that mixel. It also assumes that the signal is stationary within a crop class, as described by [95]. Hence, the mathematical expression of the mixing model is as follows:

$$NDVI_{j}(t) = \sum_{i=1}^{n} fc_{i,j} \times NDVI_{i}(t) + \varepsilon_{j}(t) \quad (1)$$

where $NDVI_{j}(t)$ is the NDVI value of the mixel $j$ at time $t$, $fc_{i,j}$ is the fraction of crop $i$ in mixel $j$, $NDVI_{i}(t)$ is the NDVI value of crop $i$ at time $t$, $n$ is the number of crops within mixel $j$, and $\varepsilon_{j}(t)$ is the error in mixel $j$ at time $t$.

To derive the crop-specific NDVI value within the mixel, we proceeded to an inversion of an $m$ equations system, where $m$ is an integer totaling the mixel and its eight neighbors. The theoretical value of $m$ is 9, but may be reduced to exclude unsuitable mixels from the system (corresponding, for example, to mixels containing none of the targeted crops). It should be mentioned that the spatial disaggregation method considers that a crop signal is also stationary among all pixels (i.e., same growing stage and condition) in the neighborhood, as intra-plot or inter-plot phenological differences within a 9-pixel window are outside the scope of our study. All details are provided in [81].

### 2.2.2. Sentinel-2 (10 m)

Sentinel-2 images were acquired at a 10 m spatial resolution with a revisit of 5 days (10 days before 2017). For each tile shown in Figure 2, annual reflectance time series, which included 97 scenes on average, were obtained at level 2 (Bottom-of-Atmosphere). Subsequently, NDVI time series were derived from the red (B4) and NIR (B8) bands as detailed by [96].

Here, data were obtained from two distribution sites. For the 2018–2020 period, data were obtained from the Copernicus Open Access Hub site (https://scihub.copernicus.eu, accessed on 10 January 2024) via the Google Earth Engine (GEE) platform (https://earthengine.google.com, accessed on 10 January 2024). This distribution site uses the Sen2Cor processing chain for obtaining level 2 products [97]. Since these products are unavailable on the GEE platform for the previous years (2016–2017), the remaining data were obtained from the French data center THEIA (https://www.theia-land.fr, accessed on 10 January 2024), which uses the MAJA processing chain for pixel quality flags and atmospheric corrections [98].

A preliminary assessment of the temporal behavior of a few pixels with a stable vegetation cover, such as forests and natural grasslands, showed that there were no significant differences between the average reflectance values of each mentioned period. This suggests that this Sentinel-2 multi-source data, without considerable impact on the signal mean values, can be used in phenological crop assessment. The data from the GEE platform were primarily chosen because we also used GEE’s computational capabilities and we thus did not have to download images.

## 2.3. Ancillary Data

### 2.3.1. Land Use Map

To generate disaggregated NDVI time series using the SD method, we need to determine the fraction of each crop within each PROBA-V (300 m) mixel. Therefore, we performed a geometric interception between the pixel grid and a land use map, which was obtained from the French Land Parcel Identification System (LPIS), available on the Institut national de l'information géographique et forestière (IGN) data distribution site (https://www.geoportail.gouv.fr, accessed on 10 January 2024). Since 2015, LPIS has provided annual maps of crop types declared by farmers at plot level [91]. We used all available plots for winter wheat, winter barley, and oilseed rape (Table 1). Here, we added
the winter barley plots to the winter wheat ones due to their similar growth patterns and indistinguishable radiometric responses observed in this work. LPIS shows that the average plot area (Table 1) is approximately half smaller than the PROBA-V pixel surface (9 ha). Together, these crops account for approximately 50% of the arable land (excluding grasslands), while the remaining area is cultivated with over 100 other crop types.

<table><caption>**Table 1.** Land use map average area (2016–2020).</caption><thead><tr><th></th><th>Total Area (ha)</th><th>Average Plot Area (ha)</th><th>% of LPIS Area</th></tr></thead><tbody><tr><td>Winter wheat</td><td>4,804,300</td><td>4.78</td><td>17.20</td></tr><tr><td>Winter barley</td><td>1,331,200</td><td>4.57</td><td>4.76</td></tr><tr><td>Oilseed rape</td><td>1,350,500</td><td>6.10</td><td>4.83</td></tr></tbody></table>

### 2.3.2. Ground Data

We used two phenology databases to evaluate the extracted phenometrics. Firstly, a comparison at regional scale between PROBA-V-based phenomaps and the Céré'Obs dataset was performed. The latter was obtained from the National Phenological Observatory for Cereals website (https://cereobs.franceagrimer.fr/cereobs-sp/#/, accessed on 15 April 2024). This dataset provides an indicator of the regional phenology progression, i.e., the percentage of area reaching a given phenological stage, as a function of date. This information is based on the analysis of field observations by agricultural advisers in *chambres d'agriculture*, further aggregated at the administrative regional level. Secondly, a comparison at the municipality level between satellite-based (PROBA-V and Sentinel-2) phenomaps and the TEMPO dataset was performed. The latter was obtained from the national phenological observation network website (https://tempo.pheno.fr, accessed on 15 April 2024). This dataset provides the date at which a given phenological stage happens at the plot level. However, the precise location of the plot is not provided, but only the coordinates of the commune to which it belongs. Table 2 presents the crop-specific phenological stages of interest, as well as the sensor from which they are estimated and their ground evaluation dataset. The selected phenological stages are described in Table S1, while Table S2 summarizes the available observations in each ground evaluation dataset. It should be noted that these phenological stages, as well as their spatial comparison level, were established according to the availability of ground data. In addition, the geographical spread of the evaluation was determined according to the sensor and the computational capacity it needed.

<table><caption>**Table 2.** Crop-specific phenological stages of interest. For each phenological stage, a phenometric (SOS or EOS) was estimated from PROBA-V (300 m) and eventually Sentinel-2 (10 m) data, depending on the spatial comparison level. The regional or municipal phenometrics median dates were compared to phenological stage median dates derived from Céré'Obs or TEMPO, respectively.</caption><thead><tr><th></th><th>BBCH Code <sup>1</sup></th><th>Phenometric <sup>2</sup></th><th>Sensor Data</th><th>Ground Data</th><th>Comparison Level</th><th>Geographic Spread</th></tr></thead><tbody><tr><td rowspan="6">Winter wheat</td><td>31</td><td>SOS</td><td>PROBA-V</td><td>Céré'Obs</td><td>region</td><td>nationwide</td></tr><tr><td>51</td><td>SOS</td><td>PROBA-V</td><td>Céré'Obs</td><td>region</td><td>nationwide</td></tr><tr><td>99</td><td>EOS</td><td>PROBA-V</td><td>Céré'Obs</td><td>region</td><td>nationwide</td></tr><tr><td>29</td><td>SOS</td><td>PROBA-V and Sentinel-2</td><td>TEMPO</td><td>municipality</td><td>inter-comparison sites</td></tr><tr><td>51</td><td>SOS</td><td>PROBA-V and Sentinel-2</td><td>TEMPO</td><td>municipality</td><td>inter-comparison sites</td></tr><tr><td>75</td><td>EOS</td><td>PROBA-V and Sentinel-2</td><td>TEMPO</td><td>municipality</td><td>inter-comparison sites</td></tr></tbody></table>
<table><caption>**Table 2. Cont.**</caption><thead><tr><th></th><th>BBCH Code <sup>1</sup></th><th>Phenometric <sup>2</sup></th><th>Sensor Data</th><th>Ground Data</th><th>Comparison Level</th><th>Geographic Spread</th></tr></thead><tbody><tr><td rowspan="3">Oilseed rape</td><td>31</td><td>SOS</td><td>PROBA-V and Sentinel-2</td><td>TEMPO</td><td>municipality</td><td>inter-comparison sites</td></tr><tr><td>65</td><td>NDVI<sub>local_min</sub> <sup>3</sup></td><td>PROBA-V and Sentinel-2</td><td>TEMPO</td><td>municipality</td><td>inter-comparison sites</td></tr><tr><td>73</td><td>EOS</td><td>PROBA-V and Sentinel-2</td><td>TEMPO</td><td>municipality</td><td>inter-comparison sites</td></tr></tbody></table>

<table><caption>**Table 3.** Crop-specific phenometrics extraction. For each phenometric, the occurrence date was estimated from fitted NDVI time series. These dates were extracted using either the threshold or the derivative method, depending on the associated phenological stage. The threshold value is calibrated on an annual basis.</caption><thead><tr><th></th><th>Phenometric <sup>1</sup></th><th>Associated Phenological Stage</th><th>BBCH Code</th><th>Method</th></tr></thead><tbody><tr><td rowspan="5">Winter wheat</td><td>SOS<sub>42-53</sub></td><td>End of tillering</td><td>29</td><td>threshold-based</td></tr><tr><td>SOS<sub>54-60</sub></td><td>Stem elongation</td><td>31</td><td>threshold-based</td></tr><tr><td>SOS<sub>95-98</sub></td><td>Heading</td><td>51</td><td>threshold-based</td></tr><tr><td>EOS<sub>60-89</sub></td><td>Development of fruit</td><td>75</td><td>threshold-based</td></tr><tr><td>EOS<sub>10-15</sub></td><td>Senescence</td><td>99</td><td>threshold-based</td></tr><tr><td rowspan="3">Oilseed rape</td><td>SOS<sub>30-45</sub></td><td>Stem elongation</td><td>31</td><td>threshold-based</td></tr><tr><td>NDVI<sub>local_min</sub> <sup>2</sup></td><td>Flowering</td><td>65</td><td>derivative-based</td></tr><tr><td>EOS<sub>97-99</sub></td><td>Development of fruit</td><td>73</td><td>threshold-based</td></tr></tbody></table>

## 2.4. Crop-Specific Phenology from Time Series

### 2.4.1. Vegetation Index and Phenological Stages

In France, the end of tillering (i.e., the increase in the horizontal footprint of the plants) in winter wheat (BBCH29) happens when post-winter vegetative activity restarts (February–March), leading to a gradual increase in NDVI (Figure 3). Development advances from stem elongation (BBCH31) in March to heading (BBCH51) in May, achieving seasonal NDVI peak value. The development of fruit (BBCH75) occurs in June, characterized by a yellowing of plants due to a reduced chlorophyll content, which leads to a decrease in NDVI. In July, at the end of the cycle, senescence (BBCH 90–99) occurs, with plants drying up completely and NDVI reaching its minimum value for the growing season.

For oilseed rape, stem elongation (BBCH31) occurs in February–March, when the NDVI starts to rise after a period of stagnation during the winter (Figure 3). Flowering in April to May results in a yellow canopy as flowers hide the chlorophyllous parts. This leads to a gradual decrease in NDVI due to an increase in red-band (0.65–0.75 µm) reflectance. Then, from full flowering (BBCH65), NDVI rises as petals fall and the green parts of the plant become visible again. This convex behavior of NDVI during flowering is in contrast to the yellowness indices, such as the Normalized Difference Yellow Index used by [78], which are typically used to monitor flower dynamics. Fruit development (BBCH73) takes place after flowering (May). Here, the plant's green pods raise NDVI up to maturation stage, which is the inflection point where the slope of the curve becomes negative. Finally, senescence (BBCH90–99) appears at the end of the cycle, in July, when the plants become dry and NDVI returns to around the minimum value for the growing season.

<footnote n="1">¹ BBCH classification establishing the main and secondary phenological stages, as detailed by [23]. ² The threshold value of each indicator varies according to the phenological stage of interest. These values are indicated in Table 3.</footnote>

<footnote n="3">³ Local minimum of the NDVI time profile, used as a phenological indicator.</footnote>

<footnote n="1">¹ The threshold value, given as a % of the amplitude, changes annually and falls within the presented range of possible values. ² Local minimum of the NDVI time profile, used as a phenological indicator.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** NDVI time profile for (A) winter wheat and (B) oilseed rape. Values represent the national average derived from disaggregated PROBA-V (300 m) data in 2019, which were fitted using the Whittaker smoother model. Vertical lines indicate observed median dates across France in 2019, of all phenological stages of interest. Vertical line labels represent winter wheat tillering (BBCH29), stem elongation (BBCH31), heading (BBCH51), development of fruits (BBCH75), and senescence (BBCH99); and oilseed rape stem elongation (BBCH31), flowering (BBCH65), and development of fruits (BBCH73). These dates were obtained from the TEMPO dataset, except for those of winter wheat stem elongation (BBCH31) and senescence (BBCH99), which were obtained from the Céré'Obs dataset. Panel (A) details the amplitude definition for each side of the curve and the calibrated threshold value for each phenometric associated with each phenological stage of interest. Curve sides are relative to the maximum of the growing season. In panel (B), calibrated threshold values are shown for phenometrics associated with stem elongation (BBCH31) and development of fruits (BBCH73), which were obtained using the same amplitude definition process detailed in panel (A). For flowering (BBCH65), the associated phenometric ($NDVI_{local\_min}$) is shown within the corresponding temporal window (gray band). A local minimum occurs when the first derivative (dashed line) is zero at time $t$, negative at $t - 1$, and positive at $t + 1$.</figcaption></figure>

### 2.4.2. Phenometrics Extraction

Phenometrics were extracted at different spatial scales, according to the sensor. This was carried out pixel by pixel for PROBA-V (300 m) and plot by plot for Sentinel-2 (10 m). The accuracy of the phenometrics varies with different interpolation and extraction methods, depending on their sensitivities to artifacts and patterns in the curve. Finding a universally applicable best method is particularly challenging [41,90]. In this context, time series were reconstructed by adjusting curves following two algorithms: HANTS (Harmonic Analysis of Time Series, [99]), using a python implementation developed by [100], for general phenometrics (SOS and EOS) detection, and the Whittaker smoother [89] for the specific oilseed rape $NDVI_{local\_min}$ detection, because of its higher sensitivity to small changes in the curve. To calibrate the HANTS parameters, we used a parameter optimization technique commonly used in machine learning, known as *GridSearch* (https://scikit-learn.org/stable/modules/generated/sklearn.model_selection GRIDSearchCV.html, accessed on 25 November 2024). In this approach, we randomly selected a sample of pixels (100 per class per year) and systematically evaluated all possible combinations of pre-defined parameter values. The combination that yielded the best-performing model was then applied to all pixels. The optimized parameters were as follows (start:end:step:selected value): Number of frequencies (2:5:1:3); Fitting error tolerance (0.01:0.15:0.01:0.05); Degree of over-determinedness (0:20:1:1); Delta (0:1:0.1:0.1). Regarding non-optimizable parameters, settings were as follows: Length of base period = length of time series; HiLo = "low"; Valid range = 0, 1.

Subsequently, phenometrics were extracted using two methods: the threshold-based method applied on the HANTS-adjusted series for general phenometrics extraction, and
the derivative-based one applied on the Whittaker-adjusted curve for the specific oilseed rape NDVI<sub>local_min</sub> extraction.

For the threshold-based method, we established a value in percentage of the amplitude, and extracted the first date on which NDVI is larger (smaller) than that value (e.g., Equation (2)) on the ascending (descending) phase, for SOS (EOS), respectively. Amplitude was calculated for each side of the curve using the maximum NDVI during the growing season to differentiate between the sides (left or right), as shown in [101]. For instance, the left-side amplitude was calculated from the maximum NDVI during the growing season and the minimum NDVI observed on this side (i.e., the lowest NDVI before the occurrence date of the maximum NDVI). The left-side amplitude was used for SOS and the right-side one for EOS.

For the derivative-based method, we identified critical points in the curve, and extracted the date on which NDVI exhibited a local minimum (NDVI<sub>local_min</sub>) or maximum (NDVI<sub>local_max</sub>). For instance, a local minimum was identified when a critical point had a slope of zero and when previous and subsequent observations had a higher NDVI. Subsequently, the local minimum with the lowest NDVI value was selected from all the identified points. This was performed within a seasonal pre-established temporal window, in which the associated phenological event (i.e., oilseed rape flowering) is expected to occur (winter, spring...).

After comparing the crop-specific time series and ground data, we observed that each of the phenological stages of interest could be associated with a given phenometric (SOS or EOS using a given threshold value). This led to the application of the threshold-based method as a general approach. The threshold value for each phenometric was determined by comparing the national-level average time profile with ground data of the associated phenological stage, which were calculated from the PROBA-V (300 m) and the TEMPO dataset, respectively (Table 3). However, due to the pronounced convex behavior of NDVI during oilseed rape flowering, the derivative-based method was employed to determine the date corresponding to the local minimum (NDVI<sub>local_min</sub>). These methods, as well as those commonly found in the literature, are detailed in [90].

We calculated a phenometric's date of occurrence as follows:

$$SOS_k = t \in [0, t_{max}] \mid \left( NDVI(t) = NDVI_{min} + \frac{k}{100} \times A \right) \quad (2)$$

SOS<sub>k</sub>: Start of Season, used as an example here.

t<sub>max</sub>: Date of occurrence of the maximum NDVI during the growing season (NDVI<sub>max</sub>).

NDVI<sub>min</sub>: The minimum NDVI recorded between the 1st January and t<sub>max</sub>.

NDVI(t): NDVI value at date t.

A: Amplitude extracted from the left side of the curve. A is calculated as NDVI<sub>max</sub> - NDVI<sub>min</sub>.

k: Threshold value adjusted for each phenological event as explained above.

## 2.5. Assessment of Phenometrics

Here, the dates of occurrence of phenometrics and their associated phenological stages are referred to as estimated and observed dates, respectively.

### 2.5.1. Regional Comparison Between PROBA-V (300 m) and the Céré'Obs Dataset

We evaluated our results at the regional scale across mainland France. For each phenometric, according to ground data availability, we compared estimated (phenometric) and observed (associated phenological stage) median dates. Here, regional ground data were available for phenological stages of winter wheat (BBCH31, BBCH51, and BBCH99) as shown in Table 2. These crop-specific phenological median dates were extracted from both PROBA-V-based phenomaps and from the Céré'Obs dataset. Subsequently, statistical indicators, such as the Pearson correlation (r), coefficient of determination (r²), root mean square difference (rmsd), and bias, were calculated using the Equations (3)–(6).
Then, regional estimated and observed phenological progressions were also compared. We evaluated winter wheat SOS$_{54–60}$ and SOS$_{95–98}$ across three regions (Nouvelle-Aquitaine, Centre-Val de Loire, and Hauts-de-France), illustrating phenological contrasts in the south-west, center and north of mainland France. Here, we fitted a logistic function to the cumulative date distribution of these two datasets before comparing them using distance statistics (i.e., rmsd and bias).

$$r = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n} (x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n} (y_i - \bar{y})^2}} \quad (3)$$

$$r^2 = 1 - \frac{\sum_{i=1}^{n} (x_i - y_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \quad (4)$$

$$rmsd = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (x_i - y_i)^2} \quad (5)$$

$$bias = \frac{1}{n} \sum_{i=1}^{n} (x_i - y_i) \quad (6)$$

where *n* is the number of RoI (region or municipality), $x_i$ and $y_i$ are estimated and observed median dates at RoI *i*, and $\bar{x}$ and $\bar{y}$ are mean of *x* and *y*. These dates are in terms of Day of Year (DoY).

## 2.5.2. Municipality Comparison Between Sensors (PROBA-V and Sentinel-2) and the TEMPO Dataset

We also assessed the sensor's performance at the municipality scale across inter-comparison sites (Figure 2). For each phenometric, estimated and observed median dates were compared. Here, municipality ground data were available for phenological stages of winter wheat (BBCH29, BBCH51, and BBCH75) and oilseed rape (BBCH31, BBCH65, and BBCH73) as shown in Table 2. These crop-specific phenology median dates were extracted from both sensor-based (PROBA-V and Sentinel-2) phenomaps and from the TEMPO dataset. Subsequently, for each sensor, we calculated the comparative statistical indicators mentioned above. Finally, we also compared estimated dates from both sensors with each other, as well as the detection rate of each one. The detection rate refers to the proportion of the number of municipalities with crop-specific phenology estimates to the total number of municipalities. When a municipality lacks phenology estimates, it indicates that no pixel within that area has accumulated a satisfactory number of cloud-free observations to extract crop-specific phenometrics.

## 3. Results

PROBA-V-based crop-specific annual phenomaps were generated at the national level. In general, we observed a spatial gradient from south-west to north for all phenometrics (e.g., Figures 4 and S1), following the latitudinal climatic characteristics in mainland France [93]. Here, we evaluated these phenomaps at the regional scale across the country (Section 2.5.1), before conducting the evaluation at the municipality scale across inter-comparison sites (Section 2.5.2). In these inter-comparison sites, we also generated Sentinel-2-based phenomaps in order to compare the performance of both sensors. An annual comparison was carried out from 2016 to 2020 separately, in order to check whether or not the inter-annual variations were captured by the remote sensing methods.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** PROBA-V-based phenomaps of winter wheat. Each column represents a phenometric associated with a phenological stage available in the Céré'Obs database, i.e., SOS<sub>54-60</sub> with stem elongation (BBCH31), SOS<sub>95-98</sub> with heading (BBCH51) and EOS<sub>10-15</sub> with senescence (BBCH99), respectively. Each row represents a year in our study period. The color palette represents the day of the year (DoY) on which the phenometric was detected.</figcaption></figure>

### 3.1. PROBA-V-Based Phenomaps Versus Céré'Obs

On the regional scale, PROBA-V-based phenomaps (SOS<sub>54-60</sub>, SOS<sub>95-98</sub>, and EOS<sub>10-15</sub>) were compared with ground data for the phenological stages available in the Céré'Obs dataset (Table 2): winter wheat stem elongation (BBCH31), heading (BBCH51), and senescence (BBCH99), respectively.
### 3.1.1. Regional Comparison of Phenological Median Dates

We estimated winter wheat SOS<sub>54-60</sub> between late January and mid-April across our study period. Phenomaps revealed an earlier phenology in meridional regions and in two occidental ones (Bretagne and Pays-de-la-Loire), as shown in Figure 4. On average over 2016–2020, this phenometric was estimated to occur on the 82th DoY ± 4.3 days. Estimated (SOS<sub>54-60</sub>) and observed (stem elongation, BBCH31) dates diverged from 4.4 to 13.5 days, depending on the year (Table S3). The highest difference was observed in 2016. However, in this year, the PROBA-V-based phenomap consistently captured the spatial gradient of the observed phenology (r² > 0.60). Overall, stem elongation was estimated 1.5 days later while 60% of its observed phenology was explained. SOS<sub>95-98</sub> was estimated between mid-April and early June, with an average date at DoY 133 ± 4.8 days. Here, the highest difference with ground data (heading, BBCH51) was observed in 2018 (rmsd = 7.7 days), while average bias was –2 days. This phenometric obtained a better correspondence with the associated phenological stage compared to SOS<sub>54-60</sub> and EOS<sub>10-15</sub> (r² higher by 13% and 35%, respectively). EOS<sub>10-15</sub> was estimated between mid-June and late July, with the average date at DoY 201 ± 4.6 days. Although phenomaps showed a difference <6 days on average with the Céré'Obs dataset (senescence, BBCH99), they only captured 37% of the observed senescence variability (9% in 2018). This phenometric exhibited the lowest bias (1.3 days on average) compared to the other ones.

Overall, PROBA-V-based phenomaps were consistent with the regional ground data provided by Céré'Obs (Figure 5), with a global (all stages and years) rmsd = 6.4 days ± 2.2 days (Table S3). For all phenometrics, estimated dates consistently captured inter-annual variations of the associated observed dates. However, there was less agreement in their spatial representation, since PROBA-V-based phenomaps exhibited higher inter-regional dynamics than the Céré'Obs dataset. This suggests that phenomaps could complement these regionally aggregated ground data by providing the missing intra-regional spatial dynamics.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Regional comparison of phenometrics (from PROBA-V) and phenological stages (from Céré'Obs) for winter wheat. (A) SOS<sub>54-60</sub> vs. stem elongation, (B) SOS<sub>95-98</sub> vs. Heading, and (C) EOS<sub>10-15</sub> vs. senescence. Each point represents a region median date and its color the year of interest. Median dates are expressed in the day of the year (DoY).</figcaption></figure>

### 3.1.2. Regional Comparison of Phenological Progression

We also evaluated intra-regional estimated (phenometric) and observed (associated phenological stage) phenological progression in terms of % of area (e.g., Figure 6). Here, winter wheat SOS<sub>54-60</sub> and SOS<sub>95-98</sub> were assessed across three regions representative of mainland France's phenological contrast (Hauts-de-France, Centre-Val de Loire, and Nouvelle-Aquitaine) for all years.

On the one hand, SOS<sub>54-60</sub> progression showed an average rmsd = 22.9% of area ± 12.8% (all regions) compared to the observed (stem elongation, BBCH31) progression (Table S4), with the highest discrepancies in the south-west (Nouvelle-Aquitaine, rmsd = 26.4% of area). On the other hand, SOS<sub>95-98</sub> progression obtained better correspondence with observed (heading, BBCH51) progression showing an average rmsd = 14.6% of area ± 8.5% (all regions). In contrast to SOS<sub>54-60</sub>, the largest discrepancies for SOS<sub>95-98</sub> were observed in the north (Hauts-de-France, rmsd = 16.3% of area). In general, phenology progression was better estimated (all stages and years) in the center of the country (Centre-Val de Loire, rmsd = 17.4% of area) and with less precision for
the year 2016 (all stages and regions). Finally, phenomaps sub-estimated the observed phenology progression for most comparisons (bias = -5.3% of area on average).

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Intra-regional comparison of phenometrics and associated phenological stages in terms of phenological progression in 2019. Phenometrics progression (in % of area) was obtained from the PROBA-V-based phenomaps (dashed lines), while phenological stage progression was obtained from Céré'Obs (solid lines). Each column represents a specific winter wheat phenometric and its associated phenological stage: (A–C) SOS<sub>54–60</sub> vs. stem elongation and (D–F) SOS<sub>95–98</sub> vs. heading, while each row represents a region of interest. The curves represent the values fitted with a logistic function.</figcaption></figure>

## 3.2. Sensor-Based Phenomaps Versus TEMPO

The sensor-based (PROBA-V and Sentinel-2) phenomaps were separately compared with ground data for the phenological stages available in the TEMPO dataset (Table 2): for winter wheat, SOS<sub>42–53</sub> vs. tillering (BBCH29), SOS<sub>95–98</sub> vs. heading (BBCH51), and EOS<sub>60–89</sub> vs. development of fruit (BBCH75); for oilseed rape, SOS<sub>30–45</sub> vs. stem elongation (BBCH31), NDVI<sub>local_min</sub> vs. flowering (BBCH65), and EOS<sub>97–99</sub> vs. development of fruit (BBCH73). Here, as TEMPO does not provide the location of the observed plots, we evaluated the median aggregated dates at the municipality scale across inter-comparison sites (Section 2.5.2). Overall, both PROBA-V- and Sentinel-2-based phenomaps demonstrated similar performances whatever phenometric, but only after the launch of Sentinel-2B in 2017. These sensor-based phenomaps were in moderate agreement with observed ground data, mainly due to their greater spatial variability compared with ground observations.
### 3.2.1. Winter Wheat Phenometrics

Globally, for winter wheat, PROBA-V- and Sentinel-2-based phenomaps were moderately consistent with the ground data provided by TEMPO (Figure 7), with global (all stages and years) rmsd values of 8.8 ± 5.8 and 9.6 ± 7.9 days, respectively (Table S5). As observed in regional comparisons, and this whatever phenometric, the estimated dates reliably reflected the inter-annual variations of the observed dates. However, on average, the PROBA-V- and Sentinel-2-based phenomaps explained only 35% and 32% of the observed intra-annual phenology (all stages and years), respectively. Indeed, the annually sensor-based phenomaps reflected higher inter-municipality variability than the TEMPO dataset. Additionally, when we compared estimates and observations within inter-comparison sites, PROBA-V was able to estimate phenology in more municipalities than Sentinel-2. Phenometrics based on central or high threshold values (SOS95–98 and EOS60–89) showed better correspondences with the associated phenological stages across our study period (2016–2020) with an average rmsd = 5.8 ± 0.9 and 6.0 ± 2.2 days, by PROBA-V and Sentinel-2, respectively. These sensors demonstrated a shared lower performance in estimating SOS42–53, with an average rmsd = 11.8 ± 2.8 days over 2017–2020 and > 25 days in 2016 (Table S5).

<figure><img src="image_11.png" /></figure>

**Figure 7.** Municipal comparison of phenometrics and associated phenological stages in terms of median dates. Phenometrics dates were obtained from both PROBA-V- and Sentinel-2-based phenomaps, while phenological stages were obtained from the TEMPO dataset. Each column represents a specific winter wheat phenometric compared to associated phenological stage: (A–D) SOS42–53 vs. tillering, (E–H) SOS95–98 vs. Heading, and (I–L) EOS60–89 vs. development of fruits. First row shows the locations where ground observations of each stage were made. Second row shows scatter-plots from all these observed municipalities, while third and fourth rows show scatter-plots from only the observed municipalities within the inter-comparison sites (dashed boxes), where
we were able to extract phenometrics from Sentinel-2: third row represents estimated dates from PROBA-V, while fourth row represents those from Sentinel-2. In these last two rows, assessed municipalities were identical for both sensors. Scatter-plots and maps share the same color code representing the year of interest. Median dates are expressed in the day of the year (DoY).

### 3.2.2. Oilseed Rape Phenometrics

For oilseed rape, the statistical interpretations below were derived from years 2018 and 2019. For the other years, oilseed rape data actually had some limitations: on the one hand, for 2016–2017, cloud-free Sentinel-2 data were not sufficient to extract metrics; on the other hand, for 2020, perhaps due to the curfew associated with the health crisis, oilseed rape ground observations were not sufficient to yield consistent comparisons.

In general, PROBA-V- and Sentinel-2-based phenomaps exhibited a lower accuracy (Figure 8) than those generated for winter wheat, with global (all stages and years) rmsd values of 10.8 ± 5.9 and 12.8 ± 7.9 days, respectively (Table S5). Here, for all phenometrics, the sensor-based estimated dates also captured the inter-annual variations of the observed dates (Figure 8). However, on average, PROBA-V- and Sentinel-2-based phenomaps explained only 22% and 16% of the observed intra-annual phenology (all stages and years), respectively. Additionally, as observed for winter wheat, when we compared estimates and observations within inter-comparison sites, PROBA-V was able to estimate phenology in more municipalities than Sentinel-2.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Municipal comparison of phenometrics and associated phenological stages in terms of median dates. Phenometrics dates were obtained from both PROBA-V and Sentinel-2-based phenomaps, while phenological stages came from the TEMPO dataset. Each column represents a specific oilseed rape phenometric compared to associated phenological stage: (A–D) $SOS_{30-45}$ vs. stem elongation,</figcaption></figure>
(E–H) NDVI<sub>local_min</sub> vs. Flowering, and (I–L) EOS<sub>97–99</sub> vs. development of fruits. First row shows the locations where ground observations of each stage were made. Second row shows scatter-plots from all these observed municipalities, while third and fourth rows show scatter-plots from only the observed municipalities within the inter-comparison sites (dashed boxes), where we were able to extract phenometrics from Sentinel-2: third row represents estimated dates from PROBA-V, while fourth row represents those from Sentinel-2. In these last two rows, assessed municipalities were identical for both sensors. Scatter-plots and maps share the same color code representing the year of interest. Median dates are expressed in the day of the year (DoY).

As already observed for winter wheat, phenometrics based on central or high thresh-
old values (NDVI<sub>local_min</sub> and EOS<sub>97-99</sub>) showed better correspondences with associated
phenological stages, with average rmsd = 7.5 ± 0.9 and 7.8 ± 1.3 days, for PROBA-V
and Sentinel-2 respectively. Differences regarding SOS<sub>30-45</sub> and stem elongation were
around 20.0 days. Finally, for the period where Sentinel-2 valid data were not suffi-
cient to extract phenological metrics (2016, 2017, and 2020), PROBA-V-based phenomaps
showed, on average for all stages and years, an acceptable discrepancy with observed dates
(rmsd < 10 days). However, in 2016–2017, estimated dates for EOS<sub>97-99</sub> did not capture the
observed spatial variations of the development of fruits.

3.2.3. Inter-Comparison of PROBA-V- and Sentinel-2-Based Phenomaps

We also compared sensor-based phenomaps with each other to assess the consistency of remotely sensed data and their sensitivity to vegetation dynamics as a function of their spatial resolution. Here, we evaluated the median aggregated dates at the municipality scale across inter-comparison sites. Overall, for all crop-specific phenometrics and years, PROBA-V-based phenomaps were correlated with those generated from Sentinel-2, but only after the launch of Sentinel-2B in 2017 (Table S6). Hence, statistical metrics below refer to the 2017–2020 period. For winter wheat, all phenometrics and years, sensor-based phenomaps showed an r = 0.8 ± 0.1 and an rmsd = 7.8 ± 2.4 days (Table S6). Among phenometrics, SOS<sub>42–53</sub> was better correlated (r = 0.8), compared to SOS<sub>95–98</sub> and EOS<sub>60–89</sub>, which showed lower discrepancy (rmsd < 8.8 days). For oilseed rape, sensor-based phenomaps showed lower spatial correspondence than winter wheat ones (r = 0.5 ± 0.2), with a temporal divergence of rmsd = 10.7 ± 7.6 days. For this crop, the sensors showed better agreement with each other when estimating NDVI<sub>local_min</sub> (r = 0.7 and rmsd = 4.3 days, on average), while SOS<sub>30–45</sub> showed the lowest correspondence (r < 0.5).

Additionally, we also assessed their phenological detection rate, i.e., the proportion
of municipalities where cloud-free observations were sufficient to extract crop-specific
phenometrics. For the 2016–2020 period and all phenometrics, we were able to extract
phenometrics from PROBA-V data in 94.4 ± 7.2% of the municipalities, and from Sentinel-2
data in 71.8 ± 18.8% (Figure 9). For 2018–2020, during which Sentinel-2 yielded the highest
number of cloud-free observations over an agricultural period, we were able to map
83.4 ± 7.0% of municipalities. However, before 2018, our ability to generate phenomaps
from Sentinel-2 decreased by 29%, especially for oilseed rape NDVI<sub>local_min</sub> detection (<50%
of municipalities).
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Phenological detection rate from both PROBA-V and Sentinel-2 data. Each row represents a specific crop, and each column a specific phenometric: (A–C) winter wheat SOS<sub>42–53</sub>, SOS<sub>95–98</sub> and EOS<sub>60–89</sub>, (D–F) oilseed rape SOS<sub>30–45</sub>, NDVI<sub>local_min</sub>, and EOS<sub>97–99</sub>. The assessment was conducted based on all municipalities in the intercomparison sites, not only on those in which ground observations were carried out. On average between 2016 and 2020, 1943 municipalities were assessed for winter wheat and 1293 for oilseed rape.</figcaption></figure>

### 3.3. Application: Analyzing and Comparing the Temporal Variations in Phenometrics for Wheat and Oilseed Rape

We analyzed the spatial and temporal variations in phenometrics for the two crops and for each year in 2016–2020. For each crop, we consider only here the phenometrics with the best performances: SOS<sub>95–98</sub> (associated with heading) for the winter wheat, and NDVI<sub>local_min</sub> (associated with flowering) for the oilseed rape.

On average over the 5 years, the spatial mean of wheat SOS<sub>95–98</sub> is the DOY 134 (14 May), with a standard deviation of 8.7 days, with the western and southern parts of the country displaying earlier dates. Moreover, at each location the date is variable from one year to the other: the multi-year range of variations is computed for each disaggregated mixel, showing that in 2016–2020, on average over France the range is 24.7 days, with a standard deviation of 10.2 days. The two standard deviation values are of similar order, showing that the temporal and the spatial variations in this phenometric are almost similar. The same property is found with the interquartile ranges that are equal to 6.3 and 6.5 days, respectively.

For oilseed rape, the standard deviation of the spatial variations is 8.7 days, and the one for the range of changes across the 5 years is 11 days. These values are calculated over a much smaller area than for wheat because of the smaller cultivation area and cannot strictly be compared to them. However, we find that oilseed rape NDVI<sub>local_min</sub> occurs 22 days before wheat SOS<sub>95–98</sub>, on average over the common growing area, with very large spatial differences. This lag between the wheat and the oil seed rape rises up to 45 days in the northern part of the country but is only 6 days in the central parts. This shows a highly diverse relationship between the two crop phenologies, that may be explained either by different cultural practice calendars inside the northern half of France, and/or various phenological response of the two crops to the seasonal climate variation.

## 4. Discussion

### 4.1. Phenomaps: Regional and Municipal Accuracy

The PROBA-V-based phenomaps for winter wheat demonstrated a strong agreement with the associated phenological stages observed at a regional level, provided by the national phenological observatory Céré'Obs. Estimated dates, when aggregated by region, reflected the phenological contrast between the south, center, and north of the country,
showing a spatial gradient consistent with the latitudinal variation of climate in mainland France [93]. However, as observed from other satellite-based phenology estimates [102–104], for some cases, phenometrics were more spatially dynamic than the associated phenological stages, resulting in a low coefficient of determination (r²), even with a low rmsd. Here, observed dates were obtained from a regional synthesis of the phenological progression, derived from an analysis by agricultural advisers from local *chambres d'agriculture*. Hence, this phenological information is an extrapolation based on the phenology observed at specific points throughout the area, which might overlook the phenological diversities, particularly in regions characterized by significant climatic and environmental contrast. Therefore, this data collection approach could explain the lack of dynamism in the spatial representation of the observed regional phenology. Consequently, this discrepancy in the spatial representation of phenology could significantly account for the low r² values. The estimated intra-regional phenology progression (% of area) was more accurate in the center of the country, perhaps because this in situ observation system is more sensitive and consistent for these regions with the largest winter wheat areas in the country.

Estimated dates, when aggregated by municipality, showed lower agreement with observed dates from the TEMPO network, compared to the observed regional phenology from the Céré'Obs network. We mainly observed a low-spatial correspondence (r² < 0.50), yet reasonable temporal differences (rmsd < 15 days), in agreement with satellite-based performances usually reported in the literature [32,80,102,105]. Here, observed dates were obtained from a visual inspection of the timing of phenological events, conducted by observers in some specific plots throughout the municipality (5 plots on average). Hence, TEMPO observations do not attempt to represent municipal phenology, unlike those of Céré'Obs, which attempt to summarize regional phenology. These observed dates were provided without the location of plots, which led to their municipal aggregation to compare them with the estimated dates. Therefore, these aggregated dates are likely not the best indicator of the municipal phenological diversity, which could also explain the moderate spatial correspondence with PROBA-V- and Sentinel-2-based phenomaps. For a robust validation of satellite-derived phenometrics, ongoing efforts to enhance the spatial representation of the observed phenology are crucial.

## 4.2. Accuracy of Crop-Specific Phenometrics

Overall, winter wheat phenometrics showed better correspondences with associated phenological stages than oilseed rape ones. This can reliably be attributed to the accuracy of the disaggregated data used to extract these phenometrics. According to [81], disaggregated data perform better for the main crops in the study area because they more significantly contribute to the signal reconstruction. The spatial disaggregation algorithm is sensitive to the crop's occurrence level within the spatial window used for extracting information. Therefore, it is likely that the reconstructed NDVI time series for winter wheat, the dominant crop in our study area, more accurately captures the phenological variations of this specific crop.

In addition, crop-specific phenometrics based on central or high threshold values (associated with heading and development of fruits) exhibited higher accuracy than those with lower threshold values (associated with tillering and stem elongation), perhaps due to their more stable radiometric response. Winter wheat SOS<sub>95-98</sub> was located around maximum NDVI of the growing season, while the EOS<sub>60-89</sub> was located at the beginning of the temporal profile decline (Figure 3). Oilseed rape flowering leads to a consistent decrease in NDVI between April and May as the flowers hide the leaves. The well-defined radiometric response of this stage was also observed by [78] when they used the Normalized Difference Yellow Index (NDYI), which consistently increased over the same period. Hence, such identifiable behaviors made NDVI<sub>local_min</sub> a reliable indicator of this phenological stage. In contrast, phenometrics associated with tillering and stem elongation showed less reliability, because these associated phenological stages demonstrated more ambiguous radiometric responses, likely due to higher variability in vegetation condition and/or
coverage among plots before and during winter, as highlighted by [80] when they used optical data. For instance, a warmer autumn may stimulate vigorous plant growth without accelerating phenology. As a result, the plot will have high NDVI values during winter, potentially enabling the earlier estimation of these stages. Finally, the detection of phenometrics associated with senescence could be influenced by several factors, including the decreasing quality of the disaggregated NDVI time series and the applied metric extraction method itself. Moreover, toward the end of the growing season, plant stress arising from water or nutritional deficiencies may act as a confounding factor, potentially resulting in an earlier-than-expected estimation.

### 4.3. Inter-Comparison of PROBA-V and Sentinel-2 Performances

Disaggregated PROBA-V data allowed us to map crop-specific phenology at the national level, providing a comprehensive view of crop vegetation dynamics in mainland France. In general, disaggregated NDVI time series, with a daily temporal resolution, were not significantly affected by the cloud cover, offering an advantage over high-spatial resolution data with a lower frequency. However, this advantage did not have an impact on the accuracy of phenometrics estimation. In contrast, in terms of phenological detection rate, PROBA-V-based phenomaps showed a strong improvement over Sentinel-2, even during the era of the two Sentinel-2A and B satellites (2018–2020) (Figure 9). Finally, generating PROBA-V-based phenomaps requires less computational resources compared to Sentinel-2. This significant advantage increases the analysis capacity regarding the spatial and temporal coverage.

### 4.4. Phenometrics Extraction Method

Threshold-based methods require extensive knowledge of local production systems, as well as ground data to facilitate time profile interpretation. This is mandatory to take into account, as far as possible, the highly dynamic and complex nature of agricultural landscapes. Here, for each phenometric associated with a given phenological stage, a specific threshold value was determined by annually comparing the national average NDVI time profile with the corresponding national average observed date. Although these thresholds are relative to the amplitude of the pixel, this approach may be less sensitive to the local phenology. For instance, the spatial variation in plant development (not phenology), mainly driven by the local weather conditions and agricultural practices (density, fertilization, irrigation, etc.), could exhibit contrasted radiometric responses between regions. Hence, these thresholds may not directly be applied to local time profiles significantly different from the national time profile used for calibration. In other words, thresholds can only be transferred between areas with similar growing conditions [106]. To address this problem, current machine learning techniques could facilitate setting extraction parameters, such as shown in [79,107]. However, although promising, these techniques demand extensive ground data that are generally both expensive and time-consuming to obtain. Moreover, these methods are not exempt from generalization challenges. Finally, our detection method is based on adjusting a single curve in the vegetation index ascending phase, and another one in the descending phase, as it is generally done. It was recently shown that adjusting separate curves in a restricted time window around each phenological event before detecting the inflexion timings improves the results, as it reduces the interdependence of the detections of each event [108]. In the future, such methods could be tested over our disaggregated time series.

### 4.5. Limitations

Sensor-based phenomaps are mainly affected by the many limitations of satellite data and signal processing techniques. In this approach, a robust and automatic phenometrics extraction requires a smooth time series with a continuous time step. However, satellite data, including composite data, often contain noise caused by cloud residues or atmospheric contamination. In addition, cloud-free observations generally show irregular
temporal distances. Consequently, gap-filling methods are usually employed in time series analysis. In this case, the accuracy of the phenometric is directly influenced by the accuracy of the crop's time trajectory reconstruction, which can vary depending on the applied interpolation method [41]. Likewise, there is no agreement on the preferred method for extracting phenometrics from gap-filled data [102]. Ref. [90] reviewed sensor-derived vegetation phenometrics extraction methods, such as thresholds (absolute or relative) and derivatives approaches, and highlighted that the accuracy of estimates varies according to the applied method. In addition, sensor-based phenomaps remain difficult to evaluate due to the lack of spatial correspondence with available ground data. PROBA-V estimates also have certain limitations inherent in the spatial disaggregation technique, which is efficient only for the most dominant crops [81], and does not attempt to provide a phenometric at the plot scale: the phenological differences of the same crop growing up in two plots within the same 3 × 3 pixel neighborhood are not retrieved. Moreover, because of this limitation, and because of the need to know all crop types within the disaggregation neighborhood, our method is not designed for providing the near-real-time phenology information necessary for agricultural advice to individual farmers. Regarding input data reliability, the accuracy of the disaggregated PROBA-V data is dependent on the quality of the crop type map used. Here, disaggregated PROBA-V data were generated using the LPIS (Section 2.3.1). Although it is positioned as the reference crop type map in France, it may be subject to misdeclaration and would not include undeclared plots. Concerning this last point, a completed LPIS version (containing undeclared plots), generated by INRAE's Rural Development Observatory (Institut national de recherche pour l'agriculture, l'alimentation et l'environnement (INRAE-ODR): https://odr.inrae.fr/intranet/carto_joomla/, accessed on 25 November 2024) team [109], could improve geographic coverage of the disaggregated PROBA-V data, for available years and solely for our study site. Concerning misdeclaration, other crop type maps can be used such as the OSO-THEIA (https://www.theia-land.fr/ceslist/ces-occupation-des-sols/, accessed on 25 November 2024) product generated from high-spatial resolution satellite data [110]. From an operational perspective, the LPIS is updated annually as part of the European Union's Common Agricultural Policy. This would allow the generalization of our method to member states' territories with LPIS maps. However, the accuracy and spatial detail of the crop type maps are still undergoing harmonization across these countries.

## 5. Conclusions

By utilizing disaggregated PROBA-V data, which provide daily time-step NDVI time series for individual crops within a 300 m pixel [81], we were able to generate comprehensive annual crop-specific phenomaps at the national level from 2016 to 2020. These phenomaps provide a coherent overview of vegetation dynamics aligned with the latitudinal climatic patterns observed in mainland France. PROBA-V-based phenomaps exhibited a strong agreement with ground observations at the regional level, while showing moderate agreement at the municipal level. This discrepancy can be attributed to the limited representativeness of the observed data at the municipal level. Overall, our findings indicate that phenomaps exhibit greater spatial variability compared to currently available ground observations, resulting in low annual values for both rmsd and r² across most cases. Notably, sensor-based phenomaps successfully captured the inter-annual variations of phenology. Phenometrics based on central and high threshold values (associated with winter wheat heading and fruit development, as well as with oilseed rape fruit development) exhibited more accurate estimations than those with lower threshold values (associated with tillage and stem elongation of crops). NDVI<sub>local_min</sub> showed to be a reliable indicator of the oilseed rape flowering. This disparity can be attributed to the higher variation in plant development and coverage before and during the winter period (i.e., tillage and stem elongation of crops' occurrence period), primarily influenced by weather conditions and local farming practices. These performance levels were comparable to those achieved by Sentinel-2-based phenomaps. Sensor-based phenomaps derived from PROBA-V and Sentinel-2 demonstrated significant correlations with each other starting from 2017, coincid-
ing with the launch of Sentinel-2B. Nevertheless, PROBA-V-based phenomaps consistently exhibited a higher phenological detection rate, even when Sentinel-2 achieved its peak of cloud-free observations (2018–2020). Altogether, these findings highlight the significant potential of disaggregated PROBA-V data for operational monitoring of crop phenology in France, for winter wheat heading and oilseed rape flowering, as the performances were as good as those obtained with Sentinel-2 after 2018 but were higher before this date. This potential was used to explore the phenology differences between crops, showing that the lag between the wheat heading and the oilseed rape flowering differs from the north to the center of France.

Furthermore, these results lay the foundation for employing SPOT-VGT or MODIS data in retrospective studies, enabling the investigation of crop phenology evolution within fragmented agricultural landscapes. The method could be applied on these medium-resolution data from as early as 1999. However, before 2015, the French LPIS maps gave the crop type as a sole fraction of each farm, thus, our method would require to be adapted to this initial version of LPIS, which may be challenging. Nevertheless, such investigations can provide valuable insights into the adaptation and resilience of species in the face of ongoing climate change, as phenological changes are one of the key fingerprints of climate change in agriculture and ecosystems.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs16234521/s1. Table S1. Crop-specific phenological stages description; Table S2. Crop-specific phenological stages. Statistical summary of available observations at regional and municipal level from the Céré'Obs and the TEMPO dataset, respectively. The number of observations is denoted by n. Values represent average between 2016 and 2020. TEMPO's summary was made at inter-comparison sites level, while for Céré'Obs it was made at national level; Figure S1. PROBA-V-based phenomaps of oilseed rape. Each column represents a phenometric associated with a phenological stage, i.e., SOS30-45 with stem elongation (BBCH31), NDVIlocal_min with flowering (BBCH65) and EOS97-99 with development of fruit (BBCH73), respectively. Each row represents a year in our study period. The color palette represents the day of the year (DoY) on which the phenometric was detected; Table S3. Regional comparison of phenometrics and associated phenological stages for winter wheat. Phenometrics dates were obtained from PROBA-V-based phenomaps, while phenological stages ones were obtained from the Céré'Obs dataset. Annual statistical metrics (rmsd and bias) are expressed in days. The bias is negative when phenometrics dates are earlier than phenological stages dates. The number of compared regions is denoted by n; Table S4. Intra-regional comparison of phenometrics and associated phenological stages in terms of phenological progression. Phenometrics and their associated phenological stages progressions were obtained from the PROBA-V-based phenomaps and Céré'Obs, respectively. Annual statistical metrics (rmsd and bias) are expressed in terms of % of area. The bias is negative when the phenometric progression is lower than the associated phenological stage. The number of compared days is denoted by n; Table S5. Municipal comparison of phenometrics and associated phenological stages in terms of median dates. Phenometrics and phenological stages dates were obtained from sensor-based phenomaps and TEMPO, respectively. Annual statistical metrics (rmsd and bias) are expressed in days. The bias is negative when phenometrics dates are earlier than phenological stages dates. The number of compared municipalities is denoted by n; Table S6. Municipal comparison of phenometrics median dates obtained from both PROBA-V- and Sentinel-2-based phenomaps. Annual statistical metrics (rmsd and bias) are expressed in days. The bias is negative when PROBA-V estimates are earlier than Sentinel-2. The number of compared municipalities is denoted by n. Assessment was conducted based on all municipalities in the inter-comparison sites, not only on those in which ground observations were carried out.

**Author Contributions:** Conceptualization, C.O., F.M. and N.D.; methodology, H.R., C.O. and F.M.; software, H.R.; validation, H.R.; formal analysis, H.R.; investigation, H.R.; resources, C.O., F.M. and N.D.; data curation, H.R.; writing—original draft preparation, H.R.; writing—review and editing, C.O., F.M., N.D. and E.V.; visualization, H.R.; supervision, N.D.; project administration, N.D.; funding acquisition, E.V. All authors have read and agreed to the published version of the manuscript.
**Funding:** This research was funded by the TOSCA-POLYPHEME project of the CNES, grant number 200769/id5917 and BECAL-Paraguay program, grant number 206/2018.

**Data Availability Statement:** The raw data supporting the conclusions of this article will be made available by the authors on request.

**Acknowledgments:** The study was initiated when H.R. and N.D. were members of the PRODIG research unit which provided initial support. H.R. benefited from a research internship at the CESBIO research unit during this work.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## References

1. Lieth, H. Purposes of a Phenology Book. In *Phenology and Seasonality Modeling*; Lieth, H., Ed.; Springer: Berlin/Heidelberg, Germany, 1974; pp. 3–19, ISBN 978-3-642-51863-8.

2. Réaumur, R.A. Observations Du Thermomètre Faites à Paris Pendant l'année 1735, Comparées Avec Celles Qui Ont Été Faites Sous La Ligne, à l'Isle de France, à Alger et Quelques Unes de Nos Iles de l'Amérique. *Mem. l'Acad. R. Sci.* **1735**, 545–576.

3. Chmielewski, F.-M.; Rötzer, T. Response of Tree Phenology to Climate Change across Europe. *Agric. For. Meteorol.* **2001**, *108*, 101–112. [CrossRef]

4. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting Plant Phenology in Response to Global Change. *Trends Ecol. Evol.* **2007**, *22*, 357–365. [CrossRef] [PubMed]

5. de Beurs, K.M.; Henebry, G.M. A Land Surface Phenology Assessment of the Northern Polar Regions Using MODIS Reflectance Time Series. *Can. J. Remote Sens.* **2010**, *36*, S87–S110. [CrossRef]

6. Menzel, A.; Sparks, T.H.; Estrella, N.; Koch, E.; Aasa, A.; Ahas, R.; Alm-Kübler, K.; Bissolli, P.; Braslavská, O.; Briede, A.; et al. European Phenological Response to Climate Change Matches the Warming Pattern. *Glob. Change Biol.* **2006**, *12*, 1969–1976. [CrossRef]

7. Myneni, R.B.; Keeling, C.D.; Tucker, C.J.; Asrar, G.; Nemani, R.R. Increased Plant Growth in the Northern High Latitudes from 1981 to 1991. *Nature* **1997**, *386*, 698–702. [CrossRef]

8. Peñuelas, J.; Rutishauser, T.; Filella, I. Phenology Feedbacks on Climate Change. *Science* **2009**, *324*, 887–888. [CrossRef]

9. Richardson, A.D.; Andy Black, T.; Ciais, P.; Delbart, N.; Friedl, M.A.; Gobron, N.; Hollinger, D.Y.; Kutsch, W.L.; Longdoz, B.; Luyssaert, S.; et al. Influence of Spring and Autumn Phenological Transitions on Forest Ecosystem Productivity. *Philos. Trans. R. Soc. B Biol. Sci.* **2010**, *365*, 3227–3246. [CrossRef]

10. White, M.A.; Thornton, P.E.; Running, S.W. A Continental Phenology Model for Monitoring Vegetation Responses to Interannual Climatic Variability. *Glob. Biogeochem. Cycles* **1997**, *11*, 217–234. [CrossRef]

11. Parmesan, C.; Hanley, M.E. Plants and Climate Change: Complexities and Surprises. *Ann. Bot.* **2015**, *116*, 849–864. [CrossRef]

12. Walther, G.-R. Community and Ecosystem Responses to Recent Climate Change. *Philos. Trans. R. Soc. B Biol. Sci.* **2010**, *365*, 2019–2024. [CrossRef] [PubMed]

13. Nemani, R.R.; Keeling, C.D.; Hashimoto, H.; Jolly, W.M.; Piper, S.C.; Tucker, C.J.; Myneni, R.B.; Running, S.W. Climate-Driven Increases in Global Terrestrial Net Primary Production from 1982 to 1999. *Science* **2003**, *300*, 1560–1563. [CrossRef] [PubMed]

14. Zhao, M.; Running, S.W. Drought-Induced Reduction in Global Terrestrial Net Primary Production from 2000 Through 2009. *Science* **2010**, *329*, 940–943. [CrossRef] [PubMed]

15. Donohue, R.J.; Roderick, M.L.; McVicar, T.R. On the Importance of Including Vegetation Dynamics in Budyko's Hydrological Model. *Hydrol. Earth Syst. Sci.* **2007**, *11*, 983–995. [CrossRef]

16. Chuine, I. Why Does Phenology Drive Species Distribution? *Philos. Trans. R. Soc. B Biol. Sci.* **2010**, *365*, 3149–3160. [CrossRef]

17. Chuine, I.; Beaubien, E.G. Phenology Is a Major Determinant of Tree Species Range. *Ecol. Lett.* **2001**, *4*, 500–510. [CrossRef]

18. Forrest, J.; Miller-Rushing, A.J. Toward a Synthetic Understanding of the Role of Phenology in Ecology and Evolution. *Philos. Trans. R. Soc. B Biol. Sci.* **2010**, *365*, 3101–3112. [CrossRef]

19. Chmielewski, F.-M. Phenology and Agriculture. In *Phenology: An Integrative Environmental Science*; Schwartz, M.D., Ed.; Springer: Dordrecht, The Netherlands, 2003; pp. 505–522, ISBN 978-94-007-0632-3.

20. Chmielewski, F.-M.; Köhn, W. Impact of Weather on Yield Components of Winter Rye over 30 Years. *Agric. For. Meteorol.* **2000**, *102*, 253–261. [CrossRef]

21. Bolton, D.K.; Friedl, M.A. Forecasting Crop Yield Using Remotely Sensed Vegetation Indices and Crop Phenology Metrics. *Agric. For. Meteorol.* **2013**, *173*, 74–84. [CrossRef]

22. Tao, F.; Yokozawa, M.; Xu, Y.; Hayashi, Y.; Zhang, Z. Climate Changes and Trends in Phenology and Yields of Field Crops in China, 1981–2000. *Agric. For. Meteorol.* **2006**, *138*, 82–92. [CrossRef]

23. Meier, U.; Bleiholder, H.; Buhr, L.; Feller, C.; Hack, H.; Heß, M.; Lancashire, P.D.; Schnock, U.; Stauß, R.; Van Den Boom, T.; et al. The BBCH System to Coding the Phenological Growth Stages of Plants—History and Publications. *J. Cultiv. Plants* **2009**, *61*, 41–52. [CrossRef]
24. Nekovar, J.; Koch, E.; Kubin, E.; Nejedlik, P.; Sparks, T.; Wielgolaski, F.E. *The History and Current Status of Plant Phenology in Europe*; COST: Brussels, Belgium, 2008; ISBN 978-951-40-2091-9.

25. Sakurai, R.; Jacobson, S.K.; Kobori, H.; Primack, R.; Oka, K.; Komatsu, N.; Machida, R. Culture and Climate Change: Japanese Cherry Blossom Festivals and Stakeholders' Knowledge and Attitudes about Global Climate Change. *Biol. Conserv.* **2011**, *144*, 654–658. [CrossRef]

26. Beaubien, E.; Hamann, A. Spring Flowering Response to Climate Change between 1936 and 2006 in Alberta, Canada. *BioScience* **2011**, *61*, 514–524. [CrossRef]

27. von Linné, C.; von Linné, C. *Philosophia Botanica: In qua Explicantur Fundamenta Botanica Cum Definitionibus Partium, Exemplis Terminorum, Observationibus Rariorum, Adjectis Figuris Aeneis*; Stockholmæ, apud Godofr. Kiesewetter: Stockholm, Sweden, 1751; pp. 1–396. Available online: https://www.biodiversitylibrary.org/item/84231#page/4/mode/1up (accessed on 1 October 2024).

28. Menzel, A. Trends in Phenological Phases in Europe Between 1951 and 1996. *Int. J. Biometeorol.* **2000**, *44*, 76–81. [CrossRef]

29. Jenkins, J.P.; Richardson, A.D.; Braswell, B.H.; Ollinger, S.V.; Hollinger, D.Y.; Smith, M.-L. Refining Light-Use Efficiency Calculations for a Deciduous Forest Canopy Using Simultaneous Tower-Based Carbon Flux and Radiometric Measurements. *Agric. For. Meteorol.* **2007**, *143*, 64–79. [CrossRef]

30. Richardson, A.D.; Braswell, B.H.; Hollinger, D.Y.; Jenkins, J.P.; Ollinger, S.V. Near-Surface Remote Sensing of Spatial and Temporal Variation in Canopy Phenology. *Ecol. Appl.* **2009**, *19*, 1417–1428. [CrossRef]

31. Soudani, K.; Hmimina, G.; Delpierre, N.; Pontailler, J.-Y.; Aubinet, M.; Bonal, D.; Caquet, B.; de Grandcourt, A.; Burban, B.; Flechard, C.; et al. Ground-Based Network of NDVI Measurements for Tracking Temporal Dynamics of Canopy Structure and Vegetation Phenology in Different Biomes. *Remote Sens. Environ.* **2012**, *123*, 234–245. [CrossRef]

32. Delbart, N.; Beaubien, E.; Kergoat, L.; Le Toan, T. Comparing Land Surface Phenology with Leafing and Flowering Observations from the PlantWatch Citizen Network. *Remote Sens. Environ.* **2015**, *160*, 273–280. [CrossRef]

33. Richardson, A.D. PhenoCam: An Evolving, Open-Source Tool to Study the Temporal and Spatial Variability of Ecosystem-Scale Phenology. *Agric. For. Meteorol.* **2023**, *342*, 109751. [CrossRef]

34. Delbart, N.; Kergoat, L.; Le Toan, T.; Lhermitte, J.; Picard, G. Determination of Phenological Dates in Boreal Regions Using Normalized Difference Water Index. *Remote Sens. Environ.* **2005**, *97*, 26–38. [CrossRef]

35. Delbart, N.; Le Toan, T.; Kergoat, L.; Fedotova, V. Remote Sensing of Spring Phenology in Boreal Regions: A Free of Snow-Effect Method Using NOAA-AVHRR and SPOT-VGT Data (1982–2004). *Remote Sens. Environ.* **2006**, *101*, 52–62. [CrossRef]

36. Maignan, F.; Bréon, F.-M.; Bacour, C.; Demarty, J.; Poirson, A. Interannual Vegetation Phenology Estimates from Global AVHRR Measurements: Comparison with In Situ Data and Applications. *Remote Sens. Environ.* **2008**, *112*, 496–505. [CrossRef]

37. Reed, B.C. Trend Analysis of Time-Series Phenology of North America Derived from Satellite Data. *GISci. Remote Sens.* **2006**, *43*, 24–38. [CrossRef]

38. Zhang, X.; Friedl, M.A.; Schaaf, C.B.; Strahler, A.H.; Hodges, J.C.F.; Gao, F.; Reed, B.C.; Huete, A. Monitoring Vegetation Phenology Using MODIS. *Remote Sens. Environ.* **2003**, *84*, 471–475. [CrossRef]

39. Zhang, X.; Friedl, M.A.; Schaaf, C.B. Global Vegetation Phenology from Moderate Resolution Imaging Spectroradiometer (MODIS): Evaluation of Global Patterns and Comparison with in Situ Measurements. *J. Geophys. Res. Biogeosci.* **2006**, *111*, G04017. [CrossRef]

40. de Beurs, K.M.; Henebry, G.M. Land Surface Phenology, Climatic Variation, and Institutional Change: Analyzing Agricultural Land Cover Change in Kazakhstan. *Remote Sens. Environ.* **2004**, *89*, 497–509. [CrossRef]

41. Helman, D. Land Surface Phenology: What Do We Really 'See' from Space? *Sci. Total Environ.* **2018**, *618*, 665–673. [CrossRef]

42. Justice, C.O.; Townshend, J.R.G.; Holben, B.N.; Tucker, C.J. Analysis of the Phenology of Global Vegetation Using Meteorological Satellite Data. *Int. J. Remote Sens.* **1985**, *6*, 1271–1318. [CrossRef]

43. Maisongrande, P.; Duchemin, B.; Dedieu, G. VEGETATION/SPOT: An Operational Mission for the Earth Monitoring; Presentation of New Standard Products. *Int. J. Remote Sens.* **2004**, *25*, 9–14. [CrossRef]

44. Justice, C.O.; Vermote, E.; Townshend, J.R.G.; Defries, R.; Roy, D.P.; Hall, D.K.; Salomonson, V.V.; Privette, J.L.; Riggs, G.; Strahler, A.; et al. The Moderate Resolution Imaging Spectroradiometer (MODIS): Land Remote Sensing for Global Change Research. *IEEE Trans. Geosci. Remote Sens.* **1998**, *36*, 1228–1249. [CrossRef]

45. Sterckx, S.; Benhadj, I.; Duhoux, G.; Livens, S.; Dierckx, W.; Goor, E.; Adriaensen, S.; Heyns, W.; Van Hoof, K.; Strackx, G.; et al. The PROBA-V Mission: Image Processing and Calibration. *Int. J. Remote Sens.* **2014**, *35*, 2565–2588. [CrossRef]

46. Caparros-Santiago, J.A.; Rodriguez-Galiano, V.; Dash, J. Land Surface Phenology as Indicator of Global Terrestrial Ecosystem Dynamics: A Systematic Review. *ISPRS J. Photogramm. Remote Sens.* **2021**, *171*, 330–347. [CrossRef]

47. Pettorelli, N.; Vik, J.O.; Mysterud, A.; Gaillard, J.-M.; Tucker, C.J.; Stenseth, N.C. Using the Satellite-Derived NDVI to Assess Ecological Responses to Environmental Change. *Trends Ecol. Evol.* **2005**, *20*, 503–510. [CrossRef] [PubMed]

48. Shabanov, N.V.; Zhou, L.; Knyazikhin, Y.; Myneni, R.B.; Tucker, C.J. Analysis of Interannual Changes in Northern Vegetation Activity Observed in AVHRR Data from 1981 to 1994. *IEEE Trans. Geosci. Remote Sens.* **2002**, *40*, 115–130. [CrossRef]

49. Picard, G.; Quegan, S.; Delbart, N.; Lomas, M.R.; Le Toan, T.; Woodward, F.I. Bud-Burst Modelling in Siberia and Its Impact on Quantifying the Carbon Budget. *Glob. Change Biol.* **2005**, *11*, 2164–2176. [CrossRef]

50. Tucker, C.J.; Fung, I.Y.; Keeling, C.D.; Gammon, R.H. Relationship between Atmospheric CO2 Variations and a Satellite-Derived Vegetation Index. *Nature* **1986**, *319*, 195–199. [CrossRef]
51. Atzberger, C.; Klisch, A.; Mattiuzzi, M.; Vuolo, F. Phenological Metrics Derived over the European Continent from NDVI3g Data and MODIS Time Series. *Remote Sens.* **2014**, *6*, 257–284. [CrossRef]

52. Rodriguez-Galiano, V.F.; Dash, J.; Atkinson, P.M. Characterising the Land Surface Phenology of Europe Using Decadal MERIS Data. *Remote Sens.* **2015**, *7*, 9390–9409. [CrossRef]

53. Zhang, X.; Liu, L.; Liu, Y.; Jayavelu, S.; Wang, J.; Moon, M.; Henebry, G.M.; Friedl, M.A.; Schaaf, C.B. Generation and Evaluation of the VIIRS Land Surface Phenology Product. *Remote Sens. Environ.* **2018**, *216*, 212–229. [CrossRef]

54. Adole, T.; Dash, J.; Atkinson, P.M. Characterising the Land Surface Phenology of Africa Using 500 m MODIS EVI. *Appl. Geogr.* **2018**, *90*, 187–199. [CrossRef]

55. Berra, E.F.; Gaulton, R. Remote Sensing of Temperate and Boreal Forest Phenology: A Review of Progress, Challenges and Opportunities in the Intercomparison of in-Situ and Satellite Phenological Metrics. *For. Ecol. Manag.* **2021**, *480*, 118663. [CrossRef]

56. Sakamoto, T.; Yokozawa, M.; Toritani, H.; Shibayama, M.; Ishitsuka, N.; Ohno, H. A Crop Phenology Detection Method Using Time-Series MODIS Data. *Remote Sens. Environ.* **2005**, *96*, 366–374. [CrossRef]

57. Chen, X.; Wang, D.; Chen, J.; Wang, C.; Shen, M. The Mixed Pixel Effect in Land Surface Phenology: A Simulation Study. *Remote Sens. Environ.* **2018**, *211*, 338–344. [CrossRef]

58. Kuchler, P.C.; Bégué, A.; Simões, M.; Gaetano, R.; Arvor, D.; Ferraz, R.P.D. Assessing the Optimal Preprocessing Steps of MODIS Time Series to Map Cropping Systems in Mato Grosso, Brazil. *Int. J. Appl. Earth Obs. Geoinf.* **2020**, *92*, 102150. [CrossRef]

59. Sakamoto, T.; Van Nguyen, N.; Ohno, H.; Ishitsuka, N.; Yokozawa, M. Spatio–Temporal Distribution of Rice Phenology and Cropping Systems in the Mekong Delta with Special Reference to the Seasonal Water Flow of the Mekong and Bassac Rivers. *Remote Sens. Environ.* **2006**, *100*, 1–16. [CrossRef]

60. Claverie, M.; Ju, J.; Masek, J.G.; Dungan, J.L.; Vermote, E.F.; Roger, J.-C.; Skakun, S.V.; Justice, C. The Harmonized Landsat and Sentinel-2 Surface Reflectance Data Set. *Remote Sens. Environ.* **2018**, *219*, 145–161. [CrossRef]

61. Kowalski, K.; Senf, C.; Hostert, P.; Pflugmacher, D. Characterizing Spring Phenology of Temperate Broadleaf Forests Using Landsat and Sentinel-2 Time Series. *Int. J. Appl. Earth Obs. Geoinf.* **2020**, *92*, 102172. [CrossRef]

62. Soudani, K.; Delpierre, N.; Berveiller, D.; Hmimina, G.; Vincent, G.; Morfin, A.; Dufrêne, É. Potential of C-Band Synthetic Aperture Radar Sentinel-1 Time-Series for the Monitoring of Phenological Cycles in a Deciduous Forest. *Int. J. Appl. Earth Obs. Geoinf.* **2021**, *104*, 102505. [CrossRef]

63. Vrieling, A.; Meroni, M.; Darvishzadeh, R.; Skidmore, A.K.; Wang, T.; Zurita-Milla, R.; Oosterbeek, K.; O'Connor, B.; Paganini, M. Vegetation Phenology from Sentinel-2 and Field Cameras for a Dutch Barrier Island. *Remote Sens. Environ.* **2018**, *215*, 517–529. [CrossRef]

64. Misra, G.; Cawkwell, F.; Wingler, A. Status of Phenological Research Using Sentinel-2 Data: A Review. *Remote Sens.* **2020**, *12*, 2760. [CrossRef]

65. Boori, M.S.; Choudhary, K.; Paringer, R.; Sharma, A.K.; Kupriyanov, A.; Corgne, S. Monitoring Crop Phenology Using NDVI Time Series from Sentinel 2 Satellite Data. In Proceedings of the 2019 5th International Conference on Frontiers of Signal Processing (ICFSP), Marseille, France, 18–20 September 2019; pp. 62–66.

66. Gao, F.; Anderson, M.; Daughtry, C.; Karnieli, A.; Hively, D.; Kustas, W. A Within-Season Approach for Detecting Early Growth Stages in Corn and Soybean Using High Temporal and Spatial Resolution Imagery. *Remote Sens. Environ.* **2020**, *242*, 111752. [CrossRef]

67. Gao, F.; Anderson, M.C.; Hively, W.D. Detecting Cover Crop End-Of-Season Using VENμS and Sentinel-2 Satellite Imagery. *Remote Sens.* **2020**, *12*, 3524. [CrossRef]

68. Tian, H.; Huang, N.; Niu, Z.; Qin, Y.; Pei, J.; Wang, J. Mapping Winter Crops in China with Multi-Source Satellite Imagery and Phenology-Based Algorithm. *Remote Sens.* **2019**, *11*, 820. [CrossRef]

69. Cai, Y.; Lin, H.; Zhang, M. Mapping Paddy Rice by the Object-Based Random Forest Method Using Time Series Sentinel-1/Sentinel-2 Data. *Adv. Space Res.* **2019**, *64*, 2233–2244. [CrossRef]

70. da Silva Junior, C.A.; Leonel-Junior, A.H.S.; Rossi, F.S.; Correia Filho, W.L.F.; de Barros Santiago, D.; de Oliveira-Júnior, J.F.; Teodoro, P.E.; Lima, M.; Capristo-Silva, G.F. Mapping Soybean Planting Area in Midwest Brazil with Remotely Sensed Images and Phenology-Based Algorithm Using the Google Earth Engine Platform. *Comput. Electron. Agric.* **2020**, *169*, 105194. [CrossRef]

71. Narin, O.G.; Abdikan, S. Monitoring of Phenological Stage and Yield Estimation of Sunflower Plant Using Sentinel-2 Satellite Images. *Geocarto Int.* **2022**, *37*, 1378–1392. [CrossRef]

72. Ayu Purnamasari, R.; Noguchi, R.; Ahamed, T. Land Suitability Assessments for Yield Prediction of Cassava Using Geospatial Fuzzy Expert Systems and Remote Sensing. *Comput. Electron. Agric.* **2019**, *166*, 105018. [CrossRef]

73. Setiyono, T.D.; Quicho, E.D.; Holecz, F.H.; Khan, N.I.; Romuga, G.; Maunahan, A.; Garcia, C.; Rala, A.; Raviz, J.; Collivignarelli, F.; et al. Rice Yield Estimation Using Synthetic Aperture Radar (SAR) and the ORYZA Crop Growth Model: Development and Application of the System in South and South-East Asian Countries. *Int. J. Remote Sens.* **2019**, *40*, 8093–8124. [CrossRef]

74. Waldhoff, G.; Curdt, C.; Hoffmeister, D.; Bareth, G. Analysis of Multitemporal and Multisensor Remote Sensing Data for Crop Rotation Mapping. *ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci.* **2012**, *I-7*, 177–182. [CrossRef]

75. Hua, A.K. Land Use Land Cover Changes in Detection of Water Quality: A Study Based on Remote Sensing and Multivariate Statistics. *J. Environ. Public Health* **2017**, *2017*, 7515130. [CrossRef]

76. Gao, F.; Zhang, X. Mapping Crop Phenology in Near Real-Time Using Satellite Remote Sensing: Challenges and Opportunities. *J. Remote Sens.* **2021**, *2021*, 8379391. [CrossRef]
77. Cheng, Y.; Vrieling, A.; Fava, F.; Meroni, M.; Marshall, M.; Gachoki, S. Phenology of Short Vegetation Cycles in a Kenyan Rangeland from PlanetScope and Sentinel-2. *Remote Sens. Environ.* **2020**, *248*, 112004. [CrossRef]

78. d'Andrimont, R.; Taymans, M.; Lemoine, G.; Ceglar, A.; Yordanov, M.; van der Velde, M. Detecting Flowering Phenology in Oil Seed Rape Parcels with Sentinel-1 and -2 Time Series. *Remote Sens. Environ.* **2020**, *239*, 111660. [CrossRef] [PubMed]

79. Mercier, A.; Betbeder, J.; Baudry, J.; Le Roux, V.; Spicher, F.; Lacoux, J.; Roger, D.; Hubert-Moy, L. Evaluation of Sentinel-1 & 2 Time Series for Predicting Wheat and Rapeseed Phenological Stages. *ISPRS J. Photogramm. Remote Sens.* **2020**, *163*, 231–256. [CrossRef]

80. Meroni, M.; d'Andrimont, R.; Vrieling, A.; Fasbender, D.; Lemoine, G.; Rembold, F.; Seguini, L.; Verhegghen, A. Comparing Land Surface Phenology of Major European Crops as Derived from SAR and Multispectral Data of Sentinel-1 and -2. *Remote Sens. Environ.* **2021**, *253*, 112232. [CrossRef]

81. Rivas, H.; Delbart, N.; Ottlé, C.; Maignan, F.; Vaudour, E. Disaggregated PROBA-V Data Allows Monitoring Individual Crop Phenology at a Higher Observation Frequency than Sentinel-2. *Int. J. Appl. Earth Obs. Geoinf.* **2021**, *104*, 102569. [CrossRef]

82. Gao, F.; Anderson, M.C.; Zhang, X.; Yang, Z.; Alfieri, J.G.; Kustas, W.P.; Mueller, R.; Johnson, D.M.; Prueger, J.H. Toward Mapping Crop Progress at Field Scales through Fusion of Landsat and MODIS Imagery. *Remote Sens. Environ.* **2017**, *188*, 9–25. [CrossRef]

83. Gevaert, C.M.; Garcia-Haro, F.J. A Comparison of STARFM and an Unmixing-Based Algorithm for Landsat and MODIS Data Fusion. *Remote Sens. Environ.* **2015**, *156*, 34–44. [CrossRef]

84. Zhao, Y.; Huang, B.; Song, H. A Robust Adaptive Spatial and Temporal Image Fusion Model for Complex Land Surface Changes. *Remote Sens. Environ.* **2018**, *208*, 42–62. [CrossRef]

85. Lobell, D.B.; Asner, G.P. Cropland Distributions from Temporal Unmixing of MODIS Data. *Remote Sens. Environ.* **2004**, *93*, 412–422. [CrossRef]

86. Haertel, V.F.; Shimabukuro, Y.E. Spectral Linear Mixing Model in Low Spatial Resolution Image Data. *IEEE Trans. Geosci. Remote Sens.* **2005**, *43*, 2555–2562. [CrossRef]

87. Busetto, L.; Meroni, M.; Colombo, R. Combining Medium and Coarse Spatial Resolution Satellite Data to Improve the Estimation of Sub-Pixel NDVI Time Series. *Remote Sens. Environ.* **2008**, *112*, 118–131. [CrossRef]

88. Chen, J.; Jönsson, P.; Tamura, M.; Gu, Z.; Matsushita, B.; Eklundh, L. A Simple Method for Reconstructing a High-Quality NDVI Time-Series Data Set Based on the Savitzky–Golay Filter. *Remote Sens. Environ.* **2004**, *91*, 332–344. [CrossRef]

89. Eilers, P.H.C. A Perfect Smoother. *Anal. Chem.* **2003**, *75*, 3631–3636. [CrossRef] [PubMed]

90. Zeng, L.; Wardlow, B.D.; Xiang, D.; Hu, S.; Li, D. A Review of Vegetation Phenological Metrics Extraction Using Time-Series, Multispectral Satellite Data. *Remote Sens. Environ.* **2020**, *237*, 111511. [CrossRef]

91. Cantelaube, P.; Carles, M. Le registre parcellaire graphique: Des données géographiques pour décrire la couverture du sol agricole. *Le Cah. Tech. L'inra* **2014**, *58–64*.

92. Peel, M.C.; Finlayson, B.L.; McMahon, T.A. Updated World Map of the Köppen-Geiger Climate Classification. *Hydrol. Earth Syst. Sci.* **2007**, *11*, 1633–1644. [CrossRef]

93. Joly, D.; Brossard, T.; Cardot, H.; Cavailhes, J.; Hilal, M.; Wavresky, P. Les types de climats en France, une construction spatiale. *Cybergeo Eur. J. Geogr.* **2010**, *501*. [CrossRef]

94. Wolters, E.; Dierckx, W.; Iordache, M.-D.; Swinnen, E. PROBA-V Products User Manual; VITO: Mol, Belgium, 2018.

95. Settle, J.J.; Drake, N.A. Linear Mixing and the Estimation of Ground Cover Proportions. *Int. J. Remote Sens.* **1993**, *14*, 1159–1177. [CrossRef]

96. Tucker, C.J. Red and Photographic Infrared Linear Combinations for Monitoring Vegetation. *Remote Sens. Environ.* **1979**, *8*, 127–150. [CrossRef]

97. Main-Knorn, M.; Pflug, B.; Louis, J.; Debaecker, V.; Müller-Wilm, U.; Gascon, F. Sen2Cor for Sentinel-2. In Proceedings of the Image and Signal Processing for Remote Sensing XXIII, Warsaw, Poland, 11–14 September 2017; SPIE: Warsaw, Poland, 2017; Volume 10427, pp. 37–48.

98. Lonjou, V.; Desjardins, C.; Hagolle, O.; Petrucci, B.; Tremas, T.; Dejus, M.; Makarau, A.; Auer, S. MACCS-ATCOR Joint Algorithm (MAJA). In Proceedings of the Remote Sensing of Clouds and the Atmosphere XXI, Edinburgh, UK, 26–29 September 2016; SPIE: London, UK, 2016; Volume 10001, pp. 25–37.

99. Zhou, J.; Jia, L.; Menenti, M. Reconstruction of Global MODIS NDVI Time Series: Performance of Harmonic ANalysis of Time Series (HANTS). *Remote Sens. Environ.* **2015**, *163*, 217–228. [CrossRef]

100. Espinoza-Dávalos, G.E.; Bastiaanssen, W.G.M.; Bett, B.; Cai, X. A Python Implementation of the Harmonic ANalysis of Time Series (HANTS) Algorithm for Geospatial Data. *IHE Delft Inst. Water Educ.* **2017**. [CrossRef]

101. Huang, X.; Liu, J.; Zhu, W.; Atzberger, C.; Liu, Q. The Optimal Threshold and Vegetation Index Time Series for Retrieving Crop Phenology Based on a Modified Dynamic Threshold Method. *Remote Sens.* **2019**, *11*, 2725. [CrossRef]

102. White, M.A.; De Beurs, K.M.; Didan, K.; Inouye, D.W.; Richardson, A.D.; Jensen, O.P.; O'keefe, J.; Zhang, G.; Nemani, R.R.; Van Leeuwen, W.J.D.; et al. Intercomparison, Interpretation, and Assessment of Spring Phenology in North America Estimated from Remote Sensing for 1982–2006. *Glob. Change Biol.* **2009**, *15*, 2335–2359. [CrossRef]

103. Boyd, D.S.; Almond, S.; Dash, J.; Curran, P.J.; Hill, R.A. Phenology of Vegetation in Southern England from Envisat MERIS Terrestrial Chlorophyll Index (MTCI) Data. *Int. J. Remote Sens.* **2011**, *32*, 8421–8447. [CrossRef]

104. Zhu, W.; Tian, H.; Xu, X.; Pan, Y.; Chen, G.; Lin, W. Extension of the Growing Season Due to Delayed Autumn over Mid and High Latitudes in North America during 1982–2006. *Glob. Ecol. Biogeogr.* **2012**, *21*, 260–271. [CrossRef]
105. Zhang, X.; Wang, J.; Henebry, G.M.; Gao, F. Development and Evaluation of a New Algorithm for Detecting 30 m Land Surface Phenology from VIIRS and HLS Time Series. *ISPRS J. Photogramm. Remote Sens.* **2020**, *161*, 37–51. [CrossRef]

106. Xu, X.; Conrad, C.; Doktor, D. Optimising Phenological Metrics Extraction for Different Crop Types in Germany Using the Moderate Resolution Imaging Spectrometer (MODIS). *Remote Sens.* **2017**, *9*, 254. [CrossRef]

107. Xin, Q.; Li, J.; Li, Z.; Li, Y.; Zhou, X. Evaluations and Comparisons of Rule-Based and Machine-Learning-Based Methods to Retrieve Satellite-Based Vegetation Phenology Using MODIS and USA National Phenology Network Data. *Int. J. Appl. Earth Obs. Geoinf.* **2020**, *93*, 102189. [CrossRef]

108. Liu, L.; Cao, R.; Chen, J.; Shen, M.; Wang, S.; Zhou, J.; He, B. Detecting Crop Phenology from Vegetation Index Time-Series Data by Improved Shape Model Fitting in Each Phenological Stage. *Remote Sens. Environ.* **2022**, *277*, 113060. [CrossRef]

109. Cantelaube, P.; Lardot, B. RPG complété 2020 Région Midi-Pyrénées. *Rech. Data Gouv V4 2022*. [CrossRef]

110. Inglada, J.; Vincent, A.; Arias, M.; Tardy, B.; Morin, D.; Rodes, I. Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series. *Remote Sens.* **2017**, *9*, 95. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.