Green Space Assessment and Management in Biscay Province, Spain using Remote Sensing Technology

Our ecosystem, particularly forest lands, contains huge amounts of carbon storage in the world today. This study estimated the above ground biomass and carbon stock in the green space of Bilbao Spain using remote sensing technology. Landsat ETM+ and OLI satellite images for year 1999, 2009 and 2019 were used to assess its land use land cover (LULC), change detection, spectral indices and model biomass based on linear regression. The result of the LULC showed that there was an increase in forest vegetation by 12.5% from 1999 to 2009 and a further increase by 2.3% in 2019. However, plantation cover had decreased by 3.5% from 1999–2009; while wetlands had also decreased by 9% within the same period. There was, however, an increase in plantation cover from 2009 to 2019 by 2.1% but a further decrease in wetlands of 4.3%. Further results revealed a positive correlation across the three decades between the widely used Normalized Differential Vegetation Index (NDVI) with other spectral indices such as Enhance Vegetation Index (EVI) and Normalized Differential Moisture Index (NDMI) for biomass were: for 1999 EVI (R2 = 0.1826), NDMI (R2 = 0.0117), for 2009 EVI (R2 = 0.2192), NDMI (R2 = 0.3322), for 2019 EVI (R2 = 0.1258), NDMI (R2 = 0.8148). A reduction in the total carbon stock from 14,221.94 megatons in 1999 to 10,342.44 megatons 2019 was observed. This study concluded that there has been a reduction in the amount of carbon which the Biscay Forest can sequester.


Introduction
The environment around us today has been impacted to varying degrees by anthropogenic activities. The largest terrestrial reservoirs of carbon on our planet are in the forest ecosystem which plays an important part in the global carbon cycle [1] since it absorbs carbon dioxide through the process of photosynthesis [2]. However, in some places around the world, this valuable green environment is being depleted as population increases in order to meet needs such as food, energy, construction, shelter, etc. The effect of this depletion is the increase in carbon dioxide concentration, another greenhouse gases in the atmosphere [3] and climate change. Nature based solutions are increasingly being promoted, including reforestation and afforestation, with the hope that they will enhance the carbon stock and mitigate the adverse impacts of greenhouse gases and climate change [4].
Also, the United Nations Framework Convention on Climate Change (UNFCCC) committed countries to reducing concentrations of greenhouse gases in the atmosphere, thereby ensuring the "stabilization of greenhouse gas concentrations" which will allow the ecosystem time to adapt naturally to climate change [5]. The Kyoto Protocol (KP), which was the main implementing instrument, sets emissions targets for developed countries which are binding for periods lasting from 2008-2012 and another from 2013-2020 for the emissions of different GHGs [6], land use, land use change and forestry activities (LULUCF) [7,8]. Some major industrialized countries did not ratify the Kyoto Protocol.
However, there has been an increase in the awareness of the place of the green environment in combating climate change, hence the need to determine the quantity of forest and how much carbon it can sequester. The Paris Agreement in 2016 recognized the major role forests play in climate change mitigation and aimed to limit global warming to less than 2°C, while pursuing efforts to limit the rise to 1.5°C compared to the pre-industrial period, and the EU regulation 2018/841 for LULUCF [9]. Also, rules and procedure were defined in Durban Climate Change conference [5] that introduced changes to the first commitment period (2008-2012) of the Kyoto Protocol to include the mandatory accounting of forest management with a new forest carbon accounting method, among others. In addition, the EU defined its own target within its Climate Policy Framework 2030 (2021-2030) to reduce greenhouse gas (GHG) emissions to 40% below 1990 levels by 2030 [10].
According to [11], the Spanish forests consist of a varied composition and complex structure of more than 150 tree species, most of which are located within mountainous regions. The forest is dominated by Pinus, Quercus, Fagus, Abies or Betula species. The report of [12] stated that in the Mediterranean area, the vegetation includes the open woodlands (dehesas), dense forests (coppices) dominated by Quercus and Fraxinus along with pinewoods, productive plantations of Populus and Eucalyptus. Near natural forest coexists with homogeneous coniferous reforestations plantations of Pinus and Eucalyptus in the Atlantic region [12]. The forest biomass of Biscay is reported to be about 66% evergreen forests throughout the province, with exceptions in Enkarterri and Durango, while in Gipuzkoa about 56% of the trees are conifers, and only 28% in Álava [13].
A previous study investigated the biomass production of predominant forest species Pinus radiata D. Don and Ecualyptus globulus Labill of Biscay using data from the Spanish Forest Inventories and reported that the Forests of Biscay stored 12.084 teragrams of biomass of dry basis in 2005 and 14.509 terahrams of biomass of dry basis in 2011 [14]. Although, the forest is a renewable resource, its accounting and management is vital for sustainability. Remote sensing technology, which involves the acquisition of information without coming into direct physical contact, has over the years been used to assess, quantify and manage forest lands. Generally, remote sensing data are recorded by sensors or imaging systems carried on platforms including: satellites, aircraft, ground, and unmanned aerial vehicles (UAV) etc. The selection of a platform is dependent on a number of factors which could include, altitude, times, cost, the desired instantaneous field of view (IFOV) etc. Its temporal, spectral, radiometric and spatial resolutions enable information to be acquired at a variety of scales which enables the modeling of forest conditions and changes under different scenarios and from a few cm to some km, with the most commonly used device being the optical form [15,16]. Previous studies have shown that remote sensing technology have been used to investigate forest growth on a spatial and temporal bases [17]; modeling forest cover attributes in a regional context [18]; phenological differences in Tasseled Cap indices improve deciduous forest classification [19]; field measurements in estimating the above ground biomass (AGB) of forests [20] and field measurement for biomass estimation [21].
Spectral indices which are a combination of spectral reflectance are indicators developed on the simple mathematical formula at given wavelengths that describes the condition of vegetation and estimate the quantity of biomass [22,23]. Several studies have used satellite based spectral indices as indicators to monitor changes to vegetation [24][25][26]. This study, therefore, used remote sensing to assess the green space land cover, spectral indices and biomass of Biscay, Spain from 1999 to 2019; evaluated spectral indices for the estimation of biomass and modeled the above ground biomass and carbon stocks based on linear regression model.

