Monitoring of Spatiotemporal Change of Green Spaces in Relation to the Land Surface Temperature: A Case Study of Belgrade, Serbia

Understanding the relationship between land use and land cover and thermal environment has recently become an emerging issue for urban planners and policy makers. We chose Belgrade, as a case study, to present a cost- and time-effective framework for monitoring spatiotemporal changes of green spaces in relation to the land surface temperature (LST). Time series analysis was performed using Landsat 5 TM and Landsat 8 OLI/TIRS imagery from 1991 to 2019 with an approximate 5-year interval (18 images in total). Spectral vegetation indices and supervised land cover classifications were used to examine changes of green spaces. The results showed a fluctuating trend of the normalized difference vegetation index (NDVI) and the normalized difference water index (NDWI). The highest values were recorded in 2019, indicating vegetation recovery in the last decade. A significant positive correlation was determined between the spectral vegetation indices and the amount of precipitation during growing season. The land cover classification showed that the share of vegetated and bare land decreased by 11.74% during the study period. The most intensive conversion of green and bare land into built-up land cover occurred in the first decade (1991–2000). To assess spatiotemporal changes in the LST, Landsat Collection 2 Surface Temperature products were used. We found a negative correlation between change in the spectral vegetation indices and change in the LST. This indicates that the reduction in vegetation was associated with an increase in the LST. The municipalities that were the most affected in each decade were also identified with our framework. The findings of this study are of great relevance for actions targeting an improvement in urban thermal comfort and climate resilience.


