Combining historical aerial photography with machine learning to map landscape change impacts on dry grasslands in the Central Alps

Striking land-use changes after WW II characterize the past century in the European Alps with impact on ecosystems and biodiversity. Documenting land-use changes is often difficult due to limited information from the past. Mapping landscape history with aerial photography can foster the understanding of human-induced changes in vulnerable ecosystems, such as the remnants of dry grasslands in the Central Alps. We aimed to assess changes in grassland vegetation and their current extent in Val Venosta (European Alps, Italy) in relation to overall landscape settings, anthropogenic drivers of change and the effectiveness of the protected areas. We performed a land-cover classification based on a mixed machine learning approach including several auxiliary classifiers in a random forest model to characterise the extent and state of (dry) grasslands. We calculated landscape metrics between 1945 and 2015 to assess shape-related changes, especially regarding their landscape embedding and the protection status of sites. Three main processes related to a changing extent in grassland habitat prevail: (i) agricultural intensification, (ii) settlement expansion at the valley bottom and (iii) forest expansion (afforestation and encroachment due to decreasing pasture activities) on the valley slopes. The remaining grassland habitat is increasingly isolated and fragmented, leaving only few core areas of dry grassland, which tended to be better conserved within protected areas. The changes in extent of dry grasslands revealed marked changes. Transformations are assumed to be predominantly caused by human impact and successional changes. Our results confirm the importance of protected area networks. The pronounced landscape changes underline the urgent need for future research with explicit focus on the changes at community level and the underlying causes. Identifying all relevant drivers of change should be a key element in targeted conservation efforts.


Introduction
Dry grasslands are intrinsically tied to the cultural landscape in Europe and have been influenced substantially by human activity for millennia (Hejcman et al. 2013). They are known to be particularly rich in vascular plants (Wilson et al. 2012) and offer indispensable habitats for a wide range of animal taxa such as birds and insects (Nagy 2014;Hochkirch et al. 2016). However, dry grasslands are amongst the most negatively affected habitats across Europe (Diekmann et al. 2019;Dengler et al. 2020), primarily due to land-use change (Tasser and Tappeiner 2002). The on-going trend of declining dry grassland area is associated with increasing habitat fragmentation. In turn, habitat fragmentation can provoke spatial and genetic isolation (Leuschner and Ellenberg 2017). A lack of connectivity within a wider landscape can result in the extinction of habitat specialists especially of those with small population size ). Consequently, this may lead to species impoverishment and homogenization in communities (Wilhalm 2018).
Dry grasslands that resemble those of the Eastern European steppe zone, occur as extra-zonal outposts in few central-alpine valleys with steppelike climatic conditions due to orographic seclusion (Braun-Blanquet 1961;Boch et al. 2020). As these inner-alpine dry grassland communities are increasingly threatened by agricultural intensification of neighboring areas, abandonment, or afforestation (Tasser et al. 2009;Zulka et al. 2014;Wilhalm 2018), there is an urgent need for conservation efforts. Therefore, identifying and understanding the main drivers of decline and degradation in grassland area is of crucial importance (van Looy et al. 2016), especially for the development of efficient protection strategies such as the Natura 2000 network in the framework of the Fauna-Flora-Habitat Directive (92/43/EEC; European Commission 2013). Since long-term monitoring data is scarce, it remains a challenge to evaluate the ecological effectiveness of conservation measures in-and outside of protected areas (EEA 2020), especially in human shaped habitats. Hence, understanding the spatio-temporal vegetation dynamics in a broader landscape context could improve the knowledge on transformation status, conservation priorities and restoration potential, which remains a research gap for dry grassland habitats of the Alps (Lübben and Erschbamer 2021).
Past changes in landscape structure significantly influenced the current structure and functioning of ecosystems and, in the same way, the current management will influence the future structure and functioning (Forman and Godron 1986). In this context, the consideration of land-use history can be a key element in explaining current patterns and resulting change processes (Aggemyr and Cousins 2012). Consequently, the quantification and characterization of landscape patterns is fundamental for identifying ongoing ecological processes, but also areas at risk of degradation and transition (Uuemaa et al. 2013). Landscape metrics derived from land cover (LC) maps generally offer high potential for capturing changes in the spatial configuration and extent of habitats by comparing remotely sensed images from different time-steps (McGarigal et al. 2012;Uuemaa et al. 2013). In this context, historical aerial photographs, as witnesses of the past, may provide valuable information for inferring former landscape structure (Nebiker et al. 2014). Simultaneously, the variable quality and the lack of calibration details of such panchromatic images constrains accurate classification of the actual spatial landscape elements (Laliberte et al. 2004). While many studies that use historical data still rely on manual classification (Aggemyr and Cousins 2012;Garbarino et al. 2014;García et al. 2014), such approaches are known to be very laborand cost intensive (Tasser et al. 2009;Liu et al. 2018). Alternatively, the application of machine learning algorithms in the context of LC classification offers a promising approach either by pixel-based moving windows analyses (e.g., Berhane et al. 2018;Liu et al. 2018) or object-based image analysis summarizing pixels with similar topological features into spatial segments (e.g., Berhane et al. 2018;Gobbi et al. 2019).
To assess the possibility of classifying the historical landscape by using historical aerial photographs dating back more than 70 years, we (i) developed a machine learning approach at two different time steps and compared differences of the historical and recent land cover, (ii) evaluated the spatio-temporal grassland characteristics regarding their shape and landscape composition and (iii) analyzed differences in transformation trends of dry grasslands in protected versus non-protected areas. With this approach, we evaluated changes between 1945 and 2015 regarding habitat properties and landscape composition with focus on remnants of dry steppe grassland in the inner-alpine dry valley of Val Venosta (European Alps, Italy).
We applied a novel mixed machine learning classification procedure (pixel-and object-based) including several auxiliary classifiers. By integrating spatial, textural and structural classifiers, we were able to take advantage of object-based image analysis translating the individual pixels into meaningful spatial units. With our combined approach, the classification was neither limited to previously generated segments nor was it limited to the information of individual pixels.
To capture landscape changes as well as changes in our focus LC class grassland at both time steps, we selected a small study-specific set of uncorrelated landscape metrics that can reveal potential habitat fragmentation and -isolation tendencies (Schindler et al. 2008;Walz 2011;McGarigal et al. 2012). Finally, we estimated the magnitude of transformation in dry grassland areas in relation to their environmental matrix to detect areas of change within-and outside protected areas.