Study Area
The Biscay Province is located in the northern part of Spain located within Latitude 42°57′40′′ N, Longitude 02°21′08′′ W and Latitude 43°34′02′′ N, Longitude 03°29′11′′ W (Fig. 1). It is bordered by the communities of Cantabria and Burgos to the west, the province of Gipuzkoa to the east, and Álava to the south, and the Bay of Biscay to the north. According to the standard climate values for Bilbao [27], the region's climate is oceanic and experiences high precipitation all year round. The average temperature in cities like Bilbao falls between 13°C in January and 26°C in August while there are more extreme temperatures in the higher lands of inner Biscay, where snow is more common during winter.

Data Acquisition
This study used Landsat TM, ETM+ and OLI covering 1999, 2009 and 2019 within the same season. They are of Level 2 processing acquired from earthexplorer.usgs. gov in the United States Geological Surveys (USGS) Global Visualization online portal. To correct for measurement and geometry computation difficulties, all images and administrative shape files were harmonized to fit into a uniform coordinate system ETRS89 UTM 30N. In addition, Google Earth Image was used as ancillary data for this study. Table 1 shows more information on the data.

Satellite Image Pre-Processing
These images were subsequently subjected to digital image processing (DIP) in ERDAS Imagine 2014 and ArcGIS Pro software; and the bands stacked. In layer stacking, the extracted individual bands for each Landsat scene were stacked into a single multispectral scene and a subset of the Area of Interest (AOI) extracted [28]. This was clipped to the administrative boundary shape file (.shp) of the Basque Country for analysis.

