Study on Factors Influencing Forest Distribution in Barcelona Metropolitan Region

1. Introduction

As one of the precious natural resources on the earth, forests have made significant contributions in protecting the ecological environment, balancing climate change, and coordinating urban development [1]. It plays an indispensable role in regulating global carbon balance, promoting energy cycle and maintaining ecosystem stability [2]. At the same time, forests are the main body of terrestrial ecosystems and play a decisive role in regional and global ecological services [3]. It is not only an important material basis for the sustainable development of a region or country, but also an indispensable renewable resource in economic construction and ecological environment construction [4]. In addition, forests are important habitats for biodiversity and help maintain the balance and stability of ecosystems. However, the forest system is also highly sensitive to changes in the external environment, and has the characteristics of adapting to the evolution of the natural environment and reflecting the intensity of human activities [5]. With the rapid development of economy and science and technology, climate change and ecological impacts caused by human activities have attracted more and more attention. Global environmental change and sustainable development are the main challenges facing mankind in modern times. The Sixth Assessment Report of the United Nations Intergovernmental Panel on Climate Change (IPCC) further confirmed the objective facts of global climate warming in the past century, and the frequency and intensity of extreme heat. And the intensity and duration of heat waves are increasing, and even if global warming stabilizes at 1.5°C, it will increase further in the future [6]. The signals of the impact of human activities on global climate warming are becoming increasingly clear. The uncertainty of climate change dynamics and the impact of greenhouse gases produced by human activities on the global ecological environment are unprecedented [7]. The impact of climate change on forest ecosystems is multifaceted. At the same time, the continuous expansion of urban construction, industrial land and infrastructure has also led to the occupation and destruction of a large amount of forest land, damaging the ecosystem, deteriorating the ecological environment, and ultimately forming a vicious cycle. A good understanding of climate change and the impacts of factors on forest ecosystems can improve resilience to a variety of potential adverse impacts. Therefore, increasing the protection of forest ecosystems and conducting long-term monitoring and evaluation of forest ecosystems are the guarantee for achieving harmonious and rapid development of economy, society and the environment, and are urgent issues to be solved in the construction of ecological civilization. In-depth study of the long-term dynamics of different forest types and their driving factors can effectively promote the understanding of the evolution of spatiotemporal patterns of forest resources and their dynamic response mechanisms, and is of great significance to the research and formulation of forest resource protection and rational utilization strategies.

With the development of remote sensing technology, traditional manual methods of surveying forest resources have been replaced, which have high precision and accuracy but are high cost, low efficiency and cannot be implemented on a large scale [9]. Monitoring the dynamic change characteristics of vegetation at different regional scales based on long-term series of vegetation indexes has become an important topic in the field of global change research [10]. The application of remote sensing technology in forest monitoring and assessment was initially based on coarse-resolution optical remote sensing images. As early as 1995, Achard et al. used optical remote sensing data from NOAA/AVHRR to generate a vegetation index, and used a supervised classification method to classify tropical rainforests in Southeast Asia [11]. Although forest categories cannot be identified in detail, it demonstrates the potential of remote sensing data in forest cover mapping in Southeast Asia, and becomes the bud for human use of remote sensing to conduct forest resource investigation and research. In 2006, Liu Aixia and others used principal component analysis (PCA) and neural network classification (NNC) methods based on AVHRR remote sensing data from 2000 to 2001 to effectively distinguish forests and non-forests in China [12]. This provides a beneficial reference for forest map drawing in other areas. Later, Jing et al. used elevation data, soil index and optical NDVI index to classify the forest in the study area into coniferous forest, broadleaf forest and shrub forest using a new forest vegetation classification method based on multi-temporal remote sensing, which greatly improved the classification accuracy [13]. With the advancement of remote sensing impact analysis methods and the gradual maturity of satellite sensing technology, high-resolution satellite data such as Landsat and Sentinel are widely used in land cover research. Datasets such as CORINE Land Cover [14] and Global Land Cover_FCS30 [15] provide long-term, high-resolution land use classifications in different countries, regions and even the world, which include more detailed distinctions between forest areas.

Currently, the most commonly used vegetation index to study vegetation coverage, phenological changes, and vegetation dynamics in time series is the Normalized Difference Vegetation Index (NDVI) [16]. It is considered to be an effective indicator that reflects the vegetation growth status and coverage, and can reflect the regional ecological environment quality and vegetation coverage changes to a certain extent, so it is widely used in the study of vegetation dynamic changes [17]. To a large extent, scholars believe that there is a positive correlation between land vegetation cover and changes in NDVI. Some scholars have studied that in the past 30 years, China's vegetation NDVI has also shown an overall improvement trend [18], and vegetation has generally shown a gradual increase trend, with a higher proportion of green space represented by forest land [19]. In addition, Wang Xiaoxia and others subdivided specific forest types and found that broad-leaved forest had the fastest improvement rate in NDVI, followed by mixed forest and coniferous forest [20]. In 2022, Zhang Xu et al. analyzed the ecological environment quality of different land use types in Europe's metropolitan areas and found that the NDVI in forest areas ranked highest among all land use classifications [21]. They then conducted another experiment in the Beijing-Tianjin-Hebei region of China, and the results were the same [22]. In 2023, Tian Rui combined meteorological factors and changes in vegetation NDVI to analyze the spatiotemporal changes in forest ecosystem coverage in the Parlung Zangbo Basin and Zayu River Basin from 1971 to 2020, revealing the forest evolution rules and main sensitive factors in the study area [23].

