Framework for Spatial and Temporal Monitoring of Urban Forest and Vegetation Conditions: Case Study Zagreb, Croatia

: Urban forest and vegetation conditions are an important variable in urban ecosystem management decision-making. However, it is difﬁcult to evaluate and monitor solely on the basis of ﬁeld measurements. Remote sensing technologies can greatly contribute to the faster extraction and mapping of vegetation health status indicators, on the basis of which agronomy and forestry experts can draw conclusions about the condition of urban vegetation in larger areas. A new remote sensing-based urban forest and vegetation cover monitoring framework is presented and applied to a case study of the city of Zagreb, Croatia. In this study, Sentinel-2 multi-temporal imagery was used to derive and analyze the current state of urban forest cover. Vegetation indices (NDVI, RVI, and GRVI) were calculated. K-means unsupervised classiﬁcation of the vegetation indices was conducted. In this way, the dimensionality of the vegetation indices was reduced, while all the data contained in it were used to represent their graded values. Vegetation that was in a poor condition stood out better that way. Finally, PCA-based change detection was performed on the vegetation indices graded values, and a map of change was produced. These results need to be interpreted and validated by foresters and agronomists in further research.


Introduction
Urban forests and vegetation are of crucial importance for modern towns, because they provide a wide range of ecosystem services that increase citizens' quality of life. They are also an important part of the urban ecological system, which plays an important role in protecting the urban ecological environment [1]. Vegetation in general has a multitude of urban ecosystem functions [2] and provides multiple benefits for city-dwellers. It makes cities more resilient to atmospheric conditions, sudden climate changes, and the occurrence of thermal islands [3][4][5][6]. The overall urban vegetation inside the city borders (individual trees, parks, and forest communities) affects the physical environment of the city. It aids the selective absorption and reflection of incident radiation and also regulates the exchange of permanent and sensible heat [7,8]. Urban vegetation reduces the potential for creating surface urban heat islands, while the removal of vegetation areas contributes to creating them [9,10]. Furthermore, the spatial distribution and density of forest and vegetation cover within city borders is recognized as a major factor influencing numerous biophysical processes in the urban environment [11,12].
In addition to conventional in situ measurement, which involves sampling by field measurement, satellite multispectral imagery and remote sensing methods can provide an effective, fast insight into the current state of urban vegetation [13,14] over a much larger area without physical contact. Identifying and monitoring the state of (urban) vegetation on land surfaces using remote sensing imagery depends on the spectral signatures of vegetation on multispectral images [15,16]. As a result, spectral variations provide fairly precise changes to isolate and display forest depletion areas in the Bosomtwe Range forest reserve. Similarly, [60] used the classical image classification algorithms, K-mean and ISODATA, for unsupervised classification and the maximum likelihood classification for supervised classification for the purpose of vegetation extraction from remote sensing imagery. In [26], a methodological framework for monitoring changes in the extent and canopy density of forest stands was presented. The supervised classification algorithm and normalized difference vegetation index (NDVI) were combined within a GIS environment to quantify the extent and density change of forest cover stands.
The objective of this paper is to present the framework for spatial and temporal monitoring of urban forest conditions, which relies exclusively on remote sensing methods and can be totally independent of in situ measurements and sensor type. The main idea is to monitor relative changes in the state of urban forest cover over time instead of determining annual changes over large areas using complicated methods [61][62][63]. In this way, trends of the situation in a certain part of the year or the phenological calendar can be monitored and highlighted. In this paper, the preliminary results of the framework for spatiotemporal analysis of urban forest distribution and its condition in the city of Zagreb, Croatia are presented.
In order to detect changes in the state of urban forest cover more effectively, a PCAbased change detection of the results of the classifications was performed. Change detection implies recognizing differences in the representation (current state) of an object or phenomenon by observing it at different times [64]. Timely change detection of features on the Earth's surface (in our case, urban vegetation) provides a foundation for understanding better relationships and interactions between humans and natural phenomena, in order to obtain additional data and improve decision-making [65].