Land Cover Classification
This processed subset image was categorized into five land cover classes for interpretation after the order of the Anderson Land Cover Classification Scheme [21,29,30]. The land cover classes are: built-up, water body, forest, plantation and wetlands. Next, a step-by-step process of training sample selection using a combination of the spectral signatures of each class and data from Google Earth was performed. The image signature editor in ERDAS Imagine 2014 software was used to facilitate the delineation and extraction of training samples and imagine signatures using the Maximum Likelihood Algorithm/Bayesian classifier [31].
The maximum likelihood decision rule is based on the probability that a pixel belongs to a specific class. The basic equation assumes that these probabilities are equal for all classes, and that the input bands have normal distributions. The equation for the maximum likelihood/Bayesian classifier is as follows [31]: where: D -weighted distance (likelihood), c -particular class, X -measurement vector of the candidate pixel, Mc -mean vector of the sample of class c, ac -percent probability that any candidate pixel is a member of class c (defaults to 1.0, or is entered from a priori knowledge), cov c -covariance matrix of the pixels in the sample of class c, |cov c| -determinant of cov c (matrix algebra), (cov c) -1 -inverse of cov c (matrix algebra), ln -natural logarithm function, T -transposition function (matrix algebra).
The pixel is assigned to the class, c, for which D is the lowest.
Once the maximum likelihood algorithm (also called the supervised classification) was performed on each of the Landsat images, they were checked for accuracy. This accuracy assessment was done using the Google Earth Imagery data as the referenced data and the classified image as the observed data.

Land Use Land Cover Change (LULCC)
Land use land cover change involves assessing the changes that have taken place between different land cover types over a period of time. This can be computed in terms of a percentage [32]. The LULCC to determine the trend from 1999 to 2019 was calculated using: where: P -percentage change of land use/land cover for a particular purpose within a specified time interval, A -area under that particular purpose of land use/land covers after the time interval, B -area under that particular purpose of land use/land covers before the time interval.

Spectral Indices
Several studies have established the use of spectral indices to assess the greenness, health status and chlorophyll content of forest vegetation which is done by using various band combination for its computation [33][34][35][36]. In this study, five spectral indices were performed using Landsat images of 1999 to 2019 for this study. These indices majorly covered bands for forest, water and built-up areas. The following relevant spectral indices were used to assess the study area:

Normalized Differential Vegetation Index (NDVI)
The NDVI measures the health/green vegetation. Vegetation reflects highest in the near infrared band than the red band. This is a combination of its normalized difference formulation and use of the highest absorption and reflectance regions of chlorophyll [37]. The value of this index ranges from −1 to 1, where −1 means no vegetation and 1 means high vegetation:

Normalized Differential Built-Up Index (NDBI)
This index enhances urban areas giving it a higher reflectance in the shortwave-infrared (SWIR) region, compared to the near-infrared (NIR) region [37]:

Modified Normalized Differential Water Index (MNDWI)
This index maximizes the reflectance of water by using green wavelengths and minimizes the low reflectance of NIR by water features [38]. Thus, the water features are enhanced and have positive values while vegetation and soil are suppressed, usually have zero or negative values: where: Green -green band such as TM band 2 (Landsat 7 and band 3 for Landsat 8), NIR -near infrared band such as TM band 4 (Landsat 7 and band 5 for Landsat 8).

Enhanced Vegetation Index (EVI)
This index is used to quantify vegetation greenness and is similar to the Normalized Difference Vegetation Index (NDVI), however, this corrects for some atmospheric conditions and canopy background noise and is more sensitive in areas with dense vegetation [39]. EVI values should range from 0 to 1 for vegetation pixels:

Normalized Differential Moisture Index (NDMI)
This index is used to determine the vegetation water content using the ratio between the NIR and SWIR. The SWIR reflects the changes in both the vegetation water content and the spongy mesophyll structure in vegetation canopies, while the NIR is affected by leaf internal structure and leaf dry matter content and not by water content [40]. The valid range is from −10,000 to 10,000 and a scale factor of 0.0001: In Landsat 4-7, NDMI = (Band 4 -Band 5)/(Band 4 + Band 5).

Biomass Computation
All the bands that reflect vegetation and biomass can be quantified; however, the near infrared band absorbs most of the chlorophyll it receives and reflects vegetation at its highest. In this study, band 4 was used in computing the biomass for TM, ETM+ and band 5 for OLI images of Landsat. This applies to 1999 and 2009 and band 5 for 2019. About 30 random points (vector) of forest biomass were selected from the Landsat images and only 29 points returned with biomass values. It is possible that one point was out of the study area (Biscay). This analysis was done in QGIS 3.4 software. Linear regression was used correlate the relationship between indices. Furthermore, the above ground biomass (AGB), below ground biomass (BGB) and carbon stock was computed. The AGB was computed by first removing pixel without data, then reclassifying and assigning classes to them. The zonal statistics was computed and used to compute the biomass based on the number of pixels (for forest, plantation and wetlands with vegetation classes) multiplied by the spatial resolution of the image and its mean. The equation is [41]: The BGB is 20% of the AGB. To compute for the total carbon stock, the total carbon density or total carbon stock was converted to CO 2 by multiplying carbon density or stock by 3.67. This is the ratio of molecular weights between carbon dioxide and carbon [42].