In the context of climate change, the structure, function and productivity of forest ecosystems face serious challenges [23]. Initially, Smith et al. used the Holdridge model to predict that the distribution of forest types would change due to global climate change, with the area of warm temperate and subtropical forests decreasing [24]. Neilson et al. also found that forest cover shifted significantly with global climate change [25]. In 2003, Cao Heqin and others used camera capture data for seven consecutive years and found that the start of the forest leaf phenology growing season was advanced and the end of the growing season was delayed, thus extending the length of the growing season [26]. In addition, there are other natural factors and man-made factors that interact with forest distribution. Based on the NDVI data from MODIS, Hu Junde et al.'s research showed that there is a high correlation between vegetation coverage changes and precipitation in the Ordos region, and their dynamic change patterns are highly consistent [27]. The contradiction between limited urban land space resources and the rapid spread of urban construction land in recent years has caused urban expansion to occupy a large amount of surrounding forests and agricultural land [28, 29]. In 2022, Arellano et al. pointed out that the obvious lack of green vegetation coverage in urban Barcelona makes it less resistant to global warming and heat waves [30].

There are many factors that affect forest distribution changes. In previous studies, there was a lack of comprehensive assessment that considered the impact of climate change, human activities, green plant distribution and other natural geographical factors on long-term forest distribution changes. In addition, affected by differences in ecological functions and physiological growth characteristics, different vegetation types have different responses and adaptation levels to influencing factors. Therefore, in-depth multi-faceted research on the dynamics of different vegetation types and their driving factors can effectively promote the understanding of the evolution of vegetation spatiotemporal patterns and their dynamic response mechanisms. From this perspective, this article will analyze the changes affecting forest land distribution and morphology in Barcelona Metropolitan Region based on the more accurate CORINE Land Cover (CLC) data set, and find the important factors affecting this change. The most important thing is to find the interaction between forest land changes and influencing factors.

2. Materials and Methods

2.1. The Field of Study

Barcelona Metropolitan Region (BMR, Figure 1) is located in the northeast of the Iberian Peninsula, in the center of the Mediterranean corridor connecting Spain to the rest of continental Europe, and has a typical Mediterranean climate. It is the largest metropolitan area in Catalonia and serves as its political, economic and cultural core. It covers an area of approximately 3,224.7km² and contains 164 municipalities. With a population of around 4.7 million inhabitants, the BMR is the most densely populated metropolitan area in the European Union. The Llobregat and Besòs rivers are the two main rivers that flow through the metropolitan region of Barcelona. The coastal and pre-coastal mountain ranges delimit the coastal and pre-coastal depressions, where the main population centers are located. The urban core area of Barcelona is about 100km², with a population of more than 1.65 million, and a population density of more than 16,500 unhabitants/km². The Collserola “serralada” (the central part of the coastal mountain range) occupies the central place of the metropolitan area, making it the "green lung" around which the built spaces are located. The Barcelona Metropolitan Territorial Plan (2010) pays special attention to the protection of the ecological environment and green spaces of the BMR. Its forest area is about 1,380 km², which represents approximately 43% of the total urban area. Based on the above, we selected the Barcelona Metropolitan Region as a research area to analyze the factors and mechanisms that affect forest distribution.

2.2. Methodology

CORINE Land Cover (CLC) provides us with relatively accurate 100-meter resolution European land cover data for 2006, 2012 and 2018, in which forest land is divided into three categories according to different growth types. This provides the possibility for us to study the long-term dynamic evolution of forest cover. Within the studyable time range, that is, between 2006 and 2018, we suppose that the forest area in the Barcelona Metropolitan Region (BMR) has decreased, and the degree of fragmentation has become more and more serious, but the change process may not be single. At the same time, there is a certain interaction between the green space index NDVI of the BMR forest area and the area and form of forest cover, and the mechanism of this interaction is complex.

In order to achieve the research purpose of this paper and verify the above hypothesis, this research will be completed through the following steps.

• First of all, we need to determine the evolution of BMR forest land distribution proportions. In addition to presenting the area change process of the overall forest land and the internal classification of the forest, we will analyze it from the following seven indicators.

(a) Landscape proportion (PLAND) reflects the proportion of various land types in the total area, with the largest area being the main landscape.

$$\mathrm{P}\mathrm{L}\mathrm{A}\mathrm{N}\mathrm{D}=\frac{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}}{\mathrm{A}} *100\%,$$

(1)

Where：

a_ij represents the area of the j-th patch in the i-th landscape type;

A is the total area of the landscape.

When the patch area percentage value is close to 0, it indicates that there are few patch types in the landscape; when the ratio is equal to 1, it indicates that the entire landscape is only composed of type 1 patches.

(b) Patch density (PD) can reflect the overall heterogeneity and fragmentation of the landscape as well as the degree of fragmentation of a certain type, and reflects the heterogeneity of the landscape unit area.

$$\mathrm{P}\mathrm{D}=\frac{{\mathrm{n}}_{\mathrm{i}}}{\mathrm{A}} *100\%,$$

(2)

Where：

n_i is the number of patches in type i landscape;

A is the total area of all landscapes.

(c) Maximum patch index (LPI) is used to determine the dominant patch type in the landscape.

$$\mathrm{L}\mathrm{P}\mathrm{I}=\frac{a_{max}}{\mathrm{A}} *100\%,$$

(3)

Where：

a_max refers to the area of the largest patch in the landscape or a certain patch type;

A is the total area of the landscape.

The size of this index value can help determine the dominant patch type in the landscape, indirectly reflecting the direction and size of interference from human activities.

(d) Compactness reflects the compactness of the landscape in the area.

$$\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}=\frac{P_i}{{\mathrm{A}}_{\mathrm{i}}},$$

(4)

