PRELIMINARY CONCERNS ABOUT AGRONOMIC INTERPRETATION OF NDVI TIME SERIES FROM SENTINEL-2 DATA: PHENOLOGY AND THERMAL EFFICIENCY OF WINTER WHEAT IN PIEMONTE (NW ITALY)

TELECER project is supported through Rural Development Programme regional action of EU CAP and is aimed at providing Precision Agriculture–devoted services for cereals monitoring in the Piemonte Region (NW-Italy) context. In this work authors explored some general and preliminary issues mainly aimed at demonstrating and formalizing those evident relationships existing between NDVI image time series and the main ordinary agronomic parameters, with special focus on phenology and thermal efficiency of crops as related to Growing Degrees Day (GDD). Winter wheat was investigated and relationships calibrated at field level, making possible to spatially characterise environmental and management effects. Two different analysis were achieved: (i) one aimed at mapping crop phenological metrics, as derivable from NDVI S2 time series; (ii) one aimed at locally modelling relationship linking GDD and NDVI to somehow test the thermal efficiency of crops in the different parts of the study area. The first analysis showed that the end of season appears to be the most constant phenological metric in the study area possibly demonstrating a time concentration of harvest operations in the area. Differently, the peak of season and the start of season metrics showed to be largely varying in the study, thus suggesting to be stronger predictors: (i) of crop development; (ii) of the effects induced by local agronomical practices. Several base temperatures were used to compute correspondent GDD. These were tested against NDVI and modelled by a parabolic model at field level. Model coefficients distribution were analysed and mapped the correspondent agronomic interpretation suggested.


Sentinel-2 Data in Agriculture
For the last years, Earth Observation (EO) data have been supporting many human activities including agriculture. In this context, the Sentinel-2 (S2) mission of the European Copernicus programme, providing high-resolution optical imagery for land monitoring purposes (Kuc and Chormański, 2019;Phiri et al., 2020), plays a key role. Starting from the S2 imagery several spectral vegetation index (VI) can be computed able to retrieve vegetation-related information at the proper spatial scale. The Normalized Difference Vegetation Index (NDVI), that combine RED (S2 -band 4) and NIR (S2 -band 8) bands, is widely used in agriculture (Huang et al., 2021). In particular, NDVI-based applications in agriculture include: crop phenology monitoring (Boori et al., 2019;Misra et al., 2020;Testa et al., 2018); crop yield estimation (Parida et al., 2021;Rahman and Robson, 2020;Zhao et al., 2020); crop nitrogen content estimation (Clevers and Gitelson, 2013;Sharifi, 2020); disaster events monitoring (Caballero et al., 2019;De Petris et al., 2021); risk management (Bacchini and Miguez, 2015); insurance activities for damage estimation Shrestha et al., 2017); CAP controls in agriculture (Kanjir et al., 2018;Sarvia et al., 2021c). Additionally, S2 VIs have been applied in precision agriculture (PA) (Epinat et al., 2001;Ma et al., 2006;Sarvia et al., 2021b, p. 2). It is worth to remind that PA is a business management strategy that uses new technologies (sensors, processing software, machinery actuators) to support farmers decisions. In particular, PA adopts remotely-sensed data for many applications like pre-growth soil fertility and moisture analyses, crop growth and yield forecasting (Sishodia et al., 2020). In this framework, S2 data well fit PA requirements due to its spectral (13 bands in the range 350-2500 nm), temporal (nominally 5 days) and geometric (10 m GSD, Ground Sampling Distance) resolutions (Segarra et al., 2020). Besides spectral approaches, several other methods can be used to retrieve vegetation-related information. Specifically, thermal data could be useful to explore physiological plant response. In fact, plants require a specific amount of heat to develop; consequently, one of the main driving phenological indicators can be found in the heat accumulated by vegetation during its growing season. The amount of heat energy that an organism accumulates can be expressed through the Growing Degree Day (GDD) (Durand et al., 1967). Thanks to the increasing adoption of meteorological sensors (temperature, rain, etc.), GDDs are currently and widely used in agrometeorological models to describe plant phenology (Ghamghami et al., 2019). Currently, several services and platforms have been developed by private companies combining S2 data with meteorological and ground data aiming at supporting the agricultural context. These platforms are expected to simplify and optimize farmers' work and to improve environmental and economic sustainability of their activities. Specifically, services concern: a) GIS-based field mapping, b) satellite-based crop monitoring, c) generation and delivery of crop disease prevention models, d) production and delivery of decision support systems; e) weather data analysis. The most of them are supplied by the Internet like OneSoil, Agrivi, FarmsLogs, Agrimap Crop Monitoring, Agricolus, Arvatec, FarmB.sat. Unfortunately, these services are not free of charge and the information they provide is uncertified. Moreover, they operate globally, with no neglecting to adopt specifications related to local agronomic conditions and crop calendars. According to Aletdinova (Aletdinova and Lenski, 2021) main existing platforms releasing services, do not consider optimization models thus limiting farmers in obtaining possible solutions/suggestions based on real-time data. Moreover, these services suffer from two main limitations: (i) phenological metrics (PMs) are often neglected in data processing, when, differently, they can significantly improve data interpretation; in fact, PMs are able to summarize complex information about crop conditions, highlighting local anomalies that can be useful for crop management; (ii) data are processed through pre-trained models whose capability of generalization is hard to be demonstrated given the high variability of local conditions that can affect crop developments. These limitations raise some doubts about the general applicability of such services in heterogeneous agricultural contexts where spatial and temporal variability of meteorological, soil and agronomical conditions are significant. Consequently, they are difficult to be easily accepted by farmers. Given these premises, the TELECER (hereinafter called TELE) project titled "Productive and qualitative improvement of cereals cultivation within an advanced supply chain based on remote sensing" was launched in 2020 January. Within this project, Farmers' Consortia (Piemonte Agricultural Consortium for Agri-Food and Cereals -CAPAC), Universities (Department of Agricultural, Forest and Food Sciences -DISAFA) and software companies were involved in developing and providing a shared, official and certified service for cereals (and maize) detection and monitoring. TELE is founded through the Regione Piemonte -PSR 2014-2020, 16.1.1 Action 2 (namely 2014IT06RDRP009: Italy -Rural Development Programme). TELE is in charge of developing a service combining satellite, meteorological and field data for crop monitoring (wheat, barley, corn and soybean) easy to be accessed by farmers. Information from the service is intended for addressing field agro-technical management practices like irrigation, fertilization and defence, with special focus on environmental conditions of the Piemonte region (NW Italy). Several local farmers were included in the project for a major sharing of problems, requirements and technical limitations that should guarantee a better compliance of the final service with robustness, reliability, effectiveness, simplicity and usability criteria.