Land Use Land Cover Classification
The result of the land use land cover classification of Biscay can be seen in Figure 2.

Land Use Land Cover Change
The area covered by various land cover classes and their percentage changes in Biscay over a three decade period can be seen in Table 2. The land use land cover shows that there has been an increase in the forest vegetation by 12.5% from 1999 to 2009 and a further increase by 2.3% in 2019. However, the plantation area decreased by 3.5% from 1999 to 2009; while wetlands also decreased by 9% within the same period. There was, however, an increase in plantation areas from 2009 to 2019 by 2.1% but a further decrease in wetlands of 4.3%. In addition, it was observed that the water body had decreased by 2.5% from 1999 to 2009 but increased by 0.4% by 2019. The built-up area increased by 2.5% from 1999 to 2009 and seems to have decreased by 0.5% by 2019. While 0.5% could be an insignificant change over a period of ten years, it is possible that some structures that were present in 2009 had been pulled down by 2019. Also, inference could be made from the type of sensors used to acquire both data in which Landsat 8 sensor can discriminate more finely than the Landsat 7. However, there was an overall increase in built-up areas from 1999 to 2019 of 4,563 ha while water bodies, plantations and wetlands had decreased by 4,534 ha, 3,181 ha, and 29,499 ha respectively.

Spectral Indices
The result of the five spectral indices is shown below in pseudo colour (Figs. 3-5). The result indicates the varying change that has taken place over the years obvious in their grey colour as the years increase. It was observed that highest values for the NDVI in 1999 was 0.94 and the lowest −0.57. The NDBI values ranges from 0.36 to −1.14, MNDWI from −0.76 to 1.12, EVI from −0. 16

Fig. 5. Spectral indices for 2019
Source: own elaboration derived from Landsat using QGIS 3.4 software It can be observed from Figure 5 that there are varying intensity changes of all the spectral indices over the three decades. From the NDVI and EVI, it is clear that the luxuriant forest biodiversity, greenness and chlorophyll content found in 1999 had diminished by 2019; however, the MNDWI and NDMI increased within the same period.

Spectral Indices and Biomass Computation
The result of the random points (vector) of forest biomass selected from the Landsat images and the spectral indices for the three decades are shown in Tables 3-5. The values are the reflectance values in which during the preprocessing, digital numbers (DN) were converted to radiance and radiance to reflectance. The NDVI which is the most common vegetation spectral indices was used to correlate its relationship to other indices using the linear regression (Figs. 6-8) for the three decades.

Biomass Computation
This was performed by removing the pixels with no data, then reclassifying the entire image of the study area, computing the zonal statistics and subsequently computing the above ground biomass. The result of the estimated above ground biomass for the three decades of the entire study area for three land cover types: forest, plantation and wetland are presented in Table 6. It can be seen from Table 6 that there is no value for the wetland land cover type for 2019. This is because there was a decrease in the wetland area from 2009 to 2019 of 4.3%, as seen in Table 2. Thus, the amount of vegetation in wetlands was drastically reduced. Although there was an increase in the forest cover from 2009 to 2019, there was no vegetation (biomass) in the wetlands to be accounted for.

Below Ground Biomass (BGB)
BGB is estimated from AGB using a non-destructive method [43]. This method uses 20% of the above ground biomass to compute the below ground biomass values for vegetation [21,43]: In 1999:

Carbon Stock
The total carbon density or total carbon stock was converted to CO 2 by multiplying carbon density or stock by 3.67. This is the ratio of molecular weights between carbon dioxide and carbon [42]. The result is shown below.

