Mapping Grassland Frequency Using Decadal MODIS 250 m Time-Series: Towards a National Inventory of Semi-Natural Grasslands

Semi-natural grasslands are perennial ecosystems and an important part of agricultural landscapes that are threatened by urbanization and agricultural intensification. However, implementing national grassland conservation policies remains challenging because their inventory, based on short-term observation, rarely discriminate semi-natural permanent from temporary grasslands. This study aims to map grassland frequency at a national scale over a long period using Moderate Resolution Imaging Spectroradiometer (MODIS) 250 m satellite time-series. A three-step method was applied to the entire area of metropolitan France (543,940 km²). First, land-use and land-cover maps—including grasslands—were produced for each year from 2006–2017 using the random forest classification of MOD13Q1 and MYD13Q1 products, which were calibrated and validated using field observations. Second, grassland frequency from 2006–2017 was calculated by combining the 12 annual maps. Third, sub-pixel analysis was performed using a reference layer with 20 m spatial resolution to quantify percentages of land-use and land-cover classes within MODIS pixels classified as grassland. Results indicate that grasslands were accurately modeled from 2006–2017 (F1-score 0.89–0.93). Nonetheless, modeling accuracy varied among biogeographical regions, with F1-score values that were very high for Continental (0.94 ± 0.01) and Atlantic (0.90 ± 0.02) regions, high for Alpine regions (0.86 ± 0.04) but moderate for Mediterranean regions (0.62 ± 0.10). The grassland frequency map for 2006–2017 at 250 m spatial resolution provides an unprecedented view of stable grassland patterns in agricultural areas compared to existing national and European GIS layers. Sub-pixel analysis showed that areas modeled as grasslands corresponded to grassland-dominant areas (60%–94%). This unique long-term and national monitoring of grasslands generates new opportunities for semi-natural grassland inventorying and agro-ecological management.


Introduction
Grasslands are one of the most extensive terrestrial ecosystems on Earth and a source of food for livestock [1]. Definitions of "grassland" differ among disciplines and include a wide variety of land-use and land-cover (LULC) types. Temporary grasslands, which are part of crop rotations, are composed of seeded vegetation [2] and mostly have a life span of less than 5 years [3]. Conversely, permanent grasslands that are not part of crop rotations can be defined as "land on which vegetation is composed

Rationale of the Approach
The approach first uses RF modeling to identify grasslands each year among the other LULC classes and then combines the annual LULC maps to characterize the spatio-temporal dynamics of grasslands from 2006-2017 ( Figure 2). Multi-class modeling (water, woods, urban, crop and grassland) was used instead of single-class modeling (grassland) to prevent arbitrarily selecting a probability threshold for discretizing the grassland class and to quantify over-detection errors [51].

Rationale of the Approach
The approach first uses RF modeling to identify grasslands each year among the other LULC classes and then combines the annual LULC maps to characterize the spatio-temporal dynamics of grasslands from 2006-2017 ( Figure 2). Multi-class modeling (water, woods, urban, crop and grassland) was used instead of single-class modeling (grassland) to prevent arbitrarily selecting a probability threshold for discretizing the grassland class and to quantify over-detection errors [51].