Where：

P_i is the perimeter of the i-th landscape;

A_i is the area of type i landscape.

The smaller the index, the more compact the landscape is, and the larger the index, the more dispersed the landscape is.

(e) Shannon entropy (ENT) is a measurement index based on information theory that can reflect landscape heterogeneity and structural complexity. It is particularly sensitive to the uneven distribution of patch types in the landscape, that is, it emphasizes the contribution of rare patch types to information, which is also different from other diversity indices.

$$\mathrm{E}\mathrm{N}\mathrm{T}=-{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{p}\mathrm{i}\right)\mathrm{*}\left({{\mathrm{l}\mathrm{o}\mathrm{g}}_2}^{\mathrm{p}\mathrm{i}}\right),$$

(5)

Where：

n refers to the total number of patch types in the landscape,

Pi refers to the area ratio of patch type i to the entire landscape.

The larger the index, the greater the degree of landscape fragmentation in the area, the more complex the structure, and the higher the negativity.

(f) Perimeter area fractal dimension (PAFRAC) refers to the non-integer dimension of the irregular geometric shape of the landscape, reflecting the complexity of the landscape shape. The value is between 1 and 2.

$$\mathrm{P}\mathrm{A}\mathrm{F}\mathrm{R}\mathrm{A}\mathrm{C}=\frac{\frac{2}{[{\mathrm{n}}_{\mathrm{i}\mathrm{j}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}\mathrm{l}\mathrm{n}{\mathrm{p}}_{\mathrm{i}\mathrm{j}}-\mathrm{l}\mathrm{n}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}]-\left[\right({\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}\mathrm{j}})\left({\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}})\right]}}{({\mathrm{n}}_{\mathrm{i}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{l}\mathrm{n}{\mathrm{p}}^2}_{\mathrm{i}\mathrm{j}})-{({\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{l}\mathrm{n}\mathrm{p}}_{\mathrm{i}\mathrm{j}})}^2},$$

(6)

Where：

a_ij refers to the area of the j-th patch in the i-th type of landscape;

p_ij represents the perimeter of the j-th patch in the i-th type of landscape;

n_i is the number of patches.

The closer the calculation result is to 1, the more regular the shape of the patch, or the simpler the patch, and the interference factor is considered to be large; conversely, the closer the calculation result is to 2, the more complex the shape of the patch is, and the interference factor is considered to be small.

(g) Cohesion reflects the aggregation and dispersion status of patches in the landscape, with values ranging from -1 to 1.

$$\mathrm{P}\mathrm{A}\mathrm{F}\mathrm{R}\mathrm{A}\mathrm{C}=[1-\frac{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}\mathrm{j}}}{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}\mathrm{j}}\sqrt{{\mathrm{a}}_{\mathrm{i}\mathrm{j}}}}]{[1-\frac{1}{\sqrt{\mathrm{A}}}]}^{-1}\mathrm{*}100\mathrm{\%},$$

(7)

Where：

a_ij refers to the area of the j-th patch in the i-th type of landscape;

p_ij represents the perimeter of the j-th patch in the i-th type of landscape;

A is the total area of the landscape.

When the index result is -1, the patches are completely dispersed, when the index result is 0, they are randomly distributed, and when the index result is 1, they are aggregated.

For the seven forest landscape indicators mentioned above, we should not only analyze the general indicators of the forest land distribution pattern of the BMR, but it is also important to evaluate each forest class analyzed.

2.

Secondly, in order to intuitively and accurately reveal the long-term change characteristics of forest NDVI, the change of each pixel is calculated to reflect the increase or degradation trends of NDVI over time.

In this study, not only the overall change trend of winter and summer NDVI in forest areas from 2006 to 2018 was analyzed, but changing trends within each time period will also be calculated to clarify the changing process of NDVI in the Barcelona forest area and its relationship with the evolution of forest morphology. NDVI data can be obtained from MODIS with a resolution of 250m.

3.

Next, we will try to use NDVI to evaluate the forest green environmental quality of BMR. The average NDVI of each type of forest land was extracted to establish a normalized green quality assessment system. The average NDVI of each type of forest land is extracted, and for each year, we establish a normalized green quality assessment system through formula (8). Set the minimum value to 0 and the maximum value to 1 as the assigned value (E_i) of the green environmental quality proportion of each type of forest species.

$${\mathrm{E}}_{\mathrm{i}}=\frac{{\mathrm{X}}_0-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{\mathrm{X}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}},$$

(8)

Where:

E_i is the dimensionless green quality weight index;

X₀ is the average NDV of each type of forest land in that year;

X_max is the maximum value of NDV of various types of forest land in that year;

X_min is the minimum value of NDV of various types of forest land in that year.

It is then evaluated according to the following formula (formula 9).

$${\mathrm{G}\mathrm{Q}\mathrm{I}}_{\mathrm{t}}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{A}}_{\mathrm{i}}\mathrm{*}\frac{{\mathrm{E}}_{\mathrm{i}}}{\mathrm{T}\mathrm{A}}$$

(9)

Where:

GQI_t is the green environmental quality index in year t;

i is the forest type;

A_i is the area of that type of land use;

E_i is the weight of the ecological quality index of this type of land use;

and TA is the total area of the forest.

4.

