Detection and counting of meadow cuts by copernicus sentinel-2 imagery in the framework of the common agricultural policy (CAP)

ABSTRACT Greening is a subsidy provided by the Common Agricultural Policy (CAP), related to mowing and designed to protect environment. National or regional paying agencies (PP) monitor and verify compliance of farmers’ declarations with CAP rules. In this work, an operational procedure is proposed aimed at supporting PPs in detecting, mapping and quantifying the number of times mowing occurred in a meadow field. In particular, 72,539 meadows fields within the Piemonte region (NW – Italy) were analysed with a time series of Sentinel-2 (S2) data. The procedure is based on the processing of filtered and regularized time series of NDVI maps. The Fast Fourier Transform (FFT) was applied at field level to decompose the local NDVI temporal profile. The frequency ($${\omega _{peak}})$$ωpeak) corresponding to the maximum amplitude ($${A_{peak}})$$Apeak) was therefore considered. $${A_{peak}}$$Apeak value was used to detect not-mowed meadows by thresholding based on the application of the non-parametric Kolmogorov-Smirnov test. Mowing counting were achieved with reference to $${\omega _{peak}}$$ωpeak and the correspondent map (called MCM) generated for the study area. MCM was, finally, tested against the validation set (285 fields). Results showed an Overall Accuracy (OA) > 87%, confirming the effectiveness of the proposed procedure in detecting, mapping and quantifying the number of times mowing occurred.


