Characterization and Mapping of Public and Private Green Areas in the Municipality of Forl ì (NE Italy) Using High-Resolution Images

: Urban Green Spaces (UGS) contribute to the sustainable development of the urban ecosystem, positively impacting quality of life and providing ecosystem services and social beneﬁts to inhabitants. For urban planning, mapping and quantiﬁcation of UGS become crucial. So far, the contribution of private green spaces to ecosystem services in urban areas has yet to be studied. At the same time, in many Italian cities, they represent a considerable part of the urban green cover. This study utilises a methodological approach and provides insights into the contribution of urban public and private green spaces by the consideration of a case study area in Northeast Italy. To achieve this goal, the main steps were: (i) NDVI extraction from very high-resolution (20 cm) orthophotos, (ii) classiﬁcation of property status and (iii) analysis of the degree of the greenness of land cover units. From our results, the total amount of the green spaces is 5.70 km 2 , of which 72.1% (4.11 km 2 ) is private, and 28.9% (1.59 km 2 ) is public. As for the land cover, three NDVI classes were identiﬁed, highlighting different degrees of homogeneity in NDVI reﬂectance response within each urban land cover unit. These results will support the planning of new green areas in the post-epidemic National Recovery and Resilience Plan.


Urban Green Spaces at the Local and European Framework
The climate crisis poses new challenges for the resilience of cities both in Europe and worldwide. In the fragile Italian urban environment, often characterised by the historical development of compact cities with high population density, it is crucial to identify sustainable solutions to complement traditional urban planning (UP) practices. At this level, a thorough knowledge of the urban territory and its resources, such as urban green spaces (UGS), is necessary to provide reliable working tools to support UP. Active collaboration between research institutes, universities and stakeholders, such as urban planners, could resolve the aforementioned issue [1] and has been encouraged by the scientific community, especially in urban green space analysis and quality of knowledge in the Mediterranean areas [2]. There is a need to implement procedures that could be replicated from the local to regional and national scale through the use of data within public administration databases and geoportals, such as cadastral and topographic datasets, and orthophotos or any highly detailed imagery product [3,4].
The importance of urban green infrastructures (UGI) is quantitatively defined [5] and recognised [6] also in terms of perception and adequate provision of ecosystem services (ES) [2,7]. This function was also recognised by the European Union (EU) in 2013, with its

Characterisation of Urban Green Spaces for Ecosystem Services Provision
Characteristics of UGSs, such as their abundance, spatial distribution, and species composition, play an essential role in UP and public health [14][15][16]. According to the World Health Organization [17], there is evidence of the beneficial effects of UGSs as a consequence of psychological relaxation and stress alleviation as well as the reduction in exposure to air pollutants. The concomitant COVID-19 pandemic situation of recent years, in particular during lockdowns, highlighted the importance of green areas inside and around cities [18], particularly concerning their socio-health and cultural and recreational value [19,20], along with the urgent need for managing these areas following integrated criteria and principles towards urban resilience [21]. Indeed, UGSs, such as parks, playgrounds, or vegetative systems in public and private spaces also provide a comprehensive source of ecosystem services, the so-called Urban Ecosystem Services (UES) [22]. These are chiefly provisioning and regulating services [23,24] and include: (i) carbon stock [25,26] and greenhouse gas (GHG) sequestration in vegetation biomass and soil, together with the mitigation of air pollution [27,28]; (ii) food production with urban agriculture [29,30]; (iii) microclimate thermoregulation by mitigation of urban heat islands [31,32]; (iv) biodiversity conservation [33]; and (v) flood regulation from the impact of extreme weather [34][35][36]. Since private green areas, such as gardens, agricultural fields and private parks are a substantial proportion of UGSs, they are essential for maintaining UES [7] and biodiversity [37]. Nevertheless, little importance is given to private green spaces, both at the management and governance levels [38]. The governance should define rules and decision-making processes for private stakeholders to support individual interests and collective benefits. What is needed regarding UGSs management is the inclusion of private green spaces and the acquisition of an integrated ecosystemic approach in UP, particularly the recognition of the city as an urban ecosystem in which people benefit from nature, either public or private [39]. Furthermore, some authors have investigated cultural UES [23,24], focusing on the differences between public and private green spaces. Analysing consumer preferences and local policies in some Western European countries, they emphasised that these green spaces are not simply substitutes for each other in terms of the living environment within an urbanized area [40]. In addition, as noted by [41], it appears that, at present, there are no valid national or European policy instruments aimed at compensating for the loss of private green space at the urban level by public green space provision. For the time being, there is only a suggestion to use models to study the welfare effects of applying for ecological compensation in cases where urban green spaces are at risk of being exploited [42].
Remote sensing and analytical techniques substantially contribute to the characterization and mapping of UGSs [14,43]. Based on recent advances such as high-resolution imagery and open access data policies [3,44], the production of detailed vegetation maps has become a crucial aid for assisting city planners in designing methods of optimising UES and climate change adaptation strategies [45,46]. Studies have demonstrated that urban vegetation dynamics within and across cities can be assessed through remote sensing by using satellite image analysis and proximal sensing in order for high-resolution assessment to be achieved [47]. For instance, the Normalized Difference Vegetation Index (NDVI) is the main indicator of photosynthetic activity [48], frequently employed in urban landscape ecology [49][50][51][52][53]. Such tools, connected to the internet of things, could become efficient monitoring tools in the future, especially since they are promoted by the guidelines in the Italian and European Recovery Plan as equipment for the Member States for monitoring and predicting possible hazards, promptly assessing their immediate impacts on systems (natural and infrastructures), and consequently defining optimal responses [54].
Novelty elements of this study include the detection and the quantification of the spatial degree of heterogeneity of urban green systems as a significant step for the ecological understanding of cities [55], also identifying areas of possible intervention. The flexibility of NDVI for this purpose is investigated as a good indicator describing land use and cover heterogeneity in urban landscapes [56].

