Mapping the Urban Atmospheric Carbon Stock by LiDAR and WorldView-3 Data

: Currently, the worsening impacts of urbanizations have been impelled to the importance of monitoring and management of existing urban trees, securing sustainable use of the available green spaces. Urban tree species identiﬁcation and evaluation of their roles in atmospheric Carbon Stock (CS) are still among the prime concerns for city planners regarding initiating a convenient and easily adaptive urban green planning and management system. A detailed methodology on the urban tree carbon stock calibration and mapping was conducted in the urban area of Brussels, Belgium. A comparative analysis of the mapping outcomes was assessed to deﬁne the convenience and efﬁciency of two different remote sensing data sources, Light Detection and Ranging (LiDAR) and WorldView-3 (WV-3), in a unique urban area. The mapping results were validated against ﬁeld estimated carbon stocks. At the initial stage, dominant tree species were identiﬁed and classiﬁed using the high-resolution WorldView3 image, leading to the ﬁnal carbon stock mapping based on the dominant species. An object-based image analysis approach was employed to attain an overall accuracy (OA) of 71% during the classiﬁcation of the dominant species. The ﬁeld estimations of carbon stock for each plot were done utilizing an allometric model based on the ﬁeld tree dendrometric data. Later based on the correlation among the ﬁeld data and the variables (i.e., Normalized Difference Vegetation Index, NDVI and Crown Height Model, CHM) extracted from the available remote sensing data, the carbon stock mapping and validation had been done in a GIS environment. The calibrated NDVI and CHM had been used to compute possible carbon stock in either case of the WV-3 image and LiDAR data, respectively. A comparative discussion has been introduced to bring out the issues, especially for the developing countries, where WV-3 data could be a better solution over the hardly available LiDAR data. This study could assist city planners in understanding and deciding the applicability of remote sensing data sources based on their availability and the level of expediency, ensuring a sustainable urban green management system.


Introduction
To date, rapid urbanization intensely poses the need for greener landscapes in many urban areas worldwide. Green spaces allow maximizing urban resilience and livability and to positively respond to climate change effects. While cities are striving for more green space, more than half of the earth's population is already living in cities, and by 2050, 66% will be city dwellers [1]. Overexploitation of environmental resources for the huge population is indeed increasing the vulnerability of the urban dwellers to natural hazards. To keep pace with rapid urbanization, efficient urban green planning could be nothing but a time being and an expeditious solution. Conservation and expansion of existing urban vegetation based on their structural and functional roles in an urban atmosphere are some of the most effective factors of green urban planning. Thus, various approaches based on advanced technologies have been implemented to assess the contributions of urban trees, especially evaluation of their roles in atmospheric Carbon Stock (CS) is being increasingly acknowledged [2,3]. Trees in city streets and parks are now being recognized as a key tool against impacts caused by the increased rate of atmospheric carbon dioxide (CO 2 ) concentrations [4][5][6][7], since they sequester atmospheric carbon during the whole growth process and at the same time delay the adverse effects of climate change contributing to the accumulation of carbon in the soil [8][9][10]. Studies found that the total yearly reduction in carbon emission can be up to 18 kg/tree in urban areas [11][12][13], which clearly brings out the importance of planting trees along with having an efficient tree management policy, especially in a complex city environment. Trees directly impact atmospheric CO 2 fixation through photosynthesis, but in urban areas, the process is quite fitful due to tree health issues. As it is well known that the well-grown trees store far more carbon than the poorly grown ones, in urban areas, it is a huge challenge to maintain and preserve mature trees and well-managed urban forests that also include tree plantations and replacements. Therefore an efficient and timewise monitoring approach is essential to introduce an adequate urban tree management system [14,15]. In addition an effective monitoring system could be ensured utilizing an accurate and convenient species-based CS mapping approach. Most of the CS calibration and predictive models are based on the estimation of Above Ground Biomass (AGB) production [13,[16][17][18][19][20], which is considered to be primarily responsible for the atmospheric CS [10,[21][22][23]. In this study, the AGB was estimated based on the tree allometric information (i.e., Height (H), Diameter at Breast Height (DBH)) collected during the field surveys.
Currently, remote sensing-based mapping has been availed as an influential approach in monitoring functional and structural urban tree features to policymakers [24][25][26][27][28][29][30][31]. In fact, spatially extracted information on tree species and habitats over large areas are significant in understanding species' roles, such as in providing ecosystem functions and services [32][33][34][35]. Over the last few decades, remote sensing-based classification of tree species has been widely utilized either in the case of mapping specific species-based ecosystem services (ES) outcomes (i.e., [36]), or growth and yield models and, etc. (e.g., [33,37,38]). Remote sensing approaches, especially hyperspectral imagery, have significantly upgraded the tree classification outcomes either in single trees or mixed populations [33,[39][40][41][42][43]. The utilization of very high spatial resolution multispectral satellite imagery (e.g., 1-m IKONOS, 0.6-m QuickBird) and aerial photos/digital imagery has been rapidly increased, especially in spatial mapping [44][45][46][47]. As a matter of fact, recently, with the advancements of remote sensing technologies, a diversified type of very high resolution remotely sensed images (such as WorldView-3, WV-3) are commercially available, certainly introducing a wave of opportunities for the accurate mapping of urban trees at a significant level [31,33,[48][49][50][51][52]. Moreover, in the case of this study, a high-resolution WV-3 image has been successfully utilized to classify the dominant tree species in Brussels, which has been found useful for further CS mapping as it was earlier in the case of Sassuolo [36].
Additionally, there are also many convincing applications of Light Detection and Ranging (LiDAR) based calibration of the tree CS utilizing the individual tree metrics (i.e., [53][54][55][56][57][58][59][60][61]). On the other hand, much less evidence is available in the CS calibration of the urban trees utilizing only the multispectral satellite data [11]. In the case of urban areas, tree species mapping is still a considerable challenge due to having spatially heterogeneous land cover types from isolated trees to the dense forest, high tree species diversity along with heavily and regularly managed trees, as well as the interruptions by buildings and their shadows [7,[62][63][64][65][66]. Considering these facts, Geospatial Object-Based Image Analysis (GEOBIA) has been utilized in this study to classify the dominant tree species. However, it is yet a crucial concern to dig out the most convenient and compatible ways to map and predict the urban tree CS in a specific urban area. A method could be considered convenient in various ways, such as its application, time consumption, and execution expenses. Even LiDAR application is the most acceptable and widely reliable, it is still expensive and hardly cost-effective for the more significant part of the world. Therefore, it would be a timely consideration to analyse the utilization of multispectral satellite data (i.e., Sentinel-2, WorldView-2/3/4) regarding CS computation possibilities of the trees in an urban area. A remote sensing-based biomass assessment has been employed in many studies [10,[67][68][69]] to obtain forest information over a large area at a reasonable cost with acceptable accuracy and minimal effort [70]. It is also evident that the method of determining relationships between field estimations and remote sensing data-derived variables and then extrapolating these relationships over large areas is very useful [10,[71][72][73][74][75][76][77][78]. Here, the main goal of this study was to map the urban tree CS based on field measurements and the application of remote sensing tools considering the following: • A comparative analysis of the application of two different remote sensing data sources (i.e., LiDAR and WV-3 image data) regarding CS mapping in the case of dominant urban tree species; • To recommend an approach in the case of CS mapping for policymakers involving urban green management.
In a word, this study has been done to provide a fundamental tool considering urban CS mapping, which is one of the most critical issues for sustainable urban green management systems and their policymakers.

Field Data
The study area, covering an area of around 49 km 2 in the eastern part of the capital region in Belgium (Figure 1), was selected considering the availability of airborne LiDAR.
As the main goal of this study was to identify only the dominant urban tree species, all other non-dominant vegetations have been excluded during sampling. The sample plots were randomly selected, covering only the streets of the whole study area. Since in the parks most of the cases of tree crowns were overlapped and or completely overshadowed by the other species. That is why overcrowded tree populations were excluded to avoid misinterpretations of the species dominance information, during the final CS mapping. During field data collection, 75 plots (yellow dots in Figure 2) of 100 m 2 (10 m × 10 m) each were selected throughout the study area following the well-known Simple Random Sampling (SRS) approach [79][80][81][82]. Only the areas with tree species dominance (i.e., woody or tall perennial plants) have been considered, excluding the ornamental herbs, shrubs, or grassland areas. The sample plots were also used during the training and validation of the tree species for GEOBIA classification. Among the 75 plots, 20 plots (red square boxes in Figure 2) were considered for the CS mapping and validation. The Diameter at Breast Height (DBH) was measured for each tree in the plot. The height (H) of trees was measured utilizing the hypsometer Nikon Forestry 550, a laser rangefinder with angle compensation technology optimized for forestry use [83], and a field computer was used to mark the plots on QGIS. Field data (H, DBH) were collected in the summer of 2019.

Remote Sensing Data
The airborne LiDAR dataset had been collected in Summer 2015 by Aerodata Surveys Nederland BV [84]. The Crown Height Model (CHM), i.e., the height of objects obtained through the difference between the digital surface model (DSM) and the digital terrain model (DTM), was produced at a spatial resolution of 0.25 m using the LAStools software [84].
The WV3 image data for this study was acquired on 17 April 2017 (Figure 2), which provides one panchromatic band of 0.3 m and eight multispectral bands at 1.2 m spatial resolution ( Table 1). The available WV-3 image data have been pan-sharpened at the initial stage. Pan-sharpening is the process of merging high-resolution panchromatic and multispectral imagery where the outcome is an image that has the high spectral resolution

Remote Sensing Data
The airborne LiDAR dataset had been collected in Summer 2015 by Aerodata Surveys Nederland BV [84]. The Crown Height Model (CHM), i.e., the height of objects obtained through the difference between the digital surface model (DSM) and the digital terrain model (DTM), was produced at a spatial resolution of 0.25 m using the LAStools software [84].
The WV3 image data for this study was acquired on 17 April 2017 (Figure 2), which provides one panchromatic band of 0.3 m and eight multispectral bands at 1.2 m spatial resolution ( Table 1). The available WV-3 image data have been pan-sharpened at the initial stage. Pan-sharpening is the process of merging high-resolution panchromatic and multispectral imagery where the outcome is an image that has the high spectral resolu- tion of the multispectral image and also the high spatial resolution of the panchromatic image [85][86][87]. The pan-sharpen process was conducted using the hyperspectral colour sharpening (HCS) algorithm that combines the high-resolution panchromatic data with lower resolution multispectral data and specifically implemented for the WV imagery [88]. The cubic convolution resampling technique was chosen to resample the multispectral image to the high-resolution image using a 4 × 4 pixel moving window. The pan image and the multispectral image have been separately orthorectified using the DSM downloaded from the GeoPunt portal [89] and the available LiDAR data. All the steps were conducted in an ERDAS Imagine environment [90]. Later, the required shapefiles (for classification in Ecognition) were educed from the UrbIS database (UrbIS P& B and UrbIS-Adm), a general GIS database of the Brussels region [91].

Geospatial Object-Based Image Analysis (GEOBIA) Classification
Tree species classification based on spectral properties is quite a matter of contention due to the high intra-class spectral heterogeneity and/or inter-class spectral similarities [7,44,[64][65][66]. The traditional pixel-based procedures, which classify each object only based on a distinct spectral signature (ignoring other spatial/contextual information) [48,92,93], are hardly capable of reaching an acceptable accuracy [94,95]. In our study, GEOBIA was applied, which is already approved as an efficient technique in the case of high-resolution image classifications [48,85,92,[96][97][98][99]. Instead of the single pixels, the GEOBIA approach typically works on: (i) Segmenting a remote sensing image into spectral similarities (i.e., segments or objects), and (ii) Evaluating the spectral, spatial, and/or context features of these segments for image classification [48,63,94,100,101]. As a result of segmentation, this approach can consider the textual and contextual information along with the spectral information during the classification [48,94,[102][103][104]. Therefore GEOBIA classification outcomes are more acceptable than those based on the existing traditional pixel-based approach [102,[105][106][107][108][109][110][111].
The Trimble eCognition Developer ® 9 platform (Trimble, Munich, Germany) [112] was utilized to classify the dominant tree species in Brussels (Figure 3). At the initial stage, the chessboard segmentation has been introduced to initiate the GEOBIA classification approach. Then the primary classification was done utilizing the available shapefile from the Urbis database [91]. After that, the contrast split segmentation algorithm was implied to define the green and non-green areas with higher accuracy. This approach is based on the Normalized Difference Vegetation Index (NDVI) layer, where the 'contrast split' segments the scene into the dark and bright image objects based on a threshold value that maximizes the contrast between them [113]. The algorithm utilizes the optimal threshold separately for each image object which initiates a chessboard segmentation of variable scale and then performs the split on each square [114,115]. It was found quite helpful to identify the shadows, pathways, and pavements between the tree crowns. scale and then performs the split on each square [114,115]. It was found quite helpful to identify the shadows, pathways, and pavements between the tree crowns. Subsequently, the multiresolution segmentation (MRS) algorithm [116] was performed to group contiguous pixels into areas (i.e., segments) geometrically and radiometrically homogenous. The MRS algorithm was set by tuning the following parameters such as the "smoothness/compactness" that determines the preferred shape of segments, and the "colour/shape" parameter that controls the weights of spectral and shape information in the calculation of segments heterogeneity [48,94,117,118]. Considering the study of Choudhury et al. [36], green areas and streets have been identified with very large objects (as the size of the shapefile polygons) and have been classified separately from the rest. In these areas, a subsequent multiresolution segmentation had been applied to have the smaller objects. For this segmentation, rather than the thematic layers, the spectral information and the geometric information of WV-3 bands have been considered. The segmentation was done several times utilizing a different number of values for each parameter. For scale, the values were within 5 to 30, whereas for the shape and compactness, the Subsequently, the multiresolution segmentation (MRS) algorithm [116] was performed to group contiguous pixels into areas (i.e., segments) geometrically and radiometrically homogenous. The MRS algorithm was set by tuning the following parameters such as the "smoothness/compactness" that determines the preferred shape of segments, and the "colour/shape" parameter that controls the weights of spectral and shape information in the calculation of segments heterogeneity [48,94,117,118]. Considering the study of Choudhury et al. [36], green areas and streets have been identified with very large objects (as the size of the shapefile polygons) and have been classified separately from the rest. In these areas, a subsequent multiresolution segmentation had been applied to have the smaller objects. For this segmentation, rather than the thematic layers, the spectral information and the geometric information of WV-3 bands have been considered. The segmentation was done several times utilizing a different number of values for each parameter. For scale, the values were within 5 to 30, whereas for the shape and compactness, the values were between 0.1 to 0.9 [48,94,117,118]. After several attempts with different values, the ideal values utilized for segmentation, were found as scale parameter = 10, shape = 0.5 and the compactness = 0.8.
Before starting the Nearest Neighbour (NN) approach, an index known as Canopy Content Chlorophyll Index (CCCI) is used to separate the grasses from the vegetation class. This index can be calibrated as follows, Then the Nearest Neighbour (NN) algorithm [120] was performed, which is a supervised classification technique that classified all objects in the entire image based on the selected samples and the defined statistics [94]. For the sample selection and algorithm training, the sample plots had been considered. Once the algorithm training had been done, classification was initiated utilizing the "Assign class" algorithm. Three dominant tree species have been classified, such as Tilia spp. L., Acer spp. L. and Aesculus hippocastanum L. [121] covering the whole study area in Brussels. In this case, the larger and intensely green parks were not considered as it was hardly possible to differentiate crowns from a mixed or overlapped tree species population. Moreover, those trees could not be considered in further CS (AGB estimation) mapping due to the larger difference between the street and the park environment concerning tree health issues.
Validation of the classification outcomes (Table 2) had also been done at the Trimble eCognition Developer ® 9 platform [112] using confusion matrices [122][123][124] which is usually applied to compare the true classes with the ones assigned by the classifier on the generated maps. During the validation, 10%-15% of the total area for each class had been chosen as "true samples", as known species or classes, to train the algorithm. The estimated producer accuracy shows the completeness of classification, and the user accuracy indicates the correctness of the classes [36,125], while the Hellden parameter is used to estimate the mean accuracy for each class. The mean accuracy for each class i can be calculated using the equation presented in [36,126]: 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. In Table 2, KIA known as Kappa Index of Agreement (or Cohen's kappa coefficient), estimates the proportion of agreements [36,127]. Moreover, the overall accuracy (OA) has been estimated, which is the ratio of the sum of diagonal values of the confusion matrix to the total number of cell counts in the matrix [36,128].

Carbon Stock (CS) Mapping
Several studies suggest that tree AGB is the most visible, dominant, dynamic, and essential pool of the terrestrial ecosystem [8,[129][130][131], constituting around 30% of the total terrestrial ecosystem carbon pool [132]. In this study, based on the relationships among the WV-3 (NDVI) data and LiDAR (CHM) data derived variables and the field data, the predictive AGB estimation and, eventually, CS mapping has been done in both cases ( Figure 4). Several studies suggest that tree AGB is the most visible, dominant, dynamic, and essential pool of the terrestrial ecosystem [8,[129][130][131], constituting around 30% of the total terrestrial ecosystem carbon pool [132]. In this study, based on the relationships among the WV-3 (NDVI) data and LiDAR (CHM) data derived variables and the field data, the predictive AGB estimation and, eventually, CS mapping has been done in both cases (Figure 4). At first, the total AGB was calculated based on the field data (i.e., DBH, H, tree species, etc.) for each of the sample plots. For this calculation, an allometric model [133] was implied to calculate the AGB for each plot. The mean AGB/plot estimation was necessary as it is recommended that the tree above ground CS is assumed to be 50% of the total AGB [134][135][136][137][138]. Then to estimate the mean CS/plot, the mean AGB/plot was multiplied by 0.5 as a conversion factor [139][140][141]. Then in QGIS utilizing the WV-3 image data (from 2017, see Section 3.2.1 for details), the NDVI (Red edge and NIR1 band) of the whole study area was computed. The NDVI layer was considered in this study for the CS prediction and mapping, as previous studies claimed to find a strong correlation between the NDVI and total AGB of the trees [11,[142][143][144]. In the case of LiDAR data, the CHM was utilized to map the CS for the dominant tree species.
Then the NDVI-derived metrics were extracted for the sample plots utilizing the "Zonal statistics" plugin [145] at the QGIS interface. The CHM-derived matrices had been computed in the case of available LiDAR data (from 2015, see Section 3.2.2 for details). After that, the linear regression models were created in a Microsoft ® Excel™ spreadsheet, calibrating the correlation between the mean CS/plot and the NDVI derived metrics to find out the best model to estimate and map the CS covering the whole study area [36]. The CHM-derived matrices have also been done to determine the best model concerning the perspective CS mapping. A fishnet of 100 m 2 (10 × 10 m) resolution (as for the sample plots) was built in QGIS for both cases (NDVI and CHM) to recognize the minimum to maximum CS zones, based on the dominant species map (exported as a shapefile in QGIS) At first, the total AGB was calculated based on the field data (i.e., DBH, H, tree species, etc.) for each of the sample plots. For this calculation, an allometric model [133] was implied to calculate the AGB for each plot. The mean AGB/plot estimation was necessary as it is recommended that the tree above ground CS is assumed to be 50% of the total AGB [134][135][136][137][138]. Then to estimate the mean CS/plot, the mean AGB/plot was multiplied by 0.5 as a conversion factor [139][140][141]. Then in QGIS utilizing the WV-3 image data (from 2017, see Section 3.2.1 for details), the NDVI (Red edge and NIR1 band) of the whole study area was computed. The NDVI layer was considered in this study for the CS prediction and mapping, as previous studies claimed to find a strong correlation between the NDVI and total AGB of the trees [11,[142][143][144]. In the case of LiDAR data, the CHM was utilized to map the CS for the dominant tree species.
Then the NDVI-derived metrics were extracted for the sample plots utilizing the "Zonal statistics" plugin [145] at the QGIS interface. The CHM-derived matrices had been computed in the case of available LiDAR data (from 2015, see Section 3.2.2 for details). After that, the linear regression models were created in a Microsoft ® Excel™ spreadsheet, calibrating the correlation between the mean CS/plot and the NDVI derived metrics to find out the best model to estimate and map the CS covering the whole study area [36]. The CHM-derived matrices have also been done to determine the best model concerning the perspective CS mapping. A fishnet of 100 m 2 (10 × 10 m) resolution (as for the sample plots) was built in QGIS for both cases (NDVI and CHM) to recognize the minimum to maximum CS zones, based on the dominant species map (exported as a shapefile in QGIS) obtained from the WV-3 image data. The classification shapefile was essential to define the regions of interest, providing QGIS to map the estimated CS values considering only the dominant tree species [36]. Otherwise, the map will show the CS values for other areas, i.e., the area covered with grass or even in an area where there is no vegetation.
Then, to validate the mapping outcomes, 20 randomly selected plots have been utilized in both cases (NDVI and CHM). This time the linear regression models were created only for the validation plots. Therefore, the differences (Tables 3 and 4) among the QGIS computed CS values, and the ground truth-values were shown to recognize the methodology's effectiveness. The validation plots were the same in both cases, which was necessary to compare the impacts and to discuss the prospects and convenience of CS mapping in urban areas (see Section 4.2 for details).

GEOBIA Classification Results and Validation
As the WV-3 image was from the early spring, trees having dead leaves/branches made it relatively harder to identify the species in Brussels. Thus, species detection was difficult due to the spectral similarities and understory issues. However, the classified map shows that most of the streets are covered with the Acer spp. and Tilia spp., while the other species, Aesculus hippocastanum has been primarily observed in the southern part of the city ( Figure 5).

GEOBIA Classification Results and Validation
As the WV-3 image was from the early spring, trees having dead leaves/branches made it relatively harder to identify the species in Brussels. Thus, species detection was difficult due to the spectral similarities and understory issues. However, the classified map shows that most of the streets are covered with the Acer spp. and Tilia spp., while the other species, Aesculus hippocastanum has been primarily observed in the southern part of the city ( Figure 5). The producer and user accuracy were around 70% for Tilia spp. and Aesculus hippocastanum, and around 80% to 100% were for the Acer spp. during the confusion matrix estimation (Table 2). Indeed, it was quite essential to define the effectiveness of the GEOBIA approach and to show the percentage of error. The producer and user accuracy were around 70% for Tilia spp. and Aesculus hippocastanum, and around 80% to 100% were for the Acer spp. during the confusion matrix estimation (Table 2). Indeed, it was quite essential to define the effectiveness of the GEOBIA approach and to show the percentage of error.

CS Mapping in Brussels with WV-3 Image Data
The computed CS for the dominant species has been mapped showing the quantity of CS based on the species in different zones ( Figure 6). As far as most of the trees were on the streets, the canopies were not overlapped by the wider crowns. Except for some plots having trees almost leafless or trimmed canopies, the computed CS values based on NDVI derived variables were significant ( Figure 6) for most of the dominant species.

CS Mapping in Brussels with WV-3 Image Data
The computed CS for the dominant species has been mapped showing the quantity of CS based on the species in different zones ( Figure 6). As far as most of the trees were on the streets, the canopies were not overlapped by the wider crowns. Except for some plots having trees almost leafless or trimmed canopies, the computed CS values based on NDVI derived variables were significant ( Figure 6) for most of the dominant species. In Figure 6, all three-square boxes (as in CS1, CS2, and CS3) are showing moderate (i.e., pink, orange, yellow, and blue coloured zones) to a higher quantity of CS for each of the species (as in TC1, TC2, and TC3). These three zones were zoomed-in considering those plots having comparatively dense canopies, as the WV-3 image was acquired in early spring, most of the plots were with lighter crowns or leaves.
In Table 3, the validation results show that the field estimation and QGIS computed CS values have apparent variations in few cases. These kinds of differences between the field estimated and the computed values were quite natural in the case of this study. The field data were collected in the summer of 2019, while the WV-3 image was acquired in the early spring of 2017. QGIS computation was done based on the NDVI (WV-3 data) extracted values where most of the trees were leafless or less green. As a result, in the case of few plots, NDVI values were lower, which indeed showed lower CS values during the In Figure 6, all three-square boxes (as in CS1, CS2, and CS3) are showing moderate (i.e., pink, orange, yellow, and blue coloured zones) to a higher quantity of CS for each of the species (as in TC1, TC2, and TC3). These three zones were zoomed-in considering those plots having comparatively dense canopies, as the WV-3 image was acquired in early spring, most of the plots were with lighter crowns or leaves.
In Table 3, the validation results show that the field estimation and QGIS computed CS values have apparent variations in few cases. These kinds of differences between the field estimated and the computed values were quite natural in the case of this study. The field data were collected in the summer of 2019, while the WV-3 image was acquired in the early spring of 2017. QGIS computation was done based on the NDVI (WV-3 data) extracted values where most of the trees were leafless or less green. As a result, in the case of few plots, NDVI values were lower, which indeed showed lower CS values during the mapping. In addition, the year gap between the field data (2019) and the WV-3 data (2017) also significantly impacted the CS mapping outcomes, since trees in urban areas usually go through the management practices (i.e., trimming, pruning, etc.) [146], which is also responsible for the larger variations among the field and QGIS computed values (Table 3). For instance, out of the 20 validation plots, only three plots (plot no. 4, 11, and 18) are showing noticeably more CS/plot than those of the field estimations. Trees in those plots are assumed to be a subject of trimming and or other management practices, explaining the reasons for having lower AGB or CS/plot in 2019 than those of 2017 for those three plots. A regression analysis had been done to understand the significance of the mapping approach. It was quite noticeable that in most cases where the tree canopies were comparatively evident, even in the early spring season, the percentage of agreement (Figure 7) was more than 80% (R 2 = 0.83). Except for a few cases, this CS mapping approach could be applied in CS mapping, even in a complex urban environment. also significantly impacted the CS mapping outcomes, since trees in urban areas usually go through the management practices (i.e., trimming, pruning, etc.) [146], which is also responsible for the larger variations among the field and QGIS computed values (Table  3). For instance, out of the 20 validation plots, only three plots (plot no. 4, 11, and 18) are showing noticeably more CS/plot than those of the field estimations. Trees in those plots are assumed to be a subject of trimming and or other management practices, explaining the reasons for having lower AGB or CS/plot in 2019 than those of 2017 for those three plots. A regression analysis had been done to understand the significance of the mapping approach. It was quite noticeable that in most cases where the tree canopies were comparatively evident, even in the early spring season, the percentage of agreement (Figure 7) was more than 80% (R 2 = 0.83). Except for a few cases, this CS mapping approach could be applied in CS mapping, even in a complex urban environment.  Figure 8 shows that the computed CS map in QGIS was quite relevant in LiDAR data for all the identified dominant species (TC1, TC2, and TC3 in Figure 8). Even the LiDAR data was from 2015, each of the three dominant species had shown evident outcomes in CS mapping (CS1, CS2, and CS3 in Figure 8) based on the field data of 2019. As the LiDAR data (from summer 2015) and field data were acquired in summer, the AGB productions per plots were significant to predict the CS during the mapping based on the CHM-derived variables in QGIS.

CS Mapping in Brussels with LiDAR Data
For the CS mapping validation, the validation plots were the same ones as earlier in the case of WV-3 image data. The outcomes of the validation plots did not show that much difference between the field estimated and QGIS computed values (Table 4). Only a few plots were showing higher differences in the quantity of computed CS, which was assumed to be a result of tree crown management practices (i.e., trimming, pruning, etc.) [146]. Out of all 20 plots, only a few plots (i.e., plot no. 2, 9, 14, 15, and 17) were showing noticeable variations between the field and QGIS estimations (Table 4). While in most of the plots, trees were showing lower differences in the case of the total atmospheric CS considering the on-plot calibrated values.  Figure 8 shows that the computed CS map in QGIS was quite relevant in LiDAR data for all the identified dominant species (TC1, TC2, and TC3 in Figure 8). Even the LiDAR data was from 2015, each of the three dominant species had shown evident outcomes in CS mapping (CS1, CS2, and CS3 in Figure 8) based on the field data of 2019. As the LiDAR data (from summer 2015) and field data were acquired in summer, the AGB productions per plots were significant to predict the CS during the mapping based on the CHM-derived variables in QGIS.

CS Mapping in Brussels with LiDAR Data
For the CS mapping validation, the validation plots were the same ones as earlier in the case of WV-3 image data. The outcomes of the validation plots did not show that much difference between the field estimated and QGIS computed values (Table 4). Only a few plots were showing higher differences in the quantity of computed CS, which was assumed to be a result of tree crown management practices (i.e., trimming, pruning, etc.) [146]. Out of all 20 plots, only a few plots (i.e., plot no. 2, 9, 14, 15, and 17) were showing noticeable variations between the field and QGIS estimations (Table 4). While in most of the plots, trees were showing lower differences in the case of the total atmospheric CS considering the on-plot calibrated values.
A regression analysis was done to understand the relevance of the applied CS mapping approach. It was noticed that the R2 was significantly higher (84%), for most of the plots ( Figure 9). As shown earlier (Section 3.2.1), the data acquisition period (2015 and 2019) and the certain tree crown management practices (i.e., Trimming, pruning, etc.) [146] have also been found responsible for having larger variations in a few cases. As a result of regular crown trimming practices, the trees, even with a wider trunk, were showing smaller crowns in CHM during the CS computation in QGIS. A regression analysis was done to understand the relevance of the applied CS mapping approach. It was noticed that the R2 was significantly higher (84%), for most of the plots ( Figure 9). As shown earlier (Section 3.2.1), the data acquisition period (2015 and 2019) and the certain tree crown management practices (i.e., Trimming, pruning, etc.) [146] have also been found responsible for having larger variations in a few cases. As a result of regular crown trimming practices, the trees, even with a wider trunk, were showing smaller crowns in CHM during the CS computation in QGIS.  A regression analysis was done to understand the relevance of the applied CS mapping approach. It was noticed that the R2 was significantly higher (84%), for most of the plots ( Figure 9). As shown earlier (Section 3.2.1), the data acquisition period (2015 and 2019) and the certain tree crown management practices (i.e., Trimming, pruning, etc.) [146] have also been found responsible for having larger variations in a few cases. As a result of regular crown trimming practices, the trees, even with a wider trunk, were showing smaller crowns in CHM during the CS computation in QGIS.

GEOBIA Classification and Urban Trees
WV-3 image data were utilized to classify three dominant tree species applying the GEOBIA approach in the case of a heterogeneous urban area like as Brussels. The Overall Accuracy (OA) of the species classification was around 71%, considering all the identified urban tree species ( Table 2). The classification accuracy could be improved (around 90%) by applying the data fusion (i.e., LiDAR and hyperspectral) approach or estimating OA considering the general urban land cover classification [147][148][149][150][151][152][153]. Since this study's primary goal was to identify and classify only the dominant tree species for further CS mapping in QGIS, we show the OA only referring to the trees class and not for the other land cover classes such as roads, grasslands, and pavements. Unfortunately, studies on urban tree species or genera-based classification (only with multispectral data) with better accuracy (more than 70%) are hardly available unless the study area is quite smaller (i.e., [154][155][156]). Recently, Fang et al. [157] have done the tree genera-based classification (from WV-3 data) with an OA of 62 to 74% for a larger urban area in Washington D.C., USA. Considering a study area covering almost 49 km 2 , this study in Brussels also could be a good example for approving the efficacy of the GEOBIA approach. The accuracy level was different for each species, where the higher was for the Acer spp. (89%) and the lower one was found for the Aesculus hippocastanum (64%). Initially, the segmentation and later the relevant samples for training the NN algorithm significantly influenced the OA of the classification outcomes, as shown in previous studies [36,158]. Training samples must be sufficient considering each species/class, which will be enough for the algorithm to verify different classes. In this study, the samples had been selected to train the NN algorithm covering the whole study area in Brussels. However, the different species' spectral characteristics made it quite challenging to distinguish due to their unique phenological ages or conditions [147,159]. In some cases, misinterpretations of desired classes had been identified, recognizing the "mixed pixel issues" where some pixels were not solely covered by one homogeneous class [147,160]. Consequently, with an OA of 71%, this study does show that the WV-3 data could be more convenient, especially in urban areas, as the WV-3 image has an increased level of radiometric, geometric, and spatial (8 bands) resolution leading to classify urban trees at the species level [147,161].
Even remote sensing-based technologies are far better acceptable, the potential barriers to adopting remote sensing techniques include but are not limited to the level of specialist knowledge, the data collection time, instruments type and parameters for the acquisitions, etc. For instance, the classification OA could minimize such issues as the image registration error, the spectral and spatial resolution limitations, and overlapping spectral signature of the classes. These facts could be considered to initiate further research to improve and utilize the GEOBIA approach to achieve a species-level classification with higher accuracy, especially for the urban trees.

CS Mapping Approach: A Comparative Analysis
In this study, the mapping approach has been found far more convincing than that of the existing traditional approaches (i.e., approaches excluding remote sensing tools). The traditional biomass assessment methods are mostly based on field measurements which are not so convenient or practical to conduct over large areas considering a broadscale assessment [8,162]. Considering the dominant species, the CS mapping outcomes in Brussels were quite explicit and feasible, where the assessment was highly influenced by the seasonal variations between the data sources (field and remote sensing data). However, no such factors or issues could hardly be sorted out to improve the mapping outcomes in both cases. For instance, it is evident that in Brussels, the CS mapping outcomes could be more significant having both data sources (i.e., field and remote sensing data) from the same season or year. The CS computation utilizing the LiDAR data could be a good example since there was not much variation with the field estimations (Table 4). However, in all cases, there were other noticeable issues such as the narrow or too wide and/or small overlapped crowns, trimmed crowns, and so on other management practices. These issues are obvious to be considered precisely for any urban landscapes [7,63,64] that are heterogeneous in space, structurally, and functionally [44,[163][164][165]. However, the proposed remote sensing-based methodology could be availed as more essential and applicable in the case of monitoring and mapping vegetation ecosystems and their services than that of the traditional ones.
To date, it is critical to developing a high-resolution map of tree biomass within the context of carbon monitoring over terrestrial ecosystems to assess the ecosystem response to climate change [166][167][168][169][170][171], which increases the utilization of RS-based technologies for the last few decades. For instance, LiDAR data are widely being used to complete highresolution surveys of vegetation structure over forested areas and cities [53,[172][173][174][175] and estimate biomass and carbon storage in urban vegetation [26,53,176,177] due to its improved accuracy. Even though it is believed that the application of RS is quite expensive in vegetation mapping, Jones et al. (2010) [33,178] found automated methods from combined hyperspectral and LiDAR data (approximately USD 6 per ha) to be competitive against traditional aerial photograph interpretation (approximately USD 12 per ha) in terms of accuracy and cost for a study area in South-western Canada. However, the availability of LiDAR data is not always cost-effective, especially for developing countries. For those cases where LiDAR data are hardly available, high-resolution image data could be a solution for the city policymakers. For instance, studies show that high-resolution commercial data (Rapid Eye, IKONOS, Geo-eye) are available with an approximate cost of 1-14 € per km 2 , which is still more affordable than the LiDAR data (USD 62-240 per km 2 ) [179,180]. According to Hummel et al., the approximate cost of aerial LiDAR data acquisition for an area of almost 128 km 2 (31,614 acres) was around USD 40,500, excluding the processing costs (additional USD 33,424) [181]. While the WV-3 image data covering an area of 128 km 2 could be around USD 2500 (USD 19 km −2 ), considering the best possible buying options for archive imagery [182]. This could be an appropriate example of LiDAR data's application of almost 15 to 20 times higher than that of the WV-3 image data. Therefore, it could be a better initiative to imply other data sources, especially in developing countries, where it is impossible to integrate or have an available LiDAR data source. Considering the cost-efficiency and the research perspectives, this study has been conducted to identify the best outcomes in each case (LiDAR and WV-3 image data) and understand the result discrepancies among similar study areas. However, it is also clear that the mapping outcomes utilizing LiDAR data were more impressive ( Table 4) than that of the WV3 data. Nevertheless, the regression analysis still showed a significant level of acceptance (Figure 7) for the WV-3 image outcomes in Brussels. On the other hand, the percentage of LiDAR data outcomes except for only a few plots was more than 80% ( Figure 9). Moreover, city planners and managers usually plan to rely on a few urban tree attributes i.e., urban forest structure, green cover, species composition and diversity, available planting spaces, and tree condition to make short-and long-term decisions about the urban forest resource [15]. However, at the same time, the costs associated with data collection and monitoring (such as remote sensing applications) have to be compensated by reducing field measurement expenses or increasing management efficiency that leads to increased income based on improved decision-making [33]. Therefore the utilization and application of the advanced remote sensing technology depend on the purpose, economic feasibility, and the prospects of the resulted outcomes. Nevertheless, the proper understanding of the urban tree species contributions in atmospheric carbon sequestration and storage is one the most pertinent issues for urban green planning and management [183,184].

Conclusions
This mapping approach could be an efficient tool in CS mapping that will assist urban planners in ensuring proper utilization of the available green spaces considering other valuable prospects of tree species mapping in a complex city environment. This study will also recognize the prospects of the applied approach based on an efficient GEOBIA classification method for urban tree species mapping. For instance, this study shows that the OA of tree species classification could be hugely influenced by the trees' positions, crown structures, and spectral attributes, where the resulting outcomes were useful for further CS mapping in Brussels. The CS mapping approach reveals that the tree stands level and temporal variations of the data acquisition period might significantly impact the total atmospheric CS for each species in urban areas. 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 could be an efficient support for the urban planners and environmental policymakers in planning further urban air quality assessments. It also illustrates the facts considering the convenience and suitability of utilizing LiDAR and WV3 image data, especially in the case of vegetation mapping along with their functional attributes. This study could assist future research in either case of tree classification or mapping ecosystem services, having a significant prospect on remote sensing applications in the case of heterogeneous urban areas. Especially for the developing countries or where the application of LiDAR data is not that cost-effective, this mapping approach will show the alters with WV-3 image data. Soon another study will be introduced to compare the CS mapping outcomes between the WV-3 data and the Sentinel 3 data in the case of urban trees.
This study will contribute to a better understanding of the methodology in mapping structural and functional properties, such as tree CS, as well as predicting the possible urban CS in typical city areas. It might be a way out for the policymakers in mapping tree species as well as their probable ecological significance in urban areas over the traditional methods. Thus, the researchers and city planners could go forward to employ and implement the advanced ways of CS mapping for the typical urban areas assessing the possible prospects of the approach against the unavoidable impacts of climate change.