Introduction
In the last years, agricultural systems are undergoing rapid transformation because of climate change, consumer demands and worldwide and European directives (Primdahl, 2014).Soil degradation processes resulting from inadequate land management are increasing (Lal, 2015).In particular, soil erosion (Lal, 1993;Montgomery, 2007), soil compaction (Tobias et al., 2018), soil consumption, loss of organic matter and structure (X.Liu et al., 2006), salinization and desertification (Santa Olalla, 2001) appear to heavily impact the agricultural system.Contemporarily, these problems have been clearly recognized (Arora, 2018) and people are now requiring actions for experiencing clean air, water, local food, green spaces and landscapes diversification.Politicians are called to properly answer these requests.Accordingly, changes in agricultural management were promoted aimed at shifting the agronomic system from intensive farming models to multifunctional ones longing for yield and quality improvement, and environmental protection (Gasteyer, 2008).
Within this framework, meadows and grasslands management plays a key role for a sustainable agriculture (Lichtfouse et al., 2009;Velten et al., 2015).These systems are composed of various herbaceous species that strongly support wildlife and biodiversity.
Permanent meadows are well known to support carbon sequestration and to protect biodiversity and habitats (Jankowska-Huflejt et al., 2011).Additionally, they lead to a significant soil fertility improvement since the main management practices are mowing and fertilisation by organic manure.Moreover, meadows are ordinarily managed without using pesticides or herbicides reducing the air, water and soil pollution compared to other crops.Within this context, farmers can be described as land operators, able to preserve and shape local landscapes, thus providing public benefits (Primdahl & Kristensen, 2011).For example, fields biodiversity can be affected by how many times and in which period mowing are performed (Tälle et al., 2016;Uchida & Ushimaru, 2014).
Recently, the European Union (EU) has started to raise farmers' attention to the importance of their role through the Common Agricultural Policy (CAP).CAP supports farming activities for environmental purposes by requiring a proper management for grassland and pasture (Pe'er et al., 2020;Roederer-Rynning, 2010).In particular, CAP "green payments" aim at supporting those farming practices that favour the achievement of environmental and climate goals that are not reflected in market prices (Matthews, 2013;Westhoek et al., 2012).To obtain greening contributions farmers have to: a) diversify crops; b) maintain permanent meadows and pastures; c) preserve ecological focus area.Additionally, farmers cannot plough or convert permanent grassland in these areas (wwwec.europa.eu).If these requirements are not respected, an administrative sanction is applied (Article 77(6) of EU Regulation N°1306/2013, Anon.; Singh et al., 2014).It is worth to remind that since 2018 farmers have to submit a Geo Spatial Aid Application (GSAA) to apply for CAP contributions (Art.17 of Reg.(EU) No 809/2014).GSAA contain structured information about the applying farming company and the fields subsidies are asked for.Payments and related controls are demanded to the EU member States (Campinas & Rosa, 2010;López-Andreu et al., 2022).In Italy, this task is assigned to national or regional paying agencies (PP).PPs manage GSAAs through a platform based on Geographic Information System (GIS) called Integrated Management and Control System (IACS).In Italy IACS is used to properly manage the process that involves administrative (AC) and spot controls (SC).AC are automatically operated by IACS, for the 100% of paperwork, with the aim of testing formal legitimacy of farmers' applications.Differently, SC are aimed at testing technical compliance of farmers' declarations with CAP requirements (Sarvia et al., 2021).These concern a subset (5%) of received applications.Subset selection is achieved on risk and randomness criteria.Ordinarily, SCs are operated by photo interpretation of high-resolution satellite or aerial images; only in particular cases ground control campaigns support the process.These are obviously limited by related costs and availability of technical staff of the PPs; weather conditions and agronomic calendars of crop are additional limiting factors (Astrand et al., 2004).In the recent years, the European Copernicus programme has provided free satellite images with high temporal, spectral and geometric resolution.These data opened new monitoring scenarios especially in agriculture where the time domain plays an important role.Specifically, the high spectral resolution of Sentinel-2 (S2) imagery permits to derive several vegetation-related spectral indices useful for crop monitoring (Bannari et al., 1995;Leprieur et al., 1994).The Normalized Difference Vegetation Index (NDVI) is certainly the most explored spectral index for vegetation monitoring (Shanmugapriya et al., 2019).It is extensively adopted also in agriculture being effective for: (i) crop yield estimation (Parida et al., 2021(Parida et al., , 2021;;Y. Zhao et al., 2020); (ii) crop phenology monitoring (Boori et al., 2019;Misra et al., 2020); (iii) crop damage estimation in the insurance policies context (Sarvia et al., 2020); (iv) tree stability assessment (Grulke et al., 2020); (v) land cover/use mapping (Phiri et al., 2020;Steinhausen et al., 2018); (vi) natural disasters (Caballero et al., 2019;Lin et al., 2016;Pei et al., 2021;Suppasri et al., 2012); (vii) ecosystems characterisation (Andrew et al., 2014;R. Liu et al., 2020;Sarvia, De Petris et al., 2021); (viii) precision agriculture (Sarvia, et al. 2021, 2;Brisco et al., 1998;Liaghat & Balasundram, 2010).Within the CAP controls framework NDVI can also play an important role (Sarvia, et al. 2021;Campos-Taberner et al., 2019;Kanjir et al., 2018;López-Andreu et al., 2022;Sarvia et al. 2022) and EU is strongly favouring (and imposing) satellitebased services for CAP controls.
Focusing on meadows, articles 13 and 14 of the national report N° 5465 and the Regulation (EU) No 1307/2013, No 809/2014 and 746/2018 of the European Parliament and the Council impose that at least one cut per year must be performed by farmers applying for greening subsidies.Additionally, the greening contribution also includes ecological focus areas, i.e. fields that are not cultivated during the growing season.Consequently, the possibility of detecting this practice by satellite along the year in meadows is highly desirable (W.Zhao et al., 2020).In particular meadow monitoring by remote sensing has increased significantly over the last two decades (Filippa et al., 2022;Reinermann et al., 2020).The main research approach concerning the meadows can be can be summarised in three categories: i) detection of meadows from other crops; ii) classification of meadows types; iii) detection and monitoring of meadows agronomic practice.The tools ordinarily adopted, when dealing with these situations, involve the joint use of remote sensing data and on-site measurements (surveys, agronomic practices, etc.).While, considering the methods to detect the meadow mowing, time series seems to be favourable than the single spectral signature analysis.Meadow management practices are related to the farmer, the environmental condition and the plant phenology and therefore the selection of a single image for mowing detection is not recommended (Komisarenko et al., 2022).Only in the last years new algorithms were developed for the meadow mowing detection, consequently the topic is still being researched.The first method was developed in the 2010 and was based on LAI and NDVI time series (Courault et al., 2010).Subsequently, especially with the Copernicus missions able to overcome some limitations (spectral and temporal resolution) of previously used satellite missions, other approaches were proposed.For example, the Sen4CAP project (ESA's Sentinels for Common Agricultural Policy) develops a meadow mowing detection product was based on these data.Adreatta proposed the detection of grassland mowing frequency using time series of vegetation indices from S2 imagery within the Province of Trento (NE -Italy; Andreatta et al., 2022).Schwieder proposed meadows mowing detection and mapping by combining S2 and Landsat 8 time series within the Federal Republic of Germany (Schwieder et al., 2022).However, during the meadow growing season optical data availability is not always guaranteed due to cloud cover, which can lead to longer or shorter gaps within satellite image time series.For this reason, new studies have chosen to combine optical data with SAR (Synthetic Aperture Radar) technology.For example, Lobert proposed mowing event detection in permanent grasslands by Senitinel-1 (S1), Sentinel-2 (S2) and Landsat 8 time series within three different sites in Germany (Lobert et al., 2021).Komisarenko exploited time series of S1 and S2 to detect grassland mowing events using deep learning in Estonia (Komisarenko et al., 2022).In other cases SAR technology was also effective in detecting meadow mowings with Artificial Neural Networks on its own (Taravat et al., 2019).
Unfortunately, most of the proposed approaches are mainly based on supervised classification algorithm, where many training samples are required to fed classifications and finally count mowings.However, especially for CAP controls purposes, no reference data can be easily collected every years raising some doubts about a real transferability of such approaches to PPs operative routines.To fill this lack, in this work an automatic procedure based on NDVI time series frequency analysis was proposed.Specifically, this procedure, once parameterised with field data, could be used in other conditions (date and location) without requiring any additional training data.This operational procedure was explored for detecting, mapping and quantifying the number of cuts in meadows within the Piemonte Region (NW -Italy) applied to all meadow fields (about 72,500) presented in the 2021 GSAA.Finally, to test the proposed approach, reference meadows were provided by the Piemonte Regional Agency for Payments in Agriculture (ARPEA) during the 2021 ground campaign and confusion matrix results assessed.

Study area
The area of interest (AOI) of this work sizes 9636 km 2 (Figure 1).It includes the provinces of Torino, Asti and Vercelli (south part) that are located in the Piemonte region (NW Italy).In the area, climate is temperate with a continental character; yearly average cumulative rainfall and temperature are 930 mm and 11.9°C, respectively.These conditions make the area extremely suitable for agricultural activity and, therefore, a significant number of meadows is present.AOI was selected as pilot area by Piemonte Regional Agency for Payments in Agriculture (ARPEA).The study concerned the 2021 agronomic season.

Satellite data
According to the local mowing calendar (Sarvia et al., 2019) the correspondent S2 images (Level 2A products) were collected, covering the period 16 April 2021-21 September 2021.This part of the year was assumed as the reference period (RP) when grass cuts can be reasonably detected and counted.L2A S2 data are supplied as 100 km x 100 km tiles calibrated "at-the-bottom of the atmosphere" reflectance and orthoprojected in the WGS84 UTM reference frame (Delwart, 2015).They were obtained from the Copernicus Open Access Hub geoportal (scihub.copernicus.eu).A total of 4 tiles (namely T32TMQ, T32TMR, T32TLQ and T32TLR) were used to cover the entire AOI and 33 S2 L2A images were downloaded for each tile.Therefore a total of 132 S2 L2A images were used in this work according to RP. L2A is supplied together with the so called Scene Classification Layer (SCL) that maps type and quality of pixels within the scene.This can be effectively used to detect and remove "bad" observations like cloudy, shadowy and fault pixels (coded as 1, 3, 8, 9, 10 within SCL) during the pre-processing step of images (Gascon et al., 2014).

Reference data
For this work, ARPEA provided the number of cuts (NC) for 285 fields (448 ha) along the 2021 agronomic season.They were obtained by interviews to farmers and selected with a size greater than 0.1 ha.The average area of the selected plots results to be about 1.6 ha.Sample fields (hereinafter called RD) were used to train and validate the proposed methodology.RD was provided by ARPEA as a vector layer orthoprojected in the WGS84 UTM 32 reference frame (Figure 2).In RD 197 plots out of 285 were mowed showing from 2 to 4 cuts per year; 88 plots outs of 285 were not mowed (as they probably refer to ecological focus areas and therefore agronomically unmanaged).It is worth to highlight that RD did not contain any field showing a single cut.This is probably due to the local management practice for meadows that, in general, includes more than two cuts per year.
RD was randomly split in a validation (VS) and training (TrS) set, corresponding to 40% (114 fields) and 60% (171 fields) out the total, respectively.Table 1 reports the main features of RD. Figure 2 shows RD spatial distribution.

Meadows map
ARPEA supplied the 2021 GSAA map containing all fields somehow involved in CAP, as a vector layer orthoprojected in the WGS84 UTM 32 reference frame.It has to be stressed that, presently, GSAA data are not public, but a prerogative of PPs.2021 GSAA map was used to select meadow and unmanaged meadow fields making possible to generate a new vector map, namely meadow map (MM), for AOI.A total of 72,539 meadow fields were found, totally sizing 47,521 ha.

NDVI Time series generation
In this work, an approach based on the analysis of NDVI time series was proposed to detect meadows cuts.The first step consisted of mosaicking the 132 S2   images, resulting in a stack of 33 images, covering the entire AOI.Consequently, an NDVI map was generated for each of the 33 mosaicked L2A S2 images covering the reference period.NDVI was computed with reference to band 8 (NIR) and band 4 (red; Rouse et al., 1974) thus determining a geometric resolution of 10 m.NDVI maps were stacked along a time series (TS) and processed through a self-developed routine implemented in IDL v 8.0.1.In particular the local (at pixel level) NDVI temporal profile was preventively filtered removing all bad observations as mapped in the SCL layer.Remaining data were slightly smoothed by the Savitzky-Golay filter (kernel = 3, derivative = 0, degree = 1) and regularized with a spline interpolation (tensor value = 10) resulting in a stack of 30 equally timely spaced NDVI maps (Mishra et al., 2006).The amount of "tension" that is applied to the curve moves the cubic spline fit (tensor close to 0) to a polynomial interpolation (large tensor value) having as many parameters as the number of involved points (William H. Press et al., 2007).The resulting filtered and regularized NDVI time series assumed a nominal time frequency of 5 days.Finally, TS was averaged at plot level by ordinary GIS zonal statistics, in order to obtain a single NDVI profile for each field.It is worth to remind that field boundaries could generate spectral mixture but it can be assumed negligible due to averaging procedure (analysed field size > 0.1 ha, i.e. about 10 NDVI pixels).Moreover, a field selection was preventively operated to exclude from the analysis the fields that were not vegetated.With reference to the field average profile, the maximum NDVI (NDVI max ) value in RP was computed and thresholded to separate vegetated from not-vegetated fields.In particular, all GSAA fields showing a NDVI max value < 0.4 were removed from the following steps assuming them as not vegetated.This was needed since not all GSAA declarations can be assumed as true.Consequently, some fields declared as meadow cannot actually be meadow.This must be accounted for before proceeding in processing.All computations concerning zonal statistics and maximum selection were computed by ordinary GIS tools available in SAGA GIS 7.5 (Conrad et al., 2015).

Frequency analysis of NDVI temporal profile
Mowing is an agronomic practice periodically performed over meadows.Ordinarily, at least two mowings are performed by farmer in temperate zones (Dovel, 1996;Gong et al., 2020) and biomass is systematically removed after a regrowth period.This phenomenon generates abrupt changes in TS resulting into cycles of decrease/increase of NDVI values.In this situation, TS shows a periodicity due to the mowing that suggests the adoption of the Fast Fourier transform (FFT) to explore frequency harmonics.FFT is widely adopted to explore and detect signal periodicity especially over vegetation (Borgogno-Mondino & Lessio, 2018;Capodici et al., 2020).FFT is in charge of decomposing the time series into several sine/cosine components having different frequencies, amplitude and phase (Shumway & Stoffer, 2006).In this work, FFT approach was used to detect and map TS periodicity in meadows, assuming cycles related to operated local mowing practices.TS processing was achieved setting some agronomic constraints.One is related to the expected NDVI signal intensity of meadows (as described by NDVI) that in RP is expected to be generally increasing (Sarvia, De Petris, Borgogno-Mondino et al., 2022).Mowings are expected to introduce periodic variations around this generally increasing trend.Consequently, before detecting cycles, a linear trend was locally modelled and removed from the correspondent TS to make stationary the local NDVI profile.In particular first order polynomial was fitted using ordinary least squares involving all observations within the reference period.The resulting linear trend was removed from TS and the periodicity explored over the residuals (i.e.detrended TS).
The Lomb-Scargle periodogram (VanderPlas, 2018) was computed for each field using R software vs 4.1.1.(R Development Core Team, R, 2013).The frequency (ω peak ) corresponding to the higher harmonic amplitude was found along the spectrum and mapped for each field.Jointly, the correspondent amplitude (A peak ) was computed by eq. 1 and mapped accordingly.
A peak ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P peak � 2σ 2 ydt q (1) where P peak is the normalized peak power value (W.H. Press, 1992) in the Lomb-Scargle spectrum and σ 2 ydt is the variance of the de-trended TS.The A peak represent the average NDVI oscillation (related to biomass removal and regrowth) in the detrended TS. ω peak and A peak were computed for each field and mapped in AOI.

Detecting number of cuts
Number of potential cuts in a field can be easily estimated once ω peak is known according to eq. 2.
The result of multiplication has to be rounded up to get an integer value as counter.
where n obs is the number of images in RP (i.e.30).
Several paradigmatic examples corresponding to four fields from RD are given in Figure 3 in order to discuss about some issues to account for while applying eq. 2. In particular: (i) Figure 3a-c show the TS profile, the de-trended TS profile for a natural development of meadows having no cut and the correspondent spectrum, respectively; (ii) Figure 3d-f show the TS profile, the de-trended TS profile for meadows having two cuts (corresponding to ω peak = 0.065) and the correspondent spectrum, respectively; (iii) Figure 3g-i show the TS profile, the de-trended TS profile for meadows having three cuts (corresponding to ω peak = 0.105) and the correspondent spectrum, respectively; (iv) Figure 3j-l show the TS profile, the de-trended TS profile for meadows having four cuts (corresponding to ω peak = 0.115) and the correspondent spectrum, respectively; Particular attention must be paid to Figure 3 b and c being a paradigmatic case to make evident that ω peak alone is not able to detect significant variations possibly relatable to cuts.In fact, with reference to ω peak , the spectrum of Figure 3c would suggest two cuts (@ω peak =0.065), shifting the meaning of the observed field.Consequently, the attention must be transferred to P peak of Figure 3c, that corresponds to a A peak < 0.05 NDVI points.This value is close to the NDVI uncertainty as reported in literature and certainly too low for justifying a cut (Agnes Bégué et al., 2010;Zhou et al., 2021).Compliance of profile variations with the expected ones related to biomass removal by cut, has therefore to be preventively tested by considering the A peak value correspondent to ω peak .
To mask out those fields where the strongest harmonic is not significant, a procedure based on the definition of a proper threshold value of A peak able to separate mowed from not-mowed fields was developed.With reference to the above mentioned training set from MM (TrS), the cumulated frequency distribution (CFD) of correspondent A peak value was computed for both the class 0 (no cuts) and for the remaining class 2, 3 and 4 (jointly considered).CFD were then compared through the two-tailed Kolmogorov-Smirnov (KS) nonparametric test (Stephens, 1970) to test significance of difference.The KS distance (KS-D) from the KS test was used as separability metric and used to find the A peak percentile that guarantee the highest separability between the compared CFDs.This method allows to define an objective threshold able to separate the two analysed classes (Adeodato & Melo, 2016).
According to the above-mentioned strategies, ω peak and A peak were mapped at field level.Not-mowed fields were masked out using the A peak threshold from the KS-distance approach.For remaining fields, ω peak was translated into the correspondent number of cuts (Eq.2).Depending on the number of mowings, five different classes were found, and the correspondent mowing counting map (MCM) generated.Table 2 reports MCM class meaning.
MCM was finally validated with respect to VS (consisting of 114 RD) by confusion matrix.User's accuracy (UA), Producer's accuracy (PA) and Overall accuracy (OA) where computed (Hay, 1988).It is worth to stress that class 1 (i.e. one mowing along RP) is not present in VS.For this class it was therefore not possible to carry out a proper validation.

MM Pre-processing
To exclude from the analysis those fields possibly being other than meadows (in spite of GSAA declaration) a field selection was preventively operated from MM based on the NDVImax value of the local field average NDVI temporal profile and all fields showing NDVImax value < 0.4 were removed from MM.Therefore, a total of 112 fields proved to not be meadows corresponding to 0.15% of 2021 declared fields.

Frequency analysis of NDVI temporal profile
ω peak and A peak were mapped at field level and the correspondent frequency histograms generated (Figure 4).
Figure 4c shows the frequency distribution of ω peak in AOI.It can be noticed that the most of values are placed between 0.066 and 0.100, corresponding to a number of cuts of 2 or 3 (per year), respectively.It should be stressed that the histograms reported in Figure 4 include all fields even comprised the not mowed ones.
Figure 4d shows that the most of A peak values are placed around 0.08 corresponding to a decrease of NDVI values about 0.16 (the entire harmonic oscillation).This value is consistent with the expected NDVI variation related to a cut (biomass removal) as also found in previous studies from authors (Agnès Bégué et al., 2018;Sarvia et al., 2019).

Thresholding A peak
The interpretation key of A peak can be summarized as it follows: higher A peak value, higher the biomass   statistical distributions of A peak and relative KS nonparametric test.
Figure 5 shows that distribution of class 0 (not mowed fields) is significantly different from other classes (mowed fields).Specifically, A peak CFD of class 0 results to be significantly different from the other classes with a KS-D value always set around 0.920.Differently, no significant differences were found between A peak CFD classes 2, 3 and 4, resulting in KS-D always lower than 0.292.Accordingly, a threshold value for A peak able to separate mowed from not-mowed meadow fields, was defined by KS-D with respect to the compared CFDs (Figure 6).The value correspondent to KS-D was found to be 0.065.This value represents the 92 nd and 2 nd percentiles with respect to the mowed and not mowed CFDs, respectively.This threshold results to separate more than 90% of not-mowed from mowed meadow fields,  suggesting to be a promising operational tool for PPs future workflow.It is worth to highlight how this thresholding is able to take care about local variability (e.g.temperature, precipitations, fertilizations) allowing to classify mowed fields in other regions.

Counting and mapping cuts
Once ω peak was computed for all monitored fields, it was translated by eq.2 into the correspondent number of cuts that were mapped in MCM.Result is reported in Figure 7. MCM makes possible to easily synthesize local situation of meadows and derive some statistics useful for territorial management of AOI (see, Table 3).
Table 3 shows that at least one mowing was classified in 54,538 meadow fields (about 75% out of total); conversely, 18,001 fields (about 25%) were detected as not-mowed.In particular, for fields classified as mowed, 7423 are in class 1 (about 10%), 23,178 are in class 2 (about 32%), 19,895 are in class 3 (about 27%) and 4042 are in class 4 (about 6%).Similar percentage values were found with reference to areas in place of number of fields.Meadow management found in this work is similar to that mapped by Schwieder.In particular, he found in Bavaria (Germany) that more than 70% of the meadows were mowed more than 2 times a year, confirming that grassland management tends to occur several times a year (Schwieder et al., 2022).
MCM was finally validated with reference to VS and the correspondent confusion matrix computed (Table 4).
Table 4 shows that all classes present a satisfying UA (> 75%) and PA (> 76%).Globally, an OA value of 87% was found.These results demonstrate the effectiveness of the proposed methodology.It is worth to remind that computation was achieved only for those meadows that entered CAP 2021 campaign.While operationally interpreting results, PPs should cautionary consider that the procedure tends to over-estimate class 0 of about 24% (UA for Class 0 = 76%).This is equivalent to say that of the 18,000 fields mapped with 0 mowing in the MCM, potentially 24% (about 4320) would result in a misclassification.Differently, at least 2 cuts can be detected with a UA of about 93%, thus including about 7% of not-mowed fields.With reference to class 0 the reliability of result (MCM) is lower (UA = 76%) but ensure that about 97% (PA) of actually not-mowed fields are correctly mapped.
In order to evaluate and discuss the results obtained in this work, a comparison with the existing literature was performed.The formal product provided by the Sen4CAP project in 2019, based on S1 and S2 time series and developed on meadows in Belgium, achieved a precision of 79-84% in the mowing detection (De Vroey et al., 2021) close to the one obtained in this work.
Andreatta detect meadows mowing frequency using the drop of NDVI derived from S2 data time series in the Province of Trento (NE -Italy) with an OA (89-93%) close to the one from this work (Andreatta et al., 2022).Similar results were obtained by Halabuk in Slovakia with a cut-uncut classification in extensively managed grasslands resulting in OA equal to 85% (Halabuk et al., 2015).While considering a time series approach Kolecka, in the region of Canton Aargau -Switzerland, starting from a NDVI temporal profile derived from S2 data, proposed a drop-detection algorithm resulting in an OA equal to 77% of correctly detected meadow mowing (Kolecka et al., 2018).Conversely, Estel combine a NDVI time series from 2000 to 2012 derived by MODIS satellite data with agricultural statistics to map grassland management intensity in Europe obtaining an OA of 77% (Estel et al., 2018).Schwieder obtained their best results in mapping meadow mowings in Germany by combining S2 and Landsat-8 time series for the years from 2017 to 2020.Specifically, using machine learning algorithms, they report an average of f-score equal to 0.6 for all years, proving an algorithm that detects mowing events in the S2 time series based on residuals from an assumed undisturbed phenology, as an indicator of meadow use intensity (Schwieder et al., 2022).Finally, Lobert combine NDVI, backscatter cross-ratio and interferometric coherence for mowing event detection in permanent meadows obtaining a F1-score equal to 84% in Germany (Lobert et al., 2021).This comparison provides a first benchmark useful to affirm that our approach has accuracy somehow similar to the  ones currently present in literature.Therefore, it can be suggested that our approach can be adopted in other areas once the agronomic calendar of mowing has been identified.Nevertheless, the great advantage of the proposed methodology is that no training phase to fed classificator is required once A peak was identified.Specifically, once it will be verified that A peak values would be remain the same over time, even if the different climatic conditions that afflict the meadows each year still persist and therefore could affect A peak values, such threshold could be used instead of using field data to identify meadows managed and not managed within the CAP framework.In particular, counting procedure allows to classify the meadow type in: unmanaged meadow or pasture (0 mowing), meadowpasture (1 mowing) and meadow (≥ 2 mowings) supporting PPs during farmer requests control.Specifically, if more than 2 mowings were detected one can be reasonably alerted about the truthfulness of farmer declaration.While some doubts could arise if only 1 mowing is detected since it could be confused with other anomalies occurred during the phenological development of the fields like hail damage, lodging or drought.Furthermore, the major benefit of counting mowings is that PP can be ensured that the meadow has been successfully managed for the correspondent year.For example, if a meadow is mowed 4 times, the probability that the count has been affected by an extreme event, and thus misclassified, will certainly be lower than for a meadow on which 1 or 2 mowings have been detected.Furthermore, in a control context such as the one required in the CAP framework, the opportunity to reduce field data, which can only be obtained through costly control campaigns, during the classification training phases of algorithms, would allow largescale monitoring of grassland conditions and management in a practical, cost-effective and time-saving way.
In spite of these promising results, some critical points of the proposed methodology can be still recognized.The first one concerns the potential cloud cover affecting S2 data in RP (Xie et al., 2017;Zhai et al., 2018), that especially in the early-phenological stages of meadows, when adverse weather conditions are most probably present (i.e. in Europe spring), could compromise the meadows mowing detection.To achieve this issue, some experiences proved how a multi-temporal analysis of SAR (Synthetic Aperture Radar) can detect mowing (Komisarenko et al., 2022;W. Zhao et al., 2020).Future developments are expected in order to improve the proposed approach.The first one will be the integration of S2 data with S1 SAR data aiming at exploiting its clouds penetration capability and therefore to cover the whole year, with particularly concern to the spring period noted for being very cloudy.Moreover, we will attempt to include meadow fields that are cut only once during future ground checks in order to include and test this class within the accuracy classification processes.Finally, a critical issue concerns the A peak role in mowing detection.Several conditions (e.g.drought, lowfertility and water stress) can affect grassland development and consequently times of farmers' activities (e.g.mowing).These unpredictable variations may change A peak , shifting deductions from the conditions that this work was calibrated on.Consequently, an investigation regarding A peak threshold's stability over time will be required before this approach can be used by PP on an ongoing manner.Furtherly, phenological metrics, cumulated rains (derived by rain gauge or from model data) and new spectral vegetation indices related to canopy water content will be explored considering climate-induced variability (i.e.local drought) on meadows with the aim to limit the misclassification overestimation and maximizing the accuracies in the meadow mowing detection.

Conclusions
In this work a prototypal approach aimed at detecting meadow mowing was presented to support CAP controls.This study was explicitly required by Piemonte Regional Agency for Payments in Agriculture (ARPEA) to verify meadow management within the Piemonte region (NW -Italy) in a simple and direct way, aiming at satisfying the national report N° 5465 applying Regulation (EU) No 1307/2013.The proposed procedure is based on the processing of time series of NDVI maps as derivable from the Copernicus Sentinel-2 data.With reference to the 2021 growing season (16 April 2021-21 September 2021) and to a study area located within the Piemonte region (NW -Italy) a processing workflow was implemented based on the decomposition of the local NDVI temporal profile by FFT.The frequency of the strongest harmonic and its amplitude proved to be able to separate mowed from not-mowed fields and giving an estimate of the number of cuts.In particular, A peak proved to be useful to detect not-mowed meadows if an opportune threshold value is found.This was obtained by comparing the statistical distributions of A peak for the mowed and not-mowed classes by the non-parametric Kolmogorov-Smirnov test.Strongest frequency resulting from FFT was used to count and map cuts in mowed fields.The approach was finally validated with reference to the available validation set.Results showed an OA > 87%, confirming that the approach is promising.Authors believe that the proposed methodology could be a valid alternative to ordinary controls carried out by PPs in the CAP controls framework.Moreover, it can support CAP control process by focusing PP checks on fields where classification (or count) is different from the one declared by farmers.This makes possible to avoid controls about those applications that are reasonably correct, thus saving time and money.

Figure 2 .
Figure 2. a) Spatial distribution of RD class type; b) Spatial distribution of RD used within the validation and training phase.Reference system is WGS84/UTM 32 N, EPSG: 32,632.