Goals
This preliminary study is aimed at testing applicability of S2 time series for phenology and crop thermal efficiency description of winter-wheat (WW) in Piemonte region (NW-Italy) trying to figure out specific spatial patterns that could be related to environmental and management factors. In particular, starting from S2 and meteorological data, local models and phenological metrics were derived and calibrated for WW. Two different analyses were performed: (i) one aimed at mapping crop phenological metrics, as derivable from NDVI S2 time series; (ii) one aimed at locally modelling the well known relationship linking growing degree days (GDD) and NDVI (Borgogno-Mondino et al., 2020) somehow testing thermal efficiency of crops in the different parts of the study area.

Study Area
Efficiency and robustness of crop monitoring by remote sensing strictly depends on the size of the observed field with respect to the image GSD (Meier et al., 2020;Vajsová et al., 2020). Consequently, only WW fields having a size > 0.1 ha (at least 10 S2 pixels @10 m GSD) were selected. A total of 47 WW fields were provided as a vector map by CAPAC local consortium (TELE partner) and used as ground reference data where crop and management type were known. Reference fields are located in the flat area of Piemonte, spreading across the entire region (Fig. 1). Piemonte is highly devoted to cereal production, producing about 11% of Italian wheat (Borri and Trione, 2020). Since Italy shows a high climatic variability depending on latitude, crop phenology management calendars change significantly region by region. For this reason, the Piemonte agronomic calendar was acquired (Sarvia et al., 2021c) showing that WW growing season starts in late October (sowing) and concludes on early July (harvesting) having its maximum phenological expression between April and May. Reference fields are represented by blue polygons. Orange/yellow/red squares border S2 tiles (Reference System: WGS84/UTM 32 N).

