Next Article in Journal
Thermal Imagery Feature Extraction Techniques and the Effects on Machine Learning Models for Smart HVAC Efficiency in Building Energy
Previous Article in Journal
Monitoring Changes to Arctic Vegetation and Glaciers at Ny-Ålesund, Svalbard, Based on Time Series Remote Sensing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Ecology, Institute for Biological Research “Siniša Stanković”—National Institute of Republic of Serbia, University of Belgrade, Bulevar despota Stefana 142, 11000 Belgrade, Serbia
2
Faculty of Environmental Sciences, Technische Universität Dresden, Helmholtzstraße 10, 01062 Dresden, Germany
3
Department for Strategic Planning and Development, Urban Planning Institute of Belgrade, Bulevar despota Stefana 56, 11000 Belgrade, Serbia
4
Department of Landscape Architecture and Horticulture, Faculty of Forestry, University of Belgrade, Kneza Višeslava 1, 11030 Belgrade, Serbia
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(19), 3846; https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193846
Submission received: 6 September 2021 / Revised: 20 September 2021 / Accepted: 21 September 2021 / Published: 26 September 2021

Abstract

:
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.

Graphical Abstract

1. 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 south-eastern 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–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:
  • How have green spaces changed in Belgrade in the last 3 decades (1991–2019)?
  • What was the relationship between vegetation indices (the NDVI and the NDWI) and climate factors?
  • 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.

2. Materials and Methods

2.1. Study Area and Population

The city of Belgrade is the capital of Serbia and one of the largest cities in south-eastern Europe, with nearly 1.7 million inhabitants and an administrative area of 3234 km2. 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 km2 (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 m2. However, the share and dynamics of green spaces associated with different land uses (aside from forests and public green spaces) are unknown.

2.2. 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.

2.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:
N I R R E D N I R + R E D
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:
N I R S W I R N I R + S W I R
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.

2.4. 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.

2.5. 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).

2.6. 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]:
LST i LSTmin LSTmax LSTmin
where LSTi 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.

2.7. 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:
Δ LST = LST GS 2 LST GS 1
where LSTGS2 is the LST in a growing season in time point 2 and LSTGS1 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:
Δ NDVI = NDVI GS 2 NDVI GS 1
Δ NDWI = NDWI GS 2 NDWI GS 1
where NDVIGS2 and NDWIGS2 are vegetation indices in a growing season in time point 2, and NDVIGS1 and NDWIGS1 are vegetation indices in a growing season in time point 1.

3. Results

3.1. Spatiotemporal Patterns of Spectral Vegetation Indices

For each growing season, the average NDVI and NDWI were calculated and are presented in Figure 4a. During the study period (1991–2019), the NDVI and the NDWI showed a fluctuating trend. The lowest values were recorded for 2000 (end of the first and beginning of the second decade of the study period) and 2011 (end of the second and beginning of the third decade). A low NDWI was also recorded in 2015. The highest values were observed at the end of the third decade, in 2019.
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 NDWI—the 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”.

3.2. 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–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.

3.3. 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 (R2 = 0.70 and R2 = 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).

3.4. 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 (R2 = 0.75, R2 = 0.61 and R2 = 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.
Between change in the NDWI and change in the LST, a significant negative relationship was recorded in the second (2011–2000) and the third decade (2019–2011) (R2 = 0.52 and R2 = 0.74, respectively, Figure 10). The municipalities with the largest decrease (or the smallest increase in the third decade) in the NDWI had the largest increase in the LST.

4. Discussion

4.1. 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, Čukarica 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].

4.2. 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 intra-urban 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.

4.3. 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.
  • A constant decrease in green spaces and elevated LST in the central city zone negatively affects the health and life quality of the inhabitants of these parts of the city. In the past two decades, frequent heat waves have caused increased heat stress in the urban population and have had an especially negative effect on the health of vulnerable groups (the elderly, children, people with diseases of the cardiovascular and respiratory system and those with mental health issues) (Climate Change Adaptation Action Plan and Vulnerability Assessment, City of Belgrade, Secretariat for Environmental Protection, 2015). The proposed measures for mitigation of the urban heat island effect suggested development of new green spaces on built-up land cover that is undergoing conversion, preservation of existing and formation of new tree lines and establishment of green roofs and walls on technically feasible surfaces (the General Urban Plan of Belgrade 2021; Action Plan and Vulnerability Assessment, City of Belgrade, Secretariat for Environmental Protection, 2015). Besides these measures, we suggest prioritizing the preservation of old, large trees in the central city zone.
  • 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].
  • The pressure of urban development on the edges of urban forests within the “inner-ring” 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 large-scale 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.