Land Use Land Cover
This study revealed that Landsat images can be used to assess the land use land cover change and carbon stock over a period of three decades. This study revealed that the land cover types of forest, plantation and wetlands in 1999 covered about 75,632 ha, 91,630 ha, and 91,280 ha respectively. However, changes occurred over the decades, with these values being replaced by 108,284 ha (forest), 88,448 ha (plantation), and 1781 ha (wetland) by 2019. According to Spanish National Forest Inventory, Spain has a total forest area of about 50,560,000 ha which is divided between 50 of its provinces and 17 of its communities [11]. Thus, the Spanish green space is rich and buoyant. However, the findings from this study revealed some dynamics in which forest lands seem to have increased from 1999 to 2019 by 14.7% while plantation areas have been reduced by 1.4%. It was observed that the wetland areas were reduced by 13.3% within the same period. The increase in forest lands could be as a result of private forest ownership. Alberdi et al. [11] reported that about 70.89% of the forest in Spain is privately owned. This makes for proper efficient management and control of anthropogenic activities (deforestation). The changes in plantation land cover type observed in this study (which include farmlands) reduced by 1.4% while built-up areas increased by 2.1% from 1999 to 2019. It is likely that certain farmland was replaced by structural development as population increased over the three decades.
In addition, there was a drastic change observed in the wetlands of Biscay. The wetlands are fast disappearing and were reduced by 13% from 1999 to 2019. Climate change could be a major contributory factor to this reduction. The water body was also affected as it fell by 2%. Another factor may be the heating up of the region as the population of Biscay increases over the years. This agrees with [44] which reported that pasture and population centres were sources of the emission of about 0.18 million Mg of CO 2 .

Spectral Indices
This study used spectral indices to assess changes in the biomass of Biscay. The indices are used to monitor the greenness and health of the vegetation. It was observed that biomass had undergone one form of change or another over the last three decades. Random points (vector) extracted from the Landsat multispectral bands were compared to the spectral index images over the years. These random points and values were extracted from the bands with the highest reflectance: band 4 of Landsat 5-7 and band 5 of Landsat 8 and used to compute the linear regression. Linear regression was used to correlate the relationship between indices. This study correlated for the four indices with the NDVI to ascertain its relationship. NDVI is widely known to measure the health/greenness of vegetation [45]. Findings revealed that NDVI, EVI, and NDMI had a direct correlation with the forest biomass of the region. NDBI and MNDWI had a negative relationship to biomass. Moisture had a positive correlation to biomass than water. It was observed, however, that the moisture content had been reduced from 1999 to 2019. This could be seen in the values of the NDMI over the years. This could also be one of the reasons for the reduction in the overall greenness of Biscay over the last three decades as observed in the land cover classification.

Computing Biomass and Carbon Stock
This study used remote sensing technology to compute the biomass of Biscay over three decades. The findings revealed that the total above ground biomass (which includes forest, plantation and wetland areas with vegetation) had fallen from 1999 to 2019. It was further revealed that by 2019, there were no wetlands with vegetation to be found in Biscay. This affected the quantity of biomass in this region. The below ground biomass had also reduced from approximately about 646,000,000 ton/ha in 1999 to about 470,000,000 ton/ha in 2019. Furthermore, the total carbon density or total carbon stock was converted to CO 2 (2012). The findings revealed that amount of carbon that the green spaces of Biscay will be able to sequester had also fallen over the decades. The amount of carbon being sequestered was 14,221.94 megatons in 1999 but had fallen to 10,347.44 megatons by 2019.

Conclusion
This study used remote sensing technology to assess forest biomass in Biscay, with field measurements not being conducted in this study. Biomass computation and assessment uses spectral indices. Therefore, a forest landmass can increase but its content including the greenness of the forest can recede. Spectral indices give information about NDVI, EVI, the chlorophyll content etc. of forest cover. The spectral reflectance of plant species of the forest in 1999 is higher than that of 2009 despite the increase in land mass. This can be seen by the intensity of the greenness in Figure 5.
In addition, the findings show that most of the recent forest cover in Biscay is not primary forest, but forest-regrowth.
Landsat images of the three periods (1999-2019) were used to assess the green space of Biscay. Although the Biscay province is still a green environment, the amount of carbon it can sequester has fallen over the years. Green space and forest play a vital role in mitigating the impact of climate change and thus should be managed adequately. It seems the existence of private forests have had a positive impact on the biomass quantity of Biscay and thus this should be promoted. Also, the planting of trees and other forest policies that enhance sustainability should be strongly encouraged.