Sentinel-2 Data
In this work 22 S2 images were obtained ranging from 2 nd November 2020 to 20 th July 2021 ( Figure 2). Images were obtained as Level 2A products (orthoprojected and at-theground reflectance calibrated) from the Scientific Open Data Hub (https://scihub.copernicus.eu/). S2 data are provided as tiles of 100 x 100 km 2 . According to WW spatial distribution, 3 tiles were needed, namely T32TMR, T32TMQ and T32TLQ (Fig 1) to cover the entire area. It is worth to remind that S2 Level 2A data are supplied together with a scene classification layer (SCL) alerting about cloudy, shadowy and fault pixels that is often used to detect and remove bad observations from the processed image time series. S2 technical features are reported in Table 1.

Meteorological data
Since one of the goals of this work was to test and describe thermal efficiency of WW by relating GDD (eq. 3) with NDVI, the daily average air temperature was needed. Consequently, meteorological data were obtained from 11 meteorological stations (figure 3) belonging to the regional meteorological network (www.arpa.piemonte.it, last access 15 th October 2021) located close to WW fields. Data were obtained for the period 2 nd November 2020 -30 th June 2021 with daily frequency for a total of 280 temperature values.

NDVI Time Series
According to its well-known capability of describing crop phenology and in order to guarantee the highest geometric resolution, in this work NDVI was considered. Consequently, starting from the native S2 tiles the correspondent image stacks of band 4 (red), band 8 (NIR) and SCL were generated. A selfdeveloped routine was implemented in IDL v4.8 (Harris Geospatial Solutions, Inc., Broomfield, CO, USA) (Chen et al., 2004) to generate NDVI times series for the reference period according to the following steps operating at pixel level: (i) bad observations were masked out according to SCL classification code; (ii) local NDVI temporal profile was interpolated by spline (tensor = 5) to obtain a regularized NDVI stack (5 day time resolution) made of 52 interpolated NDVI maps (for each processed tile). Finally, resulting NDVI stacks were mosaicked and the WW vector layer used to compute by zonal statistics the average NDVI temporal profile for WW fields.

Phenological Metrics Computation and Mapping
PM are useful to detect key moments along the growing season of vegetation (Misra et al., 2020) able to synthesize the entire cycle. For this work the following PM were computed at field level with reference to the local average NDVI temporal profile: start of season (SOS, as DOY = Day of the Year); NDVI peak value (NDVIM) and the correspondent date (POS, as DOY); end of season (EOS, as DOY), length of season (LOS, as DOY difference), growing rate (GR, as ratio) and senescence rate (SR, as ratio). NDVIM corresponds to the maximum NDVI value reached by the crop during the considered period as derivable from the average and regularized NDVI temporal profile of the considered field; POS correspond to the DOY when NDVIM is reached; SOS and EOS were computed through a thresholdbased approach looking for the DOY when the local NDVI value reaches 0.5 before and after the NDVI peak respectively (Misra et al., 2020); LOS was computed as the difference between EOS and SOS informing about the length of the phenological season. Concerning GR and SR, they were computed according to Eq. 1 and 2.
where is the maximum NDVI value along the growing season; is the NDVI value at SOS; is the NDVI value at EOS. Some statistics were finally computed to synthesize phenological behaviour of monitored fields.

Computing GDD
Plant growth is highly related to air temperature conditions. Growing Degree Days (GDD) is a metric that is known to describe biomass growth during the phenological season accounting for the temperature effects on crops, somehow related to its thermal efficiency. In this work daily GDD was computed according to eq. 3 (Mcmaster, 1997). Daily temperature values were obtained from the meteorological station closest to the considered field and used to compute GDD.
where is the daily average temperature for the i-th observation; is the nominal base temperature (crop dependent) defining the minimum temperature value activating the phenological development of vegetation (Salazar-Gutierrez et al., 2013).
is a very complex factor depending on both the plant species and the cultivar (Salazar-The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France Gutierrez et al., 2013). Moreover, while working with winter cereals is known to vary along the phenological stages. For winter wheat ranges from 0°C to 5°C (Boogar et al., 2013;Fraisse and Paula-Moraes, 2018;Salazar-Gutierrez et al., 2013;Shaykewich, 1995), making this choice not unique. Consequently, in order to assess the most suitable TBASE for WW, three TBASE were considered: (a) TBASE = 0; (b) TBASE = 5; (c) TBASE = 0 before heading while TBASE = 5 after heading until harvest. Thus, with reference to the available meteorological data, three different GDD profiles (namely GDD_0, GDD_5, GDD_05) were generated for each WW field corresponding to the three above mentioned TBASE values. Depending on the different TBASE, GDD profiles of crops show different growing rates. Specifically, lower TBASE correspond to crops that grow along the winter season; higher TBASE reflect a phenological stop in winter time and a slower overall growth.

NDVI vs GDD
Literature refers about a high correlation between GDD and NDVI (de Beurs and Henebry, 2004) being the former related to the cause and the latter to the effect of vegetation growing. In particular, a parabolic model demonstrated to well fit this relationship (eq. 4), assuming NDVI as predicted variable and the correspondent GDD as predictor one.  (Gonçalves and White, 2005)) and the mean absolute errors (MAE, (Willmott and Matsuura, 2005)).