We need to study the relationship between forest area distribution and various possible influencing factors. We will establish a 1km grid within the metropolitan area and extract the annual proportion of forest area within the grid and the average value of winter and summer NDVI. In addition, collect various types of data listed in the table. MODIS can provide land surface temperature (LST) and normalized building index NDBI. The urban heat island (UHI) effect will be represented by the urban-rural temperature difference, that is, the difference in average surface temperature between built-up areas and rural areas is expressed based on land cover data. For each year's daytime and nighttime LST, we subtract their average values ​​in the rural land to get the approximate intensity distribution of UHI. The larger the value, the higher the UHI intensity. DEM terrain data comes from SRTM with a resolution of 30m. The impervious ground data comes from GlobeLand 30, also with a resolution of 30 meters. E-OBS can provide annual and monthly European precipitation raster data, but the resolution is 1° and the scale is very large. Therefore, we used the Kriging interpolation method to re-establish the BMR precipitation map with a resolution of 1 kilometer based on the E-OBS precipitation data. After obtaining the above data, we used the proportion of forest area in each grid as the dependent variable to establish an OLS model to analyze their importance.

Table 1. List of factors potentially influencing forest spatial distribution.

Type	Factors
Natural factors	Loçngitude
	Latitude
	Distance from coastline
	Orientation
	Altitude
	NDVI
	Precipitation
	LST
Human activity	NDBI
	Urban heat island effect
	Impermeable area
	Artificial area

Table 1. List of factors potentially influencing forest spatial distribution.

Type	Factors
Natural factors	Loçngitude
	Latitude
	Distance from coastline
	Orientation
	Altitude
	NDVI
	Precipitation
	LST
Human activity	NDBI
	Urban heat island effect
	Impermeable area
	Artificial area

5.

Once it is clear that there is an obvious interaction between NDVI and forests, we must analyze the climate factors that affect NDVI to discover the potential threats that climate change may have on forests and predict possible trends in changes in BMR forest layout caused by climate change. Based on the data obtained from the grid in step 4, we used average NDVI as the dependent variable, annual precipitation and daytime and nighttime LST as independent variables to establish three OLS regression models to analyze the correlation between them and the average NDVI. And analyze the average NDVI and the distribution pattern of the found most significant independent variables on the map.

6.

It is also very important to study the impact of various factors on forest landscape indicators. Taking 2018 as an example, we calculated the landscape indicators in each grid and used them as independent variables to establish models with the factors involved in the fourth step to analyze their relationships. In order to build models, it is necessary to cut the BMR into 1km grids and the landscape index inside each grid is calculated separately. Because of that, the results calculated by PAFRAC, COHESION and some other indices seriously deviate from the actual situation, and such analysis is meaningless. In addition, the landscape proportion (PLAND) index has been analyzed as a dependent variable in step 4. Therefore, in this section we only analyze patch density (PD), compactness and Shannon entropy (ENT). In addition, we will also consider the impact of temperature changes from 2006 to 2018 on forest pattern morphology.

7.

Finally, ArcGIS was used to analyze the transformation process within the forest areas of the Barcelona metropolitan area between 2006 and 2018, and find the types of forest land that were being lost and the types that were growing significantly over the years. In addition, land transfer in the entire metropolitan area also needs to be analyzed to find the main degradation directions of forest areas.

3. Results

3.1. Analysis of Forest-Related Landscape Indicators

3.1.1. The Overall Scale of the Forest Landscape Has Been Slightly Reduced

Analyzing forest land area and occupancy ratio can intuitively reflect the importance of forests in BMR landscapes as well as the internal structure and evolution of forests. Overall, the total area and proportion of forest land decreased between 2006 and 2018, but not significantly (Figure 2). Its area decreased from 1,384km² to 1,354km², a proportion of about 1%. Among forests, coniferous forests occupy a dominant position, with an area almost twice that of broadleaf forests, but the areas of both forests are decreasing. The smallest land use is mixed forest, which accounts for only about 1.6% of forest land, but it is the only expanding forest species.

From Figure 3, we can clearly see the distribution structure of BMR forest in these three years. Overall, the forest area is widely distributed and can almost cover the study area. They can appear in large areas of contiguous distribution in the northeast, while showing a high degree of fragmentation in the southwest. The south-central area is where the central city of Barcelona is located, and along much of the coastline there is a dense distribution of forest that is difficult to see.

3.1.2. Increased Fragmentation of Forest Landscapes

Figure 4 shows the growth process of the number of patches and patch density in the BMR forest landscape during the study period. Overall, over 12 years, the number of forest land fragments increased from 487 to 494, and patch density also increased. From the perspective of classification within the forest, the total number of patches of various woodlands and the density of forest land occupies are also growing. Coniferous forest is still the most outstanding forest species, contributing nearly two-thirds of the fragments to BMR's forest landscape. The most inconspicuous one is still the mixed forest.

3.1.3. Forests Are the Most Widespread Type of Land Use

Analyzing the maximum patch index can determine the dominant patch types in the landscape and determine the importance of different types of land in an area. To identify the landscape type with the largest area, we must first classify RMB according to the land use situation displayed by CLC. As a result, there are 11 categories in total: 1- continuous built-up area; 2- discontinuous built-up area; 3- industrial land; 4- transportation land; 5- mine, dump and construction sites; 6- leisure land; 7- cropland; 8- woodland; 9- grassland; 10- barren land; 11- water bodies. Descriptions of various land use classifications can be reviewed in the Appendix A.

After the classification was completed, we calculated the total area of each land type and found that forest has the largest area among all types of land in BMR. In other words, forest land is the most advantageous land type in BMR. And based on previous results, we already know that coniferous forest is the most important forest species in the forest. The results of the maximum patch index of BMR and forests are presented in Figure 5. Obviously, more than 62% of BMR forests are coniferous forests. Forest occupies more than 40% of the BMR's land, but it is slowly being invaded and lost year y year, even if the reduction is not rapid.

3.1.4. The Structure of Forest Land Is Becoming More Complex and the Degree of Human Interference Is Deepening

