Progress in Grassland Cover Conservation in Southern European Mountains by 2020: A Transboundary Assessment in the Iberian Peninsula with Satellite Observations (2002–2019)

Abstract

Conservation and policy agendas, such as the European Biodiversity strategy, Aichi biodiversity (target 5) and Common Agriculture Policy (CAP), are overlooking the progress made in mountain grassland cover conservation by 2020, which has significant socio-ecological implications to Europe. However, because the existing data near 2020 is scarce, the shifting character of mountain grasslands remains poorly characterized, and even less is known about the conservation outcomes because of different governance regimes and map uncertainty. Our study used Landsat satellite imagery over a transboundary mountain region in the northwestern Iberian Peninsula (Peneda-Gerês) to shed light on these aspects. Supervised classifications with a multiple classifier ensemble approach (MCE) were performed, with post classification comparison of maps established and bias-corrected to identify the trajectory in grassland cover, including protected and unprotected governance regimes. By analysing class-allocation (Shannon entropy), creating 95% confidence intervals for the area estimates, and evaluating the class-allocation thematic accuracy relationship, we characterized uncertainty in the findings. The bias-corrected estimates suggest that the positive progress claimed internationally by 2020 was not achieved. Our null hypothesis to declare a positive progress (at least equality in the proportion of grassland cover of 2019 and 2002) was rejected (X2 = 1972.1, df = 1, p < 0.001). The majority of grassland cover remained stable (67.1 ± 10.1 relative to 2002), but loss (−32.8 ± 7.1% relative to 2002 grasslands cover) overcame gain areas (+11.4 ± 6.6%), indicating net loss as the prevailing pattern over the transboundary study area (−21.4%). This feature prevailed at all extents of analysis (lowlands, −22.9%; mountains, −17.9%; mountains protected, −14.4%; mountains unprotected, −19.7%). The results also evidenced that mountain protected governance regimes experienced a lower decline in grassland extent compared to unprotected. Shannon entropy values were also significantly lower in correctly classified validation sites (z = −5.69, p = 0.0001, n = 708) suggesting a relationship between the quality of pixel assignment and thematic accuracy. We therefore encourage a post-2020 conservation and policy action to safeguard mountain grasslands by enhancing the role of protected governance regimes. To reduce uncertainty, grassland gain mapping requires additional remote sensing research to find the most adequate spatial and temporal data resolution to retrieve this process.

1. Introduction

Grasslands are among the most threatened ecosystems in the world due to conversion, desertification and degradation [1]. In Europe, mountain grasslands are particularly endangered [2,3]. These high diversity and species-rich ecosystems that coevolved with centuries of livestock grazing have been declining since the 1960′s [4]. The decline of mountain grasslands has multiple impacts on the planet. Locally, the unique biodiversity [5], cultural heritage [2] and ecosystem services [6] are compromised. Regionally, the set of pressures once created by agro-pastoral activities and commodities (e.g., meat production) is moved to lowlands, and globally, very sensitive ecological regions are degraded (e.g., tropical regions) to replace the set of commodities lost in mountains [7]. Given the impacts, halting the loss of mountain grasslands has become a priority for international and European agendas. This is a priority of the Aichi Biodiversity Agenda Target 5 “Habitat loss and degradation”, European Biodiversity agenda (Target 3 and 6) [8] and agri-environment schemes and payments of the Common Agriculture Policy (CAP; [9,10]). Now, in 2021, and with a new set of working documents on the table (e.g., European Strategy for Biodiversity 2030), what is the progress made in the conservation of mountain grasslands by 2020? Understanding the progress made is necessary for policy action and policy making.

Several studies have been conducted in mountain grassland conservation in Europe in the last decades [5,6,11,12,13]. Surprisingly, very few comprised a period of analysis close enough to the international target year of 2020 (e.g., 2015–2020). The current understanding of progress made by 2020 relies on a handful of studies [14,15,16,17], with most focused on major mountain regions (e.g., Carpathians, Apennines), with smaller ranges (e.g., northwestern Iberian mountains) less considered. A lag in time and space seems to exist between scientific and policy agendas. This may be related to the strong emphasis of the scientific community on the capitalization of nature to establish safeguards (e.g., ecosystem services). Hence, the data available for the period 2015–2020 needs to be enriched. Based on a set of studies near to 2020 (2015–2020), loss or decline was the prevalent pattern in mountain grasslands. In the Apennines, grasslands decreased by 4.6% between 1987–2016 at the expense of forest expansion after land abandonment [17]. This trajectory in boreal seminatural grasslands was more dramatic. The extent of seminatural grasslands in 2015 was 49.1% less compared to 1960 due to conversion into arable land and forests [14]. Jaworek-Jakubska et al. [16] in the Polish mountains indicated a 10.4% decline of grasslands between 2000–2016 in favor of forests. In contrast, a remarkable increase of managed grasslands was found between 1999–2015 (+39%) in the Carpathians [15]. More direct payments of CAP for mountain grasslands were identified as the main reason for this trajectory. However, Halada et al. [15] emphasized that conservation of mountain grasslands remains a challenge because of their poor accessibility.