Satellite Data
The satellite data include MODIS MOD13Q1 and MYD13Q1 composites for the period 2006-2017 (12 years), providing a combined temporal resolution of 8 days [52]. MODIS data from [2003][2004][2005] were not used due to the unavailability of reference data before 2006. These MODIS data were downloaded for the whole of France in WGS84 (EPSG 4326) projection using the USGS Application for Extracting and Exploring and Extracting Analysis Ready Samples (AppEEARS) [53]. We used the red (620-670 nm) and near-infrared (841-876 nm) spectrum bands with a spatial resolution of 250 m. We did not use vegetation indices such as normalized difference vegetation index (NDVI) or enhanced vegetation index (EVI) since they do not increase modeling accuracy significantly [42]. To generate annual grassland maps, MODIS satellite data were stacked by vegetative year, ranging from November 1 of year n to October 30 of year n+1.
Although MOD13Q1 and MYD13Q1 products are cloudless, some pixels had missing values on specific dates, especially in winter, when the days are short and cloudy. To ensure continuous time-series for all pixels, a time interpolation step was required. In this study, the nearest neighbor interpolation method was applied, where missing (NA) values were replaced with the values of the closest date. While the smoothing of the time series restores vegetation phenological profiles [54], the use of raw data was preferred because (i) the dataset consists of a dense time series (8 days) with therefore less missing data than in the 16-day series, (ii) it is still a challenge to discriminate noise from natural variations in vegetation, and part of the information contained in raw time series is lost after filtering [55], and (iii) the contribution of smoothing to the accuracy of LULC classification remains marginal or even negative [56]. The range of pixel values varied with dates and spectral bands, which could lead to the dominance of certain bands in the RF model. Hence, all MODIS pixel values were standardized by setting the overall mean of each spectral band to 0 and the standard deviation to 0.25.