Materials and methods
Biogeographical and historical outline of the study area Our study area, the Val Venosta/Vinschgau, is situated in the Central Alps within the Autonomous Province of Bolzano-South Tyrol, Italy (Fig. 1). The west-east oriented valley is enclosed by high alpine mountain ranges (Ortles-, Ötztal-and Sesvenna Alps reaching up to 3905 m a.s.l.), which cause exceptional biogeographic conditions in terms of the climatic and biological setting (Rampold 1997). These high mountains isolate the valley against maritime air masses, which comes along with low precipitation partly below 500 mm a −1 at the valley bottom (Wilhalm 2018). High temperatures due to often cloudless conditions combined with strong foehn winds lead to dry Fig. 1 Map of the study area and the position of dry grassland biotope areas. The upper boundary (delineated in black) for our study follows the maximum altitude records for dry grassland areas (at 1500 m a.s.l.). GIS-data was retrieved from the catalogue of the Autonomous Province of Bozen/Bolzano (https:// geoka talog. buerg ernetz. bz. it/ geoka talog) conditions especially on steep south-facing slopes. Based on this climatic setting, the valley offers a habitat for highly specialized grassland species adapted to a dry and hot climate up to 1500 m a.s.l. (Braun-Blanquet 1961;Schwabe and Kratochwil 2004;Lübben and Erschbamer 2021). The limited water availability on the rocky, and steep south-facing slopes determines the type of vegetation that can colonize these slopes (Leuschner and Ellenberg 2017). Further, very early human intervention through grazing, fire and cutting favored the persistence and spread of open grasslands since neolithic times and throughout the Holocene (e.g. Dietre et al. 2017). Since the early twentieth century, increased afforestation efforts promoted the migration of Pinus nigra but also Robinia pseudoacacia to stabilize the steep south-facing slopes (Braun-Blanquet 1961). These afforestation efforts together with decreasing pasture activities on the valley slopes and agricultural intensification at the valley bottom have led to major LC change in the Val Venosta (Wilhalm 2018). Furthermore, the focus of agriculture in the Val Venosta has changed tremendously since the 1940s with a heavy expansion of more profitable apple cultivation (ASTAT 2018). Consequently, the valley bottom is characterized by a vegetation mosaic where homogenization processes prevail, characterized by an impoverished and trivial flora (Hilpold et al. 2018;Wilhalm 2018). Contrastingly, shrub encroachment and forest is intermingled with the original pasture-shaped continental dry grassland communities of the order Festucetalia vallesiacae on the slopes (Lübben and Erschbamer 2021).
The conservation efforts regarding the remaining dry grassland areas in the valley included the designation of five major dry grassland protected habitat complexes (Fig. 1). The protected sites are actively managed by maintaining grazing regimes or the active removal of shrubs and invasive alien species (Table 1). Hence, we considered the Val Venosta to be ideal for studying LC change related impacts on (dry) grassland habitats.

Historical aerial photographs
The historical LC was derived from a set of 18 historic aerial photographs in 1945 (Istituto Geografico Militare, Firenze IT). They belong to the oldest aerial documents for the landscape of Val Venosta. The panchromatic images were obtained on August 27th with a Fairchild camera operating at an altitude ranging from 7900 to 8200 m a.s.l. and a focal length of 137 mm. The aerial photographs were digitized with a 2500 dpi scanner and subsequently georeferenced, orthorectified, and projected (i.e., ETRS89-UTM-32N). Data were provided by the Amt für Landesplanung und Kartografie of the autonomous province of Bolzano (2021).

Recent orthophoto
For identifying the current LC, we used a recent orthophoto in the geodetic system WGS84-UTM-32N recorded from July to September in the years 2014 and 2015 (Amt für Landesplanung und Kartografie 2021). The product is a 3-band (RGB) true color orthophoto with a resolution of 0.2 × 0.2 m. The high resolution outperforms open-source satellite products (e.g. Landsat, Sentinel-2) and presents a valuable benefit for identifying changes of inner-alpine dry grassland areas, given their small remaining distribution.

Class definition and pre-processing
We selected our classes (Table 2) purposefully for the most dominant classes upon visual inspection (Chen and Stow 2002), screening for the universally defined categories according to Anderson et al. (1976). The classification scheme was developed with the aim of describing the variety of LC features in the study area. Since we focused on the change detection of grasslands, we identified main vegetation classes with clear and unambiguous class boundaries in reflectance values aiming at the best distinction between the individual classes. Generally, classes with strongly overlapping spectral signatures impede clear class distinction (Pringle et al. 2009). Hence, we only included dominant land cover classes that were present in the valley at both time steps (upon visual inspection and delineation; Table 2). Each of the chosen classes covers major parts in the valley and is considered a representative class of the valley's landscape characteristics (e.g., the dominant apple cultivation in the valley bottom). Subordinate classes (such as other farming practices) were summed up to larger categories. The class boundaries were tested statistically for significant differences in spectral reflectance values based on a randomly selected subset of each class (n = 20,000 pixel). We used pairwise Wilcoxon rank sum test (Bonferroni corrected) and additionally confirmed distinctiveness of classes visually via boxplots for historical and recent data sets (see Suppl. Mat. 1- Fig. S1). Since the historical reflectance values partly lack distinctive characteristics between the individual classes, we corrected the reflectance values, accounting for high within-class variability due to topographic effects like differing solar illumination resulting from rugged terrain characteristics (Teillet et al. 1982;Riano et al. 2003). We applied a c-correction with a terrain smoothing factor of 5 (out of 5 tested methods) as the best compromise between reduction of variability while best maintaining the class mean reflectance. We validated the performance by inspection of the effective interquartile range (IQR) reduction (Ma et al. 2021) and the coefficient of variation (COV). The COV represents the reduction of the standard deviation from its mean (Wu et al. 2016). This resulted in an IQR reduction of 12.9% and a COV of 2.3% on average (details on the illumination modelling and correction methods are provided in Suppl. Mat. 2).
To account for high inter-class similarities (e.g., similar reflectance values for objects in different LC classes like apple trees versus trees in a forest), which potentially increases misclassification rate in LC classification (Huth et al. 2012), we created a mask dividing our data set into valley bottom and valley slopes. This allowed differentiating spatially segregated LC class distributions where, e.g., the apple orchards are present on the valley bottom and the forest area is mostly restricted to the valley slopes. We subdivided the study area based on the pronounced cut in the valley's topography. The valley floor is predominantly flat, while the north-and south-facing valley slopes distinctly emerge from the valley bottom with a steep incline, frequently exceeding 30°). We manually determined the division line according to this prominent change in topography using a DEM at 1 m resolution. The approach to mask areas according to the research objective is consistent with other studies (e.g. Pringle et al. 2009;Eitzel et al. 2016). We selected training areas for each data set that best represented the average and variation in spectral reflectance of each class (see Table 2). We identified the training areas by manual delineation of polygons spread evenly across the valley to avoid missing gradual changes of class reflectance values. We continuously selected training areas until no further variation in the spectral profile was added (by monitoring the spectral class profiles in ArcGIS, version 10.7, ESRI, Redlands, CA, USA). Finally, each class consisted of at least 30 training areas, but more were added if needed. For the actual LC classification, we resampled both datasets to 1 × 1 m resolution using bilinear interpolation, which represented a compromise between computing power and best possible spatial resolution. We extracted the pixels within the designated training areas subsequently previously to the supervised classification.

Land cover classification
The general classification framework is provided in Fig. 2. It displays a schematic overview of the processing steps as explained in detail in the following sections.

Auxiliary model classifiers
LC classification of panchromatic aerial photography is acknowledged as a difficult objective because the spectral information is limited to a single spectral band, which inevitably comes with rather high intraclass but low inter-class variability (Ratajczak et al.

2019). Comparable studies with similarly limited
Classes were carefully selected visually according to the most dominant land cover classes present at both points in time (i.e., historic and recent) *All parts indicated with asterisk are incorporated during the post processing as an overlay obtained as vector layer from the catalogue of the Autonomous Province of Bolzano **Even though dry grasslands present the natural occurring grassland type in the valley, we refer to the class as "grassland" where not explicitly indicated differently (based on additional in-situ inspection, e.g. within protected areas) input data (e.g. Liu et al. 2018;Gobbi et al. 2019) identified object-based analysis as best classification method for panchromatic data. However, a disadvantage of object-based approaches especially in the context of vegetation analyses regards the omission of small landscape elements such as single trees/shrubs (De Giglio et al. 2019). In our study site, some shrubs barely cover 1 m 2 . Clustering areas of similar spectral characteristics into larger image segments could dissolve such elements. However, they are important to quantify the encroachment by woody LC classes as prevalent and recent development in our study area (Lübben and Erschbamer 2021). This highlights the advantage of using pixel-based approaches. Since both approaches generally achieved a high classification quality in comparative studies (Suppl. Mat. 3), we decided to draw from the advantages of both methods with a mixed approach, which is seen as promising approach that could enhance LC classification (Malinverni et al. 2011;Berhane et al. 2018). Therefore, we derived five additional spatial classifiers with different information content on mean object-and texture characteristics while also implementing the individual pixel reflectance values of the aerial photographs in the classification model (Table 3, extracts of the layers in Suppl. Mat. 4). Object-based classifiers are derived by a mean shift segmentation (spectral detail: 17, spatial detail: 2, minimum segment size: 10; following Lourenço et al. 2021) in ArcGIS and textural classifiers by a 7 × 7 moving window based on a grey-level co-occurrence matrix with "glcm" package (Zvoleff 2020) of the statistical programming language R (version 4.2.1, R Core Team 2021). All classifiers were calculated for each available band (red, green, blue), which resulted in a total set of 18 classifiers for the recent dataset and 6 classifiers for the historical dataset.

Model setup and goodness-of-fit-assessment
In context of landscape ecological analyses, the most frequently applied machine learning approach are random forest models (Stupariu et al. 2022). Random forest is a non-parametric machine learning algorithm that determines the local LC class (response vector) based on the original reflectance values and the auxiliary classifiers as predictor variables (Breimann 2001). Prior to all classification steps, we divided the pixels within the previously delineated training areas into 70% training data and 30% validation data. The pixels of these manually designated training areas were extracted randomly (compare section "Class definition and pre-processing"). We implemented a random forest model, exploratively determining the optimal number of trees (n tree = 300) and classifier used per decision (m try = 3) by aiming for the best classification output at the lowest computational power. To reduce collinearity among the classifiers and, further, dimensionality we applied backwards variable elimination based on their respective importance for the LC class attribution (Diaz-Uriarte 2007). All models were set up in a supervised classification manner based on the R package "caret" (Kuhn 2022).
After selection we trained our model on the remaining classifiers based on the training pixels (Hufkens et al. 2020). We randomly permuted the variables at each node decision and implemented a tenfold leave-one-out cross validation with a final majority voting over all built trees (Kuhn and Johnson 2013). We validated the model performance based on the 30% of the test pixels with classical error estimates: the out-of-bag proportion (OOB), overall accuracy (OA), Kappa statistics and LC class accuracies by confusion matrix measures (Keshtkar et al. 2017;Berhane et al. 2018;Hufkens et al. 2020). Finally, we classified the full study area applying the best performing model.
After classification, we merged the individual datasets (valley bottom, valley slopes) and added the missing classes manually (i.e. the few remaining patches of forest on valley bottom) or with the auxiliary layers, i.e. rivers and roads were added based on publicly available vector data provided by the autonomous province of Bolzano for the recent data set (compare Table 2). Due to the heterogeneous reflectance values of the LC classes agriculture (other), builtup and apple orchard historically, we limited their presence to larger manually categorized areas. We verified the correction by repeated inspection of the class accuracy measures.

Landscape change assessment
To assess the landscape changes, we calculated spatial landscape metrics specifying the landscape composition and spatial configuration at both time steps (Haines-Young and Chopping 1996). We selected a set of metrics that capture changes in grasslands in Table 3 Auxiliary classifiers, their calculation and purpose for enhancing potential classification accuracy *The authors applied the standard deviation, which we specified to the pixel deviation. The standard deviation would equal one value per segment, we were interested in the homogeneity/heterogeneity in the segments Improving the identification of spatially scattered noisy LC classes (e.g. differing illumination due to shadowing as present in forested areas) Saatchi et al. (1999), Caridade et al. (2008) Contrast (7 × 7) Difference between the highest and lowest reflectance value in the pixel neighborhood Measure of the local variation within the pixel neighborhood to distinguish systematically rough LC classes (e.g. apple orchard rows) and class borders Saatchi et al. (1999), Caridade et al. (2008) terms of its patch characteristics related to (core) area, edge, shape, and aggregation of the class's patches in the landscape. By examining the Pearson correlation coefficient, we filtered highly correlated variables that were significant for both time steps (r > 0.8).
We selected significant metrics (Table 4) based on the Variance Inflation Factor (VIF < 10) to reduce multicollinearity (Plexida et al. 2014). Additionally, we extracted the mean and standard deviation of the five selected metrics for each LC class to highlight changes on class and landscape level between the time steps.
Furthermore, we compared the tendency of dominance of an individual class expressed as the Shannon Evenness Index ( SHEI ) for both time steps distinguished in valley bottom and valley slopes (McGarigal et al. 2012). To examine significant temporal pattern variation in grassland characteristics represented by the landscape metrics, we applied the nonparametric Wilcoxon rank sum test (α = 0.05) due to our non-normally distributed data with heterogeneous variances. We tested for true effect size on the significant differences by applying Cohen's d (package "effsize"; Torchiano 2020). Finally, we assessed the changes of the dominant LC classes between the time steps as well as the overall change magnitude for grassland areas. We calculated the change magnitude based on proportional transformations of grassland considering the surroundings of each grassland pixel of 0.5 ha to capture the average changes (compare Schwabe and Kratochwil 2004). Therefore, we implemented a moving window analysis with a circular buffer around each grassland pixel to assess the relative change of grassland compared to the historical grassland distribution. This revealed grassland areas that may be in transition to another LC class or have changed completely. In this context, we particularly consider potentially different developments on protected dry grassland areas in the five major dry grassland biotopes compared to non-protected areas. Therefore, we analyzed the transformation extent (as expressed by the magnitude) separately in-and outside of protected areas. Differences could potentially result from varying management practices. By considering the transformation in a broader neighborhood, we view the grassland changes in its broader landscape context, i.e., an area will only show a high magnitude in change if the whole grassland area in its surroundings has transformed accordingly. All analyses have been conducted in R by applying the package "landscapemetrics" (Hesselbarth et al. 2019).

Model performance
In general, the four final models (i.e. two models per time step representing valley bottom and valley slope, respectively) performed well (Table 5), but there were differences in accuracy, especially for the two different points in time. The classification of the recent datasets achieved excellent overall model accuracies with 99.2% (valley slopes) and 99.4% (valley bottom) and high model reliability (κ = 0.99), whereas the In brackets the class accuracy after post-processing is indicated. True positive = correctly identified, false negative = missed to identify as respective LC class. The final choice of classifiers is ordered according to their importance for the model historical dataset achieved good accuracy with 76.8% (valley slopes) and 71.2% (valley bottom) also with substantial model reliability (κ = 0.69, and κ = 0.62, respectively). The most important model classifiers (Table 5) were the segmentation-as well as textural layers.
For the historical dataset, particularly the classes shrub and agriculture on the valley slopes showed significant error rates, which were not improved with the post-classification processing. The greatest mixup of these classes occurred with the class grassland (Suppl. Mat. 5). Classification of the class agriculture was similarly difficult on the valley bottom. Nevertheless, the post classification procedure improved the overall accuracies by 13.2% (valley slopes) and 14.2% (valley bottom). We inspected the effectiveness of our post-classification for all classes by manual comparison of the final map to the aerial photographs and also by anew extraction of the class accuracy metrics (Table 5; given in brackets). Overall, the class accuracies were good, whereby especially well distinctive classes due to clear spectral signatures such as forest profited from the post-classification.
In contrast, the recent datasets performed well for all classes. Only the class shrubs was slightly underclassified (false negative: 4%) on the valley slopes. Further errors were considered negligible (compare Suppl. Mat. 5), which was also confirmed by in-situ inspection.

Proportional landscape changes
Fundamental change in land cover occurred in the study area since 1945 that can be summarized in two main trends differentiating the developments on valley bottom from the ones on valley slopes: (i) On the valley bottom, a massive expansion of apple orchards and built-up area at the cost of other cropland and grassland habitat occurred. (ii) On valley slopes formerly mixed (semi-)natural vegetation patches mostly consisting of shrubs, grassland and forest developed into homogenous forest areas.
The most prominent LC shifts over the whole area were settlement expansion and agricultural intensification (Fig. 3). At the valley slopes built-up area also increased slightly. Further changes were observed in the LC class forest, which increased on the valley slopes while shrubs markedly decreased in area (compare details (a) and (b) in Fig. 3). Shrubs and grassland transformed mostly into dense forest. This indicates that relevant ecological processes in the study area such as the transformation of former grassland area to successional forest, were identified with the classification. The high transformation of shrubs, however, points towards a potential misclassification historically. While a significant shift in the grassland pattern was prominent (Fig. 3), an overall decline in grassland may have been obscured by historical misclassification as shrubs (compare classification performance Table 5).
Investigating the historical grassland that has changed to other LC classes today, revealed a divergent trend on valley bottom and valley slopes (Fig. 4). Only 30% (valley slopes) and 19% (valley bottom) remained unchanged grassland pointing towards substantial successional processes on the valley slopes and clear anthropogenic influences on the valley bottom. 54% of all transformed grassland areas on the valley slopes are woody classes today (forest 50%, shrubs 4%) followed by agricultural expansion (mostly to fertilized grassland areas for forage production). On the valley bottom, the majority of all transformations (61%) consisted of conversion to apple orchards, followed by settlement expansion (9%). Substantial shifts in grassland extent are apparent and confirm high landscape dynamics in the valley.

Spatial patterns of grassland changes
The spatial distribution of the changes in grassland area exhibited clusters that featured substantial differences in change magnitude (Fig. 5). In the central-and eastern parts of the study area grasslands were almost fully transformed at the valley bottom, whereas in the western study area they partly remained unchanged. Changes on the valley slopes were less equivocal and rather heterogeneous with the eastern part of the study area (i.e. Laces and Naturno) featuring a minor magnitude of change. Furthermore, the magnitude of change showed that grassland extent was tendentially better preserved under protection. Especially protected sites with grazing management (compare Table 1) appeared better preserved. These prominent differences are displayed in Fig. 5.3a (grazed areas) compared to Fig. 5.4a (non-grazed areas), which is confirmed in the average magnitude of change of 31.5% versus 72%, respectively. Overall, the transformation status of grassland effectively revealed the patterns of areas in transition to other vegetation types (yellow zones in Fig. 5). The map highlights that only a few grassland areas that were not affected or intruded by neighboring habitat types remained stable and still featured a substantial core area. The Wilcoxon rank sum test indicated significant differences in all landscape metrics between the two points in time for grassland. The landscape metrics CAI, SHAPE, FRAC and CONTIG differed significantly on valley slopes with low to medium effect sizes (Cohen's d) between the two points in time. Similar effect sizes were found at the valley bottom for SHAPE and CON-TIG (Table 6). Today's grassland patches were characterized by a mean decrease in core area (− 59.9%) and a less complex shape (− 27.1%) but also less connected patches (− 20.8%). These results point towards more isolated and more compact grassland patches with a reduced core area (i.e., patch separation). The differences between the two points in time were less pronounced on the valley bottom. Nevertheless, a decrease in shape complexity (− 26.7%) and greater spatial separation between the individual grassland patches (− 15.0%) was evident. Next to these trends on valley slopes and valley bottom we found a homogenization tendency. This was indicated by the diversity metric SHEI that decreased by 29.69% for valley slopes (historic: 0.64, recent: 0.45) and 9% for valley bottom (historic: 0.70, recent: 0.64). The dominant classes today were forest covering 73% of the valley slopes and apple orchards covering 60% of the valley bottom. Both LC classes showed a visible trend of spatial clustering (Fig. 3).

Grassland characteristics in relation to a changing cultural landscape setting
In line with other studies (Tasser et al. 2009;Wilhalm 2018;Dengler et al. 2020) our LC change analysis revealed a significant shift in the pattern of grasslands in many areas. The shift in pattern was mostly due to the spread of neighboring habitats which also led to a decline in grassland extent. Especially at the valley bottom the decline was visible, with agricultural intensification and settlement expansion as major causes for habitat loss. Land-use change on the valley bottom is mostly driven by the heavy expansion of apple cultivation areas. Nutrient input through direct fertilization and heavy pesticide application is associated with biodiversity decline in close vicinity (Hilpold et al. 2018;Diekmann et al. 2019) but also on adjacent slopes as shown by the degradation of butterfly communities in the valley (Tarmann 2009). Mountain slopes experienced a change in management regime towards abandonment of grazing and afforestation that is reflected by increased forest cover (Wilhalm 2018). These results are in line with trends in grassland habitats across Europe (Zulka et al. 2014;Deák et al. 2016;Gossner et al. 2016;Magnes et al. 2020Magnes et al. , 2021Choler et al. 2021). Historical grassland areas identified as such have changed markedly, indicating a considerable change in grassland extent. However, the classification of the historical dataset was difficult since historical data was only available as panchromatic photographs with a single spectral reflectance value. Potentially, terrain features (e.g., ridges or furrows) obscured correct classification of the LC classes and further impeded to capture overall changes in grassland extent especially in larger shadowed areas (compare Suppl. Mat. 7). Partially similar spectral characteristics (Table 2) on areas with low growing crops or fertilized pastures as well as in shaded areas caused by the terrain prevented a flawless classification of the historical data. The terrain characteristics caused shadowing effects and attenuated the reflectance values (Vanonckelen et al. 2014), especially on open LC classes such as grassland. This limited sound inferences on the overall changes in the entire valley and required to focus on areas that were visually confirmed as grassland. On the valley bottom, the high reflectance values of some objects (e.g. buildings made of concrete) in the class built-up led to confusion with the similar spectral characteristics of grassland. Goodness of fit measures (Table 4) and visual inspection of both LC maps and aerial photographs confirmed the partial misclassification of grassland on the historical LC map. Another potential cause is the spatial resolution of the historical photographs, which can cause small and frequently changing landscape features to blur into larger coherent areas (Laliberte et al. 2004). Although we resampled both datasets to equivalent spatial resolution, it seems confirmed that the differences in data quality affect the landscape metrics as a result of potential classification error and grain size (Uuemaa et al. 2005;Frazier and Kedron 2017). In contrast, the clearly superior data quality and the additional spectral bands of the recent orthophoto allowed the detection of grasslands also behind ridges or furrows. Thus, the overall results must be viewed critically due to potential errors resulting from classification inaccuracies and data quality of the historical aerial photography. Still, such historical data can well serve as base for showing trends and shifts in patterns (Eitzel et al. 2016). The trend of encroachment already existed in the past, which is reflected by the structural patch characteristics. There is a tendency that grassland patches, which still exist today, are characterized by less core area (CAI) and higher isolation and spatial discontinuity (CONTIG). This trend is contrasted by a reduction in shape complexity (SHAPE, FRAC ). Usually, spatially isolated habitat structure is associated with higher shape complexity because homogeneous patches are intruded by species of surrounding habitats close to the patch boundaries (Labadessa et al. 2017). Thus, we assume our results to be influenced by two potential factors: (i) Intensive human influence coming along with geometrization of landscape elements (Forman and Godron 1986), and (ii) the loss of ecological corridors (i.e. connection between the patches of equal LC class; McGarigal et al. 2012). Formerly heterogeneous patch structures that were still connected and once built a bigger patch are now cut by other habitats and isolated. Only in the western study area spatially clustered grassland patches remained unchanged. This could explain the less pronounced habitat discontinuity on the valley bottom. Translating these developments back to past landscape structure our results suggest: The historic grassland area was rather continuous with a larger core area and better connectivity. Still, surrounding shrubs already intruded the uniform grassland patches, which resulted in a more complex shape of the grassland areas. Transferring patch-related findings to landscape level, today a picture of a rather homogeneous landscape emerges, whereby few LC classes (valley bottom: apple orchard, valley slopes: forest) are dominant (as confirmed by SHEI). Diacon-Bolli et al. (2012) emphasized that the conservation of grasslands in Europe, as historical co-production of the cultural landscape, depend on a heterogeneous landscape with high structural diversity. In this context, our findings underline the need for counteracting further declines in the grassland's extent, which are embedded in a highly transformative landscape.

Protection needs of dry grassland fragments
The increasing isolation of dry grassland patches on the valley slopes goes hand in hand with habitat fragmentation. Only a fraction of the historical extent in dry grassland areas persisted, mostly restricted to valley slopes. For the preservation and restoration of a heterogeneous landscape mosaic (Walz 2011) they are crucial elements of the regional pastoral system that provide valuable connectivity networks for species exchange and seed dispersal (Aggemyr and Cousins 2012). The conservation of dry grasslands becomes even more important as recent research on their evolutionary history revealed that as western outposts of the continental dry grasslands they were cut from their main distribution in the Eurasian steppe zone long ago. The inner-alpine valleys, such as the Val Venosta, present a unique and genetically-independent refuge for dry grassland species 2022). Even though they are incorporated in protection frameworks, the effectiveness of current management strategies is doubtful and might not be sufficient to mitigate species decline, especially for highly specialized taxa (Diekmann et al. 2019). So far, there is no single or best solution for preserving as many vulnerable species as possible . In this context, our results point towards better preservation of the dry grassland extent within protected areas, which could also be confirmed by in-situ inspection of dry grasslands at non-protected and protected sites. The smaller magnitude of transformed dry grassland (Fig. 5) within the protection boundaries confirms the use of proper management strategies within these sites. This has also been stressed by Magnes et al. (2020Magnes et al. ( , 2021 for Austrian dry grasslands. However, we need to consider that potentially also other factors could have evoked LC changes on the non-protected areas (such as a change in climatic factors, altered water availability or other geohazards on the steep mountain slopes). We cannot fully ascertain that protected sites are better preserved solely due to their protection status. Thus, we suggest that further research is needed that also considers the changes at community and species level.
To maintain dry grassland remnants, various studies highlighted the importance of grazing to curtail encroachment processes (e.g. Zulka et al. 2014;Boch et al. 2019). Furthermore, a reduced pasture activity was addressed as potential cause for changes in floristic composition of these grasslands in Val Venosta (Lübben and Erschbamber 2021). Accordingly, our findings confirm a possible relationship between the extent of dry grasslands and the presence or absence of grazing. This is indicated by the lower magnitude of transformed grassland in protected sites with a prevailing grazing regime (Fig. 5). Our results underpin the need for preserving or reintroducing extensive grazing regimes to curtail further encroachment within the valley. Where grazing is re-introduced, the supression of woody species and litter accumulation is usually beneficial for the occurrence of specialized flora (Elias et al. 2018) and fauna (Fonderflick et al. 2014).

Benefits and limitations of the classification
Generally, the spectral information of the historical panchromatic data was limited to a single spectral band and for the recent orthophoto to three spectral bands in the visible range. Compared to more recent and frequently utilized multi-or hyperspectral Fig. 5 Land cover change magnitude map for grassland of the whole study area (upper three displays) and detail views for two dry grassland biotope areas: (3) a grazed protected site and (4) a non-grazed protected site, each with the recent aerial photograph in the background. The protection boundaries are encircled in white for all five biotope areas: (1) Tartscher Bühel, (2) Obere Leiten, (3) Kortscher Leiten, (4) Schlanderser Leiten, (5) Sonnenberg. The change magnitude is the relative share in changed grassland to the total grassland area in a 0.5 ha neighborhood. Details for 1, 2 and 5 are provided in Suppl. Mat. 6- Fig. S6 ◂ satellite imagery (e.g. Keshtkar et al. 2017;Berhane et al. 2018;Zhang et al. 2018), no secondary spatial information based on spectral reflectance in other bandwidths could be deduced for the LC classification. The lack of multiple spectral bands especially in the historical data set impeded clear LC class distinction, resulting in overlapping spectral signatures (Pringle et al. 2009). Yet, we considered a labor-intensive manual classification as impracticable because our study area spanned across more than 50 km from east to west. By applying a mixed model approach of object-and pixel-based methods we expected to bypass known limitations of the individual approaches (Gao and Mas 2008). This was based on previous evidence indicating a mixed approach could potentially enhance classification (Malinverni et al. 2011;Berhane et al. 2018). For the recent dataset our results confirmed the benefit of a mixed approach. The overall accuracy was at least similarly good to superior compared to comparative studies implementing both machine learning algorithms (see Suppl. Mat. 3). However, the classification of the historical dataset performed comparable only after post-classification, which presents a weakness of our mixed approach. The still high amount of falsely classified pixels confirms general difficulties in classifying historical aerial photographs dating far back (Lydersen and Collins 2018). The basic obstacle of similar spectral signatures between some LC classes can only be solved by complementary classifiers (Lillesand et al. 2015), which was confirmed by the results of our mixed approach. Overall, we found that the spatial-and texture-based classifiers were most important for the final model which highlights the benefit of such complementary classifiers in pixel-based approaches.
A clear limitation of using such historical data presents the time-consuming and computationally intense preparatory work: (i) The correction of the spectral signatures to enhance highly shaded-and reduce overexposed areas to minimize inter-class variability, (ii) the calculation of additional classifiers (also for the recent orthophoto), and finally (iii) the less intuitive choice of high-quality training areas due to the black and white appearance. General inaccuracies of the input data and errors made during the training process propagate and reduce the overall accuracy (Lillesand et al. 2015). In terms of change detection, geometric errors and relief displacement inhibit a one-to-one-pixel comparison, which would require identical data sets (Hussain et al. 2013). Such errors are difficult to eliminate, because historically no further ground-truthing is possible (compare e.g. Liu et al. 2018). A specific drawback in our study is also the presence of linear gaps and shifts in the historic data set as a cause of mosaicking multiple photographs without geo-located reference points historically (compare Suppl. Mat. 7). Still, the results Table 6 Results of landscape metrics for grassland represented in mean-and standard deviation (SD) values per time step and area (valley slope and valley bottom) and their change in percent compared to the past Differences in all landscape metrics between the two points in time were highly significant with p < 0.05. Cohen's d maximally reached medium effect size, for which effect size range is coded as: < 0.2 (negligible, +), ≥ 0.2 to < 0.5 (small, *), ≥ 0.5 to < 0.8 (medium, **), ≥ 0.8 (high, ***) on mean changes and the structure related metrics rendered highly useful for LC change interpretation. The lack of validation possibilities next to in-situ confirmation also applies to the recent LC map. Even though there exists recent LC information for the chosen study site (such as CORINE land cover), the coarse resolution inhibits quality assessment and verification of the landscape changes of our detailed LC map at 1 m resolution.

Conclusions
The historical cultural landscape in Europe is subject to constant changes, with anthropogenic influences often being the main cause of habitat degradation and -decline, in particular for grassland areas. By looking back into the past, we could detect changes in the landscape matrix of an inner-alpine dry valley (Val Venosta, Italy), ranging from overall LC change to detailed analysis of individual dry grassland complexes. In this context, we present a valid option for extracting LC change based on historical and recent aerial photography by integrating textural and structure-related information into a pixel-based LC classification approach. Our study underlines the great value of such historical memories, even if the quality of the historical data generally prevents pixel-wise comparison and there are no possibilities for validation going beyond visual inspection of the historical photography. Strong agricultural expansion and intensification prevail on the valley bottom, whereas potentially the cessation of grazing provoked successional changes of the former grassland area towards woody LC classes. Today's grassland patches are substantially fragmented, offering only few core refugia for dryadapted species within an increasingly homogeneous landscape setting dominated by a few LC classes. Moreover, our analysis revealed heterogeneous patterns of change between protected and non-protected areas, whereby dry grasslands appear best conserved within protected areas that still have a grazing regime. In this regard, the limited extent of the protected sites compared to the non-protected sites inhibited definite conclusions on the conservation status of these grasslands, which highlights future research needs. While our results underline the importance of protected area networks and the eligibility of grazing as conservation measure of dry grassland areas in general, it further remains unclear to what extent these areas support the persistence of individual dry grassland specialist species. The alarmingly high transformation rate of the former grassland area towards shrubs and forest in the valley threatens future persistence of these fragmented habitats. Therefore, we emphasize that further research is needed to elucidate the cause-effect relationship between the loss of dry grassland species and anthropogenic pressures, but also to consider that other drivers such as climatic changes may play a important role. It will also be crucial to understand the effect of conservation efforts in relation to a multitude of drivers of change.
Author contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by EK. The first draft of the manuscript was written by EK and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Libera Università di Bolzano within the CRUI-CARE Agreement. This research was carried out under the PhD scholarship (Mountain Environment and Agriculture, XXXV° cycle) funded by the Free University of Bozen-Bolzano.
Data availability All data that was used for this study are openly available or can be requested for scientific purposes at the Autonomous Province of Bolzano: https:// geoka talog. buerg ernetz. bz. it/ geoka talog.

Declarations
Competing interests The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.