**1. Introduction**

Landscape can be defined as a spatial configuration of patches of dimensions relevant to the phenomenon under consideration or to the selected organism, which exists only at the moment in which it is perceived by the senses. The landscape can also be considered as a portion of the real world within which we are interested in describing and interpreting processes and patterns; this context can lead to different conclusions, depending on whether we use abiotic and/or biotic factors [1]. These factors make it intrinsically dynamic, both at a temporal [2] and spatial scale [3], since it is conditioned by environmental conditions [4], the ecological processes that take place in it, changes in land use, and anthropic disturbances [5].

Landscapes with marked heterogeneity have a complex structure of habitats, which translates into a high index of diversity [6]. This heterogeneity, as a factor of organization of ecological systems, presents a permanent character of landscapes and determines the generation of differentiated environmental mosaics [7]; which makes their study at a large scale extremely difficult. To the natural dynamics perpetuated by geomorphological [8–10], hydrological [11], and biotic [12] processes, to cite some significant examples, we must add human activities, which are currently the main factor in landscape evolution worldwide [13–16].

Remote sensing imagery is widely used for land cover classification, target identification, and thematic mapping from local to global scales due to its technical advantages such

**Citation:** Sánchez Sánchez, Y.; Martínez Graña, A.; Santos-Francés, F.; Reyes Ramos, J.L.; Criado, M. Multitemporal Analysis of Land Use Changes and Their Effect on the Landscape of the Jerte Valley (Spain) by Remote Sensing. *Agronomy* **2021**, *11*, 1470. https://doi.org/10.3390/ agronomy11081470

Academic Editor: Francisco Manzano Agugliaro

Received: 9 June 2021 Accepted: 22 July 2021 Published: 24 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

as multiresolution, wide coverage, repeatable observation, and multispectral/hyperspectral– spectral records [17,18]. Several classification methods based on satellite images can be classified as supervised or unsupervised classification methods giving higher priority to supervised classification [19] because of its learning method [20,21].

Each land use class is to an image what a patch is to a landscape mosaic, so with this premise in mind, landscape metrics can be applied to measure the effect of land use changes on landscape shape [22]. Different studies calculated multiple levels of landscape metrics to measure landscape patterns in order to analyze land use, mainly focused on urban land use [23–25], and on the classification of satellite images rather than on the analysis of spatial patterns [26]. The research studies where a landscape analysis was performed used metrics based on average patch characteristics [27], ignoring the distribution, size, and changes occurring in the last decades, therefore, this study focuses on diversity indices and indices of the entire landscape.

The implementation of geographic information systems and remote sensing [28] have helped to carry out large-scale landscape studies, and with the support of different spatial pattern analysis programs [29], they have allowed a more precise study of the temporal dynamics of the landscape, which has provided a significant leap in the quality of the studies. For the analysis of spatial patterns, multiple indices have been developed that allow the study of landscape configurations at different temporal moments, evaluating their composition and configuration, the proportion of each class, or the shape of the elements [30].

The objective of this article was to know the effect of changes in land use on the landscape configured as a fluvial-structural valley during the last three decades, comparing the potential landscape with the landscape in 1994, and analyzing the spatial patterns in the subsequent years 2000, 2010, and 2019. Using Landsat 5 and Landsat 8 satellite images, and by analyzing spatial patterns and indices of diversity, dominance, shape, and fragmentation with Fragstat software, four scenarios were characterized and studying their tendency in the configuration of the environmental mosaic allowed us to understand the natural dynamics and the influences of human activities on the landscape.

Landscape metrics and the indices obtained have shown an important role in the analysis of changes in land use, so this study will enable land managers to implement appropriate measures for the maintenance of physical and functional connectivity in an anthropized environment, in order to achieve the objectives, set by the strategies of conservation and improvement of the landscape, for the sake of future sustainability.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Jerte valley is located in the north of the province of Cáceres (Spain) (Figure 1). It presents a geomorphological structure of a valley sandwiched between two mountain complexes (Sierra de Béjar, Sierra de Tormantos).