Regarding the progress in mountain grassland conservation, another two gaps seem to exist. The first is the performance of different governance regimes (e.g., protected areas, PAs). For instance, it was not possible to identify studies assessing how protected governance regimes affect grassland cover conservation. However, state-controlled PAs are a main global conservation strategy to address concerns, such as those relating to the endangered mountain grasslands of Europe [18]. They should be at the front line of conservation action and political commitment to attenuate negative tendencies [19,20]. The second gap relates to the use of maps as source data for conservation progress analysis. Most of the studies assume mapped areas as true areas [14,15,17], whereas maps have errors (e.g., commission, omission) that can affect the estimates and, therefore, the trajectories observed in international conservation targets. Uncertainty in grassland change studies remains poorly documented. Nowadays, good practices recommend the provision of uncertainty estimates. The conventional error matrix of producer’s and user’s accuracy assessments cannot indicate the level of uncertainty in areal terms. Nonetheless, this approach can be extended to create unbiased area estimates with associated confidence levels (95% confidence interval) [21,22]. Overall, knowledge in the performance of protected governance regimes is necessary and of interest to policy action and PAs to support their missions [23]. Unbiased area estimates of mapped areas, as already in use in deforestation and carbon emissions issues [24], can provide a more complete view e.g., for the two decades of European and international policy commitment with mountain grasslands protection.

Multispectral satellite data and image classification play an important role in the assessment of progress made in mountain grasslands conservation. They allow increasingly accurate land cover and vegetation change maps [25,26] that can support the quantification of spatial targets at multiple spatial [13] and temporal extents [26]. For this, the expansion of open satellite data at moderate resolution (e.g., 10 to 30 m, Sentinel-2, Landat), of the set of classifiers (e.g., random forests) and sources of actual or near-ground information (e.g., platform with free access to very high resolution imagery) in the last decade has provided a substantial contribution. Nonetheless, a review of remote sensing literature shows that grasslands have been relatively poorly studied compared to other land cover types (e.g., forests, [27]). Furthermore, the heterogeneity of grasslands in fragmented landscapes, as is common to European mountains, remains the major challenge when classifying from remote sensing data. To reduce classification issues, an important opportunity may exist in the multiple classifier ensemble (MCE) approach. That provides the opportunity to integrate independent classification results of different algorithms (e.g., those based on artificial intelligence, such as random forests, k-nearest neighbour, decision trees) through ensemble strategies (e.g., soft voting). Improvement of the accuracy of classification results compared to single classifiers has been observed [28,29]. MCE also yields class-allocation probability layers at the pixel level, which can be used to generate class-allocation uncertainty information (e.g., Shannon entropy, [30]) and support refinement operations before map comparisons, reporting pixel assignment confidence and analyses of the relationship between class-allocation uncertainty and classification accuracy.

This study conducted mapping and analysis of bias-corrected grassland cover data based on Landsat remote sensing imagery (2002–2019, 30 m resolution) over a transboundary mountain region of the northwestern Iberian Peninsula (Peneda-Gerês) to address three main questions:

(1)

What is the progress made in grassland cover conservation by 2020 according to the extent and net shift observed between 2002 and 2019?

(2)

To what extent do protected and unprotected land governance regimes differ in grassland cover conservation?

(3)

What is the accuracy level of MCE approaches in mountain grasslands? What is the relationship between class-allocation uncertainty and classification accuracy?

The main goal was to reduce the data and critical gap regarding the shifting character of European mountain grasslands towards the international target established by 2020. Grasslands of the study area are representative of ancient mountain pastoral systems of the Iberian Peninsula, which are rapidly transforming with farmland change and abandonment, and are a target of European policies [13]. Specifically, we used Landsat imagery and derivatives within a multiple classifier ensemble approach [28] to map grassland cover in 2002 and 2019. Post classification comparison of uncertainty refined grassland cover maps to identify changed areas. Accuracy and uncertainty assessments were made to characterize the confidence of pixel assignment (Shannon entropy) and create bias-corrected change areas estimates with 95% confidence intervals following the good practices protocol of Olofsson et al. [22]. Using these data, the progress made in grasslands cover conservation by 2020 was summarized at four extents: transboundary, lowlands, mountains protected and mountains unprotected. The trajectories in protected and unprotected governance regimes were compared to identify the effect, if any, of governance regime in grassland conservation. The relationship between class-allocation uncertainty and classification accuracy was analysed. We conclude this study by discussing the alignment of the progress registered with the European and international policy ambitions by 2020 and the uncertainty as a major feature of mountain grasslands classification.

2. Materials and Methods

2.1. Study Area

The study area was the Peneda-Gerês transboundary mountain region in the northwestern Iberian Peninsula (41°42′59″ N, 8°08′60″ W, 6352.92 km², see Figure 1). Here, the elevation ranges from 15 to 1513 m above sea level (a.s.l) and steep topography is a main feature. The climate is temperate oceanic sub-Mediterranean with mean annual precipitation of 1300 mm and mean annual temperature of 13 °C [13]. The study area includes a lowlands domain (375,759 ha), and a mountain domain corresponding to the transboundary Gerês-Xurés Biosphere reserve (259,533 ha). Within the mountain domain are protected (Peneda-Gerês National Park, Portugal; Baixa Limia-Serra do Xurés Natural Park, Spain) and unprotected land conservation regimes (189,932 ha). In this study, the mountain protected conservation regime used for analysis refers to the Portuguese Peneda-Gerês National Park (69,543 ha) established at the end of the authoritarian period (1972), which makes it an excellent location to understand governance regime effects. The Peneda-Gerês National Park is managed according to specific land use planning and has renown in natural and cultural heritage as the last remaining native forest ecosystems enclosed in man-made farmland mosaics with extensive human influence [31]. At the socio-economic level, this transboundary study area provides testimony to the three major socio-economic disruptions that have occurred in the Iberian Peninsula in the last half-century: the end of the authoritarian periods and establishment of democracy (1974), the Europeanization of society with adhesion to the European Union (EU, 1986) and the achievement of a Europeanized society with entrance in the EuroZone (1999). The Europeanization of society brought strong socio-economic growth supported by European funding, convergence of environmental policy with EU directives (e.g., new law for environment) and the Common Agriculture Policy (CAP) as a player in the farming systems of the region.