Subsequently, a new index called
was computed (eq.5) in order to better describe the phenological peak as predicted by temperature data.

= −
where is the maximum NDVIG as predicted by the quadratic model. Finally, the differences between NDVIGMAX and NDVIM were calculated at field level.

Phenological Metrics
Once PM were computed at field level, statistical distribution ( Fig. 4) of SOS, EOS, POS, LOS and NDVIM from WW fields were analysed in order to qualify and characterize fields within AOI. Based on boxplots following considerations can be carried out: (a) average WW SOS is placed at 6 th observation with a fairly high variability, which means that SOS starts around at 27 th November fitting expected local agronomical calendar (Sarvia et al., 2021b); (b) average WW EOS is placed at 46 th observation with a low variability, which means that EOS ends at 15 th June in almost all cases; (c) average WW LOS is equal to 39 observations with a fairly high variability, which means that LOS persist around 195 days. (d) average WW POS is placed at 35 th observation with a fairly high variability, which means that POS occurs around 21 st April. (e) average WW NDVIM is equal to 0.83 with a fairly high variability. These results appear to be consistent and supported by the local climatic, environmental and agricultural conditions. SOS high variability can be interpreted as related to farmers' technical choices and, in particular, with the sowing date of each field. Conversely, EOS low variability can be associated to the ripening and maturity stages that usually occur at similar dates (June) in AOI. Consequently, LOS high variability can be associated mostly to the one that SOS suffers from. The high variability of NDVIM and POS are probably due to WW different agricultural, meteorological and environmental local conditions that can positively, or negatively, affect crop development. It is worth to remind that the most relevant drivers that could influence NDVI (and therefor crop development) are: soil type, nutrients, pathologies, water availability and terrain aspect (Neenu et al., 2013). As far as GR and SR are concerned, the correspondent statistical distributions were computed and summarized in the boxplots of figure 5. Averagely, GR was found to be 0.12 with very low variability; this means that WW growth tends to be quite similar in the monitored fields. This result is consistent with winter crop phenological development which is characterized by slower growth (GR), probably due to winter low temperatures within AOI. Conversely, the average SR value was found to be -0.04 with fairly high variability. This result can still be related to crop reaction to environmental conditions, leading to significant variations within crop development and maturity.

GDD Computation
As previously stated, GDD profiles strongly depend from TBASE. Figure 6 shows the average GDD profiles corresponding to the three considered TBASE values, as computed with reference to the whole WW dataset. Figure 6. GDD profiles (GDD_5, GDD_0, GDD_05) as computed with reference to the three considered TBASE values. NDVI and GDD profiles are the average ones from the whole WW dataset.
Looking at figure 6 it can be noted that GDD_0 profile continues its growing even in winter time without long pause periods. Conversely, GDD_5 profile appears to stop its growing from December to early March. As far as GDD_05 profile is concerned it looks similar to the GDD_0 one in its early phase (until April); after, it slows down assuming a slope similar to the GDD_5 one. These differences highlight that a high degree of uncertainty is present depending on the selected model. A more objective way for selecting the proper model is therefore mandatory. The possibility of testing NDVI vs GDD correlations is certainly one of the possible tool for obtaining the answer to this question (see next section). Furtherly looking at figure 6, on 30 th June (last observation) the following thermal sums can be found: 1294, 2309 and 1926 for GDD_5, GDD_0 and GDD_05, respectively. These values appear to be consistent with the ones from literature. Both GDD_0 and GDD_5 show final values similar to those reported in (Salazar-Gutierrez et al., 2013) for Georgia, USA and in (Gill et al., 2014) for the Punjab region, India. Concerning GDD_05, no study was found in literature to support or contradict our results.