The Shannon entropy index reflects the heterogeneity and structural complexity of the landscape, while the perimeter area fractal dimension can explain the complexity of the landscape shape and human interference. From Figure 6, we can see that the overall Shannon entropy of BMR has increased, from 2.2 to 2.25, which shows that the land patches of RBMR are more scattered, the distribution is more uneven, and the structure is more complex. However, it experienced a small decrease between 2012 and 2018. The complexity of forest land has also been increasing, and the index increase has been higher than that of the BMR as a whole. However, the forest's perimeter area fractal dimension is decreasing and experienced its lowest point in 2012. Judging from the shape of the forest patches, they have become more regular than before. They may have experienced interference from human factors, but this process is very small. The perimeter area fractal dimension of the mixed forest is the largest among the three forest classification species, its shape distribution is the most irregular, and this index is larger than at the beginning. Broadleaf forest has also experienced growth, but it is the forest species with the lowest index, very close to 1, and it is the most disturbed. The coniferous forest with the largest area has been increasingly affected by humans, and its index has declined the fastest among forests.

3.1.5. Forest Land Is Increasingly Dispersed

According to the compactness and cohesion index, the aggregation of landscape land distribution can be analyzed. Figure 7 shows us the evolution results of the two indices. During this process, BMR forest land became increasingly fragmented, as evidenced by changes in both indices. The change process of broad-leaved forest is the most obvious, and it is also the most similar to the overall change pattern of the forest. The coniferous forest has not changed much, but it has gradually developed in a more dispersed direction. The exception was the mixed forest, which was more compactly distributed over 12 years. Judging from the aggregation index, although it is lower than the initial value, it has experienced an increase between 2012 and 2018, that is, it is experiencing aggregation distribution. However, it must be pointed out that the mixed forest is the most dispersed among the three types of forest, and its degree of aggregation is obviously much lower than the others.

3.2. NDVI Change Trend Analysis

3.2.1. The Development of NDVI in BMR and Forest Areas Is Relatively Optimistic

First, we compared each pixel of NDVI in 2006 and 2018 from the overall and forest local level, and found that the changes in NDVI during the entire study period were relatively optimistic (Figure 8). On the whole, the NDVI in most areas of the BMR is improving, especially in the northern and central areas, where the improvement is relatively large. However, some areas have experienced deterioration, which is more serious in the northeast and southwest. Figure (b) illustrates that the development of NDVI in the forest area is very good, and almost all parts have been improved, except for a small area in the northeast that is worthy of concern.

3.2.2. The Development Status of NDVI in 2012 Was Disappointing

Figure 9 shows us the development of NDVI in BMR and forest areas in two periods. The most obvious thing is that the performance in the first period is worse than that in the second period. Between 2006 and 2012, although the NDVI in many places in the BMR was improving, areas with greater deterioration were fragmented and scattered in almost every corner of the study area. Moreover, the distribution of improved areas is not concentrated enough. NDVI development in forest areas is also unsatisfactory, with deteriorating and unchanged marks occupying most of the area.

But what is surprising is the substantial improvement in NDVI between 2012 and 2018. Overall, the trend of improvement is very obvious and is concentrated in the northern region. Degraded and unchanged markers are distributed alternately, but there are fewer instances of degradation and the distribution is more dispersed. There is concentrated deterioration only in the northern end. The changes in NDVI in the southern region were significantly better than those in the previous period. For forest areas, the NDVI is also improving in a good direction, with almost no degradation, even if a small area of deterioration marks is concentrated in the northern end.

3.2.3. The Forest's NDVI Is the Most Outstanding

In addition, we also extracted and counted the annual changes in the average NDVI of each land type in BMR, and the results are presented in Figure 10(a). We clearly see that among all land uses, the highest annual NDVI contribution is forest, with a minimum of 0.62, followed by grassland. The smallest NDVI is in the continuous built-up area, but it has continued to rise in three years. The NDVI of industrial land and its development pattern are very similar to those of continuous built-up areas. All land uses have experienced a substantial increase in NDVI in the second stage, that is, from 2012 to 2018. As can be seen from Figure (b), the average NDVI of the forest remained almost unchanged from 2006 to 2012, but rapidly increased to 0.67 from 2012 to 2018.

3.3. Forest Green Environmental Quality Assessment

3.3.1. Broadleaf Forest Has the Highest Assessment Weight

The NDVI weight evaluation index system that can be applied to changes in land use types in different periods and regions can more scientifically and reasonably evaluate the green environmental quality of forest areas [23]. Table 2 shows the weight index used by different forest lands to assess environmental quality. Coniferous forest, which accounts for the largest area, is not assigned the largest weight. Broad-leaved forest has become the most important land occupation type in the environmental quality assessment process. Mixed forests have the smallest weight numbers.

3.3.2. The Quality of Forest Green Environment Is Increasing Year by Year

We show the evaluation results of the green environmental quality of the forest area in Figure 11. Satisfactorily, the green quality of the Barcelona forest is gradually improving. The environmental quality of forests is relatively high. In 2006, the index was 0.627. After two increases, it became 0.635. And the growth rate is getting faster and faster every time.

3.4. Analysis of Factors Affecting Forest Distribution