4.4. 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.

5. 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–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.

Author Contributions

Conceptualization, M.M., S.Č. and M.P.; methodology, M.M., J.C. and M.P.; formal analysis, M.M. and J.C.; investigation, J.C. and M.P.; resources, J.C. and A.T.; writing—original draft preparation, M.M., S.Č., Z.P. and J.T.-D.; writing—review and editing, M.M., Z.P. and J.T.-D.; visualization, M.M., A.T. and S.Č.; supervision, M.P. All authors have read and agreed to the published version of the manuscript.

Funding

M.M. was supported by DAAD funding for “Research Stays for University Academics and Scientists, 2021”, funding program number 57552334.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets generated and/or analyzed during the study are available from the corresponding author upon reasonable request.

Acknowledgments

The authors acknowledge the Ministry of Education, Science and Technological Development, Republic of Serbia, contract numbers 451-03-9/2021-14/ 200007 and 451-03-9/2021-14/200169.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Grimmond, S. Urbanization and global environmental change: Local effects of urban warming. Geogr. J. 2007, 173, 83–88. [Google Scholar] [CrossRef]
  2. Imhoff, M.L.; Zhang, P.; Wolfe, R.E.; Bounoua, L. Remote sensing of the urban heat island effect across biomes in the continental USA. Remote Sens. Environ. 2010, 114, 504–513. [Google Scholar] [CrossRef] [Green Version]
  3. Bolund, P.; Hunhammar, S. Ecosystem services in urban areas. Ecol. Econ. 1999, 29, 293–301. [Google Scholar] [CrossRef]
  4. Bowler, D.E.; Buyung-Ali, L.; Knight, T.M.; Pullin, A.S. Urban greening to cool towns and cities: A systematic review of the empirical evidence. Landsc. Urban Plan. 2010, 97, 147–155. [Google Scholar] [CrossRef]
  5. Gill, S.E.; Handley, J.F.; Ennos, A.R.; Pauleit, S. Adapting Cities for Climate Change: The Role of the Green Infrastructure. Built Environ. 2007, 33, 115–132. [Google Scholar] [CrossRef] [Green Version]
  6. Smid, M.; Russo, S.; Costa, A.C.; Granell, C.; Pebesma, E. Ranking European capitals by exposure to heat waves and cold waves. Urban Clim. 2019, 27, 388–402. [Google Scholar] [CrossRef]
  7. Shih, W.Y.; Ahmad, S.; Chen, Y.C.; Lin, T.P.; Mabon, L. Spatial relationship between land development pattern and intra-urban thermal variations in Taipei. Sustain. Cities Soc. 2020, 62, 102415. [Google Scholar] [CrossRef]
  8. Tran, D.X.; Pla, F.; Latorre-Carmona, P.; Myint, S.W.; Caetano, M.; Kieu, H.V. Characterizing the relationship between land use land cover change and land surface temperature. ISPRS J. Photogramm. Remote Sens. 2017, 124, 119–132. [Google Scholar] [CrossRef] [Green Version]
  9. Nikolov, P. Demographic Factors, Urbanization and Urban Sprawl in Sofia, Rome and Belgrade. In Proceedings of the 14th International Scientific Conference VSU’2014, Sofia, Bulgaria, 5–6 June 2014. [Google Scholar]
  10. Slaev, A.D.; Nedović-Budić, Z.; Krunić, N.; Petrić, J.; Daskalova, D. Suburbanization and sprawl in post-socialist Belgrade and Sofia. Eur. Plan. Stud. 2018, 26, 1389–1412. [Google Scholar] [CrossRef]
  11. Liu, Y.; Wang, Y.; Peng, J.; Du, Y.; Liu, X.; Li, S.; Zhang, D. Correlations between urbanization and vegetation degradation across the world’s metropolises using DMSP/OLS nighttime light data. Remote Sens. 2015, 7, 2067–2088. [Google Scholar] [CrossRef] [Green Version]
  12. Weng, Q. Thermal infrared remote sensing for urban climate and environmental studies: Methods, applications, and trends. ISPRS J. Photogramm. Remote Sens. 2009, 64, 335–344. [Google Scholar] [CrossRef]
  13. Zhou, D.; Xiao, J.; Bonafoni, S.; Berger, C.; Deilami, K.; Zhou, Y.; Frolking, S.; Yao, R.; Qiao, Z.; Sobrino, J.A. Satellite remote sensing of surface urban heat islands: Progress, challenges, and perspectives. Remote Sens. 2019, 11, 48. [Google Scholar] [CrossRef] [Green Version]
  14. Liu, L.; Zhang, Y. Urban heat island analysis using the landsat TM data and ASTER Data: A case study in Hong Kong. Remote Sens. 2011, 3, 1535–1552. [Google Scholar] [CrossRef] [Green Version]
  15. Ahmed, B.; Kamruzzaman, M.D.; Zhu, X.; Shahinoor Rahman, M.D.; Choi, K. Simulating land cover changes and their impacts on land surface temperature in dhaka, bangladesh. Remote Sens. 2013, 5, 5969–5998. [Google Scholar] [CrossRef] [Green Version]
  16. Huang, L.; Shen, H.; Wu, P.; Zhang, L.; Zeng, C. Relationships analysis of land surface temperature with vegetation indicators and impervious surface fraction by fusing multi-temporal and multi-sensor remotely sensed data. In Proceedings of the 2015 Joint Urban Remote Sensing Event (JURSE), Lausanne, Switzerland, 30 March–1 April 2015; Volume 2015. [Google Scholar] [CrossRef]
  17. Kumar, D.; Shekhar, S. Statistical analysis of land surface temperature-vegetation indexes relationship through thermal remote sensing. Ecotoxicol. Environ. Saf. 2015, 121, 39–44. [Google Scholar] [CrossRef]
  18. Luo, X.; Li, W. Scale effect analysis of the relationships between urban heat island and impact factors: Case study in Chongqing. J. Appl. Remote Sens. 2014, 8, 084995. [Google Scholar] [CrossRef]
  19. Fattah, M.A.; Morshed, S.R.; Morshed, S.Y. Multi-layer perceptron-Markov chain-based artificial neural network for modelling future land-specific carbon emission pattern and its influences on surface temperature. SN Appl. Sci. 2021, 3, 1–22. [Google Scholar] [CrossRef]
  20. Karnieli, A.; Agam, N.; Pinker, R.T.; Anderson, M.; Imhoff, M.L.; Gutman, G.G.; Panov, N.; Goldberg, A. Use of NDVI and land surface temperature for drought assessment: Merits and limitations. J. Clim. 2010, 23, 618–633. [Google Scholar] [CrossRef]
  21. Karnieli, A.; Ohana-Levi, N.; Silver, M.; Paz-Kagan, T.; Panov, N.; Varghese, D.; Chrysoulakis, N.; Provenzale, A. Spatial and seasonal patterns in vegetation growth-limiting factors over Europe. Remote Sens. 2019, 11, 2406. [Google Scholar] [CrossRef] [Green Version]
  22. Chi, Y.; Sun, J.; Sun, Y.; Liu, S.; Fu, Z. Multi-temporal characterization of land surface temperature and its relationships with normalized difference vegetation index and soil moisture content in the Yellow River Delta, China. Glob. Ecol. Conserv. 2020, 23, e01092. [Google Scholar] [CrossRef]
  23. Guha, S.; Govil, H.; Besoya, M. An investigation on seasonal variability between LST and NDWI in an urban environment using Landsat satellite data. Geomat. Nat. Hazards Risk 2020, 11, 1319–1345. [Google Scholar] [CrossRef]
  24. Sun, D.; Kafatos, M. Note on the NDVI-LST relationship and the use of temperature-related drought indices over North America. Geophys. Res. Lett. 2007, 34. [Google Scholar] [CrossRef] [Green Version]
  25. Pongrácz, R.; Bartholy, J.; Dezso, Z. Application of remotely sensed thermal information to urban climatology of Central European cities. Phys. Chem. Earth 2010, 35, 95–99. [Google Scholar] [CrossRef]
  26. Obradović-Arsić, D.; Filipović, D. Physical-geographic factors of development of Belgrade. In Challenges of Spatial Development of Ljubljana and Belgrade; Krevs, M., Djordjević, D., Pichler-Milanović, N., Eds.; Scientific Publishing House of the Faculty of Arts: Ljubljana, Slovenia, 2010; pp. 37–54. [Google Scholar]
  27. Justice, C.O.; Townshend, J.R.G.; Holben, A.N.; Tucker, C.J. Analysis of the phenology of global vegetation using meteorological satellite data. Int. J. Remote Sens. 1985, 6, 1271–1318. [Google Scholar] [CrossRef]
  28. Gao, B.C. NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 1996, 58, 257–266. [Google Scholar] [CrossRef]
  29. Marusig, D.; Petruzzellis, F.; Tomasella, M.; Napolitano, R.; Altobelli, A.; Nardini, A. Correlation of field-measured and remotely sensed plant water status as a tool to monitor the risk of drought-induced forest decline. Forests 2020, 11, 77. [Google Scholar] [CrossRef] [Green Version]
  30. Lu, D.; Weng, Q. A survey of image classification methods and techniques for improving classification performance. Int. J. Remote Sens. 2007, 28, 823–870. [Google Scholar] [CrossRef]
  31. Taylor, L.; Hochuli, D.F. Defining greenspace: Multiple uses across multiple disciplines. Landsc. Urban Plan. 2017, 158, 25–38. [Google Scholar] [CrossRef] [Green Version]
  32. van der Linden, S.; Rabe, A.; Held, M.; Jakimow, B.; Leitão, P.; Okujeni, A.; Schwieder, M.; Suess, S.; Hostert, P. The EnMAP-Box—A Toolbox and Application Programming Interface for EnMAP Data Processing. Remote Sens. 2015, 7, 11249–11266. [Google Scholar] [CrossRef] [Green Version]
  33. Elias, I.; Hernandez, R.; Shi, W. A Random Forests classification method for urban land-use mapping integrating spatial metrics and texture analysis. Int. J. Remote Sens. 2018, 39, 1175–1198. [Google Scholar] [CrossRef]
  34. Goldblatt, R.; You, W.; Hanson, G.; Khandelwal, A.K. Detecting the boundaries of urban areas in India: A dataset for pixel-based image classification in google earth engine. Remote Sens. 2016, 8, 634. [Google Scholar] [CrossRef] [Green Version]
  35. Lu, L.; Weng, Q.; Xiao, D.; Guo, H.; Li, Q.; Hui, W. Spatiotemporal variation of surface urban heat islands in relation to land cover composition and configuration: A multi-scale case study of Xi’an, China. Remote Sens. 2020, 12, 2713. [Google Scholar] [CrossRef]
  36. Amiri, R.; Weng, Q.; Alimohammadi, A.; Alavipanah, S.K. Spatial-temporal dynamics of land surface temperature in relation to fractional vegetation cover and land use/cover in the Tabriz urban area, Iran. Remote Sens. Environ. 2009, 113, 2606–2617. [Google Scholar] [CrossRef]
  37. Dobbs, C.; Nitschke, C.; Kendal, D. Assessing the drivers shaping global patterns of urban vegetation landscape structure. Sci. Total Environ. 2017, 592, 171–177. [Google Scholar] [CrossRef]
  38. Threlfall, C.G.; Ossola, A.; Hahs, A.K.; Williams, N.S.G.; Wilson, L.; Livesley, S.J. Variation in vegetation structure and composition across urban green space types. Front. Ecol. Evol. 2016, 4, 66. [Google Scholar] [CrossRef] [Green Version]
  39. Pataki, D.E.; Carreiro, M.M.; Cherrier, J.; Grulke, N.E.; Jennings, V.; Pincetl, S.; Pouyat, R.V.; Whitlow, T.H.; Zipperer, W.C. Coupling biogeochemical cycles in urban environments: Ecosystem services, green solutions, and misconceptions. Front. Ecol. Environ. 2011, 9, 27–36. [Google Scholar] [CrossRef]
  40. Brovkin, V. Climate-vegetation interaction. J. Phys. IV 2002, 12, 57–72. [Google Scholar] [CrossRef]
  41. Holdridge, L.R. Determination of world plant formations from simple climatic data. Science 1947, 105, 367–368. [Google Scholar] [CrossRef]
  42. Hawinkel, P.; Thiery, W.; Lhermitte, S.; Swinnen, E.; Verbist, B.; Van Orshoven, J.; Muys, B. Vegetation response to precipitation variability in East Africa controlled by biogeographical factors. J. Geophys. Res. Biogeosci. 2016, 121, 2422–2444. [Google Scholar] [CrossRef] [Green Version]
  43. Weltzin, J.F.; Loik, M.E.; Schwinning, S.; Williams, D.G.; Fay, P.A.; Haddad, B.M.; Harte, J.; Huxman, T.E.; Knapp, A.K.; Lin, G.; et al. Assessing the Response of Terrestrial Ecosystems to Potential Changes in Precipitation. Bioscience 2003, 53, 941–952. [Google Scholar] [CrossRef]
  44. Nemani, R.R.; Keeling, C.D.; Hashimoto, H.; Jolly, W.M.; Piper, S.C.; Tucker, C.J.; Myneni, R.B.; Running, S.W. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 2003, 300, 1560–1563. [Google Scholar] [CrossRef] [Green Version]
  45. Fensholt, R.; Horion, S.; Tagesson, T.; Ehammer, A.; Grogan, K.; Tian, F.; Huber, S.; Verbesselt, J.; Prince, S.D.; Tucker, C.J.; et al. Assessing drivers of vegetation changes in drylands from time series of earth observation data. In Remote Sensing Time Series. Remote Sensing and Digital Image Processing; Kuenzer, C., Dech, S., Wagner, W., Eds.; Springer: Cham, The Netherlands, 2015; Volume 22, pp. 183–202. [Google Scholar] [CrossRef] [Green Version]
  46. Wramneby, A.; Smith, B.; Samuelsson, P. Hot spots of vegetation-climate feedbacks under future greenhouse forcing in Europe. J. Geophys. Res. Atmos. 2010, 115, 1–12. [Google Scholar] [CrossRef] [Green Version]
  47. Grimmond, C.S.B.; Oke, T.R. Evapotranspiration rates in urban areas. In Impacts of Urban Growth on Surface Water and Groundwater Quality: Proceedings of an International Symposium Held during IUGG 99, the XXII General Assembly of the International Union of Geodesy and Geophysics, at Birmingham, UK 18–30 July 1999; IAHS Press, Institute of Hydrology: Wallingford, UK, 1999; pp. 235–243. [Google Scholar]
  48. Konarska, J.; Uddling, J.; Holmer, B.; Lutz, M.; Lindberg, F.; Pleijel, H.; Thorsson, S. Transpiration of urban trees and its cooling effect in a high latitude city. Int. J. Biometeorol. 2016, 60, 159–172. [Google Scholar] [CrossRef] [PubMed]
  49. Rahman, M.A.; Moser, A.; Gold, A.; Rötzer, T.; Pauleit, S. Vertical air temperature gradients under the shade of two contrasting urban tree species during different types of summer days. Sci. Total Environ. 2018, 633, 100–111. [Google Scholar] [CrossRef] [PubMed]
  50. Rahman, M.A.; Moser, A.; Rötzer, T.; Pauleit, S. Within canopy temperature differences and cooling ability of Tilia cordata trees grown in urban conditions. Build. Environ. 2017, 114, 118–128. [Google Scholar] [CrossRef]
  51. Duncan, J.M.A.; Boruff, B.; Saunders, A.; Sun, Q.; Hurley, J.; Amati, M. Turning down the heat: An enhanced understanding of the relationship between urban vegetation and surface temperature at the city scale. Sci. Total Environ. 2019, 656, 118–128. [Google Scholar] [CrossRef]
  52. Li, X.; Zhou, W.; Ouyang, Z.; Xu, W.; Zheng, H. Spatial pattern of greenspace affects land surface temperature: Evidence from the heavily urbanized Beijing metropolitan area, China. Landsc. Ecol. 2012, 27, 887–898. [Google Scholar] [CrossRef]
  53. Logan, T.M.; Zaitchik, B.; Guikema, S.; Nisbet, A. Night and day: The influence and relative importance of urban characteristics on remotely sensed land surface temperature. Remote Sens. Environ. 2020, 247, 111861. [Google Scholar] [CrossRef]
  54. Wang, Y.; Zhang, Y.; Ding, N.; Qin, K.; Yang, X. Simulating the impact of urban surface evapotranspiration on the urban heat island effect using the modified RS-PM model: A case study of Xuzhou, China. Remote Sens. 2020, 12, 578. [Google Scholar] [CrossRef] [Green Version]
  55. Gunawardena, K.R.; Wells, M.J.; Kershaw, T. Utilising green and bluespace to mitigate urban heat island intensity. Sci. Total Environ. 2017, 584–585, 1040–1055. [Google Scholar] [CrossRef]
  56. Zhang, Q.; Ban, Y. Evaluation of urban expansion and its impact on surface temperature in Beijing, China. In Proceedings of the 2011 Joint Urban Remote Sensing Event, Munich, Germany, 11–13 April 2011. [Google Scholar]
  57. Coutts, A.M.; Harris, R.J.; Phan, T.; Livesley, S.J.; Williams, N.S.G.; Tapper, N.J. Thermal infrared remote sensing of urban heat: Hotspots, vegetation, and an assessment of techniques for use in urban planning. Remote Sens. Environ. 2016, 186, 637–651. [Google Scholar] [CrossRef]
  58. Goldblatt, R.; Deininger, K.; Hanson, G. Utilizing publicly available satellite data for urban research: Mapping built-up land cover and land use in Ho Chi Minh City, Vietnam. Dev. Eng. 2018, 3, 83–99. [Google Scholar] [CrossRef]
  59. Kabisch, N.; Selsam, P.; Kirsten, T.; Lausch, A.; Bumberger, J. A multi-sensor and multi-temporal remote sensing approach to detect land cover change dynamics in heterogeneous urban landscapes. Ecol. Indic. 2019, 99, 273–282. [Google Scholar] [CrossRef]
  60. Zhou, W.; Huang, G.; Cadenasso, M.L. Does spatial configuration matter? Understanding the effects of land cover pattern on land surface temperature in urban landscapes. Landsc. Urban Plan. 2011, 102, 54–63. [Google Scholar] [CrossRef]
  61. Powell, R.L.; Roberts, D.A.; Dennison, P.E.; Hess, L.L. Sub-pixel mapping of urban land cover using multiple endmember spectral mixture analysis: Manaus, Brazil. Remote Sens. Environ. 2007, 106, 253–267. [Google Scholar] [CrossRef]
  62. Quintano, C.; Fernández-Manso, A.; Shimabukuro, Y.E.; Pereira, G. Spectral unmixing. Int. J. Remote Sens. 2012, 33, 5307–5340. [Google Scholar] [CrossRef]
  63. Yu, J.; Chen, D.; Lin, Y.; Ye, S. Comparison of linear and nonlinear spectral unmixing approaches: A case study with multispectral TM imagery. Int. J. Remote Sens. 2017, 38, 773–795. [Google Scholar] [CrossRef]
  64. Yang, C.; Wu, G.; Ding, K.; Shi, T.; Li, Q.; Wang, J. Improving land use/land cover classification by integrating pixel unmixing and decision tree methods. Remote Sens. 2017, 9, 122. [Google Scholar] [CrossRef] [Green Version]
  65. Krtalić, A.; Linardić, D.; Pernar, R. Framework for spatial and temporal monitoring of urban forest and vegetation conditions: Case study Zagreb, Croatia. Sustainability 2021, 13, 6055. [Google Scholar] [CrossRef]