Geographical Framework of the Study Area
Forlì is an Italian municipality located in Forlì-Cesena Province within the Po Valley in Emilia-Romagna Region (Northern Italy) at 34 m a.s.l. (Lat. 44.2227398, Long. 12.0407312), 5 km away from the foothills of the Tuscan-Romagnolo Pre-Apennines and about 26 km from the Adriatic coast. Its area is about 228.2 km 2 , of which 32.88 km is urbanised. The municipality counts 117,138 inhabitants [57] and a population density of 512.4 inhabitants km −2 . The territory is divided into 21 districts, grouped as eight "Comitati Territoriali" (Territorial Committees). The present study focuses on the six core districts encompassing the historic urban centre (q6) and its surroundings (q1-q5, Figure 1) with a total area of 13.6 km 2 , representing the most urbanised, and the more densely populated districts in the entire municipal territory.
Forlì has a warm temperate climate, and is stably humid, with a hot summer (Köppen-Geiger Cfa classification). The coldest month (January) has an average temperature of 3.9 • C (1991-2020 period), while the hottest month (July) has an average temperature of 24.6 • C (1991-2020 period). The average annual cumulative precipitation is around 769 mm for the same reference period [58].

Aims of the Study
The municipality of Forlì used the EU-funded Life project SaveOurSoils4Life (2015-2019, https://www.sos4life.it/en/, accessed on 10 January 2022) focused on recording the ecosystem services of urban soils, within which this research aimed to fill some knowledge gaps in the provision of UES in the framework of the Forlì Green Space Plan. Forlì has a warm temperate climate, and is stably humid, with a hot summer (Köppen-Geiger Cfa classification). The coldest month (January) has an average temperature of 3.9 °C (1991-2020 period), while the hottest month (July) has an average temperature of 24.6 °C (1991-2020 period). The average annual cumulative precipitation is around 769 mm for the same reference period [58].

Aims of the Study
The municipality of Forlì used the EU-funded Life project SaveOurSoils4Life (2015-2019, https://www.sos4life.it/en/, accessed on 10 January 2022) focused on recording the ecosystem services of urban soils, within which this research aimed to fill some knowledge gaps in the provision of UES in the framework of the Forlì Green Space Plan.
The first aim of this study is to quantify and classify the contribution of private and public green spaces to understand the structure of UGI in the Municipality of Forlì better, resorting to very high-resolution images. Second, we assessed the degree of greenness and the heterogeneity of vegetation cover in specific land use and cover classes regarding NDVI spectral response, contributing further to the characterisation of municipal green spaces in Forlì.

Urban Green Spaces Detection and Classification
To achieve the aim of the study, spatial analyses were based on remotely-sensed data and GIS modelling by using different sources: the land use and land cover (LULC) data at the regional level, very high spatial resolution orthophotos (0.2 m pixel size) in visible spectral range (Red, Green and Blue, RGB) and near infrared (NIR), integrated with the The first aim of this study is to quantify and classify the contribution of private and public green spaces to understand the structure of UGI in the Municipality of Forlì better, resorting to very high-resolution images. Second, we assessed the degree of greenness and the heterogeneity of vegetation cover in specific land use and cover classes regarding NDVI spectral response, contributing further to the characterisation of municipal green spaces in Forlì.

Urban Green Spaces Detection and Classification
To achieve the aim of the study, spatial analyses were based on remotely-sensed data and GIS modelling by using different sources: the land use and land cover (LULC) data at the regional level, very high spatial resolution orthophotos (0.2 m pixel size) in visible spectral range (Red, Green and Blue, RGB) and near infrared (NIR), integrated with the cadastral maps for the Municipality of Forlì and the tree census as described here below (Table 1).

Land Use and Land Cover Database and Thematic Maps
Since the 1970s, LULC has been one of the most requested and used regional geographic databases by local governments and professionals. To fully respond to the demands being made by urban planners, especially regarding soil sealing and other land use changes, the Emilia-Romagna Regional Land Use Working Group created the most recent LULC database, as used in this study, with new features compared to the previous 2008 and 2014 releases. The most updated product is the 2017 Land Cover Database (edition 2020, Figure 2), covering the study area. It was derived via photo-interpretation of the Italian Agriculture Remote Sensing Consortium (TeA) 2017 orthophotos (reference scale 1:10,000, position accuracy 5 m) with 0.2 m spatial resolution, which were acquired during May 2017. The minimum mapping unit is 1600 m 2 [59]; this results in acquiring greater detail compared to the CORINE Land Cover (CLC) products using a minimum mapping unit of 25,000 m 2 . For the construction of the database, the European specifications of the CLC project from the Copernicus Program [60] were taken as a reference from which the first three levels of the LULC classification were derived [61]. The fourth level represents the categories of detail largely defined by the Emilia-Romagna Land Use Working Group that has  For the construction of the database, the European specifications of the CLC project from the Copernicus Program [60] were taken as a reference from which the first three levels of the LULC classification were derived [61]. The fourth level represents the categories of detail largely defined by the Emilia-Romagna Land Use Working Group that has operated in past years within the Interregional Centre for Geographic and Statistical Information Systems (CISIS, [62]). Polygons in the coverage are all defined by a four-digit numerical code (COD TOT) derived from CLC classification. The first number in the classification code corresponds to the highest hierarchical level, as 1 describes artificial surfaces; 2, agricultural lands; 3, forests and semi-natural areas; 4, wetlands and 5, water bodies. Polygons in categories where the fourth level is absent have the numerical code with the fourth digit set equal to zero.

Multispectral Orthophotos, Municipality Cadastral Map and Tree Census
Detection and classification of UGSs were performed by adopting the NDVI, modelled in an open-source GIS environment (QGIS v.3.x). Multispectral analyses were performed on aerial images (ortho-rectified photos) provided by the Agricultural Disbursement Agency (AGEA) and distributed for exclusive use to the Emilia-Romagna Region. Aerial image acquisition was completed in 2020 and currently represents the most up-to-date survey available at regional and national levels. Orthophotos are characterised by a very high spatial resolution (0.2 m pixel size), with four bands in the visible range (RGB) and nearinfrared (NIR) using a minimum mapping unit of 1600 m 2 , and a cloudy cover of less than 5%. They are at 1:5000 nominal scale provided in the global reference system ETRS89 and then reprojected in the RND2008 reference system (EPSG: 6875). Images were shot from February to May 2020.
The cadastral layer, provided by the Forlì Municipality, was used to intersect the NDVI values according to cadastral parcels for each district. Each cell is identified by type (i.e., street, parcel, or water) but not by its ownership status (public or private). This layer was then overlaid with the urban green cadastral map, also provided by the municipality office, containing green space polygons managed by the municipal authority. This information was integrated with the official census of trees (updated to 2021), covering the entire municipality with a buffer of 1 m around the tree position point and including road trees ( Figure 3). By using this method, it was possible to obtain a key to extract differences among the privately owned areas ( Figure 4). The NDVI was calculated, via QGIS Raster Calculator, for each district considered, using the following equation: where NIR represents the spectral reflectance recorded in the near-infrared region of the spectrum (from 780 nm to 2500 nm) and RED represents the spectral reflectance recorded in the red (visible) region (around 620 to 750 nm). According to [63], the range of NDVI is between ±1, where higher values correspond Land 2023, 12, x FOR PEER REVIEW 9 of 21 doubtful ownership were cross-checked and corrected if required. In some areas, the ownership status was also directly requested from the municipality by checking with cadastral codes.

Statistical Analysis
Statistical analysis was carried out associating the rasterized NDVI values, sampled over a regular 1 m square grid, to the CLC class ( Figure 5). After vectorization of the NDVI raster layer, it was sampled with a 1 m regular point grid. Each point corresponded to the NDVI centroid value to which the x and y coordinates had been associated. The 1-m grid was selected after several tests and considering an adequate ratio between precision and computational availability. To avoid QGIS computation miscalculations due to the large data frames (for example, there are 1,874,650 sampling points in the q6 district), RStudio software v.4.x [66] was employed. NDVI sampled layers were imported and associated by coordinates through a "right join process" to the corresponding CLC classes. In each district, CLC classes were ranked based on three NDVI categories, and the median values (low ≤ 0.0, medium > 0 and < 0.15, high ≥ 0.15). The NDVI was calculated, via QGIS Raster Calculator, for each district considered, using the following equation: where NIR represents the spectral reflectance recorded in the near-infrared region of the spectrum (from 780 nm to 2500 nm) and RED represents the spectral reflectance recorded in the red (visible) region (around 620 to 750 nm). According to [63], the range of NDVI is between ±1, where higher values correspond to vegetation in better health conditions. Negative values generally indicate non-green areas such as built-up areas or water bodies. NDVI threshold values are not univocal, but according to [46,64], using thresholds provides valuable information as to the number of green areas in urban contexts. Values ≤ 0.1 represent absence or scarcity of vegetation, while moderate values identifying between 0.2 and 0.3 indicate shrubs and grassland [53]. After several tests, a threshold suitable for extracting green areas in the Forlì urban context was selected, and pixels with NDVI values greater than 0.15 were considered vegetated (similar thresholds have been selected by [46,65]). UGS which matched with the cadastral municipal urban green layer polygons ( Figure 4) were selected, and the difference was calculated concerning the total NDVI layer, called "Cadastral NDVI urban green" (Figure 4), thus extracting areas designated in the records as privately owned (Private green areas in Figure 4). Visual analysis was then performed using very high-resolution satellite images available on Google Maps and ground photos from Google Street View; cases of doubtful ownership were cross-checked and corrected if required. In some areas, the ownership status was also directly requested from the municipality by checking with cadastral codes.

Statistical Analysis
Statistical analysis was carried out associating the rasterized NDVI values, sampled over a regular 1 m square grid, to the CLC class ( Figure 5).

Statistical Analysis
Statistical analysis was carried out associating the rasterized NDVI values, samp over a regular 1 m square grid, to the CLC class ( Figure 5). After vectorization of the NDVI raster layer, it was sampled with a 1 m regular po grid. Each point corresponded to the NDVI centroid value to which the x and y coor nates had been associated. The 1-m grid was selected after several tests and consideri an adequate ratio between precision and computational availability. To avoid QGIS co putation miscalculations due to the large data frames (for example, there are 1,874,6 sampling points in the q6 district), RStudio software v.4.x [66] was employed. NDVI sa pled layers were imported and associated by coordinates through a "right join proce to the corresponding CLC classes. In each district, CLC classes were ranked based on th NDVI categories, and the median values (low ≤ 0.0, medium > 0 and < 0.15, high ≥ 0.1 After vectorization of the NDVI raster layer, it was sampled with a 1 m regular point grid. Each point corresponded to the NDVI centroid value to which the x and y coordinates had been associated. The 1-m grid was selected after several tests and considering an adequate ratio between precision and computational availability. To avoid QGIS computation miscalculations due to the large data frames (for example, there are 1,874,650 sampling points in the q6 district), RStudio software v.4.x [66] was employed. NDVI sampled layers were imported and associated by coordinates through a "right join process" to the corresponding CLC classes. In each district, CLC classes were ranked based on three NDVI categories, and the median values (low ≤ 0.0, medium > 0 and <0.15, high ≥ 0.15). Each class corresponds to a different degree of greenery: low, medium and high. A one-way analysis of variance (ANOVA) test was carried out to assess the differences in NDVI reflectance within the same CLC units in the six districts of the study area. The Tukey-Kramer HSD test for unequal N was used to assess the statistical significance of differences, as expressed in mean values.

Public and Private Green Spaces Mapping
UGSs spatial distribution after NDVI analysis is shown in Figure 6. Higher NDVI values (≥0.15) correspond to areas covered by vegetation, while lower values (<0.15) are unvegetated artificial surfaces or water bodies. According to the map legend, areas that visually show higher values of NDVI are located in the north-western part of the study area, where large areas characterise districts q4 and q5 under agricultural land use. Orange and dark orange colours correspond to the more densely urbanised areas in the eastern districts (q1, q2 and q3). Areas in red, such as those in q3 district, represent natural and artificial water bodies. Figure 7 presents the UGSs after the application of the threshold (0.15) of NDVI. Light grey parts represent urban and peri-urban spaces covered by unvegetated artificial surfaces or water. UGSs spatial distribution after NDVI analysis is shown in Figure 6. Higher NDVI values (≥0.15) correspond to areas covered by vegetation, while lower values (<0.15) are unvegetated artificial surfaces or water bodies. According to the map legend, areas that visually show higher values of NDVI are located in the north-western part of the study area, where large areas characterise districts q4 and q5 under agricultural land use. Orange and dark orange colours correspond to the more densely urbanised areas in the eastern districts (q1, q2 and q3). Areas in red, such as those in q3 district, represent natural and artificial water bodies.   Regarding the distribution of green spaces in the six selected districts, the total green area is 5.70 km 2 , corresponding to about 42% of the total area. Among the public green allotment, there are different urban parks. They include the "Resistenza" Park (0.45 km 2 ) in the q1 north-western district, a historic park of the city of Forlì, established in the 1800s, with a further expansion during the 1930s, and the "Via Dragoni" Park in the q2 district, which represents a significant public green area (0.06 km 2 ), occupying 3.2% of the total Regarding the distribution of green spaces in the six selected districts, the total green area is 5.70 km 2 , corresponding to about 42% of the total area. Among the public green allotment, there are different urban parks. They include the "Resistenza" Park (0.45 km 2 ) in the q1 north-western district, a historic park of the city of Forlì, established in the 1800s, with a further expansion during the 1930s, and the "Via Dragoni" Park in the q2 district, which represents a significant public green area (0.06 km 2 ), occupying 3.2% of the total 12.1% public green area.
The district q6 ranks among those with the lowest presence of green (21.5%) due to its extremely dense urban fabric, with many roadside green areas of limited size, district gardens and private gardens. However, more green areas were recently created by urban regeneration projects as the University Campus Park (in the south-eastern part of the district) and the Museums Garden in front of the San Domenico Museum complex (in the western part), which were not visible in the 2020 AGEA orthophoto as they were completed in 2021.
As for the classification based on ownership, 72.1% of the total green area is private, and 28.9% is public. A detailed quantitative summary for each district is reported in Table 2, and the corresponding cartographical representation is shown in Figure 8.    In q6 (Centro Storico district, Figure 8), public and private areas are dimensionally comparable, with a total public green area of 0.18 km 2 (10.3%) and a total private green area of 0.20 km 2 (11.2%). The same situation is found in q1 (Spazzoli district), where private UGSs amount to 0.37 km 2 (19.7%), whereas the public areas amount to 0.31 km 2 (16.2%). Differences between publicly owned and private green areas are found in q4 (Foro Boario), with 1.44 km 2 of private areas (37.3%) and 0.29 km 2 of public green (7.6%), and in q2 (Grandi Musicisti Italiani district) where private green spaces are much larger than public ones, with 0.41 km 2 for private areas (21.9%) and 0.22 km 2 for public green (12.1%). The Romiti district, q5, characterised by a relevant share of agricultural land use, has much more private green spaces (1.53 km 2 , 51.7%) than public green spaces (0.19 km 2 , 6.6%), most of which are located along the Montone riverbanks.
The only case, among the districts analysed, in which the share of public green spaces is greater than the private one, is in q3 (Resistenza district), with 0.16 km 2 of private green spaces (12.8%) vs. 0.38 km 2 of public green spaces (31.1%), due to the presence of a large part of the "Franco Agosto" Park, which also covers a considerable share of the entire territory of the district, amounting to about 24% of total public green area.

NDVI and Corine Land Cover Units
Mean NDVI values associated with different CLC units showed statistically significant differences (p < 0.05). For example, Figure 9 reports the "box and whiskers" plot for the q6 central district with the lowest number of CLC classes. All classes present statistically significant differences (red letters in Figure 9), with an overall mean NDVI value of 0.065. Supplementary Materials report the figures for the other five districts of the study area.
Land 2023, 12, x FOR PEER REVIEW 13 of 21 (12.1%). The Romiti district, q5, characterised by a relevant share of agricultural land use, has much more private green spaces (1.53 km 2 , 51.7%) than public green spaces (0.19 km 2 , 6.6%), most of which are located along the Montone riverbanks. The only case, among the districts analysed, in which the share of public green spaces is greater than the private one, is in q3 (Resistenza district), with 0.16 km 2 of private green spaces (12.8%) vs. 0.38 km 2 of public green spaces (31.1%), due to the presence of a large part of the "Franco Agosto" Park, which also covers a considerable share of the entire territory of the district, amounting to about 24% of total public green area.

NDVI and Corine Land Cover Units
Mean NDVI values associated with different CLC units showed statistically significant differences (p < 0.05). For example, Figure 9 reports the "box and whiskers" plot for the q6 central district with the lowest number of CLC classes. All classes present statistically significant differences (red letters in Figure 9), with an overall mean NDVI value of 0.065. Supplementary Materials report the figures for the other five districts of the study area. NDVI median values and interquartile ranges for each district are shown in Figure  10, grouped by CLC code. The lowest values are found in CLC classes with a high presence of sealed surfaces (e.g., 1112 sparse residential fabric and 1222 street networks). CLC classes 1121 (urban residential fabric) and 1214 (public and private service settlements) show different NDVI median values, ranging from −0.06 to 0.31 and from −0.08 to 0.25, respectively, resulting from the heterogenous composition of the urban landscape. As expected, the highest median NDVI values are found for CLC class 1411 (parks, >0.05). NDVI median values and interquartile ranges for each district are shown in Figure 10, grouped by CLC code. The lowest values are found in CLC classes with a high presence of sealed surfaces (e.g., 1112 sparse residential fabric and 1222 street networks). CLC classes 1121 (urban residential fabric) and 1214 (public and private service settlements) show different NDVI median values, ranging from −0.06 to 0.31 and from −0.08 to 0.25, Land 2023, 12, 660 12 of 18 respectively, resulting from the heterogenous composition of the urban landscape. As expected, the highest median NDVI values are found for CLC class 1411 (parks, >0.05).
Land 2023, 12, x FOR PEER REVIEW 14 of 21 Figure 10. Graphics of Corine Land Cover classes found in overall districts of the study area., and the NDVI classification based in classes of Table S1 in Supplementary Materials. In orange, are the NDVI median values and the interquartile range. Blue lines represent the limits of the median NDVI classes (low ≤ 0.0, medium > 0 and < 0.15, high ≥ 0.15). The legend for CLC classes is given in the caption of Figure 2.
In each district, all CLC classes were ranked according to three classes of NDVI (Figure 11) based on their median values (low ≤ 0.0, medium > 0 and < 0.15, high ≥ 0.15). The correlation between the median value and the interquartile range resulted in two opposite trends: low and medium NDVI CLC classes exhibited statistically significant direct relationships (p < 0.05) between median values and interquartile ranges (r equal to 0.563 and 0.837 respectively, for low and medium NDVI), while the opposite was observed for the high NDVI class (r = −0.422). In the first case, an increase in NDVI was associated with an increase in the variability of the spectral response due to an increase in vegetated areas within a predominantly artificialized urban fabric. However, in the second case, the increase in vegetated areas reduced the variability of the spectral response as it occurred in an urban fabric already characterised by a significant number of green areas. The resulting spatial distribution of classed NDVI values for the CLC units is shown in Figure 11.  Table S1 in Supplementary Materials. In orange, are the NDVI median values and the interquartile range. Blue lines represent the limits of the median NDVI classes (low ≤ 0.0, medium > 0 and <0.15, high ≥ 0.15). The legend for CLC classes is given in the caption of Figure 2.
In each district, all CLC classes were ranked according to three classes of NDVI (Figure 11) based on their median values (low ≤ 0.0, medium > 0 and <0.15, high ≥ 0.15). The correlation between the median value and the interquartile range resulted in two opposite trends: low and medium NDVI CLC classes exhibited statistically significant direct relationships (p < 0.05) between median values and interquartile ranges (r equal to 0.563 and 0.837 respectively, for low and medium NDVI), while the opposite was observed for the high NDVI class (r = −0.422). In the first case, an increase in NDVI was associated with an increase in the variability of the spectral response due to an increase in vegetated areas within a predominantly artificialized urban fabric. However, in the second case, the increase in vegetated areas reduced the variability of the spectral response as it occurred in an urban fabric already characterised by a significant number of green areas. The resulting spatial distribution of classed NDVI values for the CLC units is shown in Figure 11.
As expected, in all the districts, the CLC classes having the lowest levels of NDVI turn out to be: (i) sparse residential fabric (CLC class 1112, NDVI median range −0.049 to −0.085), (ii) road networks (CLC class 1222, NDVI median range −0.090 to −0.122), followed, where present, by (iii) areas with industrial, artisanal, and agricultural production settlements (CLC class 1211, NDVI median range −0.080 to −0.106), (iv) commercial settlements (CLC class 1213, NDVI median range −0.090 to −0.114), (v) highways and freeways (CLC class 1221, NDVI median range −0.117 to −0.145), (vi) rail networks and (vii) ancillary spaces (CLC class 1224, NDVI median range −0.058 to −0.101). The only district that has the CLC class 1111, which corresponds to the compact and dense residential fabric, is q6, the one with the lowest level of vegetation cover. In fact, the area is composed of a dense network of historic buildings with little free space for green areas, except for the  It is worth noting the discrepancy between the NDVI values in vineyards (CLC class 2210) found in districts q4 and q5, which reflect the timing of the orthophotos: the low level of greenery in q4 is due to the lack of plantations because as seen in the orthophoto, there is no canopy cover and the soil is almost bare or only partially covered by grass (NDVI median value −0.024); while the high level of canopy cover (NDVI median value 0.224) in q5 is due to the remarkable planting beds.
Another example of discrepancies between different districts is found in CLC class 2123 (horticultural crops in open field, greenhouse and under plastic); specifically, in q4 and q5, the areas under this class show a low value of NDVI as compared to the same class in q3 district, due to coverage with artificial materials or the temporary absence of a crop.

Discussion
The NDVI-based approach, through high-resolution orthophotos, allowed the classification of UGSs [67] and links to their property status via cadastral data, whereas ancillary spatial data made it possible to analyse NDVI responses in terms of homogeneity of the reflectance response in land cover units.
Considering the outcomes of the first step of the analysis, in terms of ownership status, the prevalence of private UGS was significantly higher (72.1% covering 4.11 km 2 ) than public ones (28.9% covering 1.59 km 2 ). The authors of [46] reported similar results for the entire city of Padua (NE Italy), with a prevalence of private green, which adds up to ca. 80% of the total greenery surveyed within the municipality. In [68], the amount of urban private green spaces in Leipzig (NE Germany), resulting from back-and front-yard green space, is around 40% (20 km 2 ) of the total amount of public green spaces, but in 25% of city districts, it prevails over public green areas. Other studies focused on quantifying private garden areas, for example, estimating between 22 and 36% of the total Dunedin urban area in New Zealand [69]. In line with these data, [35] found that about 23% (33 km 2 ) of the urban area in Sheffield (UK) is covered by home gardens. Even if it is not a standardised measurement for all cities, it has been estimated that in Forlì, the total amount of public green space per capita is 32.1 m 2 per inhabitant, slightly above the national average of 31.1 m 2 per inhabitant, but lower than in other north-eastern Italian cities with an average of 50 m 2 per capita [70]. This indicates an overall optimal UGS amount per-capita value, considering that the World Health Organization (WHO) has suggested that every city should have a minimum of 9 m 2 of green space per capita, with an optimum between 10 and 15 m 2 per capita.
At the district level, q4 and q5 are characterised by having the highest percentages in private UGSs due to the prevalence of agricultural land use and smaller dense residential units. On the contrary, in q3 district, the presence of the "Franco Agosto" Park, with a relevant extension (0.27 km 2 ), increases the amount of public green spaces (31.1%) concerning, which represents only 12.8% of the total. More densely populated areas, such as the q6 district, show lower percentages of green coverage irrespective of ownership status. According to [46,68], the large number of private green spaces should be explicitly considered in UP, developing site-specific policies that consider their contribution to urban ecosystem services provision and allowing for a more complex ecosystem services-oriented UP.
The analysis conducted on NDVI values about LULC classes must consider the presence of numerous outliers in all investigated districts. For these, very low NDVI outliers indicate the presence of artificial surfaces in predominantly vegetated areas. In contrast, outliers at the distribution's upper tail denote vegetation cover in predominantly artificial areas. More densely urbanized areas, such as q6 and adjacent areas, are generally covered by artificial surfaces with a high degree of imperviousness. On the other hand, according to [71], bare soil and agricultural land may be confused with impervious surfaces with low NDVI values when the same surfaces are dry or barely vegetated, as observed for vineyards in our case study area.

Conclusions
Mapping and characterization of UGSs is a crucial node in present and future UP at local, national and European levels for more efficient UP and green areas management. At local level, this study focuses on: (i) detecting and mapping UGS in the Municipality of Forlì, identifying private and public green surfaces, (ii) providing a reference framework for further research aimed at quantifying ecosystem services provided by urban green areas within the framework of the NRRP's interventions, and (iii) characterizing the heterogeneity of the NDVI spectral responses within and among LULC classes in the urban fabric.
Our results show that the total amount of green spaces is 5.70 km 2 , of which 72.1% (4.11 km 2 ) is private, and 28.9% (1.59 km 2 ) is classified as public. As for the land cover, three NDVI classes were identified, highlighting different degrees of homogeneity in NDVI reflectance response within each urban land cover unit.
Deep understanding of the ownership status, type and distribution of UGS can be beneficial for public authorities and policymakers in developing site-specific policies and foster incentives aimed at encouraging sustainable urban green management practices [72].
In the future steps of this research, the use of specific and appropriately calibrated models, such as i-Tree (UFORE), will enable the estimation of the ecosystem services provided by the urban forest, such as reduction of water runoff, removal of air pollutants (particulate matter), carbon sequestration and its fixation within plant biomass, all based on an explicit knowledge of the structure and composition of the urban forest. The results of this study will also provide the inputs for the implementation of a prediction model of the accessibility to green areas by the public, based on lower-resolution satellite images. Further studies could also concentrate on differentiation between various types of green spaces, such as grasses, shrubs or bushes, as suggested by [73]. In addition, our results are of interest for further scientific and applied studies for city planners, as this might guide afforestation practices foreseen by the NRRP and encourage the ecosystem services provision of the urban forest.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/land12030660/s1, Figure S1: q1 (Spazzoli district) Box and Whiskers plot with HSD Tukey of NDVI spectral response on Land Cover database; Figure S2: q2 (Musicisti Grandi Italiani district) Box and Whiskers plot with HSD Tukey of NDVI spectral response on Land Cover database; Figure S3: q3 (Resistenza district) Box and Whiskers plot with HSD Tukey of NDVI spectral response on Land Cover database; Figure S4: q4 (Foro Boario district) Box and Whiskers plot with HSD Tukey of NDVI spectral response on Land Cover database; Figure S5: q5 (Romiti district) Box and Whiskers plot with HSD Tukey of NDVI spectral response on Land Cover database; Table S1: Summary table about greenery level and the Land use and cover by Corine revised classification, according to NDVI values in different districts of the study area.