In order to explore the impact of various natural and man-made factors on forest distribution, we established a grid with a side length of 1 km, and used the proportion of forest area in each grid as the dependent variable, longitude, latitude, distance from the sea, orientation, altitude, slope, average winter and summer NDVI, annual precipitation, daytime LST, nighttime LST, NDBI, daytime, nighttime urban heat island effects, impervious ground and artificial ground were used as independent variables to establish 3 OLS regression models (Table 3) for the data in 2006, 2012 and 2018, accompanied by 7 significant variables. The most significant impact on forest distribution was NDVI (+), but in the last year NDBI (-) became the most important influencing factor. This is followed by daytime urban heat island effect (-), elevation (-), nighttime LST (+), longitude (-), artificial surface (-) and precipitation (+). In 2006 and 2018, the average NDVI and NDBI could alone explain 85.9% and 83.8% of the changes in forest distribution area, but in 2012 they only accounted for 60.2% of the explanatory power. Surprisingly, the artificially created daytime urban heat island effect has a stronger impact on forest growth than natural precipitation factors. If the capabilities of NDVI, NDBI and daytime heat island effect are taken into account, they have the ability to explain more than 88% of the patterns of forest distribution changes. The reason why daytime LST and nighttime urban heat island effect are not included in the best model is that they are highly collinear with the daytime urban heat island effect and nighttime LST variables.

3.5. Analysis of Climate Factors Affecting NDVI

Once it is confirmed that average NDVI is the most significant external factor regulating forest distribution, we must analyze the climatic factors that can affect the development of NDVI to predict possible trends in changes in BMR forest layout caused by climate change. For this purpose, we again established three OLS models, in which the dependent variable is NDVI, the annual average NDVI, and the independent variables are precipitation, daytime and nighttime LST. The results of the regression analysis parameters shown in Table 4 indicate that diurnal LST (-) is the climatic factor most strongly related to NDVI. In the 2006 model, it alone can explain 83% of the spatial variation of NDVI. And this proportion continues to rise. By 2018, it had 90% of the interpretation rights. Second is precipitation, which also has a certain influence on changes in NDVI, but it is very unstable. In 2006, its impact on NDVI was positive, but in 2012, it guided the development of NDVI in the opposite direction. However, in the last year, its influence was greatly weakened (student’s t=0.9), and the proportion of nighttime LST, which had almost no impact on the dependent variable in the first two years, rose to the second place (student’s t=4.27).

We present the changes in the average NDVI and daytime LST of the three years of BMR in Figure 12, which is consistent with the results predicted by the model, and the distribution patterns of the two are opposite and their red areas and light areas can almost overlap each other even though the LST image has a lower resolution. Areas with a higher NDVI index have lusher vegetation and lower daytime temperatures, while areas with low NDVI have relatively higher temperatures.

During the 12 years of study, the NDVI of BMR has hardly changed much, with low indices concentrated in non-forest areas and coastal areas. The red part of the surface temperature center gradually deepens in color and slightly expands toward the blue zone, which illustrates that the BMR is warming during the day.

3.6. Analysis of Influencing Factors of Forest Landscape Index

We took the three landscape indicators in 2018 as examples to analyze the effects and capabilities of various possible influencing factors on different forest landscape indices. In addition to the 13 variables considered when analyzing spatial changes in forest area proportions, we also added the difference between daytime and nighttime surface temperatures between 2006 and 2018 (LST_2018-LST_2006), which reflects the temperature changes in the BMR during these 12 years. After establishing three OLS regression models (Table 4) with various landscape indices as dependent variables, the most important results we found are consistent with the previous ones. NDVI is still the most meaningful independent variable for forest landscape changes. For these three indices, NDVI has a positive effect on them. Although the relationship between NDVI and compactness index is inverse, according to the definition of compactness index, the smaller the index, the greater the compactness of the landscape.

Table 5. Forest landscape indicator OLS model selected for analysis in 2018.

Independent variable d	Model_PD a					Model_ENT b					Model_Compactness c
	B	Beta	t	Sig.		B	Beta	t	Sig.		B	Beta	t	Sig.
Constant	-4549.12	-	-3.66	0.00		-0.08	-	-3.34	0.00		-372.52	-	-0.12	0.905
Longitude	0.00	-0.44	-3.71	0.00		0.00	-0.41	-3.38	0.00		0.00	0.02	0.07	0.943
Latitude	0.00	0.47	3.66	0.00		0.00	0.43	3.31	0.00		0.00	0.07	0.30	0.768
Distance from coastline	0.00	-0.20	-2.54	0.01		0.00	-0.18	-2.28	0.02		0.00	-0.09	-0.68	0.498
Orientation	0.00	0.00	-0.28	0.78		0.00	-0.01	-0.36	0.72		-0.01	0.00	-0.17	0.862
Altitude	-0.02	-0.10	-3.42	0.00		0.00	-0.10	-3.29	0.00		0.01	0.03	0.53	0.599
Slope	1.18	0.01	0.82	0.41		0.00	0.02	1.05	0.29		-6.09	-0.04	-1.68	0.093
NDVI_MEAN	91.58	0.34	8.20	0.00		0.00	0.36	8.47	0.00		-83.81	-0.22	-2.98	0.003
Precipitation	1.62	0.05	2.38	0.02		0.00	0.05	2.31	0.02		-0.87	-0.02	-0.51	0.611
LST_DAY	-3.63	-0.22	-7.38	0.00		0.00	-0.22	-7.23	0.00		2.85	0.12	2.31	0.021
LST_NIGHT	3.38	0.11	4.06	0.00		0.00	0.11	3.85	0.00		-2.59	-0.06	-1.24	0.216
NDBI	-70.13	-0.33	-7.87	0.00		0.00	-0.32	-7.44	0.00		3.73	0.01	0.17	0.868
Impermeable area	4.19	0.03	0.58	0.56		0.00	0.04	0.86	0.39		14.88	0.06	0.82	0.411
Artificial area	-20.75	-0.13	-2.98	0.00		0.00	-0.14	-3.12	0.00		-0.86	0.00	-0.05	0.961
Difference_LST_DAY_2018-2006	4.78	0.13	6.70	0.00		0.00	0.13	6.47	0.00		-3.82	-0.07	-2.12	0.034
Difference_LST_NIGHT_2018-2006	-2.23	-0.05	-2.15	0.03		0.00	-0.04	-2.00	0.05		5.47	0.08	2.10	0.036

