Assessing Land-Cover Changes in the Natural Park ‘Fragas do Eume’ over the Last 25 Years: Insights from Remote Sensing and Machine Learning

Citation: Díaz-García, P.; Regos, A. Assessing Land-Cover Changes in the Natural Park ‘Fragas do Eume’ over the Last 25 Years: Insights from Remote Sensing and Machine Learning. Land 2024,13, 1601. https://doi.org/10.3390/ land13101601 Academic Editor: Paul Aplin Received: 14 May 2024 Revised: 25 June 2024 Accepted: 25 September 2024 Published: 1 October 2024 Copyright: © 2024 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/). land Article Assessing Land-Cover Changes in the Natural Park ‘Fragas do Eume’ over the Last 25 Years: Insights from Remote Sensing and Machine Learning Paula Díaz-García1and Adrián Regos 2,3,4,* 1Facultade de Ciencias, Universidade da Coruña, 15071 A Coruña, Spain; [email protected] 2Misión Biológica de Galicia del Consejo Superior de Investigaciones Científicas (MBG-CSIC), 36143 Santiago de Compostela, Spain 3Centre de Ciencia i Tecnologia Forestal de Catalunya, Ctra. St. Llorenç de Morunys km 2, 25280 Solsona, Spain 4BIOPOLIS Program in Genomics, Biodiversity and Land Planning, CIBIO/InBIO (Research Centre in Biodiversity and Genetic Resources), Campus de Vairão, 4485-661 Vairão, Portugal *Correspondence: [email protected] or [email protected] Abstract: The ‘Fragas do Eume’ Natural Park includes one of the best-preserved Atlantic forests in Europe. These forests are part of the Natura 2000 Network. This scientific study focuses on analysing land-cover changes in the ‘Fragas do Eume’ Natural Park (NW Spain) over a 25-year period, from 1997 to 2022, using machine learning techniques for the classification of satellite images. Several image processing operations were carried out to correct radiometry, followed by supervised classification techniques with previously defined training areas. Five multispectral indices were used to improve classification accuracy, and their correlation was evaluated. Land-cover changes were analysed, with special attention to the transitions between eucalyptus plantations and native deciduous forests. A significant increase in eucalyptus plantations (48.2%) (Eucalyptus globulus Labill.) was observed, while native deciduous forests experienced a decrease in their extent (17.6%). This transformation of the landscape affected not only these two habitats, but also cropland and scrubland areas, both of which increased. Our results suggest that the lack of effective conservation policies and the economic interest of fast-growing tree plantations could explain the loss of native deciduous forests. The results highlight the need to implement pro-active and sustainable management measures to protect these natural forest ecosystems in the ‘Fragas do Eume’ Natural Park. Keywords: Atlantic deciduous forests; Galicia; landsat satellite mission; long-term analysis; image classification; artificial intelligence; remote sensing 1. Introduction Biodiversity is essential for the proper functioning of ecosystems and the maintenance of services such as climate regulation, cultural, and provisioning [ 1 ]. Ecosystems with higher biodiversity are more resilient to changes [ 2 ], underscoring the importance of its conservation. However, in recent decades, human activities have significantly impacted biodiversity with primary causes of loss including pollution, invasive species, climate change, overexploitation, and land use changes [3]. To address biodiversity decline, 195 countries, including Spain, ratified the Convention on Biological Diversity [ 4 ]. In Europe, biodiversity protection follows guidelines like the Natura 2000 Network, under the Habitats Directive (Directive 92/43/EEC) [ 5 ] and the Birds Directive (Directive 2009/147/EC) [ 6 ]. These aim to conserve species and habitats, including the ‘Atlantic deciduous forests’ which extend along the European coast from Norway to Portugal [ 7 ]. These forests, composed of deciduous trees, such as oaks, beeches, birches, and alders [ 8 , 9 ], face threats from deforestation and agricultural intensification [ 10 , 11 ]. In Galicia, the native ‘Fragas’ forests, featuring a variety of tree species typical of Atlantic forests and a rich diversity of mosses, lichens, and Land 2024,13, 1601. https://doi.org/10.3390/land13101601 https://www.mdpi.com/journal/land Land 2024,13, 1601 2 of 13 ferns, are especially significant. These include priority habitats like the alluvial forests of Alnus glutinosa and Fraxinus excelsior (91E0*) [ 10 ]. ‘Fragas’ are dominated by two species of Quercus:Quercus robur L., the oak, or Quercus pyrenaica Willd., the pyrenean oak [12,13]. In recent decades, Galicia has witnessed a notable surge in eucalyptus monoculture, primarily driven by Eucalyptus globulus Labill., now covering 26% of forest areas—a doubling from 2002’s 13% [ 14 , 15 ]. This expansion exacerbates habitat fragmentation issues [ 8 ]. Eucalyptus plantations, dominated by Eucalyptus globulus, thrive due to their rapid growth but deplete soil nutrients [ 16 , 17 ] and bacterial diversity [ 18 ], hindering the recovery and future cultivation in affected areas [ 19 ]. The impact extends to the ‘Fragas do Eume’ Natural Park, where E. globulus presence has significantly risen, especially before its declaration as protected area [ 20 , 21 ]. This natural park, spanning 9125 hectares [ 22 ] in the Galicia’s Eume river basin, earned natural park status in 1997, protecting one of the best-preserved Atlantic forests in Europe [ 23 ]. Classified as a SAC (special area of conservation) in the Natura 2000 Network, it harbours diverse habitats and host species dating back to the Tertiary period. These habitats are: Atlantic humid oakwoods with Erica ciliaris and Erica tetralix (4020*); tree-heath forests of Laurus nobilis (5230*); pseudo-steppes with grasses and annuals of the Thero-Brachypodietea order (6220*); nardus grasslands with numerous species on siliceous substrates in mountainous areas (and sub-mountainous areas of continental Europe) (6230*); active raised bogs (7110*); slopes, screes, and ravines with Tilio-Acerion forests (9180*); and alluvial forests with Alnus glutinosa and Fraxinus excelsior (Alno-Padion, Alnion incanae, Salicion albae) (91E0*) (Directive 92/43/EEC, 1992, Annex I) [ 5 , 13 ]. The Fragas are home to threatened or vulnerable species, such as bryophytes like Sphagnum pylaesii, invertebrates like Elona quimperiana and Margaritifera margaritifera L., amphibians like Chioglossa lusitanica, and mammals like Galemys pyrenaicus [ 24 ]. It also harbours species of ferns that have survived since the Tertiary period. These species are subtropical in origin and survive thanks to the high humidity and mild temperatures of the environment, such as Culcita macrocarpa and Woodwardia radicans [23,25]. The objective of this study is to quantify the changes in the main land uses and land covers that have occurred in the ‘Fragas do Eume’ Natural Park since its declaration as a protected area to the present day. The decline of Atlantic forest in the ‘Fragas do Eume’ Natural Park is primarily driven by the expansion of eucalyptus plantations. However, we expect that this expansion should occur at a lower rate than before the park’s designation as a protected area due to its legal protection. If the legal protections are effective, we would also expect a halt in the expansion of Eucalyptus spp. plantations and a natural recovery of the native Atlantic forest over the last 25 years. More specifically, our study aims to document land cover changes within the park and introduce an innovative methodology that leverages machine learning techniques and multispectral indices to enhance the accuracy of satellite image classification. To deal with the uncertainty associated with the classification algorithms and the need to communicate such uncertainty to decision makers, we assess land-cover changes by using a wide range of classification methods and spectral information. This methodology has the potential for application in similar studies at both global and macroregional levels. 2. Materials and Methods 2.1. Study Area The study area is the ‘Fragas do Eume’ Natural Park, located in Galicia, in the northern part of the province of A Coruña (Figure 1). It covers 9125 hectares (Decree 218, 1997, Article 2) [ 22 ] extending throughout the Eume river basin. It spans across 5 municipalities: Cabanas, Pontedeume, A Capela, Monfero, and As Pontes. The Fragas has been a natural park since 1997 (Decree 218, 1997, Article 2) [ 22 ], being one of the best-preserved Atlantic forests in Europe [ 23 ]. It is also a SAC (special area of conservation) of the Natura 2000 Network (Directive 37, 2014, Annex I) [ 26 ], since it harbours different habitats of community importance (Directive 92/43/EEC, 1992, Annex I) [5]. Land 2024,13, 1601 3 of 13 Land 2024, 13, x FOR PEER REVIEW 3 of 13 Network (Directive 37, 2014, Annex I) [26], since it harbours different habitats of community importance (Directive 92/43/EEC, 1992, Annex I) [5]. Figure1.Map showing the location of ‘Fragas do Eume’. The orange area indicates the extent of the Natural Park. 2.2.MethodologicalFramework To achieve the objective of our study, remote sensing and machine learning techniques were employed. Remote sensing allows us to obtain data from the Earth’s surface using satellites [27]. The interaction between the Earth’s surface and the atmosphere can be recorded by satellite sensors [27]. Each entity or object on the Earth’s surface has a spectral signature, i.e., different electromagnetic characteristics, which enables its identification [28]. The information is collected by airborne sensors that measure and receive radiometric information. The spectral and spatial resolution at which these sensors are capable of recording the electromagnetic energy emitted by the sun and reflected by objects will determine the level of detail that can be achieved during the classification process of these images [29,30]. Methodologically, the study explores the accuracy in the supervised classification of images of the Earth’s surface obtained by NASA’s Landsat satellites using different machine learning techniques, as well as the effect of incorporating auxiliary information to the spectral bands. Thematically, a relationship will be established between the main types of land uses and land covers that are spectrally discernible and the priority habitats identified in the natural park in order to infer the possible impacts on protected habitats since its declaration as a protected area. 2.2.1. Pre-Processing of Satellite Images To detect changes in land use in the ‘Fragas do Eume’ Natural Park between 1997 and 2022, we used 2 cloud-free images for each year. These images were obtained from Figure 1. Map showing the location of ‘Fragas do Eume’. The orange area indicates the extent of the Natural Park. 2.2. Methodological Framework To achieve the objective of our study, remote sensing and machine learning techniques were employed. Remote sensing allows us to obtain data from the Earth’s surface using satellites [ 27 ]. The interaction between the Earth’s surface and the atmosphere can be recorded by satellite sensors [ 27 ]. Each entity or object on the Earth’s surface has a spectral signature, i.e., different electromagnetic characteristics, which enables its identification [ 28 ]. The information is collected by airborne sensors that measure and receive radiometric information. The spectral and spatial resolution at which these sensors are capable of recording the electromagnetic energy emitted by the sun and reflected by objects will determine the level of detail that can be achieved during the classification process of these images [ 29 , 30 ]. Methodologically, the study explores the accuracy in the supervised classification of images of the Earth’s surface obtained by NASA’s Landsat satellites using different machine learning techniques, as well as the effect of incorporating auxiliary information to the spectral bands. Thematically, a relationship will be established between the main types of land uses and land covers that are spectrally discernible and the priority habitats identified in the natural park in order to infer the possible impacts on protected habitats since its declaration as a protected area. 2.2.1. Pre-Processing of Satellite Images To detect changes in land use in the ‘Fragas do Eume’ Natural Park between 1997 and 2022, we used 2 cloud-free images for each year. These images were obtained from NASA’s Landsat satellite mission, by the Landsat 5 TM sensor on 27 March and 2 August 1997, and by the Landsat 8 OLI sensor on 24 March and Landsat 9 OLI sensor on 22 July 2022. These images are available on the website of the Earth Resources Observation and Science (EROS) Centre of the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/, accessed on 1 October 2022) [ 31 ]. These images are projected in the WGS 84 coordinate Land 2024,13, 1601 4 of 13 system, UTM projection, Zone 29 North. We chose images from winter and summer to better differentiate between deciduous and evergreen vegetation due to the different phenology of the species, which implies differential spectral behaviour between the two types of land covers. The images were downloaded with a Level 1T processing level, meaning they underwent geometric correction considering the terrain. Next, we performed the radiometric calibration and correction to improve images affected by parameters related to the attenuation of electromagnetic radiation from the sun [ 32 ]. The apparent reflectance method (“apref”) or top-of-atmosphere reflectance was applied. This calibration was done in R Studio (https://www.r-project.org/, accessed on 1 October 2022) using the ‘RStoolbox’ package (version 1.0.0) for all layers except the July layer of Landsat 9 OLI, which was processed in the QGIS program (version 3.38.3) (https://qgis.org/es/site/, accessed on 1 October 2022 ) because the ‘RStoolbox’ package did not have an update to process Landsat 9 OLI layers. 2.2.2. Supervised Classification Procedure We classified land-cover types from the radiometric information of the spectral bands. Additionally, auxiliary information was incorporated into the classification process in order to improve the accuracy of the final land cover maps. Specifically, we computed 5 multispectral indices : Normalised Difference Water Index (NDWI2), Modified Normalised Difference Water Index (MNDWI), Normalised Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Modified Soil Adjusted Vegetation Index (MSAVI). NDWI2 and MNDWI are indices related to measuring the amount of water in vegetation or soil moisture saturation levels (Table 1). NDVI, EVI, and SAVI measure the intensity of radiation from different bands of the electromagnetic spectrum emitted by vegetation (Table 1). Table 1. List of indices with their corresponding formulas and justifications for selection: Normalised Difference Water Index (NDWI2), Modified Normalised Difference Water Index (MNDWI), Normalised Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Modified Soil Adjusted Vegetation Index (MSAVI). Index Formula Rational NDWI2 (Green −NIR)/(Green + NIR) NDWI2 is used to enhance the presence of water bodies by taking advantage of the strong absorption of shortwave infrared (SWIR) and high reflectance of green wavelengths by water. MNDWI (Green −SWIR)/(Green + SWIR) MNDWI improves the delineation of water bodies in environments where NDWI might be less effective, particularly in built-up areas NDVI (NIR −Red)/(NIR + Red) NDVI is widely used to quantify vegetation health and cover by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs). EVI G×((NIR −Red)/((NIR) + (C1 ×Red) −(C2 ×Blue) + L)) EVI is designed to optimise the vegetation signal with improved sensitivity in high biomass regions and better correction for atmospheric conditions and soil background signals. SAVI (NIR −Red) ×(1 + L)/(NIR + Red + L) MSAVI is tailored to minimise soil brightness influences, making it more effective in areas with sparse vegetation cove To avoid using redundant information in the classification procedure, we calculated the Pearson correlation among the spectral indices by using the ‘cor’ function from the ‘stats’ package in RStudio. For better visualisation of the data, we created a heatmap indicating the degree of correlation between each pair of indices. Land 2024,13, 1601 5 of 13 Supervised classification is a process in remote sensing where an algorithm categorises pixels in an image into predefined classes based on training samples provided by the user. The user selects representative samples of each class from the satellite image, which are used to train the classification algorithm. The algorithm then analyses the spectral characteristics of the pixels in the image and assigns them to the most similar class based on the training samples [ 33 ]. Training areas were defined using QGIS software (https://qgis.org/es/site/, accessed on 1 October 2022) by selecting a set of pixels coinciding with well-known areas from Landsat images for each year. The areas were chosen from the natural park and its adjacent areas. We selected and identified a minimum of 10 areas/polygons per class. Those areas were homogeneous while still capturing environmental heterogeneity to account for spectral variability within each class. The chosen classes include up to 18 habitats from Annex I of Directive 92/43/EEC [ 5 ] (see Table 2). Evergreen forests are defined as dense, forested areas where the majority of trees retain their leaves throughout the year. In our study area, this typology is mainly composed of Eucalyptus globulus. Deciduous forests are characterised by trees that shed their leaves annually, typically in the autumn. These forests experience seasonal changes, with trees growing new leaves in the spring, which are then shed in the fall to conserve water and energy during the winter. In the study area, these forests are dominated by oaks. We created different training areas for land cover type and year (see Table 3). Table 2. List of habitats of the Annex I of Directive 92/43/EEC included in the land-cover classes identified in the Fragas do Eume Natural Park. * indicates priority habitat listed in the Annexes of the Habitat Directive. Land Cover Classes Habitat Code Natural Habitats Listed in the Habitats Directive Water areas 3160 Natural dystrophic lakes and ponds 3260 Lowland to montane rivers with vegetation of Ranunculion fluitantis y de Callitricho-Batrachion 3270 Rivers with muddy banks with Chenopodion rubri p.p. and Bidention p.p. vegetation Shrubland 4030 European dry heaths * 5230 Arborescent matorral with Laurus nobilis Croplands * 6220 Pseudo-steppes with grasses and annuals of the Thero-Brachypodietea * 6230 Species-rich Nardus grasslands on siliceous substrates in mountain areas (and sub-mountain areas in Continental Europe) 6410 Molinia meadows on calcareous, peaty, clay-like, or silt-laden soils (Molinion caeruleae) 6430 Hydrophilic tall herb fringe communities of the plains and of the montane-to-alpine levels 6510 Lowland hay meadows (Alopecurus pratensis, Sanguisorba officinalis) Deciduous forest * 9180 Tilio-Acerion forests of slopes, screes and ravines * 91E0 Alluvial forests with Alnus glutinosa and Fraxinus excelsior (Alno-Padion, Alnion incanae, Salicion albae) 9230 Galicio–Portuguese oak woods with Quercus robur and Quercus pyrenaica 9260 Castanea sativa woods Evergreen forest Land 2024,13, 1601 6 of 13 Table 3. Number of training areas per habitat in the years 1997 and 2022 and the total hectares of each one. Habitat 1997 2022 1997 2022 Polygons Hectares Water 12 18 43 51 Deciduous forest 11 14 89 61 Shrubland 14 20 73 92 Croplands 16 21 7 19 Evergreen forest 11 11 27 53 The supervised classification based on these pre-defined landcover classes was performed by using 9 machine learning algorithms available from the ‘Caret’ package in R: avNNet (model averaged neural network), gbm (stochastic gradient boosting), knn (k-Nearest Neighbors), mda (mixture discriminant analysis), pls (partial least squares), rf (random forest), svmPoly (support vector machines with polynomial kernel), and svmRadial (support vector machines with radial basis function kernel). The performance of the classification algorithms was assessed with error matrices, indicating the overall accuracy (i.e., the percentage of correctly classified pixels), user’s accuracy, producer’s accuracy, and kappa index. The Kappa index is a statistical measure of observed agreement that takes into account the effect of chance [ 34 ]. We compared the results obtained using the different classification algorithms mentioned earlier. To evaluate the enhancement of the classification accuracy as additional information was integrated into the classifier, we compared the performance of each classifier using identical training areas as input data. We used images only from summer months, summer months plus winter months, and summer plus winter months with spectral indices. Only indices that were not correlated, meaning they did not provide redundant information, were finally included in the classification procedure. To evaluate differences in the accuracy of the different algorithms, we created boxplots with the overall accuracy data of all algorithms for each year and map. This was done to visualise the data and make comparisons between the information provided by each algorithm in the datasets. 2.2.3. Change Analysis For the analysis of changes in the identified thematic classes, the map resulting from the classification process with the best overall accuracy and Kappa index value was selected. Once the maps were chosen, the area occupied by each habitat type was estimated for each year. To further understand the primary land-cover conversions occurring between 1997 and 2021, we constructed a transitional matrix. 3. Results 3.1. Accuracy Assessment The results show a high correlation between several of the multispectral indices, namely EVI, SAVI, and NDVI, and NDWI2 and MNDWI (Figure 2). The variables least correlated with the other indices were NDWI2 and MSAVI. These indices were included in the classification process along with the spectral bands, allowing the incorporation of information related to ‘greenness’ and humidity, proxies for primary productivity and water balance. The results do not show significant differences in the overall accuracy of the resulting maps derived from the different machine learning classifiers. The medians have very similar values, ranging from 0.8 to 0.95 in 1997 and from 0.7 to 0.9 in 2022 (Figure 3), indicating a consistently high overall classification capability of the identified land cover types. The image classification capability varies depending on the incorporated spectral information. In 1997, the incorporation of bands from two different periods decreases the Land 2024,13, 1601 7 of 13 predictive capability, while in 2022, there is less variability around the median, indicating some consistency among classifiers. Land 2024, 13, x FOR PEER REVIEW 7 of 13 Figure2.Correlation diagram between multispectral indices. Top right: March 1997, top left: August 1997, bottom right: March 2022, and bottom left: July 2022. The results do not show significant differences in the overall accuracy of the resulting maps derived from the different machine learning classifiers. The medians have very similar values, ranging from 0.8 to 0.95 in 1997 and from 0.7 to 0.9 in 2022 (Figure 3), indicating a consistently high overall classification capability of the identified land cover types. The image classification capability varies depending on the incorporated spectral information. In 1997, the incorporation of bands from two different periods decreases the predictive capability, while in 2022, there is less variability around the median, indicating some consistency among classifiers. Figure 2. Correlation diagram between multispectral indices. Top right: March 1997, top left: August 1997, bottom right: March 2022, and bottom left: July 2022. Land 2024, 13, x FOR PEER REVIEW 7 of 13 Figure2.Correlation diagram between multispectral indices. Top right: March 1997, top left: August 1997, bottom right: March 2022, and bottom left: July 2022. The results do not show significant differences in the overall accuracy of the resulting maps derived from the different machine learning classifiers. The medians have very similar values, ranging from 0.8 to 0.95 in 1997 and from 0.7 to 0.9 in 2022 (Figure 3), indicating a consistently high overall classification capability of the identified land cover types. The image classification capability varies depending on the incorporated spectral information. In 1997, the incorporation of bands from two different periods decreases the predictive capability, while in 2022, there is less variability around the median, indicating some consistency among classifiers. Figure 3. Boxplot of the overall accuracy data of the algorithms for the summer maps (Summer), summer and winter maps (Summer + Winter), and summer, winter, and indices maps (Summer + Winter + Index) in the years 1997 and 2022. The line dividing the box represents the median, the box represents 50% of the data, the lines represent 25% of the data, and the points are outliers. Land 2024,13, 1601 8 of 13 3.2. Change Analysis Changes in the Fragas de Eume Natural Park were analysed using the maps with the highest overall accuracy in 1997 and 2022, in this case, the ‘pls’ algorithm (overall accuracy of 0.87 and 0.98, respectively). From these data, we can observe significant changes in the extent of some areas over the period 1997–2022. We found that deciduous forest decreased (17.6%) from 1997 to 2022, while shrubland (18.2%), cropland (5.9%), and evergreen forest (i.e., eucalyptus plantations) increased (48.2%) (Figures 4–6). Land 2024, 13, x FOR PEER REVIEW 8 of 13 Figure3.Boxplot of the overall accuracy data of the algorithms for the summer maps (Summer), summer and winter maps (Summer + Winter), and summer, winter, and indices maps (Summer + Winter + Index) in the years 1997 and 2022. The line dividing the box represents the median, the box represents 50% of the data, the lines represent 25% of the data, and the points are outliers. 3.2.ChangeAnalysis Changes in the Fragas de Eume Natural Park were analysed using the maps with the highest overall accuracy in 1997 and 2022, in this case, the ‘pls’ algorithm (overall accuracy of 0.87 and 0.98, respectively). From these data, we can observe significant changes in the extent of some areas over the period 1997–2022. We found that deciduous forest decreased (17.6%) from 1997 to 2022, while shrubland (18.2%), cropland (5.9%), and evergreen forest (i.e., eucalyptus plantations) increased (48.2%) (Figures 4–6). The primary land-cover transitions were from deciduous forest to evergreen forest and shrubland, and from evergreen forest to shrubland (Table 4). Table4.Number of pixels and hectares for each class per year. Habitat19972022 Hectares Water 323 337 Deciduous forest 4379 3607 Shrubland 1295 1584 Croplands 2399 2540 Evergreen forest 682 1011 Figure4.This Sankey diagram illustrates the transitions between different land use types over time. The categories include croplands (yellow), deciduous forest (red), evergreen forest (green), shrubland (light green), and water (blue). The flow lines between the source (left) and target (right) nodes represent the changes in land use, measured in hectares, highlighting the dynamic nature of land cover transitions within the study area. The thickness of each line corresponds to the magnitude of change, facilitating the visualisation of how land use categories have evolved. Figure 4. This Sankey diagram illustrates the transitions between different land use types over time. The categories include croplands (yellow), deciduous forest (red), evergreen forest (green), shrubland (light green), and water (blue). The flow lines between the source (left) and target (right) nodes represent the changes in land use, measured in hectares, highlighting the dynamic nature of land cover transitions within the study area. The thickness of each line corresponds to the magnitude of change, facilitating the visualisation of how land use categories have evolved. Land 2024, 13, x FOR PEER REVIEW 9 of 13 Figure5.Bar plot of land data by land type and year. Figure6.Maps of land cover habitat classes for 1997 (top) and 2022 (bottom) in the ‘Fragas do Eume’ Natural Park. Coordinates are in the UTM coordinate system, WGS84 Zone 29N.  Figure 5. Bar plot of land data by land type and year. Land 2024,13, 1601 9 of 13 Land 2024, 13, x FOR PEER REVIEW 9 of 13 Figure5.Bar plot of land data by land type and year. Figure6.Maps of land cover habitat classes for 1997 (top) and 2022 (bottom) in the ‘Fragas do Eume’ Natural Park. Coordinates are in the UTM coordinate system, WGS84 Zone 29N.  Figure 6. Maps of land cover habitat classes for 1997 (top) and 2022 (bottom) in the ‘Fragas do Eume’ Natural Park. Coordinates are in the UTM coordinate system, WGS84 Zone 29N. The primary land-cover transitions were from deciduous forest to evergreen forest and shrubland, and from evergreen forest to shrubland (Table 4). Table 4. Number of pixels and hectares for each class per year. Habitat 1997 2022 Hectares Water 323 337 Deciduous forest 4379 3607 Shrubland 1295 1584 Croplands 2399 2540 Evergreen forest 682 1011 4. Discussion The resulting maps derived from a supervised classification procedure indicate significant changes in the main land-cover types identified in the ‘Fragas do Eume’ Natural Park. There has been an increase in the area occupied by Eucalyptus plantations, a species that