The geology of the study area is mainly formed by granitoids with some quartz outcrops, in the dykes of the Alentejo-Plasencia fault. On the southwestern edge of the study area there are outcrops of metamorphic rocks.

The geomorphological component was obtained based on the mapping of geomorphological units [31], synthesized in a mapping of geomorphological domains. The main geomorphological domains are as follows: summit surfaces and fluvial divides, embedded fluvial-glacial valleys, slopes and colluvial slopes, polygenic surfaces, hills and hillocks, glaciers, and fluvial terraces, alluvial and valley bottoms.

As for the vegetation that develops in the Jerte Valley, it is conditioned by the slopes of the hillsides, the differences in altitude, and the climate. In the areas located to the west, areas of low slope and glacial geomorphology, open formations such as pastures have developed. As one ascends in latitude and altitude there are large extensions of tree crops (cherry trees) very typical of this area. Halfway up the slope, both to the south and north of the river, there are large areas of wooded formations (pine forest, chestnut grove, oak

(cherry trees) very typical of this area. Halfway up the slope, both to the south and north of the river, there are large areas of wooded formations (pine forest, chestnut grove, oak

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 3 of 15

forest, ...). On reaching the timberline, there are different scrublands, and later, ascending in altitude, the summit pastures. Some summits have a steep slope which hinders the development of vegetation while the rock remains on the surface. forest, ...). On reaching the timberline, there are different scrublands, and later, ascending in altitude, the summit pastures. Some summits have a steep slope which hinders the development of vegetation while the rock remains on the surface.

**Figure 1.** Study area.

## *2.2. Methods*

**Figure 1.** Study area.

*2.2. Methods* 

The methodology followed (Figure 2) for landscape mapping was that described by Martínez-Graña [32], using a supervised classification of Landsat images as vegetation The methodology followed (Figure 2) for landscape mapping was that described by Martínez-Graña [32], using a supervised classification of Landsat images as vegetation mapping for the years 1994, 2000, 2010, 2019.

mapping for the years 1994, 2000, 2010, 2019. With ArcGIS 10.8 software, the mapping of homogeneous units was carried out, from the union between lithology and geomorphological domain mapping. Once the union of these mappings was carried out, 29 homogeneous units were obtained. To simplify these mappings, all units smaller than 2 hectares were filtered out, as they were not representative. The units with similar landscape development behavior, such as quartz slopes and colluviated slopes, basic rock slopes and colluviated slopes, and schist slopes and colluviated slopes, were unified. Finally, the 18 homogeneous units, most representative of the With ArcGIS 10.8 software, the mapping of homogeneous units was carried out, from the union between lithology and geomorphological domain mapping. Once the union of these mappings was carried out, 29 homogeneous units were obtained. To simplify these mappings, all units smaller than 2 hectares were filtered out, as they were not representative. The units with similar landscape development behavior, such as quartz slopes and colluviated slopes, basic rock slopes and colluviated slopes, and schist slopes and colluviated slopes, were unified. Finally, the 18 homogeneous units, most representative of the study area, were mapped.