^a The dependent variable of this model is patch density (PD). ^b The dependent variable of this model is Shannon entropy (ENT). ^c The dependent variable of this model is Compactness. ^d UHIE_DAY and UHIE_NIGHT became excluded independent variables during the regression analysis.

Table 5. Forest landscape indicator OLS model selected for analysis in 2018.

Independent variable d	Model_PD a					Model_ENT b					Model_Compactness c
	B	Beta	t	Sig.		B	Beta	t	Sig.		B	Beta	t	Sig.
Constant	-4549.12	-	-3.66	0.00		-0.08	-	-3.34	0.00		-372.52	-	-0.12	0.905
Longitude	0.00	-0.44	-3.71	0.00		0.00	-0.41	-3.38	0.00		0.00	0.02	0.07	0.943
Latitude	0.00	0.47	3.66	0.00		0.00	0.43	3.31	0.00		0.00	0.07	0.30	0.768
Distance from coastline	0.00	-0.20	-2.54	0.01		0.00	-0.18	-2.28	0.02		0.00	-0.09	-0.68	0.498
Orientation	0.00	0.00	-0.28	0.78		0.00	-0.01	-0.36	0.72		-0.01	0.00	-0.17	0.862
Altitude	-0.02	-0.10	-3.42	0.00		0.00	-0.10	-3.29	0.00		0.01	0.03	0.53	0.599
Slope	1.18	0.01	0.82	0.41		0.00	0.02	1.05	0.29		-6.09	-0.04	-1.68	0.093
NDVI_MEAN	91.58	0.34	8.20	0.00		0.00	0.36	8.47	0.00		-83.81	-0.22	-2.98	0.003
Precipitation	1.62	0.05	2.38	0.02		0.00	0.05	2.31	0.02		-0.87	-0.02	-0.51	0.611
LST_DAY	-3.63	-0.22	-7.38	0.00		0.00	-0.22	-7.23	0.00		2.85	0.12	2.31	0.021
LST_NIGHT	3.38	0.11	4.06	0.00		0.00	0.11	3.85	0.00		-2.59	-0.06	-1.24	0.216
NDBI	-70.13	-0.33	-7.87	0.00		0.00	-0.32	-7.44	0.00		3.73	0.01	0.17	0.868
Impermeable area	4.19	0.03	0.58	0.56		0.00	0.04	0.86	0.39		14.88	0.06	0.82	0.411
Artificial area	-20.75	-0.13	-2.98	0.00		0.00	-0.14	-3.12	0.00		-0.86	0.00	-0.05	0.961
Difference_LST_DAY_2018-2006	4.78	0.13	6.70	0.00		0.00	0.13	6.47	0.00		-3.82	-0.07	-2.12	0.034
Difference_LST_NIGHT_2018-2006	-2.23	-0.05	-2.15	0.03		0.00	-0.04	-2.00	0.05		5.47	0.08	2.10	0.036

^a The dependent variable of this model is patch density (PD). ^b The dependent variable of this model is Shannon entropy (ENT). ^c The dependent variable of this model is Compactness. ^d UHIE_DAY and UHIE_NIGHT became excluded independent variables during the regression analysis.

In addition to NDVI (+), NDBI (-) also has a very prominent correlation with PD and ENT, and their explanatory power together can reach more than 79% and 90% respectively. Secondly, the daytime surface temperature also has a negative impact that cannot be ignored on the spatial distribution of the two indices. Daytime temperature changes also have a greater ability to positively regulate these two dependent variables. These four independent variables explained about 84% and 95% of the distribution changes of these two variables in 2018 respectively.

For compactness, daytime LST(+) is the second most influential factor, followed by diurnal surface temperature changes. But temperature differences at different times work in different ways. Daytime temperature changes positively control the spatial distribution of forest compactness, while nighttime differences reversely guide its development. Together with NDVI, these four most significant variables have an impact on forest compactness of more than 76%.

3.7. Land Use Transfer Matrix Analysis

3.7.1. Mixed Forests Are the Only Expanding Woodlands

In order to find the changes in various cover types within the forest from 2006 to 2018, we first statistically calculated the ground cover transfer matrix of the forest in Table 6, and summarized the specific increases and decreases in the area of various types of forests in Figure 13. Approximately 1,350.56km² of forests in the BMR in 2006 and 2018 overlapped. Among them, the area of mixed forests and broadleaf forests increased slightly. This is because part of the coniferous forests converted to them. But from the perspective of BMR, the area of broadleaf forests has decreased much more than the area that has increased. Coniferous forest has shrunk the most among the three types of forest land, with a decrease of about 20km². Only the area of mixed forest increased, but the amplitude was very small, only 0.32km².

3.7.2. Forests Are Mainly Transferred Like Cultivated Land and Grassland

In addition to analyzing the transfer changes of different forests, we also compiled the transfer process of all BMR land in Table 7, and the change information is presented in Figure 14. In the previous analysis, we knew that the BMR forest range was shrinking. In fact, from the figure, we can see that its change range is the largest among all types of land cover in the BMR, with an area of approximately 28km² lost in 12 years. The lost forest was mainly used for agricultural production, about 13.4km², and the rest was turned into grassland, covering 10km2. In addition, discontinuous urban land also takes up part of the forest. Sadly, very little land has been converted to forest. Agricultural land, on the other hand, takes up a fair amount of forest, but overall it is noticeably smaller as industry, traffic, broken cities and grasslands encroach on it.