NDVI as a function of GDDC
To find which TBASE was more suitable for the AOI fields, the three GDD profiles were tested against NDVI. According to eq. 4, a parabolic model was calibrated for each WW field (47) and for all the three considered GDD profiles. A total of 47 × 3=141 models were calibrated. Correspondent SEE were computed and separately analysed for the three GDD profiles (Fig 7). Figure 7. SEE statistics as computed from the calibrated parabolic models relating local NDVI and GDD. They refer to the GDD_5, GDD_0 and GDD_05 profiles. Bottom-up: 5 th , 25 th , 50 th , 75 th and 95 th percentiles. Cross corresponds to the mean value. Figure 7 shows that the lowest SEE mean value (0.043) corresponds to the model that uses the GDD_05 profile as NDVI predictor. A similar performance in terms of SEE mean value was achieved with reference to the GDD_0 profile (0.047). Differently, GDD_5 proved to be a worse predictor determining a SEE mean value of 0.079. As far as representativeness of SEE mean value is concerned, it was possible to note that the one associated to GDD_0 was the higher one (lower standard deviation = 0.009855) even if weakly better of the one from GDD_05 (standard deviation = 0.010056); additionally, GDD_0 SEE distribution appeared the be less skewed (skewness = -0.148) than the one from GDD_05 (Skewness = 0.309), having a favourable skewness compacting higher SEE values against its mean value. On the other side, unreliability of NDVI prediction by GDD_5 is confirmed by the correspondent highest standard deviation (0.023). Taking care of these evidences, it can be stated that, in AOI, a TBASE of 5 °C has to be absolutely, avoided. Some doubts persist about the performances given by GDD_0 an GDD_05 profiles, that globally appear to be very similar. Nevertheless, a preliminary evaluation based on the above mentioned criteria, addressed us to consider GDD_0 as the best choice, in spite of the slightly higher SEE mean value. A final and decisive selection could however come by comparing the performance (Mean Absolute Error -MAE, figure 8) of the GDD_05 and GDD_0 in the last part of the growing season where the two profiles diverge (16 th April -30 th June). Figure 8. MAE of the GDD_0 and GDD_05 models (MAE_0 and MAE_05 respectively) in the period 6 th April -30 th June. Figure 8 showed that GDD_0 model leads to lower mean and standard deviation MAE values (0.052 and 0.015 respectively) than GDD_05 (0.056 and 0.021 respectively), denoting a slight difference between the two models; consequently, TBASE=0 was chosen as the reference one for this work. Once the optimal TBASE was defined, GDD_0-based model coefficients distributions (a, b and c) were analysed (Fig. 9). It is worth to remind that: c coefficient defines the model offset; b is mainly related to function maximum/minimum and slope in the variable space; a majorly conditions both slope and concavity (higher absolute values correspond to a narrower concavity). As previously reported in Eq. (5), a and b concur in determining the model maximum/minimum. Figure 9. Boxplots representing the statistical distribution of the three GDD_0-based model coefficients (a, b, c) and the NDVIGMAX predicted by the model. From bottom-up: 5 th , 25 th , 50 th , 75 th and 95 th percentiles. Cross is the mean value. Dots are outliers. All coefficients are multiplied by 10000. Figure 9 shows that mean values were 3500, 8.7 and -0.004 for c, b and a respectively; it is worth to remind that all coefficients are multiplied by 10000.
In particular, c value is consistent with the NDVI value corresponding to the begin of WW development (0.3); b values result to be quite low due to long growth period and low winter temperatures; a values result to be always negative confirming that NDVI decreases after reaching its peak (Sarvia et al., 2021a). Concerning the possible agronomic meaning of model coefficients, the following interpretation can be given: (i) c, being related to the initial NDVI values could depend on several agronomical drivers like crop seeding and emergence time, local pedological features and genetic differences of cultivars; (ii) a and b, defining the increase of NDVI as a function of GDD, could be used to summarize the strength of WW development and, consequently, its thermal efficiency.
Additionally, a refers about length of vegetation activity before harvesting. Summarizing, it can be stated that the c coefficient is completely related to the early stages of development and in particular, mainly influenced by the sowing period. Conversely, a and b seem to be related to the growing season and in particular to the growing rate. To give an operational demonstration of the generated information a small focus area (FA) was chosen and mapped. In particular the spatial variability of coefficients ( Figure 10), NDVIGMAX, NDVIM and their difference ( Figure 11) were mapped. Distribution of coefficients showed a great variability within FA, suggesting that they can be adopted for WW characterization according to the previously provided interpretation keys. For example, concerning the c coefficient of figure 10a, it can be noted that some fields showed values lower than 1000 suggesting that in November no vegetative activity was on and, therefore, that WW was not still in its emergence phase. Conversely, other fields showed values greater than 5000, proving that a weak vegetative activity was on and, therefore, that WW had already started its emergence phase. Concerning other coefficients, in FA high b and a absolute values were located where WW showed low c values. This suggest that a rapid growth rate occurred in these fields trying to compensate emergence delays. Additionally, some fields showed intermediate values of b and a suggesting that crop development was more constant. As far as the relationship between PM and GDD-related estimates is concerned, NDVI peak is known to be a good predictor of crop yield (Sultana et al., 2014). Therefore, interesting outcomes can come by comparing NDVIGMAX and NDVIM. Surprisingly, the two peaks showed different values and a different pattern in FA (Fig. 10). It is worth to remind that NDVIGMAX is estimated through the GDD-based model, and it somehow refers about WW expected thermal efficiency. Differently, NDVIM is the NDVI maximum as seen by satellite and therefore it absorbs all environmental and agronomic factors affecting crop growth (soil, cultivar, water availability, agronomic management). Since differences between NDVIGMAX and NDVIM ranged between 0.04 and 0.09 NDVI points (Fig. 11c), they can be retained significant (Borgogno-Mondino et al., 2016)). Consequently, they can be used to somehow recognize anomalies in thermal efficiency of WW that can be related to possible differences in yield. Due to the preliminary nature of this work, future developments will be conducted in order to relate these deductions to the yields and over other crops pertaining TELE (corn, barley, soybean, etc). The possibility of locally map this values, could drive to interesting information useful for a proper and more sustainable management of crops in the area.