Study Area and Data Sources
Zagreb (45 • 49 N, 15 • 59 E) is a central European city located in the central part of the Republic of Croatia on the plain bordered by the River Sava to the South and Mt. Medvednica to the North (Figure 1). Its urban area extends along the plain to Medvednica in the South and its slopes in the North. There are individual trees, avenues, parks (in the lowland areas), and forests (Medvednica). The gently sloping hills of Medvednica are covered in forests, and the flat land between Medvednica and the Sava are covered with urban vegetation (grassed areas, shrubs, single trees, and parks). The climate of Zagreb is classified as oceanic with significant continental influences, closely bordering on both a humid continental climate and a humid subtropical climate, with four distinct seasons. Zagreb summers are warm and wet, winters are cold, springs are mild and pleasant with unpredictable weather changes, and autumns are temperate but rainy and foggy during the latter half of the season (Weather Atlas). The city is a densely populated urban area with 792,875 citizens (according to the last population census in 2011), covering an area of 641.36 km 2 . The extent of the study area covers the limits of the city of Zagreb. and September in those three years (Table 1) were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/ accessed on 22 May 2021). Only images with clear sky views of the study area (no or very limited cloud cover) were included in the analysis. As no cloud-free satellite imagery was found for May 2017 and June in all three years of the survey, the data for those months were completely excluded from the analyses in this paper. The exception was April 2017. There was no cloud-free Sentinel-2 sic at the time, so the one from May was used.