Figure 3 .
Figure 3. TS profile, de-trended TS profile (field average) and relative Lomb-scale spectrum for four meadow fields.Dotted line shows the values of W peak and A peak considered (a) TS profile expected from meadows not mowed; (b) De-trended TS profile expected from meadows not mowed; (c) Lomb-scale spectrum expected from meadows not mowed; (d) TS profile expected from mowed two times; (e) De-trended TS profile expected from meadows mowed two times; (f) Lomb-scale spectrum expected from meadows mowed two times; (g) TS profile expected from mowed three times; (h) De-trended TS profile expected from meadows mowed three times; (i) Lomb-scale spectrum expected from meadows mowed three times; (j) TS profile expected from mowed four times; (k) De-trended TS profile expected from meadows mowed four times; (l) Lomb-scale spectrum expected from meadows mowed four times.

Figure 5 .
Figure 5. Boxplots of A peak distributions and related results about the Kolmogorov-Smirnov nonparametric test for Training set fields.*** correspond to KS p-value lower than 0.001.Black lines defines (bottom-up) the 5th, 25th, 50th, 75th and 95th percentiles.

Figure 6 .
Figure 6.Cumulated frequency distribution (CFD) of A peak of RD classes.Red line indicates the threshold value defined by KS-D method and used in this work to separate the mowed from unmowed meadow fields.

Table 1 .
Number and size and relative split of surveyed fields used as reference in this work.

Table 2 .
Class code meaning according to the number of mowings classified during RP.
Bégué et al., 2018)the cut (Agnès AgnèsBégué et al., 2018); differently, low A peak values may indicate mowing on sparse vegetation or simple noise due to environmental conditions.To better investigate this issue, A peak CFDs corresponding to the TrS classes were compared.Figure5shows boxplots summarizing

Table 3 .
Number of plots and corresponding area of the classes mapped in MCM.

Table 4 .
Confusion matrix involving the detected classes (number of cuts).Matrix elements are reported in number of fields.Grey cells contain correctly classified fields for each class.