CONCLUSIONS
This work was developed within the TELECER project that is aimed at improving and supporting the cereal sector within the Piemonte region (NW Italy), through the adoption of Copernicus Sentinel-2 data. Presented results are preliminary and concern the characterization of winter wheat fields in the area managed by the CAPAC. To characterize winter wheat development, two different procedures were proposed: a) one based on phenological metrics derived by NDVI time series; in particular SOS, EOS, LOS, POS, NDVIM, GR and SR were computed at field level for winter wheat; b) one based on modelling of local relationships between growing degree days and NDVI to investigate local thermal efficiency of WW. The first analysis proved that EOS is the less variable phenological metric in AOI, suggesting a common WW behaviour during its maturity/harvesting. Differently, POS and SOS appeared to be the most variable phenological metrics in AOI, able to characterize WW development and agronomical practices (sowing) respectively. According to the second analysis, several TBASE values were tested demonstrating that a TBASE of 0°C is the most proper in AOI and that a parabolic model well fits the relationship between GDD and NDVI. It was also proved that model coefficients can be used to summarize local thermal efficiency of WW. Some agronomic interpretations were proposed about model coefficients; moreover, since they can be locally computed at field level, they can also be mapped giving new spatial information for a more proper sustainable management of WW in the area. This is a remarkable consideration especially when WW fields are managed/monitored unitarily through agronomic strategies based on Consortia, like CAPAC.