2.2. Workflow

This study used imagery from the Landsat 7 and 8 satellite sensors over a transboundary mountain region in the Iberian Peninsula to detect grassland cover change, characterize the progress made in mountain grasslands conservation by 2020, and evaluate its variation among protected and unprotected land governance regimes. The overall methodology is shown in Figure 2, and includes four main components. (a) Mapping of grassland cover in 2002 and 2019. Two Landsat satellite scenes were selected and a supervised MCE approach [28] that combined three different artificial intelligence algorithms (random forests, decision trees and K-nearest neighbours) was applied to discriminate grassland cover. This resulted in maps of grasslands cover (n = 2) and probability layers (n = 2) for each class (grasslands, no grasslands) and for each date. (b) Cleaning and change analysis, a two-step procedure that first removed low probability grassland pixels from the class using the U uncertainty metric (U > 0.2) and later performed a post classification comparison of uncertainty-refined maps to create the grassland cover change map (2019–2002). (c) Accuracy and uncertainty assessment using the Shannon entropy and the bias-corrected error matrix approach of Olofsson et al. [22]. Shannon evaluated the uncertainty in the class assignment, while the error matrix assessed the accuracy of the grassland cover change map and created bias-corrected area estimates for grassland cover and change areas using reference data. (d) Statistical analysis to evaluate the progress made in the conservation of grasslands cover, the effect of land governance regime and the relationship between cumulated class-allocation uncertainty and classification accuracy.

2.2.1. Landsat Imagery and Ancillary Derivatives

Two Landsat satellite scenes covering the study area were obtained from the platforms Landsat 7 (26 March 2002) and Landsat 8 (13 February 2019). The scenes were downloaded from USGS as surface reflectance, geometrically corrected and stored in Geotiff format. They were selected based on the following criteria: (i) cloud-free scenes taken in the (ii) late winter-early spring season, and (iii) in all years of the analysis (anniversary-date scenes), in order to optimize mapping and change detection and to reduce errors related to seasonal differences [32]. In the study area, grasslands have a strong reflectance in the late winter-early spring season that clearly distinguish them from other vegetation types in winter dormancy, namely from early succession types (e.g., ferns). Due to the temporal lag between the scenes, radiometric normalization using the pseudo-invariant features technique was performed [33]. Normalization used the 2019 scene as the reference image to adjust the 2002 scene. For each scene, we also computed the Tasseled Cap components greenness (TCG), wetness (TCW) and fifth (TCF). TC components for each date were derived using the equations given in [21,34].

2.2.2. Supervised Classification of Grassland Cover

To map grassland cover, we performed supervised classification with an MCE approach [28] that through soft voting combines the predicted probabilities of three classifiers: random forests (RF), decision trees (DT) and k-nearest neighbors (KNN). MCE approaches have previously been used in land cover mapping and have generally increased the accuracy of classifications with respect to single classifier approaches [28,29]. As input data, we used the Landsat bands (Green, NIR, SWIR) and the Tasseled cap components (TCG, TCW and TCF) created above. TCG was chosen because it maximised capture of the high reflectance of grasslands in the late winter-early spring season, which provided good separation during this period. TCW and TCF were chosen based on the improvements demonstrated during grassland mapping in European landscapes [35]. Regarding the classification scheme, two classes were considered (grasslands, no grasslands). Grasslands included natural and managed grasslands that can be mown, grazed or ploughed. The training data consisted of a total of 131 training sites recognizable as grasslands (n = 82) or no grasslands (n = 49) in 2002 and 2019 (see Figure 1,

Table 1. Parameterization of the Multiple Classifier Ensemble approach and composition of the reference training data.

Classification		Reference Training Data
Method	Classifier Parameterization	Class	Training Sites	N° of Pixels
Multiple Classifier Ensemble (MCE)	Random Forests (n_estimators = 500, criterion = ‘gini’, max_depth = 4, min_samples_split= 2, min_samples_leaf = 1, max_features = ‘auto’, bootstrap = True, oob_score = True); Decision Trees (criterion = ‘entropy’, max_depth = 4, min_samples_leaf = 1, min_samples_split = 2, random_state = 100); K-nearest neighbors (n_neighbors = 4, weights = ‘distance’, leaf_size = 1);	Grasslands	82	272
		No Grasslands	49	579
		Total	131	851

). Each training site was a homogeneous polygon of grasslands or no grasslands, with these selected using previous field data and knowledge of the study area combined, and through visual inspection of very high resolution imagery available in Google Earth [36]. Polygons were rasterized at 30 m resolution, and each training site was composed of at least two pixels. The training site collection was a time-consuming procedure due to the small size, high degree of fragmentation and heterogeneity of grasslands parcels that is typical of Mediterranean mountain environments. Before classification, the parameters of each classifier were tuned using a grid-search procedure with a 10-fold cross validation of the training samples (see Appendix A, [37]). Classifiers were parameterized as presented in Table 1. This classification procedure generated grassland cover and probability layers for each class and for each date. All operations were performed in Python using the packages scikit-learn, gdal and numpy [38,39].