Reference Data
Reference data came from LPIS and Copernicus HRL GIS files (Table 1). The LPIS vector data were collected for each year from 2006-2017 to calibrate and validate modeling of the crop and grassland classes. Unambiguous crop labels 1-15, 20-23 and 25 in the LPIS were assigned to the crop class, while labels 18 (permanent grasslands, more than 5 years old [3]) and 19 (temporary grasslands, 5 years old or less) were assigned to the grassland class (Table S1). In addition, the Copernicus HRL raster data "Water & Wetness", "Tree Cover Density" and "Imperviousness" were collected to calibrate and validate modeling of the water, woods and urban classes, respectively. These layers, which show the dynamics of LULC in the EU since the 2000s, are produced at a 20 × 20 m spatial resolution from satellite data, such as Landsat-8 or Sentinel, with an overall accuracy above 80% [25]. The "Tree Cover Density" and "Imperviousness" HRLs are expressed as a density of occupancy ranging from 0-1, while the "Water & Wetness" HRL is composed of four classes: (1) permanent water, (2) temporary water, (3) permanent wetness and (4) temporary wetness. The sampling step aimed to identify and select "pure" MODIS pixels for use as calibration and validation samples for modeling grasslands. First, the footprint of each MODIS pixel was vectorized and overlaid on the reference data. Then, samples were selected for each year from MODIS pixels covering a homogeneous (i.e., single) LULC class (Figure 1). More precisely, the following rules defined the homogeneity criterion: (i) A sample of the grassland or crop class is a MODIS pixel strictly included within a parcel block that contains only grassland or one crop type (i.e., it covers >80% of the pixel's area), respectively [57]; (ii) A sample of the woods or urban class is a MODIS pixel with >85% of its area covered by a density of "Tree Cover Density" or "Imperviousness" HRL >0.8, respectively. Indeed, since "Tree Cover Density" and "Imperviousness" HRL products are expressed in density values from 0 to 1, a threshold value (0.8) was set to discriminate between wooded (or urban) areas and non-wooded (or non-urban) areas; (iii) A sample of the water class is a MODIS pixel with >90% of its area covered by the (1) permanent water class of the "Water & Wetness" HRL. Temporary water (2), permanent wetness (3) and temporary wetness (4) classes were discarded because they can characterize either water areas or grasslands.
The sample spatial density was very heterogeneous throughout France, which could lead to over-representation of some regions in RF modeling. To overcome this issue, sub-sampling was used, with only one sample being kept per mesh of 5 × 5 km per class and per year. For each sample of year n, spectral values were extracted from the MODIS pixels of year n. The number of samples used per year ranged from 21,703-24,565 for all LULC classes, of which 4394-5974 were of the grassland class (Table S2). Then, for each year, the samples were divided randomly and equally per class into a calibration and a validation dataset.

Random Forest Modeling
The RF modeling step identified areas covered by grassland each year. To this end, the five LULC classes (water, urban, woods, crop and grassland) were modeled annually from 2006-2017. For each year, a RF model was applied to the annual MODIS time-series and calibrated with the calibration dataset using a 10-fold approach with three repeats [58]. The maximum number of trees (ntree) was set to 500. Tuning was applied to define the optimal number of variables used for each branch of the RF model tree (mtry) based on the highest Kappa value [59]. Independent validation was performed for each year by applying the RF model to the validation set and calculating three metrics: overall accuracy, the Kappa index and the F1-score of the grassland class. To highlight spatial variation in modeling accuracy, the F1-score of the grassland class was also calculated by biogeographical region.

Grassland Dynamics Analysis
The next step analyzed the spatio-temporal dynamics of grasslands by combining the annual LULC maps derived from the RF modeling. To identify and correct illogical LULC transitions, we applied a temporal filter using the approach developed by Clark et al. [60]. For example, it is impossible for a pixel classified as crop one year to be classified as urban the next year and then to return to the crop class two years later. Conversely, a pixel can change from grassland to crop and then back to grassland, depending on the agricultural practices. To this end, we applied transition rules in a 3-year moving window (year n, n+1, and n+2) for each pixel starting in 2006 and repeating it annually until 2017 ( Table 2). In detail, the filter inspected each pixel to determine whether it had the same class in years n and n+2; if so, and if the class in year n+1 created an illogical transition, then the class in year n+1 was corrected to equal the class in year n. The filter was then advanced by one year in the time-series and the corrected classes of previous years were considered when applying the transition rules. A pixel's grassland frequency was calculated by summing the number of years that it was classified as grassland from 2006-2017, and the frequency values were then scaled to a [0-1] range to facilitate result interpretation (12). To visualize this layer's thematic contribution, its grassland classification was compared to those of national, European and global LULC layers for the lower Couesnon marshes, a Natura 2000 site near Mont-Saint-Michel (HER 12). This site is dominated by natural grassland habitats [43]. In total, five GIS layers were collected: (i) The components "231-Pastures", "242-Complex cultivation patterns", "321-Natural grasslands" and "411-Inland marshes" of the 2018 CORINE Land Cover layer; (ii) The "grassland" HRL of the 2015 Copernicus layers; (iii) The "18-Permanent grasslands" and "19-Temporary grasslands" components of the 2016 LPIS layer; (iv) The "211-Grasslands" component of the 2018 French national LULC layer ("OSO") [23]; (v) The grassland frequency maps, calculated for the period 2006-2017 for each of the six MODIS MCD12Q1 v6 products at 500 m spatial resolution ("IGBP","UDM", "Annual LAI", "Annual BGC", "Annual PFT", and "LCCS 3") [40].

Sub-Pixel Analysis
Since MODIS data have a coarse spatial resolution (250 m), the pixels classified as grassland have, at a finer scale, LULC diversity that varies according to the landscape. The objective of the sub-pixel analysis step was to quantify the LULC composition of MODIS pixels classified as grassland from a fine-scale reference LULC layer. This reference layer was generated in raster format at 20 m resolution Remote Sens. 2019, 11, 3041 9 of 21 by combining (i) the parcels identified as crop and grassland (i.e., permanent and temporary grasslands) classes in the 2016 LPIS, (ii) the areas classified as permanent and temporary water in the 2015 "Water & Wetness" HRL, (iii) the urban areas classified with a density >0.8 in the 2015 "Imperviousness" HRL and (iv) the woody areas classified with a density >0.8 in the 2015 "Tree Cover Density" HRL. Areas without information in the reference layer (e.g., a woods density <0.8, agricultural parcels not described by the LPIS) were described as mixed. For each pixel of the 20 m reference map, the majority LULC class relative to the HRL Copernicus and LPIS layer has been assigned. The reference layer was then intersected with the grassland layer modeled from the 2016 MODIS data corrected by the temporal filter. MODIS 250 m pixels modeled as grassland with a mixed percentage >20% were discarded. The LULC composition of MODIS pixels classified as grassland was measured as a percentage of MODIS pixel area per HER.

Identification of Grasslands
The accuracy of RF modeling of the five LULC classes from 2006-2017 was high, regardless of the year considered (overall accuracy ranging from 0.93 to 0.96, Kappa index ranging from 0.91 to 0.94), especially for the grassland class (F1-score ranging from 0.89 to 0.93) ( Table 3). Temporal filtering increased the accuracy of the five LULC classes slightly (+ 0-0.02 for overall accuracy, + 0-0.02 for the Kappa index), specifically of the grassland class (+ 0-0.02 for the F1-score). Analysis of the 12 annual confusion matrices highlighted confusion among the grassland, crop and woods classes. More precisely, most under-detections of the grassland class (producer's accuracy 87.8% ± 2.3) were due to confusion with the crop (6.5% ± 1.6) and woods (5.1% ± 1.0) classes, as were most of the over-detections (user's accuracy 93.8% ± 1.0) (crop 3.2% ± 0.6, woods 1.0% ± 0.2) (Table S7)   Table 3. Overall accuracy, Kappa index and F1-score for the grassland class per year before and after temporal filtering (the highest value for each year and each accuracy metric is in bold).

Characterization of Grassland Frequency
The map of grassland frequency in France derived from annual LULC maps from 2006-2017 illustrates areas where grassland frequency was high, such as hedged landscapes in the Armorican massif (HER 12), the Ardennes (HER 22) and the Central massif (HER 21), and mountain pastures in the Pyrenees (HER 2) and the Alps (HER 1) ( Figure 3). Conversely, grassland frequency was low on large cereal plains, such as Alsace (HER 18), the Paris basin (HER 9) and the Aquitaine basin (HER 14). In detail, the map also indicated spatial imbrication of permanent (frequency >0.8) and temporary (frequency <0.5) grasslands in the Normandy bocage (HER 12, Figure 3A) and the residual occurrence of permanent grasslands along valleys in the open-field landscape of the Coteaux calcaires de l'Est (HER 10, Figure 3B). 2006  878  198  840  173  2007  826  170  847  181  2008  837  167  823  156  2009  818  156  837  142  2010  841  154  798  132  2011  853  151  809  147  2012  835  146  829  144  2013  851  149  834  139  2014  817  130  830  131  2015  844  121  800  107  2016  758  88  751  100  2017  796  94  738  87   Table S7. Annual confusion matrices between modeling the five LULC derived from the MODIS time-series (rows) and the validation samples (columns) for the whole of France from 2006-2017. The grassland frequency map for the lower Couesnon Natura 2000 site consistently showed that it had a high frequency of grassland, although edges of the marshes were unclear due to the low spatial resolution of the MODIS data ( Figure 4). In comparison, the CORINE Land Cover layer adequately identified the natural grasslands of the lower Couesnon, but ambiguously classified them as "231-Pastures", "242-Complex cultivation patterns" or "411-Inland marshes" and not as "321-Natural grasslands". The grassland HRL under-detected most natural grasslands, especially those that are flooded several months each year, and did not discriminate between temporary and natural grasslands. The LPIS correctly classified natural grasslands as "18-Permanent grasslands", but with incomplete mapping, since the western end of the marshes (48.54 • N, 1.54 • W)-not under the EU CAP-is missing. Finally, the national land cover layer "OSO" best identified grasslands, but its classification ("211-Grasslands") could not discriminate between semi-natural and temporary grasslands. The comparison of the grassland frequency maps derived from the MCD12Q1 v6 products with the Google Earth image and CORINE Land Cover map (Figure 4) points out-beyond the coarse spatial resolution of the MCD12Q1 v6 products-a very strong under-detection of grasslands of the IGBP, UMD, Annual FTP and LCCS 3 maps, and a strong under-detection and strong over-detection of grasslands for the Annual LAI and Annual BGC maps.

Year Calibration Points Validation Points All LULC Classes Grassland Class All LULC Classes Grassland Class
grasslands of the IGBP, UMD, Annual FTP and LCCS 3 maps, and a strong under-detection and strong over-detection of grasslands for the Annual LAI and Annual BGC maps.

Land Cover Percentages in MODIS Pixels Classified as Grassland
Comparing the 2016 LULC classification at 250 m spatial resolution to the reference layer at 20 m spatial resolution resulted in the selection of 1,612,318 MODIS pixels-ca. 18% of the total area of

Land Cover Percentages in MODIS Pixels Classified as Grassland
Comparing the 2016 LULC classification at 250 m spatial resolution to the reference layer at 20 m spatial resolution resulted in the selection of 1,612,318 MODIS pixels-ca. 18% of the total area of France ( Figure 5). In the MODIS pixels classified as grassland throughout France in 2016, the grassland class was dominant on all HERs (60%-94%), followed by the crop (0%-25%), mixed (4%-10%) and woods classes (0%-4%). The water and urban classes were negligible (<1%) (Figure 6). Nonetheless, analysis by biogeographical region identified differences in grassland area among HERs. The Alpine region had the highest percentages of grassland area, both in the Pyrenees (HER 1) and the Alps (HER 2) (77% and 94%, respectively). The percentage of crop area was low (8%) in the Pyrenees and negligible (<1%) in the Alps. The Atlantic region had less grassland dominance and a relatively high percentage of crop area. Specifically, the Landes (HER 13) and Aquitaine coteaux (HER 14) had the lowest percentage of grassland area (62% and 60%, respectively) but the highest percentage of crop and wooded areas (25% and 4%, respectively). Conversely, the percentage of grassland area was higher in the Armorican massif (HER 12) and Sologne (HER 20) (78% each). The Continental region had a high percentage of grassland area, especially in the Central massif (HER 3, 17 and 21; 82%, 81% and 85%, respectively) and the Ardennes (HER 22; 80%). The percentage of crop area was low (8%-15%), except in Alsace (HER 18) and the Saône plain (HER 15) (20% and 22%, respectively). The percentage of wooded area was generally low (0%-2%), except in the Jura (HER 5) and Vosges (HER 4) mountains (4% and 3%, respectively). MODIS pixels classified as grassland in the Mediterranean region had a high percentage of grassland area (76%-83%), while the percentage of crop area was moderate (9%-15%). The wooded area was smallest (<1%), excluding the Cévennes massif (HER 8) (4%). Nonetheless, analysis by biogeographical region identified differences in grassland area among HERs. The Alpine region had the highest percentages of grassland area, both in the Pyrenees (HER 1) and the Alps (HER 2) (77% and 94%, respectively). The percentage of crop area was low (8%) in the Pyrenees and negligible (<1%) in the Alps. The Atlantic region had less grassland dominance and a relatively high percentage of crop area. Specifically, the Landes (HER 13) and Aquitaine coteaux (HER 14) had the lowest percentage of grassland area (62% and 60%, respectively) but the highest percentage of crop and wooded areas (25% and 4%, respectively). Conversely, the percentage of grassland area was higher in the Armorican massif (HER 12) and Sologne (HER 20) (78% each). The Continental region had a high percentage of grassland area, especially in the Central massif (HER 3, 17 and 21; 82%, 81% and 85%, respectively) and the Ardennes (HER 22; 80%). The percentage of crop area was low (8%-15%), except in Alsace (HER 18) and the Saône plain (HER 15) (20% and 22%, respectively). The percentage of wooded area was generally low (0%-2%), except in the Jura (HER 5) and Vosges (HER 4) mountains (4% and 3%, respectively). MODIS pixels classified as grassland in the Mediterranean region had a high percentage of grassland area (76%-83%), while the percentage of crop area was moderate (9%-15%). The wooded area was smallest (<1%), excluding the Cévennes massif (HER 8) (4%). Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 21

Can MODIS 250 m Time-Series Combined with the RF Classifier Discriminate Grasslands from Other LULC Types at the National Scale?
Overall, the RF classification of MODIS 250 m time-series accurately mapped grasslands each year from 2006-2017 at the scale of France (F1-score 0.89-0.93). These results confirm the advantage of using the RF classifier to discriminate between grasslands [1], which in France are of different types (temporary, wet, mesic, dry or alpine) and thus have different spectral patterns, as well as to manage high-dimensional data (2 spectral bands × 46 dates per year) [46]. For comparison, the overall accuracy of the RF classification applied at the national scale (0.93-0.96) is similar to that obtained using annual MODIS time-series at the regional scale for the Midlands Region of Ireland (0.92) [38] or the Grenoble Alpine region (0.88) [39]. Due to the very high temporal resolution of MODIS images (46 dates/year), the general error rate of classification between grassland and crop classes was low (± 5%) (Table S7). However, the confusion between grassland and crop classes, which was the main source of error at the national scale, reflects the strong similarity in spectral patterns between certain crops and temporary grasslands sown throughout the year [26]. Moreover, the dates of the reference data for the Water & Wetness HRL (2015) and Tree Cover Density HRL (2012 and 2015) differ up to 10 years from the dates of MOD13Q1 and MYD13Q1 products (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). Although forests, rivers and water bodies are relatively stable over time, changes in LULCs could have occurred between 2006 and 2015, which may lead to classification errors.
The use of "pure" pixels for validation may lead to a slight overestimation of the classification accuracy. However, the use of "mixed" pixels would probably underestimate the classification accuracy. In any case, given the large number of samples used per year (± 12,000), the sample quality (farmers' declarations), sample independence (from training samples) and their homogeneous distribution throughout the territory (sub-sampling 5 × 5 km), the estimated classification accuracy is considered reliable. At the regional scale, while grasslands were classified with great accuracy in the Atlantic and Continental biogeographical regions (F1-score 0.90 and 0.94, respectively)-which cover most of France-more errors were observed in the Alpine (F1-score 0.86) and Mediterranean (F1-score 0.62) regions. In the Alpine region, confusion was observed mainly between the grassland and woods classes (± 18%). This could be because many parcels described as "grassland" in this region in the

Can MODIS 250 m Time-Series Combined with the RF Classifier Discriminate Grasslands from Other LULC Types at the National Scale?
Overall, the RF classification of MODIS 250 m time-series accurately mapped grasslands each year from 2006-2017 at the scale of France (F1-score 0.89-0.93). These results confirm the advantage of using the RF classifier to discriminate between grasslands [1], which in France are of different types (temporary, wet, mesic, dry or alpine) and thus have different spectral patterns, as well as to manage high-dimensional data (2 spectral bands × 46 dates per year) [46]. For comparison, the overall accuracy of the RF classification applied at the national scale (0.93-0.96) is similar to that obtained using annual MODIS time-series at the regional scale for the Midlands Region of Ireland (0.92) [38] or the Grenoble Alpine region (0.88) [39]. Due to the very high temporal resolution of MODIS images (46 dates/year), the general error rate of classification between grassland and crop classes was low (± 5%) (Table S7). However, the confusion between grassland and crop classes, which was the main source of error at the national scale, reflects the strong similarity in spectral patterns between certain crops and temporary grasslands sown throughout the year [26]. Moreover, the dates of the reference data for the Water & Wetness HRL (2015) and Tree Cover Density HRL (2012 and 2015) differ up to 10 years from the dates of MOD13Q1 and MYD13Q1 products (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). Although forests, rivers and water bodies are relatively stable over time, changes in LULCs could have occurred between 2006 and 2015, which may lead to classification errors.
The use of "pure" pixels for validation may lead to a slight overestimation of the classification accuracy. However, the use of "mixed" pixels would probably underestimate the classification accuracy. In any case, given the large number of samples used per year (± 12,000), the sample quality (farmers' declarations), sample independence (from training samples) and their homogeneous distribution throughout the territory (sub-sampling 5 × 5 km), the estimated classification accuracy is considered reliable. At the regional scale, while grasslands were classified with great accuracy in the Atlantic and Continental biogeographical regions (F1-score 0.90 and 0.94, respectively)-which cover most of France-more errors were observed in the Alpine (F1-score 0.86) and Mediterranean (F1-score 0.62) regions. In the Alpine region, confusion was observed mainly between the grassland and woods classes (± 18%). This could be because many parcels described as "grassland" in this region in the LPIS include moorland facies with a strong mixture of grassland and woods, which are difficult to discriminate using remote sensing [23]. The Mediterranean region was the least accurately modeled, with a significant under-detection of grasslands. The facies of many Mediterranean grasslands are similar to those of scrubland with a low percentage of ground cover and sometimes the presence of many shrubs [64]. This would explain why some Mediterranean grasslands were under-estimated due to confusion with the urban and woods classes (Table S11). The confusion between crops and grasslands in the Mediterranean region can be explained by the high proportion of vineyards and orchards that have herbaceous vegetation between the rows of vines or fruit trees. The highly distinctive landscape of the Mediterranean region requires a specific RF model with an adapted nomenclature [65].
The quality of MODIS products used in this study could also explain some of the errors in the grassland classification. The modeling was based on the red and NIR bands of the MOD13Q1 and MYD13Q1 products, which may lead a lower modeling accuracy of grasslands in the Alpine biogeographical region compared to those achieved in the Atlantic and continental lowland regions, since these two spectral bands are sensitive to terrain distortions. However, we chose to use the red and NIR bands rather than the NDVI and EVI of the MOD13Q1 product, because these vegetation indices have some limitations: (i) they reduce spectral dimension from two features (red and NIR bands) to one (NDVI or EVI), which can reduce class separability; (ii) the EVI is computed using the 500 m blue band, which can generate artifacts in images at 250 m spatial resolution; (iii) the NDVI index saturates in high values, while grasslands have high reflectance values in spring and summer. Besides, MOD13Q1 and MYD13Q1 products-which are 16-day composite products-contain inconsistencies in pixel values related to differences in incidence angles and albedo values between the images used to create the composite products. These spectral biases can generate errors in LUCL classification-especially in winter, when atmospheric disturbances are more important. More consistent spectral products with observations every 8 days-time series of globally distributed spectral MODIS NBARs (nadir bidirectional reflectance distribution function adjusted surface reflectance)-are available [66]. However, the low spatial resolution of this product (500 m) would not be appropriate for the detection of the majority of grasslands that were detected using the MOD13Q1 and MYD13Q1 products (250 m).
Sampling can also explain the regional differences in the accuracy of the RF model. Despite the application of 5 × 5 km grid sub-sampling, the Atlantic and Continental regions had more samples than the Alpine and Mediterranean regions since they cover the largest area in France (Tables S3, S4, S5 and  S6). As a result, the RF model, which is calibrated with the highest Kappa index [59], discriminated between grasslands in the Atlantic and Continental regions better, since their larger number influenced overall accuracy more than the smaller number of grasslands in the Alpine and Mediterranean regions. The accuracy of RF models is also strongly influenced by the quality of the sampling [46,67]. One of the challenges when analyzing data with moderate spatial resolution, such as MODIS, is selecting enough samples that correspond to "pure" pixels. Samples selected from "mixed" MODIS pixels have a negative influence on classification accuracy [68]. For this reason, we carefully selected only MODIS pixels covering >80% of the same LULC class. This strong sample selection was possible because the LPIS contains a large amount of reference data (i.e., more than 200,000 polygons across France for each year since 2006). Applying this approach outside the EU would be more difficult due to the scarcity of availability of field databases [1,47]. However, using expert-rule classification methods based on annual NDVI time profiles could be an alternative solution to discriminate grasslands from other LULC types in regions where field data are not widely available [47,69].
The grassland maps included in the global MCD12Q1 v6 products were obviously less accurate than the one obtained from MOD13Q1 and MYD13Q1 composites. While the six LULC maps of the MCD12Q1 v6 product were generated by a single RF model applied throughout the planet using about 3,000 reference samples selected by visual interpretation of very high spatial resolution images [40], the MODIS 250 m LULC map was computed at a national scale with one RF model calibrated using about 12,000 reference samples validated in the field. We visually compared these LULC maps on a Natura 2000 site. A thorough statistical comparison performed at a national scale would deserve a dedicated study. The reference data-i.e., pure MODIS pixels at 250 m spatial resolution-produced over the period 2006-2017 from LPIS are publicly available [57] and can be integrated for future global LULC classification. To our knowledge, this is the first study conducted to identify "permanent" grasslands using satellite remote sensing for such a long period (12 years); in comparison, Garcia-Feced et al. [6] and Lasseur et al. [39] monitored grasslands over a 5 or 4-year period, respectively. Given the temporal depth of the analysis (12 years), grasslands identified as permanent (i.e., a frequency close to 1) have a high probability of being semi-natural [2]. Nevertheless, certain 12-year permanent grasslands may be long-term "temporary" grasslands that have not yet restored the floristic composition of semi-natural grasslands [3]. Adding new MODIS images to the time-series (>2017) would yield a longer time-series (20 years) that could address this uncertainty.
The two major limitations of current national and international LULC maps for monitoring semi-natural grasslands are their (i) coarse spatial resolution and (ii) low thematic resolution, which hinder clear correspondence between semi-natural and temporary grasslands [13]. The national grassland frequency map addresses these limitations and improves knowledge about semi-natural grasslands in two ways. First, its 250 m resolution identifies semi-natural grasslands more accurately than maps produced from LPIS or CORINE Land Cover at a spatial resolution of 5 km [22] or from the EVA database at a spatial resolution of 10 km [19]; second, it uses long-term monitoring to discriminate between semi-natural and temporary grasslands, which is something that LULC maps derived from annual time-series analysis of remote sensing data do not do [23,25,30].
However, this grassland frequency map should be read carefully. While interpreting the extreme values (0 and 1) is straightforward, this is less true for intermediate values. First, the analysis does not consider damage to semi-natural grasslands from 2006-2017 that would affect the grassland frequency. For example, a semi-natural grassland converted into an urban area in 2011 would be classified as grassland only from 2006-2010 and would have a frequency of ca. 0.5. Second, although annual grassland modeling was highly successful (F1-score 0.89-0.93), the few errors generated each year may have accumulated when calculating grassland frequency as a aggregation of annual LULC maps [39,70]. Given these issues, the grassland frequency map should be used to help pre-locate semi-natural grasslands rather than as an inventory per se. For example, pre-locating semi-natural grasslands could improve implementation of restoration strategies for agricultural systems and biodiversity in Europe [6].

Is the 250 m Spatial Resolution of MODIS Data Adequate for Identifying Grasslands in Fragmented Landscapes?
Results of the sub-pixel analysis for 2016 reveal that MODIS images at 250 m spatial resolution successfully identified grassland-dominated areas. These results confirm that MODIS images can detect grasslands in homogeneous landscapes [36,39], such as mountain pastures, and in more fragmented landscapes such as the Armorican massif. Although grassland percentage predominated in all HERs (60%-94%), it may be more appropriate to use the term "grassland landscape" rather than "grassland". A MODIS pixel classified as grassland generally included 15% of crop area. Conversely, it is likely that some "isolated" and small semi-natural grasslands were classified as crop or woods. Although these results are based on a large sample (ca. 18% of the area of France), they may be biased slightly by (i) the quality and spatial resolution of the reference layers (HRL, LPIS), which do not represent small elements of LULC (e.g. hedges, small water bodies, roads), and (ii) the low-but existing-rates of modeling error (producer's accuracy 87.8%, user's accuracy 97.8%).

Conclusions
Analysis of MODIS 250 m time-series from 2006-2017 (12 years) for the whole of France revealed the presence of permanent grasslands. This map is a new aid for pre-locating semi-natural grasslands at the national scale. Results demonstrate that a single RF model correctly discriminates between different types of grassland in the Atlantic, Continental and Alpine biogeographical regions, but identifying semi-natural grasslands in the Mediterranean region requires a specific RF model and nomenclature. Although the spatial resolution of MODIS data is coarse compared to parcels sizes, sub-pixel analysis highlights that areas modeled as grassland correspond to grassland-dominant areas. Perspectives of this approach include (i) characterizing grasslands by temporally analyzing their inter-annual spectral profiles using deep-learning approach, (ii) identifying grasslands at a finer scale using Sentinel-1/2 time-series and (iii) monitoring grasslands over a longer period and for the whole of Europe.