Nation-wide mapping of tree-level aboveground carbon stocks in Rwanda
Aerial imagesWe use publicly available aerial images of Rwanda at 0.25 × 0.25 m2 resolution, collected in June–August of 2008 and 2009. The images were acquired from 3,000 m altitude above ground level, originally with a mean ground resolution of 0.22 × 0.22 m2 pixel size then resampled to 0.25 × 0.25 m2, using a Vexcel UltraCam-X aerial digital photography camera34. The images exhibit a red, green and blue band stored under 8 bit unsigned integer format. The aerial images cover 96% of the country and the remaining 4% was filled with satellite images from WorldView-2, Ikonos, Spot and QuickBird satellite sensors which are part of the publicly available dataset.Environmental dataWe use locally available climate data: mean annual rainfall, mean annual temperature and elevation data (10 × 10 m2 resolution) to assess relationships between tree density, crown cover and environmental gradients. We also use land cover data to extract the spatial extent of plantations, forest, farmland, and urban and built-up areas for our landscape stratification. Climate data were obtained from the Rwanda Meteorological Agency as daily records from 1971 to 2017. The national forest map was manually created in 2012 using on-screen digitizing techniques over the 2008 aerial images35. A forest was defined as ‘a group of trees higher than 7 m and a tree cover of more than 10% or trees able to reach these thresholds in situ on a land of about 0.25 ha or more’51. A shrub was defined as ‘a group of perennial trees smaller than 7 m at maturity and a canopy cover of more than 10% on a land of about 0.25 ha or more’. The forest dataset was composed of 105,690 forest polygons, classified as either natural forest (closed natural forest, degraded natural forest, bamboo stand, wooded savanna and shrubland) or ‘forest plantations’ (Eucalyptus spp., eucalyptus; Pinus spp., pine; Callitris spp., callitris; Cupressus spp., cypress; Acacia mearnsii, black wattle; Acacia melanoxylon, melanoxylon; Grevillea robusta, grevillea; Maesopsis eminii, maesopsis; Alnus acuminata, alnus; Jacaranda mimosifolia, jacaranda; mixed species, mixed; and others) (Extended Data Fig. 7i). We separate shrubland from natural forest and merged it with savanna into the class ‘savannas and shrublands’. We further separated tree plantations and grouped them into Eucalyptus and non-Eucalyptus plantations. Then, a farmland map was acquired from the Rwanda Land Management and Use Authority (RLMUA)52 and overlaid with the 2012 forest cover map as a reference to clean the overlapping parts, under an assumption that the overlap is due to land use dynamics. Finally, a layer marking urban and built-up areas was acquired from RLMUA as well and the same preprocessing step as done for farmlands was applied. The combination of the land cover datasets resulted in our stratification scheme with six classes: natural forests, savannas and shrublands, Eucalyptus plantations, non-Eucalyptus plantations, farmland and urban and built-up.Mapping of individual trees using deep learningWe used the open-source framework developed by ref. 17 to map individual tree crowns. The framework uses a deep neural network based on the U-Net architecture53,54. We trained the network using 97,574 manually delineated tree crowns spread over 103 areas/bounding boxes representing the full range of biogeographical conditions found across Rwanda. To cope with the challenge of separating touching tree crowns, we used a higher weight for boundary areas between crowns, as suggested in refs. 17,53. Crown sizes in the predictions were found to be 27% smaller as compared to the manual delineations within the 103 training areas, due to the applied boundary weight that emphasizes gaps between tree crowns. Therefore, to calculate the real canopy cover, we extended each predicted tree crown by 27% and dissolved the touching crowns into continuous features. We counted single tree crowns for each hectare presented here as tree density and the percentage of each hectare covered by the extended tree crowns as canopy cover.We developed a postprocessing method that separates clumped tree crowns and fills any gap inside a single crown (Extended Data Fig. 2). Our postprocessing method, which we refer to as detect centre and relabel (DCR), determines the crown centres in the model predictions assuming that tree crowns have a round shape and then relabels the model predictions on the basis of weighted distances to the identified crown centres. First, DCR performs a distance transform, computing for each pixel the Euclidean distance to the nearest pixel predicted as background. Let the transformed image be distance-transformed (DT). Then an m × m maximum filter is applied to DT, where m depends on the size of the smallest object to be separated. We store all pixels for which the original DT value is the same before and after max-filtering. These pixels are the instance centres as they are furthest away from the boundary and have the highest distance values within the area defined by m. In the case of several connected instance centres in regions where multiple connected pixels have the same distance from the background, only a single instance centre is kept. Finally, each pixel x predicted as a crown in the original image is assigned to its nearest instance centre, where the distance function penalizes background pixels on the connecting line between the instance centre and x.Allometry for biomass and carbon stock estimationGenerally, allometric equations define a statistical relationship between structural properties of a tree and its biomass55,56. In our case, we assume a relationship between the crown area and aboveground biomass (AGB), which varies between biomes36. Since destructive AGB measurements are rare, we established biome-specific relationships between crown diameter (CD) derived from the crown area (CD = 2√(crown area/π)) and stem diameter at breast height (DBH) (equations (3) and (6)). DBH has been shown to be highly correlated with AGB36,37,38,39,40. We then used established relationships from literature to derive AGB from DBH for savannas and shrublands (equation (4)), tree plantations (equation (5)) and natural forests (equation (7)). AGB was predicted for each tree and summed for 1 ha grids to derive AGB in the unit Mg per ha. Values were multiplied by 0.47 (refs. 57,58) to derive aboveground carbon (AGC). Summed numbers over land cover classes are considered as carbon stocks. The bias as reported here was calculated following the approach from ref. 36 reporting the relative systematic error in per cent:$$mathrm {bias} = frac{1}{N}mathop {sum}limits_{i = 1}^N {frac{{(Y_{mathrm {obs}} – Y_{mathrm {pred}})}}{{Y_{mathrm {obs}}}}}times 100$$
(1)
The error for the evaluation with NFI data was defined by:$$mathrm{bias} = frac{{left| {mathop {sum}nolimits_N {(Y_{mathrm{obs}} – Y_{mathrm{pred}})} } right|}}{{left| {mathop {sum}nolimits_N {Y_{mathrm{obs}}} } right|}}$$
(2)
For trees outside natural forests, we used the database from ref. 36 including 10,591 field-measured trees from woodlands and savanna plus 952 samples from agroforestry landscapes in Kenya37 to establish a linear relationship between CD and DBH (Extended Data Fig. 3a). The Kenyan dataset is compatible with the trees in Rwanda. To ensure compatibility, the Kenya data contained open-grown trees most of which are of the same families or genus as in Rwanda grown under the same conditions, the latter factor shown to be important for generalizing37.A major axis regression (average of four runs each 50% of the data) led to equation (3):$${{{mathrm{DBH}}}}_{{{{mathrm{predicted}}}}},{{{mathrm{in}}}},{{{mathrm{cm}}}} = – 4.665 + 5.102 times {{{mathrm{CD}}}}$$
(3)
Equation (3) showed a reasonable performance with a very low bias (average of four runs on the 50% not used to establish the equation (3)): r² = 0.71; slope = 0.95; root mean square error (RMSE) = 6.2 cm; relative RMSE (rRMSE) = 42%; bias = 1%). We tested equation (3) on an independent dataset from Kenya consisting of 93 trees where AGB was destructively measured (Fig. 3b). The Kenyan database provides an uncommon opportunity to use destructive samples in which the carbon mass is not estimated indirectly and the relationship between crown area and carbon is direct: we do not need to invoke a second allometry to derive the dependent variable. All trees were open-grown trees in the same growing conditions as the agricultural areas of Rwanda. On these 93 trees, DBH can be predicted reasonably well from CD using equation (3) (r² = 0.84; slope = 0.86; RMSE = 8 cm; rRMSE = 25%; bias = 6%). We then applied an allometric equation from literature37 established for non-forest trees in East Africa to estimate AGB from DBHpredicted and compared the predicted AGB with the destructively measured AGB (r² = 0.81; RMSE = 511 kg; rRMSE = 55%; bias = 25%) showing an acceptable performance (Extended Data Fig. 3c) but indicating a systematic bias, which will be further tested with biome-specific field data (next section). We apply equation (4) to estimate AGB for trees outside forests in Rwanda in savannas and shrublands:$${{{mathrm{AGB}}}}_{{{{mathrm{predicted}}}}},{{{mathrm{in}}}},{{{mathrm{kg}}}} = 0.091 times {{mathrm{DBH}}_{{mathrm{predicted}}}}^{2.472}$$
(4)
Given the different structure of trees in farmlands, urban and built-up areas and plantations as compared to trees in natural forests and in natural non-forest areas, we used a different equation for trees in these areas. It was established in Rwanda using destructive samples from tree plantations39:$${{{mathrm{AGB}}}}_{{{{mathrm{predicted}}}}},{{{mathrm{in}}}},{{{mathrm{kg}}}} = 0.202 times {{mathrm{DBH}}_{{mathrm{predicted}}}}^{2.447}$$
(5)
A different CD–DBH relationship was established for natural forests. Here, we conducted a field campaign in December 2021 sampling 793 overstory trees in Rwanda’s protected natural forest. We measured both CD and DBH and established a logarithmic major axis regression model with a Baskerville correction59 between the two variables to predict DBH from CD (Extended Data Fig. 3d). We did four runs each using 50% of the data to establish equation (6) (average of the four runs) and the other 50% to test the performance also averaged over the four runs (r² = 0.71; slope = 0.99; RMSE = 13 cm; rRMSE = 45%; bias = 19%). Note that CD is extended by 27% to account for underestimations of touching crowns in dense forests (see previous section):$$begin{array}{l}{mathrm{DBH}}_{{mathrm{predicted}}},{mathrm{in}},{mathrm{cm}} = left({mathrm{exp}}left(1.154 + 1.248 times {mathrm{ln}}({mathrm{CD}} times 1.27) right)right.\left. times left({mathrm{exp}}(0.3315^2/2) right) right)end{array}$$
(6)
We then used a state-of-the-art allometric equation established for tropical forests38 to predict AGB from DBH for natural forests in Rwanda:$$begin{array}{l}{{{mathrm{AGB}}}}_{{{{mathrm{predicted}}}}},{{{mathrm{in}}}},{{{mathrm{kg}}}} = {{{mathrm{exp}}}}Big[ {1.803 – 0.976{{{E}}} + 0.976,{{{mathrm{ln}}}}left( rho right)}\+ 2.673;{{{mathrm{ln}}}}left( {{{{mathrm{DBH}}}}} right) – 0.0299left[ {{{{mathrm{ln}}}}left( {{{mathrm{DBH}}}} right)} right]^2 Big]end{array}$$
(7)
where E measures the environmental stress38 (a gridded layer is accessible via https://chave.ups-tlse.fr/pantropical_allometry.htm) and ρ is the wood density. Here, we used a fixed number (0.54), which is the average wood density for 6,161 trees from ref. 40, weighted according to the abundance of the species in the plots. The relative error was calculated by the quadratic mean of the intraplot and interplot variations, which is 18.2% (Extended Data Table 1b). No destructive AGB measurements were found that showed a similar CD–DBH relationship as we measured during the field trip in Rwanda’s forest. We could thus not evaluate the performance for natural forests at tree level but had to rely on plot-level comparisons (next section).Evaluation and uncertainties of the allometryBiomass estimations without direct measurements of height or DBH inevitably include a relatively high level of uncertainty at tree level38,60. Uncertainty does not only originate from the CD to DBH conversion but also the equation converting DBH to AGB. As shown in the previous section, no strong systematic bias could be detected for the CD to DBH conversion but the evaluation of the CD-based AGB prediction with an independent dataset from destructively measured AGB revealed a bias of 25%. However, this comparison (Extended Data Fig. 3c) may not be representative for an entire country having a variety of landscapes and tree species, so a systematic propagation is unlikely. We also did not have sufficient field data to evaluate the conversions in natural forests. Here, we used data from 15 natural forest plots with 6,161 trees published by ref. 40 and ref. 41 and directly compared the summed biomass of the trees we predicted over their plots. The median measured biomass for the plots is 121 MgC ha−1 and we predict a median biomass of 81 MgC ha−1 (plot-based rRMSE = 54%; bias = 11%; bias on summed plots = 26%). The overall underestimation by our prediction is not necessarily a model bias but may be partly explained by the contribution of the understory trees, which cannot be captured by aerial images. Interestingly, our C stock estimates are in the same range of magnitude as global biomass products43,44,45,61 (Extended Data Fig. 4), indicating that overstory tree-level carbon stock assessments are possible from optical very high resolution images, even in tropical forests. Several global products overestimated biomass for non-forest areas like savannas or croplands, which is probably because they are calibrated in denser forests. The most recent products of ref. 42 and ref. 61 are much closer to the estimates from our results and the NFI. This is also seen in the grid-based correlation matrix where ref. 42 correlates best with our map, followed by ref. 61.We further use NFI data from 2014 to measure the uncertainty of the final carbon stock estimates and evaluate if systematic differences between AGB predictions and field assessments can be found for different land cover classes (Extended Data Table 1). For the NFI data, a total of 373 plots with 2,415 trees were measured and species-specific allometric equations applied62. To identify systematic errors at landscape scale, we extracted averaged values for areas around the plots from our predictions and calculated statistics on averages over all plots. Interestingly, our predictions for farmlands only show a bias of 5.9%: we estimate on average 2.46 MgC ha−1 and the inventories measure 2.37 MgC ha−1 on their 150 plots. For savanna and shrublands, we estimate 4.16 MgC ha−1 while inventories measure 3.31 MgC ha−1 (bias = 18.9%). For plantations, we estimate lower values (8.16 compared to 16.79 MgC ha−1; bias = 52.6%). To calculate the total uncertainty on country-wide C stock estimates, we weighted the bias from the different classes according to their relative area. We estimate a total uncertainty on the carbon stock predictions of 16.9% at the national scale (Extended Data Table 1).We found a very low bias for estimated C density in farmlands (5.9% bias) which make up most of the areas outside natural forests in Rwanda (Extended Data Table 1, Extended Data Fig. 6). The high bias for plantations can be explained by three factors: large bare areas considered part of plantations by the manual delineation of plantation areas (Extended Data Fig. 1); regular harvesting and continual thinning which keep many plantation trees young and small; and the fact that our aerial images are from 2008 while plantation trees have grown until 2014 with a few new NFI plots initiated after 2008. The bias in savannas and shrublands can be explained by the following factors: the presence of multistemed trees with large crowns such as Acacia spp. and Ficus spp. among others; the fact that a crown-based method overestimates C stocks of shrubs with a small height; and presence of shrub trees with both small height and small (multiple) stems. If tree-level based carbon stock assessments derived from crown diameter as presented here should become standard to complement national inventories, a database with sufficient samples to evaluate for systematic errors needs to be established for each biome and inventory and satellite/aerial image-based methods need to be further harmonized.To further quantify the error propagation of the CD to DBH conversion for our application, we established four equations each randomly using 50% of the dataset and predicted the carbon stock for each tree in Rwanda with each equation. We did this separately for natural forests and trees outside natural forests. We calculated the rRMSE between the aggregated carbon stocks for each hectare. We averaged the rRMSE for each land cover class and show that the uncertainty for all classes does not exceed 5% (Extended Data Table 2a).Evaluation and uncertainties of tree crown mappingWe created an independent test dataset, which was never seen during training and was also not used to optimize hyperparameters. The test set consists of 6,591 manually labelled trees located in 15 random 1 ha plots (Extended Data Fig. 5). Thanks to the size of the country, the plots represent all rainfall zones and three major landscapes of the country. The plot-level comparison yielded very high correlations between the predictions and the labels and is shown in Extended Data Fig. 5. We also calculated a confusion matrix showing an overall per pixel accuracy of 96.2%, a true positive rate of 79.6% and a false positive rate of 6.8% (Extended Data Table 2b). Trees outside natural forests are easy to spot and count for the human eye, so we have confidence in the plot-based evaluation. However, it is often challenging in natural forests. Here, we used again the field measurements from 15 plots with 6,161 trees40,41. We find that we underestimate the total tree count by 22.6%, which may, at least partly, be explained by understory trees hidden by overstory trees and which are, therefore, not visible in our images. New field campaigns are needed to better understand and calibrate our results and possibly correct for systematic bias.Application and evaluation beyond RwandaWe acquired 83 Skysat scenes at 80 cm for Tanzania, Burundi, Uganda, Rwanda and Kenya. The model trained on the 25 cm resolution aerial images of Rwanda from 2008 was directly applied on the Skysat images. Forest and non-forest areas were manually delineated to decide which allometric equation to use for the carbon stock conversion. We randomly selected 150 1 × 1 km2 patches and aggregated the predicted carbon density per patch and compared the results with previously published maps42,43,44,45. Results show that the model can directly be applied to comparable landscapes on different datasets. Note, however, that accurate carbon stock predictions need local adjustments with field data. We then tested the tree crown model transferability on aerial images from California (NAIP; 60 cm) and France (20 cm) and found that the model delivers realistic results without any local training or calibration (Extended Data Figure 8).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More