2.2.3. Cleaning Phase and Map of Grasslands Cover Change

To detect changes in grasslands, we performed a two-step approach. First, we implemented a cleaning procedure that removed low probability grassland pixels from the grasslands cover maps obtained above. Figure 3 illustrates the procedure performed. The U uncertainty measure (U = 1–pmax) was used to identify low probability grassland pixels [30]. U indicates the strength of the class assignment. A high U value indicates a confusion between the grasslands and no grasslands classes. The estimation of U relies on the degree of probability that a certain belongs to the grasslands class (pmax), which was obtained from the modelled probability vector created for grasslands (pu) during the ensemble classification. We set a limit that all pixels with U > 0.2 would be removed from the grassland class and assigned to the no grassland class. This threshold follows Loosvelt et al. [30] that found the majority of incorrectly classified pixels were associated with U values larger than 0.2. Later, we performed a post classification comparison with the “cleaned” maps to create grassland cover change maps (2019–2002). Each pixel was identified as either having changed or not changed. This process resulted in a four-class change map (stable grassland, no grassland, grassland loss, grassland gain).

2.2.4. Accuracy and Uncertainty Assessment: Pixel Class Assignment and Estimation of Bias-Corrected Grasslands Cover Change Areas

Reporting the progress made in grassland cover conservation and its uncertainty, as described for forest cover, is still not a practice [40]. Since we aimed to provide a more complete view for the almost two decades of European and international policy commitment with mountain grassland protection (2002–2019), we analysed accuracy and uncertainty in the classification and change detection phases. In the classification phase, we assessed the confidence of pixel class assignment using the Shannon entropy metric (H, Equation (1), [41]). H is a weighted uncertainty measure that differs from the U uncertainty measure by weighting the probabilities obtained by each class in a certain pixel to provide an overview of pixel reliability. H is derived from the probability vector (pu) created by the MCE classifier [28] that holds the probability p(i) of being classified into class i (pu = (p(1), p(2), …, p(c)), with c being the total number of land cover classes (here c = 2). A pixel has minimal entropy (H = 0) with a maximum probability (pmax = 1), while maximal entropy (H = 1) occurs when all classes have the same probability, e.g., p(i) = 1/c, and no specific classis preferred (grasslands or no grasslands). The uncertainty about the pixel’s true class is maximal. The Shannon entropy uncertainty metric is a per-pixel layer.

$$H=-{{\displaystyle \sum }}_{i=1}^c\left[p\left(i\right)·log\left(p\left(i\right)\right)\right],$$

(1)

In the change detection phase, we assessed the accuracy of the grassland cover change map and generated bias-corrected area estimates with 95% confidence intervals [22]. For that purpose, we followed the good practices protocol of Olofsson et al. [22]. This approach uses stratified random sampling to define validation sites where the confrontation between classified and reference classes generates an error matrix. The error matrix is then used to obtain accuracy estimates in terms of proportions of area, which are further used to compute the bias-corrected area estimates for each of the four-class changes at a 95% confidence level. The bias-corrected area estimate is obtained by including the area of map omission error and by leaving out the area of map commission error [22]. Consequently, for the grassland change map, an error matrix was produced in terms of proportions of area. When map categories are in the rows (i) and the reference categories in the columns (j), Atot refers to the total area of the map, Am,I is the mapped area (ha) of category i in the map and Wi = Am,iAtot is the proportion of the mapped area as category I and pij is then as given in Equation (2).

$${\widehat{p}}_{ij}=W_i\frac{n_{ij}}{n_{i }},$$

(2)

The error-adjusted area estimate (${\widehat{A}}_j$) of a category is obtained as in Equation (3).

$${\widehat{p}}_{ij}=W_i\frac{n_{ij}}{n_{i }},$$

(3)

The estimated standard error of the estimated proportion of area $S\left({\widehat{p}}_{ j}\right)$ is given by Equation (4).

$$S\left({\widehat{p}}_{ j}\right)=\sqrt{{\displaystyle \sum }_{i=1}^q W_i^2\frac{\frac{n_{ij}}{n_{i }}\left(1-\frac{n_{ij}}{n_{i }}\right)}{n_{i }-1}.}$$

(4)

Later, the standard error of the stratified area estimate, and the 95% confidence interval for ${\widehat{A}}_j$ are derived as in Equations (5) and (6), respectively.

$$S\left({\widehat{A}}_j\right)=A_{tot}\times S\left({\widehat{p}}_{ j}\right),$$

(5)

$${\widehat{A}}_j\pm 2\times S\left({\widehat{A}}_j\right),$$

(6)

A set of 708 validation sites (see Figure 1) was defined using a stratified random sampling with unequal sample (stable grasslands (n = 43), stable no-grasslands (n = 608), grassland loss (n = 32) and grassland gain (n = 25)) informed by the grassland change maps. For each validation site, the class was assigned through visual inspection of very high resolution imagery available in Google Earth [36]. This procedure was initially applied to the entire study area. Later, however, in order to create bias-corrected area estimates also for the remaining extent of analysis (lowlands and mountains domain, protected and unprotected mountain extents), we made the assumption that the pattern of map omission and commission error observed across the study area would be valid also for each scale of analysis. Therefore, using the error matrix of the study area (n = 708), we created area estimations for each extent of analysis by updating only the mapped areas.

2.3. Statistical Analysis

To assess the progress made in the conservation of mountain grasslands in the transboundary study area, the grassland covers of 2019 and 2002 were compared using a two-proportion Z-test with continuity correction that produced a two-tailed probability value that evaluated the null hypothesis of equality of proportions in the two times. The two-proportion Z-tests were separately made for the entire study area (635,292 ha), lowlands (375,759 ha) and mountain domains (259,533 ha), protected (69,543 ha) and unprotected governance extents (189,932 ha). If differences were significant (p < 0.05), positive progress had been achieved in the case where the proportion of grassland cover was at least maintained or expanded between 2002 and 2019. Conversely, progress was negative if the proportion of grassland cover in 2019 was lower than in 2002.

To assess the impact, if any, of protected and unprotected governance regimes in the progress made, we also used a two-proportion Z-test that determined if the net loss of grasslands fitted a null hypothesis of equal distribution between protected and unprotected governance extents. In the study area, the current protected extent is larger than that used in this statistical analysis. The protected extent of this analysis refers only to the Peneda-Gerês National Park (69,543 ha) in the Portuguese territory. We focused only on this protected extent because it has been protected since 1971 or over the entire period of analysis. The protected extent in the Spanish territory was only created in 2008. The function prop.test of R the statistical software was used to perform the two-proportion tests.

To analyse the relationship between class-allocation uncertainty and classification accuracy, we grouped the results of validation (n = 708) into correct classified and misclassified, and estimated the cumulated Shannon entropy (Shannon entropy 2002 + Shannon entropy 2019) for each site using the Shannon entropy layers created in Section 2.2.4. The comparison of cumulated Shannon entropy values of groups was made using the nonparametric Wilcoxon-Mann-Whitney test with 10,000 sample Monte Carlo distribution approximations implemented in the R coin package [42].

3. Results

3.1. Grassland Cover Change in the Peneda-Gerês Transboundary Mountain Region

This study used bias-corrected change estimates to understand grassland cover change, which are estimates obtained from the adjustment of mapped areas for the omission and commission errors found during the accuracy assessment (Figure 4). In the period 2002–2019, the bias-corrected estimates revealed considerable changes over the transboundary mountain region of Peneda-Gerês in the northwestern Iberian Peninsula (Figure 5). Grasslands experienced changes at all extents of analysis. Across the entire study area, grasslands have always been a non dominant cover type (2019–2002; no grasslands = 558,345 ± 7829 ha, 88.7% of the study area,

Table 2. Total estimated bias-corrected areas for the stable grasslands and no grasslands classes with margin of error at the 95% confidence level (ha) at all extents of analysis: study area, lowlands and mountain domains, protected and unprotected mountain governance extents. The standard error and the 95% bias-corrected areas (ha) were derived from the error matrix according to the good practices protocol of Olofsson et al. [22].

Period (2019–2002)	Mapped Area (ha)		Bias-Corrected (ha)		Standard Error of Bias-Corrected Area Estimate (ha)		Confidence Interval (95%; ha)
	Stable Grasslands	No Grasslands	Stable Grasslands	No Grasslands	Stable Grasslands	No Grasslands	Stable Grasslands	No Grasslands
Study area	37713	546995	46378	558345	3866	3995	7577	7829
Lowlands	26959	310315	32849	319598	2694	2756	5280	5402
Mountains	10752.9	236680	13259	238747	1228	1297	2406	2542
Mountain protected	885	67345	1603	67173	264	288	518	564
Mountain unprotected	9869	169335	11926	171574	995	1041	1950	2040

) and 46,378 ± 7577 ha of grasslands remained stable (67.1 ± 10.1% relative to 2002 grassland cover). Grassland loss occurred in 22,695 ± 4945 ha (−32.8 ± 7.1% relative to 2002 grassland cover) and grassland gain occurred in 7875 ± 4587 ha (+11.4 ± 6.6%), resulting in a net loss of grassland cover of 14,820 ha (−21.4%) over the study area. When narrowing the analysis into lowland and mountain domains, we observed a similar trajectory. Grassland loss in lowlands reached 17,428 ± 3706 ha (−34.6 ± 7.3% relative to 2002 grassland cover in lowlands) and grassland gain occurred in 5884 ± 3435 ha (+11.7 ± 6.8%), resulting in a net loss of grassland cover of 11,544 ha (−22.9%) in the lowlands of the study area. Here, 13,529 ± 2406 ha of grasslands remained stable (65.3 ± 10.5% relative to 2002 grassland cover in lowlands). In the mountain domain, grassland loss occurred in 5267 ± 1302 ha (−28.0 ± 6.9% relative to 2002 grassland cover in the mountain domain) and grassland gain occurred in 1991 ± 1156 ha (+10.5 ± 6.1%), resulting in a net loss of grassland cover of 3276 ha (−17.4%) over the mountain domain. Stable grassland amounted to 71.9 ± 12.8% relative to 2002 grasslands cover in the mountain domain. Across the governance regimes (mountain unprotected and protected), the net loss pattern prevailed. In the unprotected extent, grasslands loss occurred in 4697 ± 1094 ha (−28.2 ± 6.5% relative to 2002 grasslands cover in the unprotected extent) and grasslands gain occurred in 1725 ± 1008 ha (+10.4 ± 6.1%), resulting in a net loss of grassland cover of 2962 ha (−19.7%). Stable grasslands amounted to 11,926 ± 1950 ha (71.7 ± 11.7% relative to 2002 grassland cover in the unprotected extent). In the protected area, grassland loss reached 570 ± 242 ha (−26.2 ± 11.1% relative to 2002 grasslands cover) and grassland gain occurred in 256 ± 150 ha (+11.7 ± 6.8%), resulting in a net loss of grassland cover of 314 ha (−14.4%). Stable grassland amounted to 1603 ± 518 ha (73.7 ± 22.8% relative to 2002 grasslands cover in the protected area.

3.2. The Progress in Mountain Grasslands Conservation by 2020 and the Difference among Land Governance Regimes

To understand the progress made in mountain grasslands conservation by 2020, two-proportion Z-tests compared the bias-corrected grasslands cover estimates of 2019 and 2002 in each extent of analysis (

Table 3. Results of the two proportion Z-tests that compared: the proportions of grassland cover in 2019 and 2002 across the different extents of analysis to analyse the progress made in mountain grasslands conservation by 2020 (upper table); and the grasslands net loss between mountains protected and unprotected land governance regimes to understand to what extent does the protected governance regime affect grassland cover change and attenuate negative conservation tendencies (lower table).

Extent of Analysis	Grasslands Cover (ha)		Sample Size (ha)	Pearson’s Chi-Squared (X2)	P (Two-Tailed)
	Year 2019	Year 2002
Study area	58159	72401	635292	1972.1	p < 0.001
Lowlands	38733	50277	375759	1698	p < 0.001
Mountains	15519	18796	259533	334.9	p < 0.001
Mountains protected	1859	2173	69602	25	p < 0.001
Mountains unprotected	13661	16623	189932	314.6	p < 0.001
	Grasslands Net Loss (ha)		Sample Size (ha)	Pearson’s Chi-Squared (X2)	P (Two-Tailed)
Mountains protected	314		2173	14.9	p < 0.001
Mountains unprotected	2962		16623

; e.g., study area, lowlands, mountains, unprotected and protected governance). This suggested that the positive progress claimed internationally by 2020 was not achieved. The differences in grasslands cover were significant (p < 0.05), but decline was the prevalent pattern and, therefore, negative progress was observed. A positive progress would be achieved where the proportion of grassland cover of 2019 was at least equal to 2002 (null hypothesis). The equality of proportions hypothesized was rejected at all extents of analysis. In 2019, grassland cover was significantly lower with respect to 2002 in the study area (X² = 1972.1, df = 1, p < 0.001), lowlands (X² = 1698.0, df = 1, p < 0.001), mountains (X² = 334.9, df = 1, p < 0.001), protected (X² = 25.0, df = 1, p < 0.001) and unprotected governance extents (X² = 314.6, df = 1, p < 0.001). To assess the impact, if any, of protected and unprotected governance regimes in the progress made, we evaluated if the net loss of grasslands fitted a null hypothesis of equal distribution between protected and unprotected governance extents (Table 3). This showed that proportions of net loss of grasslands were significantly lower in the protected governance regime (X² = 501.2, df = 1, p < 0.001) suggesting that protected land areas are attenuating the negative tendency of mountain grasslands decline.

3.3. The Confidence of Grasslands Pixel Assignment and Uncertainty in Grasslands Cover Change Estimates

The Shannon entropy (H) examined the pixel class assignment performed by the MCE approach that resulted in uncertainty maps that locate the strengths and weaknesses in the grasslands class assignment in 2019 and 2002 (Figure 6). A grasslands pixel with a low Shannon values indicates confidence in that assignment. Our Shannon uncertainty map, and the shape of frequency distributions of values, indicated that the majority of grassland pixels values were below the threshold of 0.2 (mean value in the dashed line). The 0.2 threshold was associated with accurate pixel assignment [30]. In 2019, 67.3% of the pixels assigned to the grasslands class registered values below this threshold (mean value = 0.18). In 2002, 79.2% of grassland pixels were below this threshold (mean value = 0.16). Therefore, in 2019 and 2002 there were 32 and 28% of grassland pixels where the grasslands class was more uncertain based on the Shannon entropy.

While the Shannon entropy stated the confidence in the class assignment in each map, the error matrix assessment following the good practices protocol of Olofsson et al. [22] confronted the mapped grasslands cover changes in relation to reference data. Across the study area, the grassland cover change map had an overall accuracy of 0.95 ± 0.02 (

Table 4. Estimated error matrix for each extent of analysis (study area, lowlands, mountains, mountain protected and mountain unprotected) with cell entries expressed as the estimated proportion (p_ij) of area in cel ij. Overall, user’s and producer’s accuracies with 95% confidence intervals. Appendix B gives the error matrix as sample counts (n_ij, left) and estimated proportions (p_ij, right).

	Study Area	Reference categories
Map categories		No grasslands	Stable grasslands	Grasslands gain	Grasslands loss	Total	Wi	User’s	Producer’s	Overall
	No grasslands	0.853	0.007	0.000	0.001	607	0.86	0.99	0.97	0.95 ± 0.02
	Stable grasslands	0.008	0.051	0.000	0.000	42	0.06	0.86	0.70
	Grasslands gain	0.014	0.010	0.010	0.001	25	0.03	0.28	0.79
	Grasslands loss	0.004	0.005	0.003	0.033	34	0.04	0.74	0.92
	Total	0.879	0.073	0.012	0.036	708	1.00
	Lowlands	Reference categories
Map categories		No grasslands	Stable grasslands	Grasslands gain	Grasslands loss	Total	Wi	User’s	Producer’s	Overall
	No grasslands	0.818	0.007	0.000	0.001	607	0.826	0.99	0.96	0.93 ± 0.02
	Stable grasslands	0.010	0.061	0.000	0.000	42	0.072	0.86	0.70
	Grasslands gain	0.017	0.012	0.012	0.002	25	0.044	0.28	0.78
	Grasslands loss	0.005	0.007	0.003	0.043	34	0.059	0.74	0.93
	Total	0.851	0.087	0.016	0.046	708	1.00
	Mountains	Reference categories
Map categories		No grasslands	Stable grasslands	Grasslands gain	Grasslands loss	Total	Wi	User’s	Producer’s	Overall
	No grasslands	0.903	0.008	0.000	0.002	607	0.912	0.98	0.98	0.96 ± 0.01
	Stable grasslands	0.006	0.036	0.000	0.000	42	0.041	0.86	0.68
	Grasslands gain	0.009	0.006	0.006	0.001	25	0.022	0.28	0.81
	Grasslands loss	0.002	0.003	0.001	0.018	34	0.024	0.74	0.88
	Total	0.92	0.052	0.008	0.002	708	1.00
	Mountain Protected	Reference categories
Map categories		No grasslands	Stable grasslands	Grasslands gain	Grasslands loss	Total	Wi	User’s	Producer’s	Overall
	No grasslands	0.958	0.008	0.000	0.002	607	0.968	0.99	0.99	0.97 ± 0.01
	Stable grasslands	0.002	0.011	0.000	0.000	42	0.013	0.86	0.47
	Grasslands gain	0.005	0.003	0.003	0.000	25	0.011	0.28	0.87
	Grasslands loss	0.001	0.001	0.001	0.006	34	0.008	0.74	0.75
	Total	0.965	0.023	0.004	0.008	708	1.00
	Mountain Unprotected	Reference categories
Map categories		No grasslands	Stable grasslands	Grasslands gain	Grasslands loss	Total	Wi	User’s	Producer’s	Overall
	No grasslands	0.883	0.007	0.000	0.001	607	0.892	0.99	0.99	0.96 ± 0.01
	Stable grasslands	0.007	0.045	0.000	0.000	42	0.052	0.86	0.71
	Grasslands gain	0.011	0.007	0.007	0.000	25	0.026	0.28	0.81
	Grasslands loss	0.003	0.004	0.002	0.022	34	0.030	0.74	0.90
	Total	0.903	0.063	0.009	0.025	708	1.00

), but the accuracy of individual classes varied considerably. Most of the confusion occurred between the grassland gain class, the no grasslands (pij = 0.014, n = 10) and stable grasslands (pij = 0.01, n = 7) classes (Appendix B gives the error matrix as sample counts). The producer’s and user’s accuracy revealed that the change map tended to omit the stable grasslands class (user’s: 0.87, producer’s: 0.70) and overestimate grassland gain (user’s: 0.28, producer’s: 0.79) and grasslands loss (user’s: 0.74, producer’s: 0.92). Based on the error matrix, the mapped areas were biased low for the stable grasslands (−18.6%) and biased high for the grassland gain (+281%) and grassland loss (+25.2%) relative to the estimated areas. The accuracy assessment supported the need to adjust the area obtained from pixel counting by taking into account the information contained in the error matrix and the use of these bias-corrected estimates with 95% confidence interval to assess the progress in mountain grasslands conservation by 2020. This accuracy pattern prevailed at each extent of analysis (lowlands, mountains, mountain protected and unprotected; Table 4). This was expected, since the error matrix used as the input for bias-corrected area estimates at these extents was the same as the study area, varying only the mapped area.

Regarding the relationship between class-allocation uncertainty in the classification stage and accuracy in the validation stage, the nonparametric Wilcoxon-Mann-Whitney test with 10,000 sample Monte Carlo distribution approximations showed that median values of cumulated Shannon entropy were significantly different between corrected classified and misclassified reference areas (z = −5.6947, p = 0.0001, n = 708). Areas correctly classified have low cumulated Shannon entropy values (Figure 7).

4. Discussion

4.1. Progress in the Conservation of Mountain Grasslands in Southern Europe by 2020

Our grassland cover change analysis based on bias-corrected area estimates over the transboundary mountain region of Peneda-Gerês in the Iberian Peninsula, a representative of ancient and traditional pastoral systems of Southern Europe, revealed progress that is inconsistent with the European and international targets for mountain grassland conservation by 2020. From the point of view of being strictly based on cover, conservation has negatively progressed in mountain grasslands. Decline was the prevalent pattern between 2002 and 2019 in the study area, lowlands, mountains, mountains protected and unprotected extents. The results are congruent with the trajectory of decline documented in the European mountains since the 1960s [11,43,44], in the neighbourhoods [45], the set of studies close to the 2020 target year [14,15,16,17], and in the study area between 2000–2010 by Regos et al. [13]. Over in the area corresponding to our mountain and mountain protected extents, Regos et al. [13] found a decline of 33.8% and 43.1% in grassland cover respectively, using mapped area estimates (not corrected for classification errors). By 2020, our net loss in grasslands over these two extents (mountains, −17.4%; mountains protected extent, −14.4%) was lower than previously registered. This suggests that the rates of loss are slowing down and conservation efforts performed in the region by different policies in the last 10 years (e.g., decoupled aid of CAP since 2010) are at least smoothing negative tendencies in mountain grasslands. The decoupled aid of CAP (since 2010) increased direct payments to farmers (30%) based on compliance with “greening measures” (e.g., maintaining of grasslands and biodiversity [46]. These payments convinced farmers to maintain and convert annual crops to grasslands and, where possible, reopen grasslands. This is plausible given the strong socio-economic crisis that affected Portugal (2008–2016). People that remained in the region used the reopening of farmlands as one of the strategies to cope with the crisis. In the Carpathians, Halada et al. [15] identified a remarkable increase of managed grasslands (1999–2015; +39%) and CAP efficiency as the main booster for this trajectory.

Despite the deceleration observed, the net loss of grasslands over our study area by 2020 (−21.4%) was fourfold higher than in other southern European mountains (Apennines: −4.6%; [17]) and twice that of registered in Central European mountains (−10.4%, [16]). Large declines were only registered in boreal seminatural grasslands (1960–2015, −49.1%; [14]. While we have not specifically assessed the environmental implications of grassland decline for the Iberian mountains, the literature shows that biodiversity can decrease [47] either through abandonment or intensification of practices in the remaining grasslands, increased wildfire proneness and extreme events [48], natural encroachment of cultural heritage [49], alternations in the biomass balance affecting cattle diet [50] and reductions in national food security [51].

A smoothing that moderates the losses of grasslands seems to be exerted by the protected land governance regime. However, this requires careful interpretation. The differences found between protected and unprotected governance regimes requires further counterfactual investigation to obtain an evidence-based effect of governance regime [52]. The net grassland loss in the protected was −5.3% lower than in the unprotected regime. Ridding et al. [20] also found that sites receiving protection retained more grasslands compared to unprotected sites. Unfortunately, in Europe, grasslands are underrepresented among nationally protected areas [53]. Conservation policy has to address this issue. The lack of studies makes further discussion difficult in this regard.

4.2. Accuracy and Uncertainty

In fragmented landscapes such as the European mountains where grasslands can be smaller than 1 ha, grassland mapping with earth observations as denoted by Lopes et al. is still a challenge [27,54]. With a two-class classification scheme embedded in an MCE approach, our accuracy and uncertainty at moderate resolution support this idea. Despite the efforts to minimize sources of noise (e.g., phenological differences, chosen classifier, quality of pixel assignment), the accuracy and uncertainty registered suggest that grasslands cover in mountains should always be estimated with a confidence interval. Assuming the mapped area to represent the true grassland cover in such internationally sensitive issues as conservation and policy targets has risks. By providing unbiased area estimates with 95% confidence intervals and Shannon entropy at the pixel level, we made progress with respect to previous studies [13,55]. Readers, conservationists, policy makers are now aware about uncertainty in the area estimates and where doubts persist spatially.

Our high overall classification accuracy, and the decrease observed in change classes (grassland loss, grassland gain), followed previous grassland studies with Landsat data in mountains [35,56,57]. The producer’s and user’s accuracy of grasslands loss classes was surprising in a positive way, since they represent an improvement with respect to previous findings [35,56]. The accuracy issues with the grassland gain class were not new. Griffiths [56] also reported low producer’s (0.36) and user’s (0.61) accuracy for this type of change. Our low user’s accuracy represents a weakness of the present study. Remote sensing and conservation communities must work together to find a solution for this classification issue. Despite that, bias-corrected data showed a predominance of grassland loss with respect to gain. Even considering the lower bias-corrected estimate (ha, lowest 95% confidence interval), we still found a trajectory of decline in grasslands between 2002 and 2019.

Loosvelt et al. [30] found the majority of incorrectly classified pixels were associated with large Shannon entropy values. Our results, based on the comparison of the cumulated Shannon entropy values of the correctly and misclassified validation areas, supports this finding. This uncertainty pixel layer may thus support intermediate refinement of classification outputs and inform the spatial variation of classification accuracy [58].

5. Conclusions

This study attempts to provide a characterization of the progress made by 2020 in mountain grassland conservation using a satellite remote sensing approach that copes with the frequently unconsidered uncertainty issue. The results presented here fill a data gap around the year 2020, and a critical gap in understanding the shifting character of mountain grasslands due to two decades of European and international policy commitment over this sensitive socio-ecological domain. It illustrates also how uncertainty can be incorporated in mapping studies and improves policy assessment and reporting over mountain grasslands. The conservation progress described here is consistent with the notion that European and international policy commitments over mountains are slowing negative tendencies in mountain grasslands but not bringing these close to zero or reversing them. The slightly lower negative tendency of the mountain protected governance regime, compared to unprotected, is consistent with the alleviation effect that protected areas can exert in systems under pressure, such as mountain grasslands. Expanding protection of grasslands by supporting farmers to maintain and reopen pastoral activities, and regulating land conversions in lowlands, seems necessary to safeguard mountain grasslands. Despite this, the reduction of uncertainty in grassland change areas, namely grassland gain, remains an issue that future studies should address in view of more accurate perspectives. Class-allocation uncertainty metrics at the pixel level can be of support in the refinement of classifications, and benefit accuracy. For a more complete perspective in conservation and policy efficiency in mountains, a characterization of the ecological condition of grasslands is necessary, namely in key components such as degradation, quality and productivity, where satellite remote sensing can also contribute.