Figure 1. Administrative unit of Belgrade (dashed line) and study area defined by the General Urban Plan of Belgrade 2021 (solid line). Map credit: Google Maps, 2021.
Figure 1. Administrative unit of Belgrade (dashed line) and study area defined by the General Urban Plan of Belgrade 2021 (solid line). Map credit: Google Maps, 2021.
Remotesensing 13 03846 g001
Figure 2. Changes in population density for 12 municipalities, according to censuses from 1991, 2002 and 2011.
Figure 2. Changes in population density for 12 municipalities, according to censuses from 1991, 2002 and 2011.
Remotesensing 13 03846 g002
Figure 3. Methodological workflow of the study.
Figure 3. Methodological workflow of the study.
Remotesensing 13 03846 g003
Figure 4. Temporal changes (a) and seasonal variation (b) in the NDVI and NDWI for Belgrade. * Data for 2006 is based only on one image per growing season that met the criteria of clear atmospheric conditions. ** Data for May is based only on one image during the study period that met the criteria of clear atmospheric conditions.
Figure 4. Temporal changes (a) and seasonal variation (b) in the NDVI and NDWI for Belgrade. * Data for 2006 is based only on one image per growing season that met the criteria of clear atmospheric conditions. ** Data for May is based only on one image during the study period that met the criteria of clear atmospheric conditions.
Remotesensing 13 03846 g004
Figure 5. Municipalities (a) and spatial zones according to the General Urban Plan of Belgrade 2021 (b). Spatiotemporal changes of the annual NDVI and NDWI in 12 municipalities of Belgrade. (c) * Data for 2006 is based only on one image per growing season that met the criteria of clear atmospheric conditions.
Figure 5. Municipalities (a) and spatial zones according to the General Urban Plan of Belgrade 2021 (b). Spatiotemporal changes of the annual NDVI and NDWI in 12 municipalities of Belgrade. (c) * Data for 2006 is based only on one image per growing season that met the criteria of clear atmospheric conditions.
Remotesensing 13 03846 g005
Figure 6. Land cover classification for 1991, 2000, 2011 and 2019.
Figure 6. Land cover classification for 1991, 2000, 2011 and 2019.
Remotesensing 13 03846 g006
Figure 7. Belgrade’s land cover change maps for the last three decades.
Figure 7. Belgrade’s land cover change maps for the last three decades.
Remotesensing 13 03846 g007
Figure 8. Relationship between the NDVI (a) and NDWI (b) and precipitation in the last three decades. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 8. Relationship between the NDVI (a) and NDWI (b) and precipitation in the last three decades. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001.
Remotesensing 13 03846 g008
Figure 9. Relationship between the change in the NDVI and the change in the LST in the last three decades. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 9. Relationship between the change in the NDVI and the change in the LST in the last three decades. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001.
Remotesensing 13 03846 g009
Figure 10. Relationship between the change in the NDWI and the change in the LST in the last three decades. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001.
Figure 10. Relationship between the change in the NDWI and the change in the LST in the last three decades. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001.
Remotesensing 13 03846 g010
Table 1. Remote sensing and climate data used in the study.
Table 1. Remote sensing and climate data used in the study.
Image IDObservation Date
[DD/MM/YY]
SensorLST 1
[°C]
AT 2
[°C]
P 2
[mm]
LT05_L2SP_186029_19910610_20200915_02_T110/06/1991
8:45 a.m.
Landsat 5 TM28.119.2338.8
LT05_L2SP_186029_19910712_20200915_02_T112/07/1991
8:45 a.m.
Landsat 5 TM33.5
LT05_L2SP_186029_19910930_20200915_02_T130/09/1991
8:45 a.m.
Landsat 5 TM26.6
LT05_L2SP_186029_19950723_20200912_02_T123/07/1995
8:26 a.m.
Landsat 5 TM32.819.7357.6
LT05_L2SP_186029_19950808_20200912_02_T108/08/1995
8:25 a.m.
Landsat 5 TM32.7
LT05_L2SP_186029_19950824_20200912_02_T124/08/1995
8:24 a.m.
Landsat 5 TM28.8
LT05_L2SP_186029_20000517_20200907_02_T117/05/2000
8:57 a.m.
Landsat 5 TM37.121.8157.3
LT05_L2SP_186029_20000602_20200907_02_T102/06/2000
8:57 a.m.
Landsat 5 TM31.4
LT05_L2SP_186029_20060721_20200831_02_T121/07/2006
9:14 a.m.
Landsat 5 TM37.120.1353.6
LT05_L2SP_186029_20110601_20200822_02_T101/06/2011
9:11 a.m.
Landsat 5 TM26.421.7230.3
LT05_L2SP_186029_20110719_20200822_02_T119/07/2011
9:10 a.m.
Landsat 5 TM41.6
LT05_L2SP_186029_20110820_20200820_02_T120/08/2011
9:10 a.m.
Landsat 5 TM33.1
LT05_L2SP_186029_20110905_20200820_02_T105/09/2011
9:10 a.m.
Landsat 5 TM37.8
LC08_L2SP_186029_20150612_20200909_02_T112/06/2015
9:20 a.m.
Landsat 8 OLI/TIRS36.322.3247.3
LC08_L2SP_186029_20150815_20200908_02_T115/08/2015
9:21 a.m.
Landsat 8 OLI/TIRS41.9
LC08_L2SP_186029_20150831_20200908_02_T131/08/2015
9:21 a.m.
Landsat 8 OLI/TIRS39.0
LC08_L2SP_186029_20190725_20200827_02_T125/07/2019
9:21 a.m.
Landsat 8 OLI/TIRS36.021.6376.0
LC08_L2SP_186029_20190810_20200827_02_T110/08/2019
9:21 a.m.
Landsat 8 OLI/TIRS36.5
1 Land surface temperature in Belgrade estimated based on Collection 2 Surface Temperature product. 2 Climate data obtained from Republic Hydrometeorological Service of Serbia (2020) based on measurements from three meteorological stations: Observatory, Košutnjak and Surčin airport. AT—average air temperature in growing season; P—sum of precipitation in growing season.
Table 2. Input data for the calculation of NDVI and NDWI.
Table 2. Input data for the calculation of NDVI and NDWI.
SensorREDNIRSWIRSpatial
Resolution
[m]
Landsat 5 TMBand 3Band 4Band 530
Landsat 8 OLI/TIRSBand 4Band 5Band 630
Table 3. Land cover class areas and accuracy assessment for 1991, 2000, 2011 and 2019.
Table 3. Land cover class areas and accuracy assessment for 1991, 2000, 2011 and 2019.
YearClassAreaAccuracy
Absolute [ha]Relative [%]User’sProducer’sOverallKappa
1991Green space39,779.1951.2099.299.7196.950.95
Bare land23,558.6730.3294.6691.38
Water3445.654.4399.3399.31
Built-up10,913.9414.0582.0186.6
2000Green space26,552.4334.1798.2799.093.390.89
Bare land30,079.838.7192.5367.41
Water3308.764.2699.8599.95
Built-up17,756.4622.8567.592.5
2011Green space34,259.2244.0999.5399.6395.100.92
Bare land22,563.7229.0496.2382.18
Water3252.874.1999.9799.99
Built-up17,621.6422.6877.8195.59
2019Green space38,522.9749.5899.7299.7795.430.93
Bare land15,694.5620.2096.7373.85
Water3144.244.0599.8799.96
Built-up20,334.3326.1780.8897.76
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Marković, M.; Cheema, J.; Teofilović, A.; Čepić, S.; Popović, Z.; Tomićević-Dubljević, J.; Pause, M. Monitoring of Spatiotemporal Change of Green Spaces in Relation to the Land Surface Temperature: A Case Study of Belgrade, Serbia. Remote Sens. 2021, 13, 3846. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193846

AMA Style

Marković M, Cheema J, Teofilović A, Čepić S, Popović Z, Tomićević-Dubljević J, Pause M. Monitoring of Spatiotemporal Change of Green Spaces in Relation to the Land Surface Temperature: A Case Study of Belgrade, Serbia. Remote Sensing. 2021; 13(19):3846. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193846

Chicago/Turabian Style

Marković, Milena, Jasmin Cheema, Anica Teofilović, Slavica Čepić, Zorica Popović, Jelena Tomićević-Dubljević, and Marion Pause. 2021. "Monitoring of Spatiotemporal Change of Green Spaces in Relation to the Land Surface Temperature: A Case Study of Belgrade, Serbia" Remote Sensing 13, no. 19: 3846. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13193846

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop