Urban Tree Species Identification and Carbon Stock Mapping for Urban Green Planning and Management

Recently, the severe intensification of atmospheric carbon has highlighted the importance of urban tree contributions in atmospheric carbon mitigations in city areas considering sustainable urban green planning and management systems. Explicit and timely information on urban trees and their roles in the atmospheric Carbon Stock (CS) are essential for policymakers to take immediate actions to ameliorate the effects of deforestation and their worsening outcomes. In this study, a detailed methodology for urban tree CS calibration and mapping was developed for the small urban area of Sassuolo in Italy. For dominant tree species classification, a remote sensing approach was applied, utilizing a high-resolution WV3 image. Five dominant species were identified and classified by applying the Object-Based Image Analysis (OBIA) approach with an overall accuracy of 78%. The CS calibration was done by utilizing an allometric model based on the field data of tree dendrometry—i.e., Height (H) and Diameter at Breast Height (DBH). For geometric measurements, a terrestrial photogrammetric approach known as Structure-from-Motion (SfM) was utilized. Out of 22 randomly selected sample plots of 100 square meters (10 m × 10 m) each, seven plots were utilized to validate the results of the CS calibration and mapping. In this study, CS mapping was done in an efficient and convenient way, highlighting higher CS and lower CS zones while recognizing the dominant tree species contributions. This study will help city planners initiate CS mapping and predict the possible CS for larger urban regions to ensure a sustainable urban green management system.


Introduction
In recent decades, ceaseless urbanization and intensifying global warming have led to increased levels of atmospheric carbon in urban areas. For city planners and researchers, understanding the contribution of trees in mitigating atmospheric carbon in urban areas has become one of the paramount concerns. The planning and management of urban parks and trees to decelerate the worsening impacts of severe carbon emissions have been explored in many different studies [1][2][3][4][5][6][7]. For instance, empirical research found that the mean aboveground urban tree carbon stock (dry biomass) across the Seattle urbanizing region was 89 ± 22 Mg C/ha, which is a significant quantity considering the severe atmospheric carbon in city areas [8,9]. Thus, there is now a growing need to study urban green This study was conducted in the northern part of Italy focusing on the small city area of Sassuolo, which is also known as one of the most severely polluted urban areas in Europe. An allometric model [63] developed specifically for the trees in Italy was utilised for the CS calibration. Five dominant species were identified covering the whole study area, and the CS/plot was calibrated and mapped based on each species.
Herein, the main goal was to understand the potential of this methodology in mapping the tree CS in an urban area considering the role of dominant urban tree species in atmospheric carbon mitigation. This method might be an efficient tool to assist urban planners in CS mapping to ensure the proper utilisation of the available green space. Thus, this study highlights that accurate tree CS mapping is crucial for estimating and identifying the dominant species contributing a significant level of atmospheric CS, which will prove to be useful support for urban planners and environmental policymakers in planning further urban air quality assessments. This study will contribute to a better understanding of this methodology in mapping structural and functional properties, such as tree CS, as well as predicting the possible urban CS in typical city areas.

Data Set
The study area is the city of Sassuolo, province of Modena, located in the Po Valley in the northern part of Italy ( Figure 1). The Po valley experiences strong anthropic pressure due to its large urban areas, intensive agriculture (among the most productive agricultural areas within Europe), and large manufacturing districts, along with topographic and meteorological conditions unfavourable to pollutant dispersion [64][65][66]. This study was conducted in the northern part of Italy focusing on the small city area of Sassuolo, which is also known as one of the most severely polluted urban areas in Europe. An allometric model [63] developed specifically for the trees in Italy was utilised for the CS calibration. Five dominant species were identified covering the whole study area, and the CS/plot was calibrated and mapped based on each species.
Herein, the main goal was to understand the potential of this methodology in mapping the tree CS in an urban area considering the role of dominant urban tree species in atmospheric carbon mitigation. This method might be an efficient tool to assist urban planners in CS mapping to ensure the proper utilisation of the available green space. Thus, this study highlights that accurate tree CS mapping is crucial for estimating and identifying the dominant species contributing a significant level of atmospheric CS, which will prove to be useful support for urban planners and environmental policymakers in planning further urban air quality assessments. This study will contribute to a better understanding of this methodology in mapping structural and functional properties, such as tree CS, as well as predicting the possible urban CS in typical city areas.

Data Set
The study area is the city of Sassuolo, province of Modena, located in the Po Valley in the northern part of Italy ( Figure 1). The Po valley experiences strong anthropic pressure due to its large urban areas, intensive agriculture (among the most productive agricultural areas within Europe), and large manufacturing districts, along with topographic and meteorological conditions unfavourable to pollutant dispersion [64][65][66]. The total investigation area is an urban area of around 273 ha ( Figure 1). Five dominant tree species (Quercus spp., Acer campestre, Populus nigra, Platanus spp., and Tilia platyphyllos) were identified, covering all the main streets and parks.
The key data for this study were WorldView 3(WV3) images, shapefiles available from the Geoportal of the Emilia Romagna Region [67], and data collected from the sample plots. Shapefiles were extracted from the CORINE Land Cover (CLC) map [68], a database of land use in the territory. In this study, artificial surfaces were extracted from the map, such as streets and public green areas. The WV3 image data for this study were acquired on 31 July 2018 ( Figure 2). The WV3 includes one panchromatic band with a 0.3 m resolution and eight multispectral bands with a 1.2 m resolution (Table 1). The total investigation area is an urban area of around 273 ha ( Figure 1). Five dominant tree species (Quercus spp., Acer campestre, Populus nigra, Platanus spp., and Tilia platyphyllos) were identified, covering all the main streets and parks.
The key data for this study were WorldView 3(WV3) images, shapefiles available from the Geoportal of the Emilia Romagna Region [67], and data collected from the sample plots. Shapefiles were extracted from the CORINE Land Cover (CLC) map [68], a database of land use in the territory. In this study, artificial surfaces were extracted from the map, such as streets and public green areas. The WV3 image data for this study were acquired on 31 July 2018 ( Figure 2). The WV3 includes one panchromatic band with a 0.3 m resolution and eight multispectral bands with a 1.2 m resolution (Table 1). The image also includes eight SWIR bands with a spatial resolution of 7.5 m. As this spatial resolution is not very suitable for tree species classification, these bands were not considered in this study. However, the additional four bands-i.e., coastal, yellow, red-edge and NIR2-did provide some significant advantages for vegetation identification and classification [47,69].

Data Collection
For the field data collection, 22 plots of 100 m 2 (10 m × 10 m) each ( Figure 3) were randomly selected considering the dominant tree species covering the whole study area of Sassuolo. Among them, seven plots were considered for the CS mapping and validation. The fundamental considerations were to measure the DBH and H of each tree, along with each tree's geographical coordinates, and to determine the names of the species. DBH measurements were done at 1.3 m above ground level. To obtain the H of the trees, a photogrammetric approach was applied. For this purpose, 2D pictures were taken at each plot utilizing a handheld digital camera while considering 3D model development through the SfM photogrammetric approach [23]. Several pictures were acquired from different positions, thereby differentiating between short and tall trees. Field data (H, DBH) were collected in October 2018 and were essential for the CS calibration. The sample plots were also used during the training and validation of the tree species classification. The time gap between the acquisition dates of the WV3 images (July) and the field data collection (October) did not affect the identification of the different species and their classification. Indeed, 3D reconstruction of the trees to calculate their H is more important when there are only a few leaves (better vision of the tops of the trees). Moreover, for tree species classification, the image from July was especially convenient, as the CS computation (In QGIS; Section 2.2.4) was based on the NDVI, which is strictly correlated with tree biomass.  The image also includes eight SWIR bands with a spatial resolution of 7.5 m. As this spatial resolution is not very suitable for tree species classification, these bands were not considered in this study. However, the additional four bands-i.e., coastal, yellow, red-edge and NIR2-did provide some significant advantages for vegetation identification and classification [47,69].

Data Collection
For the field data collection, 22 plots of 100 m 2 (10 m × 10 m) each ( Figure 3) were randomly selected considering the dominant tree species covering the whole study area of Sassuolo. Among them, seven plots were considered for the CS mapping and validation. The fundamental considerations were to measure the DBH and H of each tree, along with each tree's geographical coordinates, and to determine the names of the species. DBH measurements were done at 1.3 m above ground level. To obtain the H of the trees, a photogrammetric approach was applied. For this purpose, 2D pictures were taken at each plot utilizing a handheld digital camera while considering 3D model development through the SfM photogrammetric approach [23]. Several pictures were acquired from different positions, thereby differentiating between short and tall trees. Field data (H, DBH) were collected in October 2018 and were essential for the CS calibration. The sample plots were also used during the training and validation of the tree species classification. The time gap between the acquisition dates of the WV3 images (July) and the field data collection (October) did not affect the identification of the different species and their classification. Indeed, 3D reconstruction of the trees to calculate their H is more important when there are only a few leaves (better vision of the tops of the trees). Moreover, for tree species classification, the image from July was especially convenient, as the CS computation (In QGIS; Section 2.2.4) was based on the NDVI, which is strictly correlated with tree biomass.

Pre-Processing
In the first stage, the WV3 multispectral image was processed to convert Digital Numbers into the Top of Atmosphere (TOA) radiance. Then atmospheric corrections were performed using the FLAASH plugin of the Environment for Visualizing Images (ENVI) software [70]. In this way, it was possible to retrieve the Surface Reflectance (SR). Then, some of the required features-i.e., the Normalized Difference Vegetation Index (NDVI), the Grey Level Co-occurrence Matrix (GLCM), and the Principal Component Analysis (PCA) for OBIA classification [47]-were computed.
NDVI is a very common vegetation index used in several remote sensing studies [71]. For evaluating a vegetative area, the equation for calculation of the NDVI is as follows: = reflectance in the near-infrared band; = reflectance in the red band.
For this work, several NDVI calculations were performed using both the NIR1 and NIR2 Near Infrared bands along with the red edge band.
GLCM is a tabulation of the frequency of different combinations of grey levels at a specified distance and orientation in an image object [72]. For principal component analysis, only the first and second principal components were used. In this way, it was possible to retrieve the less-correlated information from the original eight multispectral bands [73].

Object Based Classification
In urban areas, the traditional pixel-based image classification method usually provides low classification accuracy due to the high spectral variability within the land cover classes affected by the sun angle, gaps in the tree canopies, and shadows [74,75]. Indeed, urban tree classification with higher accuracy is still a considerable challenge, and most studies recommend the application of the Object Based Image Analysis (OBIA) approach to improve classification accuracy in urban areas [49][50][51][52]. Herein, the OBIA approach was also applied to improve classification accuracy [49][50][51][52]76]. The OBIA method includes not only spectral information but also other added information, such as context, texture, geometry, and spatial features [36,77], which can minimise the number of units to be considered for the classification [78]. Segmentation is the key procedure used to divide an image into different significant objects in which the spectral and spatial features will be computed. The segmentation procedure divides the image into spatially continuous and homogeneous regions [79]

Pre-Processing
In the first stage, the WV3 multispectral image was processed to convert Digital Numbers into the Top of Atmosphere (TOA) radiance. Then atmospheric corrections were performed using the FLAASH plugin of the Environment for Visualizing Images (ENVI) software [70]. In this way, it was possible to retrieve the Surface Reflectance (SR). Then, some of the required features-i.e., the Normalized Difference Vegetation Index (NDVI), the Grey Level Co-occurrence Matrix (GLCM), and the Principal Component Analysis (PCA) for OBIA classification [47]-were computed.
NDVI is a very common vegetation index used in several remote sensing studies [71]. For evaluating a vegetative area, the equation for calculation of the NDVI is as follows: where ρ NIR = reflectance in the near-infrared band; ρ RED = reflectance in the red band.
For this work, several NDVI calculations were performed using both the NIR1 and NIR2 Near Infrared bands along with the red edge band.
GLCM is a tabulation of the frequency of different combinations of grey levels at a specified distance and orientation in an image object [72]. For principal component analysis, only the first and second principal components were used. In this way, it was possible to retrieve the less-correlated information from the original eight multispectral bands [73].

Object Based Classification
In urban areas, the traditional pixel-based image classification method usually provides low classification accuracy due to the high spectral variability within the land cover classes affected by the sun angle, gaps in the tree canopies, and shadows [74,75]. Indeed, urban tree classification with higher accuracy is still a considerable challenge, and most studies recommend the application of the Object Based Image Analysis (OBIA) approach to improve classification accuracy in urban areas [49][50][51][52]. Herein, the OBIA approach was also applied to improve classification accuracy [49][50][51][52]76]. The OBIA method includes not only spectral information but also other added information, such as context, texture, geometry, and spatial features [36,77], which can minimise the number of units to be considered for the classification [78]. Segmentation is the key procedure used to divide an image into different Forests 2020, 11, 1226 6 of 22 significant objects in which the spectral and spatial features will be computed. The segmentation procedure divides the image into spatially continuous and homogeneous regions [79] and limits local spectral variation [49,50]. For this study, the Trimble eCognition Developer ® 9 platform (Trimble, Munich, Germany) [80] was utilised for the OBIA approach.
In this study, both segmentation and classification were implemented in successive steps to obtain the best results. At first, a chessboard segmentation was carried out using the shapefiles of the Municipality of Sassuolo extracted from the Corine Land Cover (CLC) map of the Emilia-Romagna Region. In particular, the shapefiles of streets and public green areas (i.e., parks) were utilized.
Green areas and streets were identified with very large objects (used as the sizes of the shapefile polygons) and classified separately from the other objects. In these areas, subsequent multiresolution segmentation was applied to the smaller objects. For this segmentation, instead of thematic layers, the spectral information and geometric information of the WV3 bands were considered. After that, based on the study of Li et al. [47], the spectral and textural features (NDVIs, GLCM, and PCA) were included to obtain better results in classification. A total of 79 features were added to the eight bands of the image.
The result after the second segmentation was a series of smaller objects with a good degree of homogeneity among themselves. The resulting medium-sized objects were suitable for the study of tree crowns. The parameters used for segmentation, after numerous tests, were determined as follows: scale parameter = 10, shape = 0.2, and compactness = 0.5.
By using the threshold values for some of the features such as the Brightness and NDVI, the initial classes such as shadows, pathways and grass were identified. Then, relying on the sample plots, several ground-truth samples (ROIs, Regions of Interest) were chosen for the different tree species. Dominant species were identified distinctly for parks and roads. This was done to obtain classification with better accuracy, as the tree species diversity in parks is usually higher than that on the streets. In the case of parks, the main classes were identified as follows: In the case of streets, the main classes were identified as follows: The sample plots were used for the classification algorithm training, where the chosen method used the Nearest Neighbour (NN) algorithm [81].
After the training step, the NN algorithm was implemented for the assignment of objects to the different classes. Features like the mean and the standard deviation for all bands, the average values of the different NDVI indexes, the Grey level Co-occurrence Matrix, and the Principal Components were included in the feature space [59]. In this study, these features were chosen to provide improved classification under the OBIA approach. For instance, the NN algorithm was used separately for trees in the parks and on the streets. This algorithm was applied to select the different dominant species and thereby obtain more accurate tree species classification. Thus, in the classification map, the suffix -street or -park was used even for the same tree species (see the "Results" section). In this way, the accuracy of the classification was improved because the sampling points were chosen separately for parks and streets. Once the classification was done, other samples were used for the validation phase. In this way, it was possible to calculate the accuracy of the proposed methodology in this area.

Photogrammetric Approach
For the tree H measurement, the SfM photogrammetry approach was introduced ( Figure 4). This approach usually allows for the 3D reconstructions of objects through the acquisition of 2D images for estimating the required features [23]. In this study, a 3D model was developed for each plot to measure the individual tree H to calibrate the CS via the allometric model [63]. The Agisoft Metashape (AM) software was utilized to produce 3D models applying the SfM algorithm [23,59,82]. In each plot, all the trees were photographed from different directions to ensure maximum data redundancy and the best accuracies. The images were captured with a handheld digital SLR camera from photo-points at regular intervals around the perimeter of each plot. A semi-circular pathway around each plot from the inner to outer part of the plot was followed to capture each image successively, thereby covering all the trees ( Figure 5). Images were acquired with different camera inclinations to avoid occlusions. About 120-150 photos/plot were taken depending on the sizes of the trees.
Forests 2020, 11, x FOR PEER REVIEW 7 of 21 measure the individual tree H to calibrate the CS via the allometric model [63]. The Agisoft Metashape (AM) software was utilized to produce 3D models applying the SfM algorithm [23,59,82].
In each plot, all the trees were photographed from different directions to ensure maximum data redundancy and the best accuracies. The images were captured with a handheld digital SLR camera from photo-points at regular intervals around the perimeter of each plot. A semi-circular pathway around each plot from the inner to outer part of the plot was followed to capture each image successively, thereby covering all the trees ( Figure 5). Images were acquired with different camera inclinations to avoid occlusions. About 120-150 photos/plot were taken depending on the sizes of the trees.  The inner circle photos were taken from a distance to include the tree-top, middle, and bottom parts of the tree, and the outer photos were oriented to include the whole tree within the frame. This approach ensured a high redundancy of images (each portion of the area of interest must be detected to have 9-10 images at least) [23]. Moreover, a varying view of the geometry must be guaranteed to  measure the individual tree H to calibrate the CS via the allometric model [63]. The Agisoft Metashape (AM) software was utilized to produce 3D models applying the SfM algorithm [23,59,82].
In each plot, all the trees were photographed from different directions to ensure maximum data redundancy and the best accuracies. The images were captured with a handheld digital SLR camera from photo-points at regular intervals around the perimeter of each plot. A semi-circular pathway around each plot from the inner to outer part of the plot was followed to capture each image successively, thereby covering all the trees ( Figure 5). Images were acquired with different camera inclinations to avoid occlusions. About 120-150 photos/plot were taken depending on the sizes of the trees.  The inner circle photos were taken from a distance to include the tree-top, middle, and bottom parts of the tree, and the outer photos were oriented to include the whole tree within the frame. This approach ensured a high redundancy of images (each portion of the area of interest must be detected to have 9-10 images at least) [23]. Moreover, a varying view of the geometry must be guaranteed to  The inner circle photos were taken from a distance to include the tree-top, middle, and bottom parts of the tree, and the outer photos were oriented to include the whole tree within the frame. This approach ensured a high redundancy of images (each portion of the area of interest must be detected to have 9-10 images at least) [23]. Moreover, a varying view of the geometry must be guaranteed to obtain accurate results. If no information about the camera position or points with known coordinates are added to the project, the reconstructed 3D point cloud will be generated within an arbitrary reference system lacking georeferencing and scale.
In this case study, to obtain a 3D scaled point cloud and measure the tree H in the 3D model, a 4 m graduated bar ( Figure 6) was added to the area of interest as a reference and depicted in various images. As the positioning of an object with known dimensions is an external constraint used for reconstruction that allows the generation of scaled products, the applied methodology is feasible for the measurement purposes of this study [83]. Tree H estimation was done by performing measurements on the reconstructed 3D model. are added to the project, the reconstructed 3D point cloud will be generated within an arbitrary reference system lacking georeferencing and scale.
In this case study, to obtain a 3D scaled point cloud and measure the tree H in the 3D model, a 4 m graduated bar ( Figure 6) was added to the area of interest as a reference and depicted in various images. As the positioning of an object with known dimensions is an external constraint used for reconstruction that allows the generation of scaled products, the applied methodology is feasible for the measurement purposes of this study [83]. Tree H estimation was done by performing measurements on the reconstructed 3D model. The AM software was used to generate the 3D model for each plot, and in each model, the reference scale bar was visible enough to constrain the reconstruction ( Figure 6). From the 3D model, H estimation was done using the scale creation tool available in the AM software. Once the markers were placed on some of the aligned images, the corresponding 3D marker was automatically added to the 3D model. Then, the lengths of the segments among the points of interest were determined. The identified lengths provided the estimated tree H utilizing the SfM approach.
Validation of the results was mandatory to assess the metric accuracy of the performed measurements. The photogrammetric SfM approach was thus compared via laser scanning. A Leica C10 scan station was used for this purpose, and the resolution and accuracy of the resulting point clouds were about a few mm (<5 mm) in size. Laser scanning and photogrammetry were independent technological methodologies used for 3D model generation.
The results obtained with the proposed photogrammetric methodology were compared with those of the laser scanning technique to validate the results and guarantee the metric accuracy of the H measurements. The reconstruction of a small portion of a building was performed using both The AM software was used to generate the 3D model for each plot, and in each model, the reference scale bar was visible enough to constrain the reconstruction ( Figure 6). From the 3D model, H estimation was done using the scale creation tool available in the AM software. Once the markers were placed on some of the aligned images, the corresponding 3D marker was automatically added to the 3D model. Then, the lengths of the segments among the points of interest were determined. The identified lengths provided the estimated tree H utilizing the SfM approach.
Validation of the results was mandatory to assess the metric accuracy of the performed measurements. The photogrammetric SfM approach was thus compared via laser scanning. A Leica C10 scan station was used for this purpose, and the resolution and accuracy of the resulting point clouds were about a few mm (<5 mm) in size. Laser scanning and photogrammetry were independent technological methodologies used for 3D model generation.
The results obtained with the proposed photogrammetric methodology were compared with those of the laser scanning technique to validate the results and guarantee the metric accuracy of the H measurements. The reconstruction of a small portion of a building was performed using both technologies. The validation area was chosen to provide better identification of homologous points. A direct comparison of the generated point clouds was not suitable for this study, so the lengths of various segments were verified (Figure 7 and Table 2) to enable a proper investigation of the accuracy of the applied methodology. During the H estimation, the errors were measured while considering the estimated values in both cases. The identified differences in the two datasets were about 4 cm; if the longer segments (the ones exceeding 10 cm) are taken into consideration, the difference is 6 cm. A direct comparison of the generated point clouds was not suitable for this study, so the lengths of various segments were verified (Figure 7 and Table 2) to enable a proper investigation of the accuracy of the applied methodology. During the H estimation, the errors were measured while considering the estimated values in both cases. The identified differences in the two datasets were about 4 cm; if the longer segments (the ones exceeding 10 cm) are taken into consideration, the difference is 6 cm.

. CS Mapping in QGIS
During CS mapping in QGIS, the shapefile of the classified trees was utilized to determine the CS/plot for each identified species. First, the total Above Ground Biomass (AGB) was calculated considering the field data (i.e., the DBH, H, tree species, etc.) for each of the sample plots (Figure 4). For this calculation, an allometric model [63] was applied to calculate the AGB for each plot. The mean AGB/plot estimation was necessary, as it is recommended that the tree above ground CS be 50% of the total AGB [84][85][86][87][88]. Then, to estimate the mean CS/plot, the mean AGB/plot was multiplied by 0.5 as a conversion factor [89][90][91]. Next, in QGIS, the NDVI (Utilizing red edge and NIR1 bands) of the whole study area was computed utilizing the WV3 image data. The NDVI layer was considered in this study for CS predictions and mapping, as several studies claimed to find a strong correlation between the NDVI and total AGB of the trees [92][93][94][95]. The NDVI-derived metrics were extracted for the sample plots utilizing the "Zonal statistics" plugin [96] in the QGIS interface. After that, linear regression models were created in a Microsoft ® Excel™ spreadsheet considering the mean CS/plot and the NDVI-derived metrics to determine the best model to estimate and map CS covering the whole study area. A fishnet of a 100 m 2 (10 m × 10 m) resolution was built in QGIS to recognize the minimum to maximum CS zones based on the dominant species classification shapefile obtained

CS Mapping in QGIS
During CS mapping in QGIS, the shapefile of the classified trees was utilized to determine the CS/plot for each identified species. First, the total Above Ground Biomass (AGB) was calculated considering the field data (i.e., the DBH, H, tree species, etc.) for each of the sample plots (Figure 4). For this calculation, an allometric model [63] was applied to calculate the AGB for each plot. The mean AGB/plot estimation was necessary, as it is recommended that the tree above ground CS be 50% of the total AGB [84][85][86][87][88]. Then, to estimate the mean CS/plot, the mean AGB/plot was multiplied by 0.5 as a conversion factor [89][90][91]. Next, in QGIS, the NDVI (Utilizing red edge and NIR1 bands) of the whole study area was computed utilizing the WV3 image data. The NDVI layer was considered in this study for CS predictions and mapping, as several studies claimed to find a strong correlation between the NDVI and total AGB of the trees [92][93][94][95]. The NDVI-derived metrics were extracted for the sample plots utilizing the "Zonal statistics" plugin [96] in the QGIS interface. After that, linear regression models were created in a Microsoft ® Excel™ spreadsheet considering the mean CS/plot and the NDVI-derived metrics to determine the best model to estimate and map CS covering the whole study area. A fishnet of a 100 m 2 (10 m × 10 m) resolution was built in QGIS to recognize the minimum to maximum CS zones based on the dominant species classification shapefile obtained from the WV3 image data. The classification shapefile was essential to define the required zones that provided QGIS data to map the estimated CS values considering only the dominant tree species. Otherwise, the map would show the CS values for other areas-e.g., the areas covered with grass or even areas with no vegetation.

OBIA Classification Results and Validation
In this study, the OBIA approach was applied based on the spectral and textural attributes of the image objects for generating rule sets during classification. Although it was challenging to detect some of the classes due to their spectral similarities, 4 species in the parks and 3 species in the streets were found to be dominant in Sassuolo (Figure 8). Otherwise, the map would show the CS values for other areas-e.g., the areas covered with grass or even areas with no vegetation.

OBIA Classification Results and Validation
In this study, the OBIA approach was applied based on the spectral and textural attributes of the image objects for generating rule sets during classification. Although it was challenging to detect some of the classes due to their spectral similarities, 4 species in the parks and 3 species in the streets were found to be dominant in Sassuolo (Figure 8). The classification validation was done through the confusion matrix [97][98][99], which allows one to determine the overall accuracy of the classification in addition to some parameters related to user and producer accuracy. The presence of validation is necessary to define the effectiveness of the classification and determine the error percentage.
Validation was carried out only for tree species that were representative of the classes of interest. The areas for validation were chosen with a percentage equal to 5% of the total area for the tree species based on the sample plots and a deep knowledge of the territory. The validation results are shown in Table 3.  The classification validation was done through the confusion matrix [97][98][99], which allows one to determine the overall accuracy of the classification in addition to some parameters related to user and producer accuracy. The presence of validation is necessary to define the effectiveness of the classification and determine the error percentage.
Validation was carried out only for tree species that were representative of the classes of interest. The areas for validation were chosen with a percentage equal to 5% of the total area for the tree species based on the sample plots and a deep knowledge of the territory. The validation results are shown in Table 3.
The producer accuracy indicates the completeness of classification, while the user accuracy indicates the classification's correctness [100]. The Hellden parameter represents the mean accuracy of individual classes. The mean accuracy for class i is calculated using the equation presented in [101]: where A is the number of correctly classified reference points for class i, B is the total number of reference points in class i in the reference data, and C is the total number of reference points classified into class i [102]. KIA is the well-known Kappa Index of Agreement (or Cohen's kappa coefficient), which measures the proportion of agreement after chance agreements have been removed from consideration. The kappa increases to one as the chance agreement decreases and becomes negative the less often that chance agreement occurs [103].
The overall accuracy is the ratio of the sum of diagonal values of the confusion matrix to the total number of cell-counts in the matrix. This value gives an idea of the accuracy for all the considered classes [104].

SfM Approach for H Estimation
The results showed that the average H was very significant for each of the species (Table 4). The average discrepancies between the photogrammetric and laser scanning estimated lengths were in the order of a few centimetres, which proved the applicability of the SfM approach for tree H estimation. In most cases, marker placement on the treetop was largely responsible for the increased average error. For instance, for taller trees-e.g., Populus nigra-the estimated mean H was 29.15 m, with an average error of 0.3 to 0.4 cm. In other cases, the average errors were lower, except for the overlapping treetops. Unfortunately, there are few studies on urban tree H estimation utilising SfM approaches with AM software, so the general European Forest Tree (EFT) database published in 2016 [105] was considered for this study to provide a possible comparison. Even though this database was released for forestry areas, comparing the results with this database was still interesting. For instance, for the tree species such as Acer campestre, Quercus spp., and Populus nigra, the estimated average heights were 10.14, 16.15, and 29.15 m, respectively (Table 4). However, according to the EFT database, the average height for Acer campestre is typically 15 m and 30 and 40 m, respectively, for Quercus spp. and Populus nigra [106][107][108]. Other recent studies also estimated the average H at 18.83 m for Quercus spp., and, for the Populus nigra, the H was determined as 23 ± 5 m in Spain and Sweden [109,110]. The differences between the H values of the trees in this study and the H values in the other studies are possibly due to the following reasons: • These other studies were done in forestry areas under different environmental conditions; • In urban areas, the sizes and shapes of the tree crowns are modified at regular intervals; • In this study, most of the trees were located in the streets, where frequent tree growth does not occur.

Urban CS Mapping and Validation
The CS for the dominant species was calculated and mapped in QGIS showing the intensity of CS based on the species in different zones ( Figure 9). Figure 9 shows that the CS map estimated based on the classification shapefile was relevant for all the identified dominant species. For the trees that were not shaded or overlapped by the other dominant species, the computed NDVI and CS mapping was significant. Among the five dominant species, the species Quercus spp., Tilia platyphyllos, and Platanus spp. were found to be most responsible for the atmospheric CS covering the areas with higher CS intensity (Red, green, and brown zones at TC in Figure 10). Only in the case of Populus nigra (the blue and pink zones at CS in Figure 10) did the values indicate moderate to low CS zones. Among the five dominant species, the species Quercus spp., Tilia platyphyllos, and Platanus spp. were found to be most responsible for the atmospheric CS covering the areas with higher CS intensity (Red, green, and brown zones at TC in Figure 10). Only in the case of Populus nigra (the blue and pink zones at CS in Figure 10) did the values indicate moderate to low CS zones.
The outcome of the CS mapping for the dominant species proved the reliability of the applied methodology, even for complex urban areas ( Figure 11). Figure 11 shows the outcome of CS mapping for the dominant species in the whole urban area represented as a box plot graph. Generally, higher values (see the median values for each species) were found for the species Acer campestre, Platanus spp., Populus nigra, and Quercus spp. in the parks. Maximum values were generally found for Acer campestre (park) and Platanus spp. (street), while lower values were found for Quercus spp. (street), as in previous studies [111]. Among the five dominant species, the species Quercus spp., Tilia platyphyllos, and Platanus spp. were found to be most responsible for the atmospheric CS covering the areas with higher CS intensity (Red, green, and brown zones at TC in Figure 10). Only in the case of Populus nigra (the blue and pink zones at CS in Figure 10) did the values indicate moderate to low CS zones. The outcome of the CS mapping for the dominant species proved the reliability of the applied methodology, even for complex urban areas ( Figure 11). Figure 11 shows the outcome of CS mapping for the dominant species in the whole urban area represented as a box plot graph. Generally, higher values (see the median values for each species) were found for the species Acer campestre, Platanus   For the CS mapping validation, seven randomly selected plots were considered. Considering all the dominant species, four plots from the streets and another three plots were selected from the parks. The results show that the CS estimations in all cases had a smaller difference with the QGIS-computed CS values, except for plots 5 and 7, which showed higher differences in their estimations ( Table 5). The NDVI-based QGIS computations were found to be more applicable for street trees, which were more evident in the WV3 image than the park trees. Smaller species found in parks (i.e., Acer campestre) showed higher CS values and were overlapped by the other dominant species. On the other hand, Populus nigra, a taller tree species, presented lower CS values and was also found to be an isolated species with higher CS values during the field estimations.  For the CS mapping validation, seven randomly selected plots were considered. Considering all the dominant species, four plots from the streets and another three plots were selected from the parks. The results show that the CS estimations in all cases had a smaller difference with the QGIS-computed CS values, except for plots 5 and 7, which showed higher differences in their estimations ( Table 5). The NDVI-based QGIS computations were found to be more applicable for street trees, which were more evident in the WV3 image than the park trees. Smaller species found in parks (i.e., Acer campestre) showed higher CS values and were overlapped by the other dominant species. On the other hand, Populus nigra, a taller tree species, presented lower CS values and was also found to be an isolated species with higher CS values during the field estimations.
A regression analysis was also performed to determine the CS computational efficiency during the mapping. In the case of the street trees, the coefficient of determination was greater than 80% (R 2 = 0.85), and the overall R 2 value was 0.42 when considering all the validation plots, except for the plot with "Populus nigra" (Figure 12). The Populus nigra showed a large difference in CS estimations because of its typical fastigiate structure and status as a true "Lombardy poplar" with a very narrow crown [112]. Due to its narrow crown, this species was not as evident and provided a lower NDVI value, leading to a lower computation of CS. However, in other cases, the park tree species (e.g., Acer campestre) were mixed and adjacent to the other dominant species and had significantly higher computed NDVI values. Due to being surrounded by other larger crowns, the small-structured species presented higher NDVI values, which is why Acer campestre also showed a larger difference in the CS estimation, as Acer campestre is a shorter tree species than the others. On the other hand, the species on the streets (i.e., Platanus spp. and Tilia platyphyllos), were found to be isolated and planted with adequate spacing. They were also more evident in the WV3 image, showing smaller differences in CS estimations ( Table 5). As the validation plots were selected at random, only two plots out of the seven plots yielded lower R 2 values, which was true only for those plots because, for the overall CS computation, both species (Acer campestre and Populus nigra) were found to be efficient in atmospheric CS like all the other dominant species (Figure 11). Thus, the proposed methodology is more efficient (R 2 = 0.85) for dominant urban tree species (Figure 11), except for those with overlapping/narrow-crowned (R 2 = 0.42) structures in a complex urban environment ( Figure 12). seven plots yielded lower R 2 values, which was true only for those plots because, for the overall CS computation, both species (Acer campestre and Populus nigra) were found to be efficient in atmospheric CS like all the other dominant species (Figure 11). Thus, the proposed methodology is more efficient (R 2 = 0.85) for dominant urban tree species (Figure 11), except for those with overlapping/narrowcrowned (R 2 = 0.42) structures in a complex urban environment ( Figure 12).

Discussion
The OBIA approach applied in this study was found to be an efficient approach for tree species classification. The overall accuracy of the thematic mapping was 78%, with a kappa index of 0.74. For the different classes of species, the accuracy level varied from high (89% for the Platanus spp.) to low (57% for the Tilia platyphyllos). Usually, the selection of training samples for use with the NN algorithm had a significant influence on the accuracy of the classification [113]. The quantity of the selected samples must be adequate to allow the algorithm to distinguish the various classes. This is why we chose to distinguish between park trees and street trees. Ad hoc samples were chosen for the dominant species for both streets and parks. Then, two separate classifications were performed. In this way, the algorithm was able to show the best outcomes in both areas with higher accuracy. The spectral characteristics of the different species of trees highlighted some common traits of spectral responses that were too difficult to differentiate. That is why classification features such as NDVI, PCA, and GLCM were included in the feature space to obtain a classification with better results [47].

Discussion
The OBIA approach applied in this study was found to be an efficient approach for tree species classification. The overall accuracy of the thematic mapping was 78%, with a kappa index of 0.74. For the different classes of species, the accuracy level varied from high (89% for the Platanus spp.) to low (57% for the Tilia platyphyllos). Usually, the selection of training samples for use with the NN algorithm had a significant influence on the accuracy of the classification [113]. The quantity of the selected samples must be adequate to allow the algorithm to distinguish the various classes. This is why we chose to distinguish between park trees and street trees. Ad hoc samples were chosen for the dominant species for both streets and parks. Then, two separate classifications were performed. In this way, the algorithm was able to show the best outcomes in both areas with higher accuracy. The spectral characteristics of the different species of trees highlighted some common traits of spectral responses that were too difficult to differentiate. That is why classification features such as NDVI, PCA, and GLCM were included in the feature space to obtain a classification with better results [47]. Consequently, the classification accuracy was almost 80%, which was also able to satisfy the requirements for this study.
The use of remote sensing techniques for urban tree species mapping ensures significant time savings compared to the traditional species mapping of the territory [114,115]. To increase classification performance, images acquired at different time periods could be used to exploit the distinction between evergreen and non-evergreen plants.
In the case of the urban areas, such as Sassuolo, as in all Italian cities, there are numerous species of trees, historic gardens, and vegetation of various types. The OBIA approach was able to distinguish trees from other types of vegetation, but differentiating between many diverse species was not be possible. The common spectral characteristics of many species will make it more difficult to produce classification with higher accuracy. It is, therefore, imperative to identify dominant species to reduce errors in classification and to effectively focus on the tree species with higher density in the area.
The SfM methodology is efficient and feasible, especially for urban tree metrics, such as dendrometry [53,59]. The applied instruments, moreover, are consumer-grade and easy to use. The graduated bar ensures the real dimensions of the detected object with minimal effort on the field. Field operations are also simple and convenient and do not require any special skills for execution. The generated model is not georeferenced but was also found to be applicable for the estimation of dendrometry in tree CS calibration.
The results obtained from the applied methodology were validated for cases with similar settings. The laser scanning technique was used as a reference for validation, as this technique provides a direct method for 3D reconstruction, allowing direct measurements of each point's position. The choice of an anthropic object for validation of the proposed SfM approach facilitated identification of the homologous segments to enable comparison between the two methodologies ( Table 2). This analysis revealed discrepancies of about 4 cm, which means that the accuracy of the applied SfM approach is in the order of a few centimeters, which is suitable for the purposes of the proposed investigations. However, the quality of the 3D reconstruction could be improved by using small UAVs (Unmanned Aerial Vehicles) during the acquisition of images.
The methodology implied in this study will have significant impacts on CS calculations [116] and mapping for the dominant species in city areas. The mapping results in this study show that the species that are most evident (i.e., located on the streets), and not shaded or overlapping, have more significant CS values based on the field estimations. The application of remote sensing for dominant tree species identification made this methodology more convenient in terms of time investment and expenses than traditional tree identification ( Figure 11). The proposed methodology will not only provide an efficient way to understand the contributions of different species in atmospheric CS but will also provide other significant information, such as the following: • An estimation of the total benefits of individual dominant tree species considering the total predicted carbon sequestration; • Understanding the consequences of ensuring adequate plant spacing when planning urban green areas; • Identifying the most efficient dominant species based on their roles in urban ecology to better utilize the available urban space.
However, further studies should be done with more sample plots for field sampling and validation. CS mapping largely depends on the accuracy of tree classification, which could be challenging in densely vegetated urban areas. This is why the present study only considered the dominant species. This study was performed to understand the perspectives and opportunities of the applied methodology.
It would be interesting to determine the CS mapping efficiency in larger city areas with a vast quantity of diverse tree species.
However, proper understanding of the urban tree species contributions in atmospheric carbon sequestration and storage remains one the most pertinent issues for urban green planning and management [117,118]. A convenient mapping method is a prerequisite not only for tree CS prediction but also to understand the roles of the dominant tree species to improve the urban environment while considering other ecological benefits. The CS mapping considered in this study could also be used for larger city areas where it remains a significant challenge to develop a map without visiting and marking each tree manually. This is why the mapping efforts utilized here provide a more efficient way to implement proper urban green management and tree plantation systems to provide an improved urban atmosphere for city dwellers.

Conclusions
While CS mapping remains time-consuming and entangled with the traditional tree by tree method, this study provides a strong baseline from which to proceed with tree species classification and CS mapping in urban areas. This study will be helpful for future studies of applied approaches like SfM for geometric parameter estimation and the OBIA classification approach. For instance, WV3 images with eight bands (i.e., the coastal band, yellow band, red edge, and NIR2 band) are far better for dominant species classification in urban areas than traditional four-band images (usually using the blue band, green band, red band, and NIR band). For estimations of the accuracy level, urban tree species classification outcomes largely depend on the positions, crown structures, and spectral attributes of the trees. In this study, tree species classification was obtained with an overall accuracy value of 78%. The distinction between the park and street trees, moreover, improved classification accuracy. Furthermore, the CS map reveals that the tree stand levels and species distribution might have a significant influence on the total atmospheric CS for each species in urban areas. Aside from a few cases, the proposed methodology would provide an adequate solution for the efficient CS mapping of dominant urban tree species. This study was designed to assess and explore the feasibility of the proposed CS mapping approach in a typical urban area to demonstrate this method's efficient application compared to traditional approaches. The present method will help researchers and city planners implement advanced methods of CS mapping for typical urban areas to ameliorate the unavoidable impacts of climate change. This study, however, is not conclusive. This approach will next be implemented in a larger city area. Indeed, an experiment is currently underway in the capital city of Belgium that will compare its results with those of another approach (LiDAR data-derived CS mapping outcomes) to develop a method that is more convenient and less resource-intensive in all cases. That will enhance the feasibility of the approach even in the case of a significantly different context, such as in the case of a larger Asian metropolis.