Satellite Data and Preprocessing
Sentinel-2A and 2B (Level-1C, SUHET, 2015) satellite data covering the study area were used to calculate the vegetation indices during the active vegetation period in 2016, 2017, and 2018. The images were acquired at 12 bits in 13 spectral bands at different spatial resolutions. Four were at 10 m (bands: 2 15,15,20,20,90, and 180 nm, respectively), and three bands were at 60 m spatial resolution (bands: 1, 9, and 10, with central wavelengths of 443.0, 945.0, and 1375.0 nm, with bandwidths of 20, 20, and 30 nm, respectively) [66]. The downloaded Sentinel-2 satellite images were in the WGS84 and UTM (Universal Transvers Mercator) projection. The data obtained in April, July, August, and September in those three years (Table 1) were downloaded from the Copernicus Open Access Hub (https://scihub.copernicus.eu/ accessed on 22 May 2021). Only images with clear sky views of the study area (no or very limited cloud cover) were included in the analysis. As no cloud-free satellite imagery was found for May 2017 and June in all three years of the survey, the data for those months were completely excluded from the analyses in this paper. The exception was April 2017. There was no cloud-free Sentinel-2 sic at the time, so the one from May was used.
Band pre-processing procedures were performed: atmospheric correction, sampling of bands with a resolution of 60 m and 20 m to a resolution of 10 m, and final cutting according to the boundaries of the area of interest. Atmospheric correction of Sentinel-2 bands was performed within the Sentinel Application Platform (SNAP) software package, with the SNAP (Sentinel-2) Toolbox plugin, Sen2Cor, for Atmospheric Correction, using the associated Level-1C product as input. They were atmospherically corrected from Top-   with the SNAP (Sentinel-2) Toolbox plugin, Sen2Cor, for Atmospheric Correction, using the associated Level-1C product as input. They were atmospherically corrected from Top-Of-Atmosphere (TOA) Level-1C to generate Level-2A (version SNAP 6.0).

Methods
The methodological framework for spatial and temporal monitoring of changes in urban forest cover is presented in Figure 2. The workflow in this framework consisted of the following steps ( Figure 2

Vegetation Indices
Spectral (vegetation) indices are mathematical combinations of different spectral bands of the electromagnetic spectrum, and they are more sensitive than individual bands to vegetation parameters [67][68][69]. They are measures of vegetation activity and are designed to find a functional relationship between vegetation characteristics and remote sensing data [21,27,70]. The main purpose of spectral indices is to enhance the information contained in spectral reflectance data by extracting the variability due to vegetation characteristics, and to minimize soil and atmospheric effects [71]. In this way, the useful 'satellite' indicators, which help in determining vegetation health, can be extracted from Sentinel-2 imagery (as in this study). Vegetation indices measure vegetation activity and show the seasonal and spatial variations of green leaves. The main objective of vegetation indices is to enhance and extract the information contained in spectral reflectance remotely- After calculating the vegetation indices, a supervised classification was carried out for the purpose of separating the urban vegetation cover and the urban non-vegetation cover. Based on the classification results, a mask was made (the non-vegetation cover was masked), and only the urban vegetation cover was used for further consideration ( Figure 3). Furthermore, the unsupervised classification of vegetation indices (in this case, they are NDVI, Ratio Vegetation Index (RVI), Green-Red Vegetation Index (GRVI)) was used to gradate their values for the purpose of better insight into their dynamics and position on the scene. Finally, the PCA method was implemented on the results of the unsupervised classification of vegetation indices, which further reduced the dimensionality of the data in order to focus attention on the essentials. In this case, this was on problematic areas with the lowest values of vegetation indices. The results were displayed on maps of changes.
Sustainability 2021, 13, x FOR PEER REVIEW 9 of 23 and for each date (a total of 12 classifications). The accuracy of vegetation cover classification ranged from 95 to 97%. Selected time series Sentinel-2 imagery was then clipped to the city administrative boundary ( Figure 3). The urban vegetation cover was extracted by merging all vegetation classes in one class, and all other areas were merged in a second class (urban non-vegetation) with no data value attached. All further processing was carried out only on the surface covered by urban vegetation. The results of the classifications were averaged over the years and showed the decrease in the urban vegetation and increase in the non-vegetation cover in the city of Zagreb from year to year (Table 3).

Vegetation Indices Analysis
The value of NDVI ranged from −1 to 1; however, because the value ranges for RVI and GRVI depended on the pixel values in the channels used, there were no fixed ranges.
The NDVI ( Figure 4, Table 4) scale varied within a range from 0.183 to 0.938, where the standard deviations were of the smallest scale in comparison to the three indices used; therefore, the smallest variable of the vegetative coverage at the scene is shown. Although the NDVI detected water bodies and bare soil (non-vegetation areas), it did not distinguish forest and domestic vegetation types, so it was not suitable for classifying areas with dense green vegetation. By analyzing the averaged mean values (Figure 5a,b), it was confirmed that the NDVI had very similar seasonal variations along with the middle values of the RVI, because both depend on the values of the same spectral channels, i.e., red and near-infrared. By observing the mean values, it was noticeable that they continuously decreased in 2017 (May to September), whereas in 2016, the value was almost constant until August, when it radically decreased. In 2018, there was a significant decline in the value between July and August ( Figure 5c). Scrutinizing the absolute minimum, variation on a much bigger scale was noticeable, which points to the enlarged capability of the NDVI

Vegetation Indices
Spectral (vegetation) indices are mathematical combinations of different spectral bands of the electromagnetic spectrum, and they are more sensitive than individual bands to vegetation parameters [67][68][69]. They are measures of vegetation activity and are designed to find a functional relationship between vegetation characteristics and remote sensing data [21,27,70]. The main purpose of spectral indices is to enhance the information contained in spectral reflectance data by extracting the variability due to vegetation characteristics, and to minimize soil and atmospheric effects [71]. In this way, the useful 'satellite' indicators, which help in determining vegetation health, can be extracted from Sentinel-2 imagery (as in this study). Vegetation indices measure vegetation activity and show the seasonal and spatial variations of green leaves. The main objective of vegetation indices is to enhance and extract the information contained in spectral reflectance remotely-sensed data, by extracting the variability due to vegetation characteristics [71]. Therefore, they are suitable for detecting within-field spatial variability [70].
In this study, three vegetation indices in the process of monitoring the urban forest condition were used: NDVI, RVI, and GRVI. NDVI [40,72] is the most common vegetation index, and it represents the ratio of the NIR and red band ( Table 2). It is used for calculation because the absorption by chlorophyll of these two bands of the electromagnetic spectrum is the highest [73]. NDVI was calculated using the (Sentinel-2) red (B4) and near-infrared (B8) bands. A common practice in remote sensing is the use of band ratios to eliminate various albedo effects. If RVI is used [40,74] (Table 2), the vegetation isoline converges at the origin. It is calculated using the (Sentinel-2) red (B4) and near-infrared (B8) bands. GRVI [75] (Table 2) was chosen because it is sensitive to photosynthetic rates in forest canopies, as green and red reflectance are strongly influenced by changes in leaf pigments. It was calculated using the (Sentinel-2) green (B3) and near-infrared (B8) channels.

Urban Vegetation Supervised Classification and Mask
Urban vegetation cover classification was conducted using Sentinel-2 imagery from April, July, August, and September 2016, 2017, and 2018 using the maximum likelihood classification algorithm in QGIS with the semi-automatic classification plugin. The maximum likelihood [76] algorithm calculates the probability distributions for the classes, related to Bayes' theorem, estimating if a pixel belongs to a land cover class. In particular, the probability distributions for the classes assumed the form of multivariate normal models [77]. The sole purpose of this classification was to determine the area covered by vegetation within the administrative boundaries of the city of Zagreb.

Unsupervised Classification of Vegetation Indices
Unsupervised classification [77] attempts to find clusters in an n-dimensional feature space and assigns those clusters to a group. Samples are not required, and it is an easy way of segmenting and understanding an image. In this case, representations of the distribution of the values of the vegetation indexes for the selected dates were made (Table 1) based on K-means unsupervised classification [78]. K-means unsupervised classification calculates initial class means evenly distributed in the data (feature) space, then iteratively clusters the pixels into the nearest class using a minimum distance technique [79]. This algorithm aims at minimizing a squared error function [78], and it can be expressed by Equation (1) [78]: where:

Principal Components Analysis-Based Change Detection
Principal components analysis (PCA)-based urban forest change detection and analysis using the results of unsupervised classification of vegetation indices derived from multitemporal satellite data was performed and is described in this paper. The PCA method has been commonly used for many years to determine the basic dimension of satellite images [80][81][82], land-cover change detection [83,84], or multitemporal satellite data [81,85,86]. Due to its simplicity and capability of enhancing the information on change, it has become one of the most popular change detection techniques [87,88]. PCA is used when it is assumed that multitemporal data are highly related and that information about changes can be highlighted in new components [65]. Compared to projecting data onto the original coordinates, projecting onto the PCs has the advantage of allowing it to detect changes in data correlation which cannot be detected in the original individual variables. Furthermore, this can ensure that any changes to the source variables will be reflected in the PC projections [82]. In the case described in this paper, the input data for PCA were the results of unsupervised classifications of three (GRVI, NDVI, RVI) vegetation indices for the months of May, July, August, and September 2017. The main idea of using PCA was to identify a set of orthogonal principal components from the original feature space, which would better represent the data in the new space. In other words, changes in the different health conditions of urban forest areas were highlighted. This provided the basic prerequisite for detecting changes in the area of interest: the availability of imagery from two or more different dates [85].

Results and Discussion
The main idea of this framework is to fully rely on remote sensing methods for the purpose of monitoring urban forest cover, the dynamics of changes in it, and detecting areas with the lowest values of vegetation indices. In this way, trends in the condition of this cover can be followed over time. Compared to the in situ trait approach, remote sensing sensors are unable to detect all traits and trait variations [89] except for some, such as biochemical-biophysical, physiognomic, morphological, structural, and functional traits at all levels of biotic organization, based on the principles of image spectroscopy across the electro-magnetic spectrum [90,91]. Remote sensing sensors can record the status, stress, disturbances, or resource limitations of vegetation over the short and long term on local and global scales [91,92]. These facts were used in creating this framework. The relative relations of the condition of the same species in the same time period are the main indicators for the work of experts in the field of agronomy, forestry, and urban planning. Namely, noticing bad trends in the state of urban forests is a trigger for experts in terms of taking action to prevent such a trend. After noticing such trends, they can go out into the field to perform in situ measurements and additional analyzes to gain insight into the absolute state of vegetation cover and take actions to mitigate or eliminate the problem.
The paper presents the unsupervised classification of vegetation indices as a novel way of presenting the values that can serve as an indicator for assessing vegetation health status. Since each of the spectral indices presents vegetation in a different way, it might be described as a kind of data fusion for improving their usage and information. This paper focuses on representations of remote sensing methods and their preliminary analysis. Forestry experts' interpretations of the results were not included in this review and are necessary in further research.

Supervised Classification Outputs
Class samples were visually selected urban vegetation (all types: grassed areas, shrubs, trees, agricultural land) and urban non-vegetation polygons (built-up areas, bare land, and water). The pixels inside the sample polygons were used to perform the classification of each Sentinel-2 set of bands (only bands 2, 3, 4, 5, 6, 7, 8, and 8a, with 10 and 20 m spatial resolution and all vegetation indices were used for supervised classification) and for each date (a total of 12 classifications). The accuracy of vegetation cover classification ranged from 95 to 97%. Selected time series Sentinel-2 imagery was then clipped to the city administrative boundary (Figure 3). The urban vegetation cover was extracted by merging all vegetation classes in one class, and all other areas were merged in a second class (urban non-vegetation) with no data value attached. All further processing was carried out only on the surface covered by urban vegetation.
The results of the classifications were averaged over the years and showed the decrease in the urban vegetation and increase in the non-vegetation cover in the city of Zagreb from year to year (Table 3).

Vegetation Indices Analysis
The value of NDVI ranged from −1 to 1; however, because the value ranges for RVI and GRVI depended on the pixel values in the channels used, there were no fixed ranges.
The NDVI (Figure 4, Table 4) scale varied within a range from 0.183 to 0.938, where the standard deviations were of the smallest scale in comparison to the three indices used; therefore, the smallest variable of the vegetative coverage at the scene is shown. Although the NDVI detected water bodies and bare soil (non-vegetation areas), it did not distinguish forest and domestic vegetation types, so it was not suitable for classifying areas with dense green vegetation. By analyzing the averaged mean values (Figure 5a,b), it was confirmed that the NDVI had very similar seasonal variations along with the middle values of the RVI, because both depend on the values of the same spectral channels, i.e., red and near-infrared. By observing the mean values, it was noticeable that they continuously decreased in 2017 (May to September), whereas in 2016, the value was almost constant until August, when it radically decreased. In 2018, there was a significant decline in the value between July and August (Figure 5c). Scrutinizing the absolute minimum, variation on a much bigger scale was noticeable, which points to the enlarged capability of the NDVI regarding the recognition and gradation of poor vegetation. In contrast, analyzing the maximum showed much less variation, which points to the familiar trait of satiation of the NDVI, i.e., its reduced ability to extract details from areas of thriving vegetation.
The span of RVI values ( Figure 4, Table 4 . April, on average, represented the period of local maximum vegetation. As the months of May and June were not included in this survey, the exact appearance of the average curve could not be determined for the indicated period. However, there was a slight decrease in vegetation intensity between April and July. The annual decline in mean RVI began in August and was further compounded by the climatological arrival of autumn in September. The value range of GRVI ( Figure 4, Table 4) was slightly smaller than for RVI (1.459-21.262). The middle and extreme GRVI values were generally lower than the corresponding RVI values. This was a consequence of the fact that the leaves of healthy vegetation reflect green radiation more strongly than red. Accordingly, the ratio between reflected infrared and green radiation was lower than the ratio of reflected infrared and red radiation. With this index, mean values increased continuously from April to August 2016 and then rapidly decreased. The same was true in 2017, with growth from May to July and then a continuous decrease until September. In 2018, the decreases occurred in steps (Figure 5e). At the maximum, greater variation was observed than at the minimum. The GRVI effectively diagnosed the distributions of agricultural vegetation and forest areas, but often failed to discern cropland (at the edges of the city's southern borders). A more detailed analysis was conducted at two locations of special importance for the residents of the City of Zagreb: Maksimir Park (Figure 1, ocher dot) and Dotrščina forest (Figure 1, yellow dot). Maksimir Park is the most popular, oldest park in Zagreb. It is located in the Eastern part, not far from the city center. In addition to five lakes, it contains many plant species, among which are century-old oak trees, which were necessarily prominent. Maksimir Park is known as the 'lungs of the city' because it produces oxygen and acts as a natural air purifier. The dominant deciduous species in Maksimir Park are pedunculate oak (Quercus robur), sessile oak (Quercus petraea), and chestnut (Castanea sativa). The dominant coniferous species are white and black pine (Pinus sylvestris and Pinus nigra) and spruce (Picea excelsa) within the oak forests. Dotršćina forest is located north of Maksimir Park on the slopes of Mt. Medvednica. Thus, the analysis focused on predominantly forest vegetation in the lowland, densely populated city center and the hilly, sparsely populated part of the city.   The objects of the analysis were the vegetation indices values of GRVI, NDVI, and RVI for the forest areas in Maksimir Park and Dotrščina forest for May, July, August, and September 2017 ( Figure 6). The hypothesis was that RVI has the greatest variability in the presentation of vegetation conditions (largest standard deviations) and that NDVI has the smallest, and this was confirmed in these areas ( Figure 6). The highest values of GRVI were measured in July in Maksimir Park, and in May in Dotrščina forest. The highest values of RVI and NDVI were measured in May 2017 in both areas ( Figure 6).
These and similar analyzes can be found in the published articles mentioned in the Introduction. In contrast to such analyzes, in this framework, these results were used to produce new results, with existing tools (unsupervised classification and PCA method), which would indicate potential problems in urban vegetation cover. Unsupervised classification of vegetation indices provides insight into the spatial distribution of their values and conducts their gradation.  A more detailed analysis was conducted at two locations of special importance for the residents of the City of Zagreb: Maksimir Park (Figure 1, ocher dot) and Dotrščina forest (Figure 1, yellow dot). Maksimir Park is the most popular, oldest park in Zagreb. It is located in the Eastern part, not far from the city center. In addition to five lakes, it contains many plant species, among which are century-old oak trees, which were necessarily prominent. Maksimir Park is known as the 'lungs of the city' because it produces oxygen and acts as a natural air purifier. The dominant deciduous species in Maksimir Park are pedunculate oak (Quercus robur), sessile oak (Quercus petraea), and chestnut (Castanea sativa). The dominant coniferous species are white and black pine (Pinus sylvestris and Pinus nigra) and spruce (Picea excelsa) within the oak forests. Dotršćina forest is located

Unsupervised Classification Outputs
The input data for the unsupervised classification were vegetation index sets (NDVI, RVI, and GRVI) for each of the dates listed in Table 1 (Figure 6). The hypothesis was that RVI has the greatest variability in the presentation of vegetation conditions (largest standard deviations) and that NDVI has the smallest, and this was confirmed in these areas ( Figure 6). The highest values of GRVI were measured in July in Maksimir Park, and in May in Dotrščina forest. The highest values of RVI and NDVI were measured in May 2017 in both areas ( Figure 6). These and similar analyzes can be found in the published articles mentioned in the Introduction. In contrast to such analyzes, in this framework, these results were used to produce new results, with existing tools (unsupervised classification and PCA method), which would indicate potential problems in urban vegetation cover. Unsupervised classification of vegetation indices provides insight into the spatial distribution of their values and conducts their gradation.

Unsupervised Classification Outputs
The input data for the unsupervised classification were vegetation index sets (NDVI, RVI, and GRVI) for each of the dates listed in Table 1. A total of 12 classifications were conducted. The urban vegetation cover in the city of Zagreb was graded into four classes, whereby the class indicated by the lightest green (vegetation 1) represented the maximum  Figure 9). Furthermore, at the Dotrščina location, significant values (apart from in May) classified into class 1 (maximum index values) appeared in each observed month, while those at the Maksimir location were insignificant, with the exception of May (Table 5). The trend of the mean movement was very similar, and the standard deviations were inversely proportional to the mean values for both locations ( Figure 10). The standard deviations were similar in May and September, and differed significantly in July and August. While the standard deviations were uniform for the Dotrščina location, the standard deviations for the Maksimir location were significantly higher in July and August. Therefore, greater variability of the vegetation cover condition was detected in the Maksimir location ( Figure 10).  (Figure 7). The unsupervised classification procedure was performed separately for each set of vegetation indices respectively. Therefore, the four classes within the vegetation were not absolutely defined, but rather represented the relative gradation of the index values for each individual set.

Principal Components Analysis of Unsupervised Classification Outputs
PCA was used to extract useful information from unsupervised classifications of vegetation indices and detect temporal changes in the current state of forest cover and creation of map changes. The principle of extracting information in this way is described for Maksimir Park and Dotrščina forest for May, July, August, and September 2017.
After transformation of the original images (unsupervised classifications), four principal components of each area were created. The first component showed the largest possible differences in variance (Table 6). However, the other components also showed significant percentages of variance ( Table 6). The vegetation indices value change map derived by PCA is shown in Figure 11. In the change maps, changes of vegetation indices values were detected and identified at the Maksimir and Dotrščina locations. These maps were created by merging the first three main components in an RGB display: R = PC3, G = PC2, B = PC1. Thus, the changes that occurred from May to September 2017 were highlighted ( Figure 11). Areas where changes occurred were marked from light pink to darker shades in both cases ( Figure 11). The observed changes are related, among other things, to significant damage to oak trees. Namely, during 2017, oak lace bug (Corythucha arcuata) spread throughout the entire area pedunculate oak and sessile oak ( Figure 11).
The data processed and presented by the maps of changes are intended for experts in forestry, agronomy, and urbanism for the purpose of monitoring the state of change in the urban forest/vegetation cover and making adequate decisions in the management of these fragile urban resources.

Principal Components Analysis of Unsupervised Classification Outputs
PCA was used to extract useful information from unsupervised classifications of vegetation indices and detect temporal changes in the current state of forest cover and creation of map changes. The principle of extracting information in this way is described for Maksimir Park and Dotrščina forest for May, July, August, and September 2017.
After transformation of the original images (unsupervised classifications), four principal components of each area were created. The first component showed the largest possible differences in variance (Table 6). However, the other components also showed significant percentages of variance ( Table 6). The vegetation indices value change map derived by PCA is shown in Figure 11. In the change maps, changes of vegetation indices values were detected and identified at the Maksimir and Dotrščina locations. These maps were created by merging the first three main components in an RGB display: R = PC3, G = PC2, B = PC1. Thus, the changes that occurred from May to September 2017 were highlighted ( Figure 11). Areas where changes occurred were marked from light pink to darker shades in both cases ( Figure 11). The observed changes are related, among other things, to significant damage to oak trees. Namely, during 2017, oak lace bug (Corythucha arcuata) spread throughout the entire area pedunculate oak and sessile oak ( Figure 11).

Conclusions
Urban and forest health assessments are very important in the formulation of ecosystem management strategies as well as managing the quality of life in cities. However, the lack of rapid, cost-effective, and quantitative analyses of the health of urban forests and vegetation hinders timely insights into the relative health status of urban vegetation. It was the motivation to create this framework entirely based on remote sensing. The described framework is a closed cost-effective system that uses free satellite images (funded by European states and free delivered by the European agency), free software (open source), and provides insight into the relative trends of the value of vegetation indices and The data processed and presented by the maps of changes are intended for experts in forestry, agronomy, and urbanism for the purpose of monitoring the state of change in the urban forest/vegetation cover and making adequate decisions in the management of these fragile urban resources.

Conclusions
Urban and forest health assessments are very important in the formulation of ecosystem management strategies as well as managing the quality of life in cities. However, the lack of rapid, cost-effective, and quantitative analyses of the health of urban forests and vegetation hinders timely insights into the relative health status of urban vegetation. It was the motivation to create this framework entirely based on remote sensing. The described framework is a closed cost-effective system that uses free satellite images (funded by European states and free delivered by the European agency), free software (open source), and provides insight into the relative trends of the value of vegetation indices and shows their spatial distribution. In this way, spatial and temporal monitoring of changes within the urban forest and vegetation cover was performed, and insight into the trends of these changes was gained. Furthermore, the areas with the highest and lowest values of vegetation indices were also highlighted. The results of these analyses and assessments of these changes need to be delivered to the experts, and they can determine their impact on environmental and human health and take action based on them. Thus, the results of the methodology within this framework are the input data for decision support systems for the purpose of urban forest and vegetation management.
In this case study, Sentinel-2 time series data were used to derive combined vegetation indices for an urban forest and vegetation analysis. Variations in the values of three vegetation indices, NDVI, RVI, and GRVI, in the months of April, July, August, and September from 2016 to 2018, were analyzed. An analysis of the mean values revealed that the vegetation condition depended on the month of observation. All the vegetation indices had the highest mean values in April and increased from year to year. The RVI had the highest perennial average in April (11.4), and during July (10.7), August (9.2) and September (7.4), it recorded a gradual fall in mean values. The NDVI behaved similarly to the RVI, with mean values peaking in April (0.807) and declining in July (0.794), August (0.764), and September (0.735). The GRVI took into account the records of the green spectral channel and, therefore, provided a somewhat different insight into the annual course of vegetation. The GRVI perennial average rose between April (7.5) and July (8.1) and fell during August (7.1) and September (6.0). These data are numerical and represent the result of measurements, but they are still estimates because they were obtained from remote sensing data, without in situ measurements that would confirm or refute them. Therefore, visualizing these assessments gives a better insight into the situation on the scene and indicates problem areas. Temporal data, on the other hand, provide insight into the trends of these problems. In this way, changes in vegetation cover can be more easily detected and visualized. After that, the attention of experts can be focused on problem areas where in situ measurements can be performed and decisions can be made based on them.
The analysis of the unsupervised classification results of vegetation indices showed that the condition of urban vegetation in the city of Zagreb depended on vegetation type, altitude, and the location of the vegetation within the city. The values of 36 spectral indices were combined to yield 12 new results by unsupervised classification. Thus, the dimensionality of the input data was reduced, and all the data contained in the spectral indices were used to display their gradation values. These reduced data were further analyzed by the PCA method in order to determine the trend of changes in the urban vegetation cover. Therefore, this research presented a remote sensing-based monitoring source with the potential to analyze continuous changes in urban forests and vegetation, and assist expert decision-makers in agronomy and forestry to develop ecosystem services management plans. This framework of spectral indices provides indicators for the rapid, cost-effective, and quantitative assessment of forest health. For the operational assessment of the proposed framework of combining spectral indices by unsupervised classification, for the purpose of obtaining indicators for assessing the condition and health of urban forests and vegetation, the results obtained should be compared with reference field data, and their interpretation should involve agronomy and forestry. Further research should be directed in the direction of using other vegetation indices and professional (forester, agronomist) interpretation of the obtained results.