study area, were mapped. Landsat images were downloaded from the U.S. Geological Survey (https://earthexplorer. usgs.gov/, accessed on 21 July 2021). These images were taken between the months of April, May, or June of the years 1994, 2000, 2010, and 2019. It was decided to use images from Landsat satellites since they were the free satellites with the highest resolution available during the entire study period. The use of images from other satellites, such as Sentinel 2, was considered, but only images from 2015 were available, so we would only have had the latest image to study and with a different resolution, so the surfaces to be compared could have varied. In the end, this option was discarded to homogenize the data as much as possible.

**Figure 2.** General methodology. **Figure 2.** General methodology.

Landsat images were downloaded from the U.S. Geological Survey (https://earthexplorer.usgs.gov/). These images were taken between the months of April, May, or June of the years 1994, 2000, 2010, and 2019. It was decided to use images from Landsat satellites since they were the free satellites with the highest resolution available during the entire study period. The use of images from other satellites, such as Sentinel 2, was considered, but only images from 2015 were available, so we would only have had the latest image to study and with a different resolution, so the surfaces to be compared could have varied. In the end, this option was discarded to homogenize the data as much as possible. Each image was subjected to radiometric calibration preprocessing, atmospheric corrections, and topographic corrections. Once the images were obtained, training areaswere selected after a field reconnaissance and with the help of orthophotos from previous years. Supervised classifications were performed on each image to obtain the land cover. It was decided to carry out the land cover mapping based on supervised classification since high-resolution orthophotos of the study years and field sampling of the last year were available and thus training areas could be delimited with great precision and a highly accurate mapping could be obtained.

Each image was subjected to radiometric calibration preprocessing, atmospheric corrections, and topographic corrections. Once the images were obtained, training areas were selected after a field reconnaissance and with the help of orthophotos from previous years. Supervised classifications were performed on each image to obtain the land cover. It was decided to carry out the land cover mapping based on supervised classification since highresolution orthophotos of the study years and field sampling of the last year were available and thus training areas could be delimited with great precision and a highly accurate Once the study areas were delimited, we proceeded to calculate the spectral signatures of each defined class, in order to extrapolate each pixel value of each class of the training area from the rest of the image. The supervised classification was performed by the maximum likelihood process, this method makes a statistical study (mean and standard deviation) of the pixel values of the training areas and calculates the probability of the values of the external indices to the training areas of belonging to one class or another, while the class with the highest probability is the one assigned to it.

mapping could be obtained. Once the study areas were delimited, we proceeded to calculate the spectral signatures of each defined class, in order to extrapolate each pixel value of each class of the After the supervised classification, a Majority filter was performed, thus filtering the neighboring contiguous cells of the larger size classes and avoiding the salt and pepper effect of supervised classification.

training area from the rest of the image. The supervised classification was performed by the maximum likelihood process, this method makes a statistical study (mean and standard deviation) of the pixel values of the training areas and calculates the probability of the To validate the historical land use mapping, the Kappa index of each supervised classification was calculated. For the mapping to be accepted for further study, the Kappa index had to have been greater than 0.75 [33,34].

values of the external indices to the training areas of belonging to one class or another, while the class with the highest probability is the one assigned to it. After the supervised classification, a Majority filter was performed, thus filtering the Land uses were classified into 10 main classes (water, forest, forest, cherry crops, treeless, scrub, snow, rocky, urban, burned area, clouds) to make the data homogeneous and to be able to analyze their parameters in the landscape as a whole.

neighboring contiguous cells of the larger size classes and avoiding the salt and pepper effect of supervised classification. To validate the historical land use mapping, the Kappa index of each supervised classification was calculated. For the mapping to be accepted for further study, the Kappa index had to have been greater than 0.75 [33,34]. The vegetation series mapping used was obtained from the Vegetation Series Map developed by Rivas Martínez in 1981 and revised in 1987 [35]. In the study area there are 5 vegetation series of the 37 existing in Spain. The vegetation series present in the study area have been reclassified (Table 1) in the denomination of the vegetation classes used in the present work in order to make a comparison of potential vegetation and real vegetation.

Land uses were classified into 10 main classes (water, forest, forest, cherry crops, treeless, scrub, snow, rocky, urban, burned area, clouds) to make the data homogeneous

and to be able to analyze their parameters in the landscape as a whole.


**Table 1.** Reclassification of the vegetation series to the legend of the vegetation mapping used for the study. Index (CON-*Pi* = proportion of the landscape occupied by patch type; *gik* = number of adjacencies between pixels of between classes, in relation to the maximum possible consid-

ೖ <sup>∑</sup> ೖ ೖసభ

ଶ ୪୬() ] × (100)

ቇ]

It measures the percentage of adjacency

tesserae of the same

types in a given area.

CONTAG=

∑ ൬

ೖ <sup>∑</sup> ೖ ೖస ] సభ × [୪୬ቈ ×

Then, the corresponding unions of the homogeneous unit mappings with the different vegetation mappings obtained from the vegetation series map and supervised classifications from Landsat images of the years 1994, 2000, 2010, 2020 were carried out, giving rise to the natural unit mappings. tance Index (CON-NECT) *cijk* = joining between patch *j* and *k* of the same patch type, based on a user-specified threshold distance; *ni* = number of patches in the landscape of each patch type (*i*). the total tesserae or of a class connected according to a threshold distance.

మ <sup>൰</sup> సభ

In the mapping of natural units, 53 units were obtained and simplified in the following way: Proximity PROX = <sup>∑</sup> ೕ ೕ మ ୀଵ Sum of the areas of

1. Very small units, smaller than 2 hectares, were eliminated. Index *aijs* = area (m2) of patch *ijs* within specified neighbor-

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 7 of 15

∑ ∑ [ × ೖసభ

Contagion

[1 +


**Figure 3.** 3D modeling of landscape units with reconnaissance field images of the most representative landscape units. **Figure 3.** 3D modeling of landscape units with reconnaissance field images of the most representative landscape units.


**Table 2.** Ending landscape unit.

Finally, with the data collected, with the software Fragstat v4.2.1 [36] the analyses of the spatial and ecological patterns of the Jerte Valley were carried out. Fragstat is a software that allows the quantification of landscape patterns, which is a prerequisite for the study of landscape relationships with the processes that degrade it. Fragstat calculates a set of indices and variables that quantitatively describe the level of fragmentation and spatial distribution of land use and land cover.

These analyses were carried out with the help of indices that describe the different classes or patches of the units or allow a description of the landscape as a whole. The indices used are those described in Table 3.


**Table 3.** Index Fragstat.

**Table 3.** *Cont.* **Index Observations** Patch Cohesión Index (COHESION) COHESION = [1 − ∑ *n <sup>j</sup>*=*<sup>i</sup> P* ∗ *ij* ∑ *n <sup>j</sup>*=<sup>1</sup> *P* ∗ *ij*√*<sup>a</sup>* ∗ *ij* ] × h 1 − <sup>√</sup> 1 *Z* i−<sup>1</sup> × (100) *Pij \** = perimeter of patch *ij* in terms of number of cell surfaces; *aij \** = area of patch *ij* in terms of number of cells. *Z* = total number of cells in the landscape. It measures the physical connectivity of the analysed category. Landscape Shape Index (LSI) LSI = <sup>25</sup>*<sup>E</sup>* ∗ √ *A E*\* = total length (m) of edge in landscape; includes the entire landscape boundary and some or all background edge segments; *A* = total landscape area (m<sup>2</sup> ). Provides a standardised measure of total edge or edge density to suit the size of the landscape. Contagion Index (CONTAG) CONTAG = [1 + ∑ *m <sup>i</sup>*=<sup>1</sup> ∑ *m k*=1 [*Pi*<sup>×</sup> *g ik* ∑ *m k*=*l g ik* ]×[ln (*Pi*× *g ik* ∑ *m k*=1 *g ik* )] 2 ln(*m*) ] × (100) *P<sup>i</sup>* = proportion of the landscape occupied by patch type; *gik* = number of adjacencies between pixels of patch types *i* and *k* based on the double-count method; *m* = number of patch types present in the landscape, including the landscape border if present. It measures the percentage of adjacency between classes, in relation to the maximum possible considering the frequency of these. Landscape Division Index (DIVISION) DIVISION = [1 − *m* ∑ *i*=1 *n* ∑ *j*=1 *aij A* 2 ] *aij* = area (m<sup>2</sup> ) of patch *ij*. *A* = total landscape area (m<sup>2</sup> ). Probability that two areas of the landscape are not located in the same habitat fragment. Connectance Index (CONNECT) CONNECT = [ ∑ *m <sup>i</sup>*=<sup>1</sup> ∑ *n j*=*k cijk* ∑ *m i*=1 *ni*(*ni*−<sup>1</sup>) 2 ](100) *cijk* = joining between patch *j* and *k* of the same patch type, based on a user-specified threshold distance; *n<sup>i</sup>* = number of patches in the landscape of each patch type (*i*). It is the percentage of the total tesserae or of a class connected according to a threshold distance. Proximity Index (PROX) PROX = *n* ∑ *g*=1 *aijg h* 2 *ijg aijs* = area (m<sup>2</sup> ) of patch *ijs* within specified neighborhood (m) of patch *ij*; *hijs* = distance (m) between patch *ijs* and patch *ijs*, based on patch edge-to-edge distance, computed from cell center to cell center. Sum of the areas of tesserae of the same class whose edges are at a specific radius. Fragmentation (F) F = *TA NP*×2×*ENNMN*×( *PD π* ) Spatial disaggregation of patches or habitat types in a given area.

## **3. Results**

The supervised rankings performed for the years 1994, 2000, 2010, and 2019 obtained a Kappa index higher than 0.75 (Table 4), so the mapping was accepted for study:



The changes in land use in the Jerte Valley can be seen in Figure 4. The land use changes obtained in the Jerte Valley can be seen in Figure 3. There is a clear increase in the cultivation of cherry trees, a minimal loss in the extent of forests, with a decrease in scrubland in the years 2000 and 2020 where wildfires were notable.

**3. Results** 

The supervised rankings performed for the years 1994, 2000, 2010, and 2019 obtained

**Table 4.** Kappa index and overall accuracy for supervised classification (material complementary).

**Year Parámetros de Imágenes Kappa Index Overall Accuracy** 

The changes in land use in the Jerte Valley can be seen in Figure 4. The land use changes obtained in the Jerte Valley can be seen in Figure 3. There is a clear increase in the

The development of human activities in the environment has greatly conditioned the

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 9 of 15

a Kappa index higher than 0.75 (Table 4), so the mapping was accepted for study:

1994 LT05\_L1TP\_202032\_19940314\_20180217\_01\_T1 0.92 0.93 2010 LT05\_L1TP\_202032\_20100310\_20161016\_01\_T1 0.76 0.80 2000 LT05\_L1TP\_202032\_20000314\_20180312\_01\_T1 0.88 0.91 2019 LC08\_L1TP\_202032\_20190303\_20190309\_01\_T1 0.86 0.88

land in the years 2000 and 2020 where wildfires were notable.

**Figure 4.** Evolution of land use in the last 4 decades in the Jerte Valley. **Figure 4.** Evolution of land use in the last 4 decades in the Jerte Valley. According to the values obtained in the diversity (SHDI and SIDI) and evenness

Based on the maps of landscape units for the study years: 1994, 2000, 2010, and 2019, the data obtained were analyzed (Figure 5). The statistical study began in 1994. Taking as a reference the vegetation series of Rivas Martínez, as potential vegetation, it can be observed that the proportion of land occupied by tree formations is much lower than would be expected (Table 5). Based on the maps of landscape units for the study years: 1994, 2000, 2010, and 2019, the data obtained were analyzed (Figure 5). The statistical study began in 1994. Taking as a reference the vegetation series of Rivas Martínez, as potential vegetation, it can be observed that the proportion of land occupied by tree formations is much lower than would be expected (Table 5). lost in the landscape pattern. The number of patches is decreasing, becoming larger, more compact, and with simpler shapes. The grouping of tesserae causes a decrease in ecotones and has direct repercussions on the richness of species present in the environment. This evolution is largely due to the increase in agroforestry practices, forest management, and disturbances such as forest fires.

(SHEI) indices, in the last three decades (Figure 6), heterogeneity and complexity are being

**Figure 5.** Multi-temporal spatial distribution of the natural units in the study area. **Figure 5.** Multi-temporal spatial distribution of the natural units in the study area.


**Table 5.** Land cover indices in the Jerte Valley, in the years 1994, 2000, 2010, 2019 and in the vegetation series.

Table 5 also shows how, in 1994, 28% of the surface area of the Jerte Valley was occupied by wooded formations (unit 2), such as oak (*Quercus pyrenaica*), chestnut (*Castanea sativa*), and pine (*Pinus sylvestris*), as opposed to 82% according to the optimum development of the vegetation series. It can also be observed that the presence of grassland areas (unit 3) is higher (13%) due to the traditional cultural uses of the land and its use for grazing, initial stages of recolonization of areas affected by forest fires, deforestation, etc. This stage of landscape development was taken as the initial situation, year 1994, and from there, the changes in the landscape over the last decades were studied.

The development of human activities in the environment has greatly conditioned the landscape in recent decades, occupying the optimal potential areas for the development of tree formations with cherry plantations (28.95%). On the other hand, there is an increase of 57% in the cultivation of cherry trees (unit 2) to the detriment of pasture and scrubland areas, 9% (unit 4 and 5).

The analysis of landscape structure was based on territorial changes in vegetation cover, geomorphology, and land use in recent decades. It was necessary to analyze land cover together with the models and management systems adopted in order to correctly assess the significance of the changes detected.

According to the values obtained in the diversity (SHDI and SIDI) and evenness (SHEI) indices, in the last three decades (Figure 6), heterogeneity and complexity are being lost in the landscape pattern. The number of patches is decreasing, becoming larger, more compact, and with simpler shapes. The grouping of tesserae causes a decrease in ecotones and has direct repercussions on the richness of species present in the environment. This evolution is largely due to the increase in agroforestry practices, forest management, and disturbances such as forest fires.

Once the spatial patterns and metrics of area, density, and landscape aggregation were analyzed (Table 6), it could be said that there is an increase in homogeneity at the landscape level in the study area, which corroborates the results obtained in reference to diversity.

The number of patches has almost halved from 16,023 in 1994 to 8183 in 2019. This has caused the patch density to decrease by 50%. Their average area of occupancy has doubled, indicating that we have fewer and larger patches. Likewise, the cohesion index, Landscape Shape Index, and the Interspersion and Juxtaposition Index have increased, which means that the masses maintain a moderately high degree of intermixing despite their tendency to homogenization.

**Figure 6.** Diversity indices, the thick lines are the values of the indices in the different years, the thinner lines are the trend of each index. (**A**) Graph of Simpson's Diversity Index, (**B**) Graph of Shannon's Evenness Index, (**C**) Graph of Shannon's Diversity Index. **Figure 6.** Diversity indices, the thick lines are the values of the indices in the different years, the thinner lines are the trend of each index. (**A**) Graph of Simpson's Diversity Index, (**B**) Graph of Shannon's Evenness Index, (**C**) Graph of Shannon's Diversity Index.

**Table 6.** Table of area, density and landscape aggregation metrics obtained by the FRAGSTAT software.


Total area (TA) 37,598 37,598 37,598 37,598 Number of Patches (NP) 16,023 15,939 9820 8183 The isolation of the spots increases, with a 74% increase of the proximity index and a 17% increase of the Euclidean Nearest-Neighbor Distance Mean. While the degree of fragmentation decreases by 200% (0.15, remember that it is inverse).

**Area, Density and Aggregation Metrics 1994 2000 2010 2019** 

Patch Density (PD) 40.48 40.31 25.48 20.69 Radius of Gyration (GYRATE\_MN) 47.21 48.90 49.29 57.87 Landscape Shape Index (LSI) 83.32 86.71 50.54 55.25 The study of the structure of the landscape at the class level makes it possible to discern the role of each tesserae in the environmental mosaic, assigning to each group a diversifying or homogenizing function in the whole.

Mean Patch Area (AREA\_MN) 2.47 2.4805 3.92 4.83 Largest Patch Index (LPI) 16.72 10.53 27.25 23.33 Euclidean Nearest-Neighbor Distance Mean (EMM\_MN) 141.70 142.01 149.61 166.43 Contagion Index (CONTAG) 30.20 29.10 39.3 40.85 Interspersion and Juxtaposition Index (IJI) 73.55 73.05 70.18 73.00 The most outstanding results of the study in this area are the increase in the area under cultivation (+7000 ha), dedicated almost exclusively to cherry (Table 7). The tesserae that make up this class increase in number, but maintain the same patterns of isolation, shape, and area. Wooded areas decrease slightly but show changes in their arrangement indicating that the wooded tesserae are larger, less fragmented, but more isolated. On the other hand, units 4 and 5 decrease (3000 ha and 1000 ha respectively) and present very small, isolated, and fragmented tesserae. The open formations (unit 7), due to their singular character, being pasture areas subjected to a strong anthropic influence, require a

Connectance Index (CONNECT) 0.36 0.35 0.48 0.48

Patch Cohesion Index (COHESION) 96.94 96.20 98.38 98.02

The study of the structure of the landscape at the class level makes it possible to dis-

The most outstanding results of the study in this area are the increase in the area

cern the role of each tesserae in the environmental mosaic, assigning to each group a di-

under cultivation (+7000 ha), dedicated almost exclusively to cherry (Table 7). The tesserae that make up this class increase in number, but maintain the same patterns of isolation,

Proximity Index (PROX) 169 124 551 295

Fragmentation Index (F) 0.07 0.07 0.16 0.22

versifying or homogenizing function in the whole.

detailed analysis and study. Large variations in fragmentation are observed, which can be explained by the agrosilvopastoral treatments to which this habitat is subjected.

**Table 7.** Table of area, density, and aggregation metrics of the patches obtained by the FRAGSTAT software. Where: TA: Total area; NP: Number of Patches; PD: Patch Density; AREA\_MN: Mean Patch Area; IJI: Interspersion and Juxtaposition Index; Frag: Fragmentation Index.


Valuable assessments can be obtained from the yearly detailed analysis. The year 2010 deserves a special comment since the initial data used for the mapping of landscape units contained parts covered by snow (unit 6). Snow represents an important and common disturbance in vegetation and, in general, in mountain areas. It highlights changes in landscape heterogeneity in annual and sub-annual space. In this particular case, it causes

notably high levels in the cohesion and proximity index, and a decrease in the interdispersion and juxtaposition index. In addition, the area occupied by bare soil (snow) increases significantly to the detriment of areas of scrub and grassland, typical of mountain areas where snow accumulates.

Finally, it should be noted that the evolution of urban infrastructures shows a stable, moderate growth and does not represent a problem for the landscape, with diffuse anthropic pressure being practically non-existent, apart from some isolated agricultural and livestock constructions.

#### **4. Discussion**

The identification of changes in land use over the last 50 years makes it possible to identify the environmental impact that has occurred in the territory. The decrease in the area of snow, the increase in crop and urban areas, demonstrate an impact generated by anthropogenic activity in the area, mainly by the cherry tree cultivation activity. In other landscape evolution simulation studies, remote sensing techniques were used [33], with classification metrics between 80 and 90% of the images studied for the classification of the extent of tropical forest fragmentation [3] and a supervised classification methodology and landscape indices for landscape fragmentation simulations [37].

The analysis of the areas of occupation of the different types of land use allows a preliminary evaluation of the evolution of the territory. As previously mentioned, the landscape is a dynamic entity in continuous spatial and temporal evolution [38]. The environmental mosaic present in the Jerte Valley is quite heterogeneous and, therefore, the agents that cause these changes are diverse. The most important are: geomorphological agents in the areas of rock and bare soil, where they are the main modelers of the landscape together with the climate; and, in the areas dominated by vegetation, (geobotanical landscape) it is the plants that modify the edaphoclimatic conditions of the area. In addition, we are in a territory with a growing human influence, which determines a much more accelerated time scale of changes due to the disturbances caused by their activities.

Curiously, as a consequence of the change in climatic conditions, the optimal areas for cherry cultivation are no longer the valley bottoms and slopes, but now occupy the higher areas (scrub and pasture). It is a phenomenon marked in other territories with other fruit crops and especially for grapevine [39]. Ultimately it is a consequence of climate change. Presumably, this trend will continue until the lithological and edaphological conditions mean that the land cannot be cultivated, even with the application of the current techniques of farming [40].

Another of the recurrent landscape disturbances in the Jerte Valley and in other areas with a Mediterranean climate are forest fires [41]. At an ecological level, they can be beneficial for the maintenance of the forest structure if they do not occur too frequently, in limited areas, or in areas with a subarid climate. In these cases, they cause simplification of the environmental mosaics [42]. In the study area, forest fires in recent years have mainly affected scrub and grassland areas (see the decrease in the area of occupancy for the year 2020), minimizing their ecological effects, compared to the foreseeable consequences of affecting complex tree formations.

The naturalization of pastures, together with the growth of cherry orchards, has caused the landscape mosaic to change and present more homogeneity and less diversity. It is an unstable and singular situation, since at the same time the natural and cultural landscapes are advancing, giving rise to a common result. This situation can lead to the generation of a cyclical dynamic based on the recurrent occurrence of fires, which is not typical of this territory and can lead to the loss of biological and environmental resources, increased erosion [43], ecological imbalances, etc.

The expansion of cherry trees towards the highlands could be a partial mitigation of the effects described above, enhancing the coexistence of valuable cultural and natural elements. This would favor the maintenance of anthropized areas with lower fire risk and greater economic and ecological potential. Another favorable measure would be

the maintenance of extensive livestock farming, both in pasture areas and in mountain pastures, since it plays a fundamental role in the conservation of these habitats, as well as the geobotanical, genetic, and landscape diversity found in them [44].

Dehesa grassland systems have been found to present a high and changing fragmentation of the stippled or mottled type due to the type of management. They are very heterogeneous but harbor a greater diversity than other equivalent potential systems and are therefore a priority for conservation.

#### **5. Conclusions**

Previous work has been carried out combining landscape metrics and remote sensing, however, changes in land use have not been studied with this methodology. In the present study, the importance of land use changes in landscape dynamics is highlighted by systematically analyzing landscape evolution using spatial patterns and indices of diversity, dominance, form, and fragmentation. The scenarios analyzed reveal a trend towards homogenization of the territory with loss of diversity and changes in natural dynamics, largely caused by human influence (cherry tree cultivation) and major disturbances created by forest fires.

The research has shown that the analysis of landscape patterns in a multi-temporal study allows both the analysis of changes in past land uses and the promotion of actions for the future conservation of land uses that are positive for the landscape, allowing compliance with conservation strategies, as well as the evaluation of the causes and consequences of large-scale actions or disturbances in the landscape. In addition, the use of geographic information systems, remote sensing, and spatial analysis software allows the evaluation and integration of many geo-environmental parameters that form the landscape, enabling the development of a base cartography of superior detail and quality; and enables the spatiotemporal and statistical analysis of landscape units and their effect on the environment. This methodology can be implemented in different regions with similar characteristics, large areas with a large number of classes, remote regions, difficult to access, etc. For the future, we intend to apply this methodology for shorter time intervals and to see the gradual modification of the landscape for annual periods to observe whether the landscape modifications are gradual or abrupt. It is a fundamental tool for the proper management of land use, land use planning, and environmental conservation. Institutions and territorial managers will be able to adapt the policies and programs of each region to each present and future scenario, adapting this generic methodology to particular situations, thus achieving objective and comparable parameters, extremely useful for decision making.

**Author Contributions:** Conceptualization, Y.S.S. and J.L.R.R.; methodology, Y.S.S., J.L.R.R., and A.M.G.; software, Y.S.S. and J.L.R.R.; validation, Y.S.S. and F.S.-F.; formal analysis, A.M.G.; investigation, Y.S.S. and M.C.; resources, A.M.G.; data curation, Y.S.S. and A.M.G.; writing—original draft preparation, Y.S.S.; writing—review and editing, A.M.G.; visualization, M.C.; supervision, A.M.G. and F.S.-F.; project administration, A.M.G.; funding acquisition, A.M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the project SA044G18 of the Regional Government of Castilla y Leon, and the GEAPAGE research group (Environmental Geomorphology) of the University of Sala-manca.

**Acknowledgments:** This research was funded by the project SA044G18 of the Regional Government of Castilla y Leon, and the GEAPAGE research group (Environmental Geomorphology) of the University of Salamanca.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