Introduction
Urban areas are characterized by higher air temperatures compared to rural areas mainly due to an altered surface cover [1,2]. Artificial materials absorb and store a significant amount of solar heat during the day, which is slowly realized during the night due to their high thermal inertia. Vegetation mitigates urban heat through evapotranspiration, shading and absorption of air pollutants [3][4][5]. Consequently, the role of vegetation and water bodies to regulate local climate in cities is of crucial importance for the well-being of the city residents.
In the last few decades, the number of heat waves in Europe significantly increased and metropolitan areas are regarded as particularly vulnerable [6]. Therefore, understanding the relationship between land use and land cover (LULC) and thermal environment has recently become an important task for urban planners and policy makers [7,8]. We chose the city of Belgrade as a case study, as it is one of the largest metropolitan regions in southeastern Europe that is experiencing urban growth and expansion. During the transition from socialism to capitalism (over the last 3 decades), suburbanisation took place rapidly without a coherent planning policy. Massive transformations of non-urban into urban land has mainly been driven by the rural-to-urban migration of economically disadvantaged citizens in search of better economic conditions [9,10]. Although it is well known that vegetation is the first to suffer during urban growth and expansion [11], changes in green spaces in Belgrade during the transitional period have not been sufficiently explored.
Thermal infrared remote sensing has been widely used for the estimation of land surface temperature (LST) and the assessment of the surface urban heat island [12,13]. Landsat data products have provided medium resolution observations of land surface since the early 1980s and are available free of charge, making them a great choice for long-term monitoring of the LST and land cover characteristics. Numerous studies have investigated the relationship between LST and vegetation abundance or health, using the normalized difference vegetation index (NDVI) [8,[14][15][16][17] and the normalized difference water index (NDWI) [18,19]. The inverse relationship of LST with the NDVI and the NDWI determined in these studies suggests that vegetation has a mitigative effect on the surface urban heat island. However, studies dealing with large spatial and/or temporal scales have reported contradictive results regarding the relationship between LST and land cover characteristics, depending on the year, season, time of the day or geographic location [20][21][22][23][24]. This indicates the complexity of this relationship and the importance of its continuous monitoring at a local level.
Spatial and temporal patterns of the LST retrieved from satellite records for Belgrade were the subject of a study conducted by Pongrácz et al. [25], but only for a short time period (2001)(2002)(2003) and in a 1 km spatial resolution. To the best of the authors' knowledge, the relationship between LST and land cover characteristics for the observational area of Belgrade has not been investigated so far. With the General Urban Plan of Belgrade 2021 (accepted in 2003/last updated in 2016, Official Gazette of the City of Belgrade No. 11,2016), the problem of the loss and degradation of green spaces was addressed and guidelines for sustainable urban development were proposed. The main goal of this study was to support ecologically based urban planning, with a novel framework for the monitoring of spatiotemporal changes of green spaces and the LST. This study addresses the following questions:
What was the relationship between vegetation indices (the NDVI and the NDWI) and climate factors? 3.
What was the relationship between the change in vegetation indices (the NDVI and the NDWI) and change in the LST in different municipalities?
The findings presented in this paper could help urban planners and ecologists to gain new insights into vegetation change and its impact, and it could guide decision and policy making towards the establishment of a high-quality, multifunctional and interconnected system of green and blue spaces, thus maintaining and enhancing the benefits delivered by those spaces.

Study Area and Population
The city of Belgrade is the capital of Serbia and one of the largest cities in southeastern Europe, with nearly 1.7 million inhabitants and an administrative area of 3234 km 2 . The urban growth (increase in population) of Belgrade was accompanied by an urban expansion of 23.56% for the period 1991-2011 [9,10]. The study area includes the territory of Belgrade within borders defined by the General Urban Plan of Belgrade 2021 and covers 778.51 km 2 ( Figure 1). This area is a part of the wider administrative unit of Belgrade and is divided into 12 municipalities. Changes in population density for municipalities according to the censuses from 1991, 2002 and 2011 are visualized on maps using QGIS version 3.16 (QGIS Geographic Information System, Open Source Geospatial Foundation Project, http://qgis.org, accessed on 22 November 2020) ( Figure 2).  Belgrade is located on the confluence of the Sava and Danube rivers and in the contact zone of the southern ridge of the Pannonian basin and the northern border of the Balkan Peninsula. North of the rivers, the terrain consists of alluvial plains and loess plateaus. The terrain rises gradually from north to south, and the southern part of Belgrade lies on hilly terrain, with the highest point at Avala Mountain (511 m a.s.l.). The territory of Belgrade is characterized by a temperate-continental climate with local variations. The mean annual air temperature is 11.7 • C and the average annual precipitation is 666.9 mm. The conditions for vegetation growth are favorable because most of the annual precipitation falls in May and June (when plants need it most). June is the month with the highest precipitation (an average of 86.6 mm), while July is the warmest month with a mean temperature of 22.1 • C [26]. According to the General Urban Plan of Belgrade 2021, forests cover 7.2%, whereas public green spaces (parks, squares, green spaces along the banks of the Danube and Sava, protective green belts, green corridors, wetlands and nurseries) occupy 1.8% of the total urban area. The estimated green space per capita (forests and public green spaces) is 52 m 2 . However, the share and dynamics of green spaces associated with different land uses (aside from forests and public green spaces) are unknown.

Satellite Imagery
Landsat time-series data from 1991 to 2019 with (approximate) 5-year intervals was selected for extracting information on land cover and LST in Belgrade. Collection 2 Level-2 Landsat products were acquired from the US Geological Survey through the EarthExplorer data portal (https://earthexplorer.usgs.gov/, accessed on 4 August 2021). In total, 18 images were downloaded for the selected years captured under clear atmospheric conditions during the growing season (May-September). Important details of Landsat data products used in this study, as well as climate data, are presented in Table 1. The methodological workflow and steps in data analysis are described in Figure 3.

Calculation of Spectral Vegetation Indices
For the remote sensing of vegetation, the most useful spectral radiances are in the red and near-infrared (NIR) regions. Healthy vegetation strongly absorbs visible light for photosynthesis and has high reflectance in the NIR region. Unhealthy or sparse vegetation, on the other hand, will generally reflect more visible light and less NIR light. As the NIR, red and short-wave infrared (SWIR) spectral regions are the most sensitive for green vegetation and vegetation changes, these bands are employed to calculate the NDVI [27] and the NDWI [28].
The NDVI is the difference between NIR and visible red reflectance values normalized over total reflectance. It is calculated thus: The NDVI values range from −1.0 to 1.0, where negative values generally represent cloud, snow or water. A higher NDVI indicates a greater level of photosynthetic activity and a greater quantity of green biomass.
The NDWI is used for remote sensing of water content of vegetation, and it is considered to be a good proxy for plant water stress [29]. The NDWI is derived from the reflectance radiated in the NIR and SWIR bands: The NDWI values vary from −1.0 to 1.0, depending on the water content, as well as the type of vegetation and cover. The high NDWI values correspond to high vegetation water content and to high vegetation fraction cover.
The calculation of vegetation indices was performed in QGIS. The Landsat band designations used for the calculation of vegetation indices are presented in Table 2. First, the NDVI and the NDWI were calculated for each image from Table 1. Then, the average values for pixels were estimated based on all available images within one growing season. For a summary of temporal changes in Belgrade and in the linear regression analysis, medians were used because distributions of the NDVI and the NDWI deviate from normal. For a summary of seasonal variations, an average value for each month was calculated based on all available images taken in that month.

Land Cover Classification and Accuracy Assessment
In order to visualize and estimate changes of green spaces in Belgrade during the last 3 decades, land cover classification was performed for selected images from 1991, 2000, 2011 and 2019 (nearly the beginning and end of each decade).
The selection of land cover classes mainly depends on the aim of the study, the study area and the resolution of the input data [30]. Given that the focus of this study is vegetated areas, and the spatial resolution was constrained to 30 m, the following land cover classes were selected: "green space", "bare land", "water" and "built-up land".
For the future relevance of this research, it is very important to specifically define the land class "green space". The authors consider green space a land covered with some form of vegetation-including forests, parks, grasslands and crops. This is the most common definition of green space [31]. One uncertainty that arises from this classification is related to agricultural areas under rotation systems used for annually harvested crops and fallow lands. If such an area is under crop at the time of observation, it will be classified as "green space", otherwise it will be classified as "bare land". Built-up land cover is characterized as an artificial and/or impervious cover, and it includes residential, commercial, industrial, transportation and communication utilities.
Training samples were selected directly from the Landsat images using a priori knowledge of the land use and Google Earth historical imagery. The samples of different land covers were randomly chosen using the polygon of pixels sampling procedure. Around 250 training polygons were collected for each Landsat image.
Supervised classification was performed in QGIS, employing the EnMAP-box Plugin [32]. A random forest was used as a classifier, and the number of trees was set to 100. The random forests are tree-based classifiers that include k decision trees. Each individual tree is trained over an independent random sample from the training data set and casts a unit vote for the most popular class. Previous studies dealing with classification in a heterogenous urban environment found that the performance of a random forest is superior to other classifiers [33,34].
In the post-processing of the classification map, the "majority filter" tool was applied in QGIS. It generalizes the map and reduces single pixel misclassifications. For an accuracy assessment of the classification, a cross-validation with three-folds was employed (2/3 of the sample was assigned as training and 1/3 as a test sample).
Based on land cover classification maps, land cover change maps were produced for three decades: 2000-1991, 2011-2000 and 2019-2011. Land cover change maps are a very useful instrument to assess and visualize where green spaces are under pressure or where additional green spaces are needed in order to create an interconnected system. Of interest to this study, there were changes from green and bare land into built-up land cover and vice versa.

Relationship between Spectral Vegetation Indices and Climate Factors
The relationship between spectral vegetation indices and climate factors was examined using linear regression analysis in RStudio version 1.3.1093 (The R Foundation for Statistical Computing, https://www.r-project.org/, accessed on 20 August 2021). Regression analysis was performed based on medians of the vegetation indices, average air temperature and sum of precipitation in the growing season (Table 1).

LST Retrieval
In order to assess spatiotemporal patterns of LST in Belgrade, Landsat Collection 2 Surface Temperature products were utilized. These products are atmospherically corrected and resampled to a 30 m resolution. Thermal data (Landsat 5 TM band 6 and Landsat 8 OLI band 10) was converted to radiance and temperature in Celsius. Data for around 20% of pixels was missing from each image.
First, the LST was retrieved for each image from Table 1. Then, the average values for pixels were estimated based on all available images within one growing season. Before the calculation of differences between growing seasons (2000 and 1991, 2011 and 2000, 2019 and 2011), the LST was normalized following recommendations of previous studies [35,36]: where LST i is the LST at pixel i, and LSTmax and LSTmin refer to the maximum and minimum LST in the study area. The LST retrieval was performed in QGIS.

Relationship between Change in the LST and Change in Spectral Vegetation Indices
The relationship between change in the LST and change in vegetation indices was examined using linear regression analysis in RStudio. In this analysis, different municipalities were compared between growing seasons (2000 and 1991, 2011 and 2000, and 2019 and 2011). The change in the LST was estimated based on the difference between two growing seasons: where LST GS2 is the LST in a growing season in time point 2 and LST GS1 is the LST in a growing season in time point 1.
The change in vegetation indices was estimated based on the difference between two growing seasons: where NDVI GS2 and NDWI GS2 are vegetation indices in a growing season in time point 2, and NDVI GS1 and NDWI GS1 are vegetation indices in a growing season in time point 1.

Spatiotemporal Patterns of Spectral Vegetation Indices
For each growing season, the average NDVI and NDWI were calculated and are presented in Figure   Seasonal variation was presented based on the average values of all calculated vegetation indices within each month. During the growing season, the NDVI was the same in May, June and August, and it was slightly higher in July. The NDWI was more variable and had the peak in June. A drop in both indices was recorded in September (Figure 4b).
For the comparison of different municipalities over time, spectral vegetation indices were summarized and are presented in Figure 5. The highest NDVI were recorded for the municipalities of "Voždovac", "Čukarica", "Grocka", "Rakovica" and "Palilula". These municipalities are mainly situated in the peripheral zone. The municipalities with the lowest NDVI were "Stari Grad", "Vračar" and "Novi Beograd", which are located in the central and middle zone of Belgrade. A similar pattern was recorded for the NDWIthe highest values appeared in the peripheral municipalities of "Čukarica", "Voždovac", "Grocka" and "Palilula", and the lowest values were again in "Stari Grad", "Vračar" and "Novi Beograd".

Land Cover Change
Classification maps were generated for 1991, 2000, 2011 and 2019, practically the beginning and end of each of the three decades ( Figure 6). The classification resulted in an overall accuracy of between 93.57% and 96.95% and a Kappa coefficient from 0.89 to 0.95 ( Table 3). The lowest accuracy was determined for the differentiation between the classes "bare land" and "built-up", especially for 2000, indicating that these two classes have the most similar spectral signature.

Land Cover Change
Classification maps were generated for 1991, 2000, 2011 and 2019, practically the beginning and end of each of the three decades ( Figure 6). The classification resulted in an overall accuracy of between 93.57% and 96.95% and a Kappa coefficient from 0.89 to 0.95 ( Table 3). The lowest accuracy was determined for the differentiation between the classes "bare land" and "built-up", especially for 2000, indicating that these two classes have the most similar spectral signature.  Absolute and relative areas of each class were also calculated and are presented in Table 3. A general increase in built-up land cover is noticeable during the entire study period (from 14.05% in 1991 to 26.17% in 2019). Urban expansion toward peripheral areas and the further densification of built-up land cover in the central city zone can be clearly seen in Figure 6. The growth of the city was accompanied by the conversion of land cover classes, mainly "green space" and "bare land" into "built-up". These changes are represented in Figure 7. The share of vegetated and bare land decreased by 11.74% from 1991 to 2019. The most pronounced conversions of the classes "green space" and "bare land" into built-up land cover occurred during the first decade (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000), while the smallest change was recorded during the third decade (2011-2019). Reverse land cover changes (from built-up land cover into "green space" and "bare land") could also be observed, but to a much lesser extent. They mainly represent abended properties undergoing natural succession. In the spatial context, significant changes occurred on the edges of urban forests in the middle city zone (the so-called "inner ring" of green spaces), as well as in the southern part of Belgrade within suburbs where green spaces and agricultural areas were converted into built-up area. In the western and northern part of Belgrade, urban expansion took place along main traffic routes, mainly at the expense of agricultural land.
In addition, it can be noticed that many smaller green spaces within the central city zone were converted into built-up areas over the whole study period.

Relationship between Vegetation Indices and Climate Factors
A significant positive linear relationship was found between both the NDVI and the NDWI, and precipitation during the growing season (R 2 = 0.70 and R 2 = 0.55, respectively, Figure 8). Years with a lower amount of precipitation during the growing season (2000, 2011 and 2015) were associated with a lower NDVI and NDWI, whereas 2019 was the year with the highest amount of precipitation and the highest values of vegetation indices.
Between vegetation indices and average air temperature during the growing season, a linear relationship was not determined (p > 0.05).

Spatiotemporal Change of the LST in Relation to Spectral Vegetation Indices
Between the change in the NDVI and the change in the LST, a significant negative relationship was determined in all three decades (R 2 = 0.75, R 2 = 0.61 and R 2 = 0.60, respectively, Figure 9). In the first two decades, the municipalities with the largest decrease in the NDVI also had the largest increase in the LST. In the third decade (2019-2011), the NDVI increased in all municipalities. Again, the municipalities with the smallest increase in the NDVI had the largest increase in the LST.

Changes of Green Spaces in Belgrade during the Last Three Decades
The diversity, composition and structure of urban vegetation primarily depend on the ecoregion in which the city is situated, being constantly modified by the urban growth and expansion in a given spatial and temporal frame [37,38]. The urban expansion of Belgrade was at its most intensive in the first decade of the studied period (1991-2000), but also continued in the following decades. It was driven by the migration of Serbian refugees from former Yugoslavia countries in the 1990s and the rural-to-urban migration of economically disadvantaged citizens. Urban expansion took place in areas of small forests, forest patches, areas under natural succession, agricultural and bare land-which are all spaces of great importance for biogeochemical cycles [39]. Considered particularly harmful are the changes in the "inner-ring" of green spaces in the middle city zone, which included the edges of the urban forests "Košutnjak" and "Zvezdara", as well as the lakeside of "Veliko blato". The high value of these spaces is reflected in the enhancement of living conditions, the contribution to adaptation to climate change and biodiversity conservation (The General Regulation Plan of the System of Green Areas of Belgrade, Official Gazette of the City of Belgrade No. 110, 2019). Another issue to be addressed is the continuous densification of the central city zone, despite a decline in population density. The central city municipalities of "Stari Grad", "Vračar", "Novi Beograd" and "Savski venac" were areas with the lowest NDVI and NDWI throughout the whole study period. The peripheral municipalities were characterized by an increase in population density during the first decade, and in some of them it continued in the following years (Surčin and Grocka). For others (Zemun, Palilula, Cukarica and Voždovac), this trend was reversed in the second decade. It is evident that peripheral areas of Belgrade differ in their demographic and urbanization dynamics, but they still represent the largest reservoir of greenness in Belgrade.
Seasonal dynamics of vegetation are dependent of ecoregion and should be considered in studies that are using remote sensing data to estimate vegetation change. Over the past three decades, we found that the NDVI and the NDWI had its peaks in July and June, respectively, and started decreasing in September. This indicates that maximal photosynthetic activity was reached in July, and maximal water content was reached in June.
Changes in green spaces in Belgrade were also driven by fluctuations in climate factors. Climate, i.e., the seasonal course of solar radiation, temperature and precipitation, determines the type of predominant vegetation on a continental, regional and local level, as well as the biogeochemical conditions that also contribute to general vegetation features (CO2 flux, carbon storage in biomass and soil) [40,41]. We found a significant positive linear relationship between the amount of precipitation during the growing season and the NDVI and the NDWI in Belgrade. Vegetation responds to variations of precipitation on a long-term and short-term time scale, which is confirmed by both on-ground assessments (numerous parameters at biological and ecological levels) and remote sensing (spectral vegetation indices) [42]. Vegetation growth is usually positively correlated with the amount of precipitation, monitored on inter-annual, intra-annual and seasonal scales [43]. The effects of climate factors tend to be co-limiting for vegetation in a particular biogeographic region, i.e., lower temperatures can be coupled with higher precipitation and limited solar radiation. However, the global bioclimatic indices for Earth's vegetated surface show that water availability is the major limiting factor for vegetation growth (40%), followed by temperature (33%) and solar radiation (27%) [44].
Maximal values of vegetation indices and an increase in greenness in 2019 may have different explanations, including driving mechanisms of both climate (an increased amount of precipitation) and human induced changes in LULC. The implementation of new regulations regarding the conservation of natural areas in Belgrade appears likely to have played a significant role in vegetation recovery.
A real challenge in monitoring vegetation distribution and health with remote sensing lies in the possibility to distinguish temporary changes caused by fluctuations in climate factors from long-term changes. Only long-term time series analysis in finer temporal scales, together with the analysis of multiple driving factors, could resolve this uncertainty and provide reliable information on vegetation change [45].

Changes of the LST in Relation to Green Spaces
The inter-dependence and dynamics of the vegetation-climate relationship are recognized through the vegetation feedbacks, which contribute to the global energy and mass fluxes and patterns of climate. In the central and eastern European domain, the evapotranspiration feedback dominates, resulting in overall negative feedback on temperature on an annual average [46]. Evapotranspiration is the sum of transpiration (water transpired by plants through stomata) and evaporation (water evaporated from the soil surface, water bodies and canopy interception). Being a crucial part of the water cycle, evapotranspiration plays a significant role in atmospheric climatic parameters [47]. The cooling effect of transpiration is achieved by the contribution of two processes: i) energy consumed for maintaining this system is heat energy drawn from the air, and ii) evaporation of water through stomata cools leaf surfaces and the surrounding air. The role of transpiration in regulating the urban thermal environment has been acknowledged, based on field data aimed at quantifying its effect on atmospheric temperature [48][49][50]. The intensity of cooling by transpiration depends mainly on the composition (i.e., percent cover) of the green space and its spatial distribution or configuration [51,52].
The assessment of LST patterns and their relationship with land cover on an intraurban scale is of great importance for the improvement of urban thermal comfort and urban climate resilience. Land covers associated with lower LST are water bodies, green areas and bare land [53,54]. However, Gunawardena et al. [55] emphasized the benefit provided by green spaces, because water bodies and bare land may also provide a warming effect during summer.
We found that an increase in the LST in Belgrade was associated with a decrease in the NDVI (in the last three decades) and the NDWI (in the last two decades). These results showed that a reduction in vegetation had a negative impact on the urban thermal environment. A previous study also found a negative relationship between change in the NDVI and the LST between two summers in Beijing, China [56]. Based on the long-term continuous data, Huang et al. [16] recorded a strong negative relationship between the LST and the NDVI and vegetation fraction in Wuhan, China. They also pointed out the importance of a fine-scale temporal resolution for the assessment of the LST-vegetation relationship. Sun and Kafatos [24] found a strong negative correlation between the LST and the NDVI over North America, but only during the warm season (from May to October). For the temperate biome in which Belgrade is situated, the relationship between LST and NDVI shifted from positive in spring (March-May) to negative in summer (June-August) [21]. In the same study, the authors inferred that water availability is the main limiting factor for vegetation growth in a temperate region.
Our framework also identified the administrative spatial units-municipalities-that were most threatened by warming in each decade. This is very useful information for targeted actions against urban heat.

Implications of Study on Urban Planning, Management and Policy
The presented study and many similar studies [57][58][59] demonstrate how freely available remote sensing observations could be applied in urban planning and management. Each study has its own local specificities (regarding climate and LULC change dynamics) but contributes to a global understanding of the urban climate-LULC relationship. There has been a large geographic bias of the cities being studied in favor of Asia and North America [13]. In Europe, the change dynamics of south-eastern cities is the least explored by remotely sensed data and methods.
This study highlights several important issues that need to be considered in urban planning of Belgrade and other cities with similar developmental challenges. In the planning and development of new green spaces, it is important to consider that water availability is the main limiting factor in temperate biome. Therefore, drought-tolerant autochthonous species should be selected for planting. Also, shrubs and trees should be favored, because these forms have a larger cooling potential in comparison to ground vegetation via both transpiration and shading [51].

3.
The pressure of urban development on the edges of urban forests within the "innerring" of green spaces is regarded as particularly harmful. The reduction of forest edges is threatening to diminish the energy flow between the forest and the surrounding areas and reduce shade for adjacent surfaces, thus increasing LST [60]. Although the protection of urban forests as natural resources and common goods is marked as a top priority according to the General Urban Plan of Belgrade 2021, this study showed an inconsistent application of this regulation in practice.
Finally, when it comes to urban planning and management, the main challenge is implementation of regulations regarding sustainable development goals. Detailed largescale urban plans are not always harmonized with higher-order plans and do not strictly comply with their strategic goals for the protection of green spaces. Instead, the protection of green spaces is often compromised to meet the interest of private stakeholders. The important issue is also a lack of measures and instruments for the implementation of principles of green infrastructure, which are mainly discussed in documents only at the conceptual level.

Limitations of the Study
The results of this study are of great relevance for better, more ecologically based urban planning. However, there are a few limitations of the proposed framework that should be considered.
The moderate spatial resolution of Landsat images (30 m) cannot entirely capture the complex spatial patterns of the urban environment. As a result of the presence of "mixed pixels", i.e., the spectral reflectance mixture of different land covers within a pixel, part of the information on land cover variability is lost [61]. In the past decades, a variety of unmixing approaches has been introduced, which can be divided into two main categories: linear and non-linear [62,63]. On multispectral images, the most widely used are linear approaches because of their simplicity and flexibility in different applications. Recently, Yang et al. [64] developed an integrated framework for mixed pixel decomposition and decision tree classification.
Besides the presence of mixed pixels, the differentiation between bare and built-up land cover was challenging in our study because of the similarity in the spectral signatures of these two classes. This led to a misclassification of certain parts of pixels, especially in 2000, and had further implications on land cover change detection. However, we believe that the overall accuracy of the supervised classification is very high and in accordance with similar studies [35,65].
Because the LST is highly variable, the analysis of temporal trend based on Landsat data is very limited. Landsat satellites have a revisited cycle of 16 days, but not all images are usable due to the effects of clouds. It is also important to mention that data for around 20% of pixels was missing for Belgrade from Landsat Collection 2 Surface Temperature products. Therefore, multi-sensor fusion methods should be employed to obtain temporally and spatially continuous LST data [16]. This represents a promising direction for further research on the thermal environment of Belgrade and urban areas in general.

Conclusions
In this study, a cost-and time-effective framework for monitoring spatiotemporal changes of green spaces in relation to the LST was presented. The framework is based on freely available satellite data, open-source software and is easily applied. The city of Belgrade was used as a case study, as one of the largest cities in south-eastern Europe that is experiencing urban growth and expansion.
Landsat 5 TM and Landsat 8 OLI/TIRS imagery was used to assess spatiotemporal changes of green spaces and the LST from 1991 to 2019 with an approximately 5-year interval. Changes of green spaces were examined based on spectral vegetation indices (the NDVI and the NDWI) and supervised land cover classification. The NDVI and the NDWI showed a fluctuating trend during the study period. The lowest values were recorded in 2000 and 2011, while the highest values were observed in 2019, indicating vegetation recovery in the third decade. The central city zone was characterized by the lowest NDVI and NDWI, throughout the whole study period, while the peripheral zone represented the largest reservoir of greenness in Belgrade. Land cover classification was performed for 1991, 2000, 2011 and 2019 and resulted in a very high overall accuracy (from 93.57% to 96.95% and a Kappa coefficient from 0.89 to 0.95). The lowest accuracy was recorded for the classes "bare land" and "built-up", meaning that differentiation between these classes can be challenging in urban areas. During the study period, the share of vegetated and bare land decreased by 11.74%. These conversions from green and bare land into built-up land cover mainly occurred in the first decade (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). Spectral vegetation indices were positively correlated with the amount of precipitation during the growing season. Vegetation recovery at the end of the third decade may have different driving mechanisms, including an increased amount of precipitation and the implementation of new regulations regarding the conservation of natural areas.
The LST was retrieved from Landsat Collection 2 Surface Temperature products. The relationship between the change in the spectral vegetation indices and the change in the LST was examined, and we found that the reduction in vegetation had a negative impact on the thermal environment of Belgrade. Our framework also identified the municipalities that were most affected by an increase in the LST in each decade. This information can be very useful for the actions targeting an improvement of urban thermal comfort and climate resilience. To conclude, we suggest the implementation of our framework (in addition to other freely available remote sensing data) in urban planning and management.  Data Availability Statement: The datasets generated and/or analyzed during the study are available from the corresponding author upon reasonable request.