4. Discussion and Conclusion

In this study, we comprehensively demonstrated the spatial changes of forests in the Barcelona metropolitan region from 2006 to 2018 through various landscape indicators, and analyzed the external factors and modes of action that may affect forest changes. Finally, we summarize the transfer process between forests and other land uses. We found that the overall area of the forest, which is the most extensive among all land types in BMR, has decreased, but the degree of fragmentation has deepened and the internal structure has become more complex. Significant amounts of forest were converted to agricultural land and grassland. And forest areas provide the highest NDVI of all areas, and it is increasing. Finally, we found that NDVI is the variable that has the greatest impact on the spatial structure of the forest and promotes each other. Daytime surface temperature is a factor that interacts significantly in the opposite direction with the development of NDVI. In addition, the diurnal surface temperature changes during the study period have also become a condition that cannot be ignored. In general, the BMR forest area is decreasing and the morphology is becoming more and more fragmented. NDVI and daytime surface temperature are the two most important factors affecting its spatial changes.

Basically, the forest in the BMR has been protected satisfactorily because it has always been the largest land use type, even though it has shrunk slightly over the past 12 years. Thanks to the approval of the Barcelona Metropolitan Territorial Plan in 2010 (https://territori.gencat.cat/ca/01_departament/05_plans/01_planificacio_territorial/plans_territorials_nou/territorials_parcials/ptp_metropolita_de_barcelona/index.html#googtrans(ca|en)), there has been a significant protection of rural spaces, especially forests, in the BMR. However, Figure 3 and Figure 5 illustrate that the morphology of the forest is becoming more and more fragmented, and the degree of human disturbance is increasing, resulting in the forest becoming more dispersed and complex in structure. Although the scale is being reduced, the forest provides the most NDVI and contributes the most to the green ecological development of the metropolitan area, and its green quality index increases year after year. The NDVI has improved overall, but it is worth noting that the precipitation in 2018 was the largest during this period, so the changes in NDVI may be much more complicated than shown in our study. At the same time, Barcelona has experienced a well-known drought in recent years. If this drought trend is consolidated, the gradual reduction of precipitation may affect the green quality of vegetation in forest areas and even the entire BMR. What is surprising is that the coniferous forest with the largest area does not have the largest green quality evaluation weight like the broadleaf forest. This may be related to the distribution pattern. Compared with coniferous forests, broad-leaved forests have a more concentrated spatial form and therefore have more advantages in providing green quality. Among the three types of forest, the mixed forest is the most dispersed and has the smallest green weight, which can also be used as evidence to prove this conjecture.

We listed a total of 14 natural and man-made independent variables that may affect the spatial distribution of BMR forests. After modeling analysis, we found that NDVI is the most closely related positive correlation factor with forest distribution, which is consistent with the results of most scholars’ studies [18-23]. The climate factor most closely related to NDVI is daytime surface temperature. The higher the NDVI index, the lower the temperature, and it can be clearly seen from the reverse distribution map of the two that they are almost coincident. Therefore, daytime LST indirectly controls forest cover patterns. Secondly, the two human factors of NDBI and daytime urban heat island effect also have a negative impact that cannot be ignored on the development of forest systems. Regarding the structure of forest distribution, taking 2018 as an example, as time goes by, NDVI increases, patch density becomes higher, the distribution becomes more compact, and Shannon entropy also increases. In fact, we believe that the main reason for higher patch density and higher entropy is not the improvement of NDVI, but the time evolution. This may imply that the younger the time, the greater the forest area in areas with higher NDVI index, but the complexity is also increased, and therefore the degree of fragmentation is also greater. NDBI also has a strong inverse relationship with PD and ENT, but has little impact on compactness. Daytime LST also clearly controls them, with higher temperatures causing forest patch density and compactness to decrease and morphological structures to become simpler. In addition, temperature changes during these 12 years also have a certain impact on the distribution structure.

Forest land in BMR has been losing over the past 12 years. Coniferous forest is the forest species with the largest decrease, but mixed forest is an exception and is the only forest that has increased in area. A large amount of forest has been converted into agricultural land and grassland. The main direction of agricultural land reduction is industry, transportation, and fragmented built-up areas.

According to the results of this study, undoubtedly, there is a very significant interaction between daytime LST and the distribution of NDVI, which is most obviously related to forests. Several studies showed that between 1971 and 2022, the temperature in the main cities of the Spanish peninsula increased an average of 3.54°C, making it one of the most evident urban climate anomalies in the world [31,32]. Based on this, we can put forward the conjecture that in the upcoming period (2018-2024), NDVI may be degraded by this impact, and the spatial distribution of the forest will also be negatively affected. This is also one of the main objectives of our further work.

This paper comprehensively considers a comprehensive assessment of the impact of climate change, human activities, and other physical and geographical factors on long-term forest distribution changes, and discovers direct and indirect influencing factors that affect forest spatial distribution and structural changes. And based on forest types, the evolution pattern within the forest and the contribution to green ecological quality were analyzed. Finally, the transfer direction of land use types in the urban area was judged. This study provides a comprehensive analysis of the forest system from a global to internal perspective. However, due to technical limitations, the resolution of various types of remote sensing impact data is not uniform, and due to the particularity of the precipitation data, the resolution is too large, so we adopt the method of re-difference, which will affect the accuracy of the experimental results. In future research, we will try our best to use or build a more accurate database with a uniform resolution.