S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

E C O LO G Y Copyright © 2023 The

Authors, some
The overlooked contribution of trees outside forests to rights reserved;
exclusive licensee
tree cover and woody biomass across Europe American Association
for the Advancement
of Science. No claim to
Siyu Liu1*, Martin Brandt1*, Thomas Nord-Larsen1, Jerome Chave2, Florian Reiner1, Nico Lang3, original U.S. Government
Xiaoye Tong1, Philippe Ciais4, Christian Igel3, Adrian Pascual5, Juan Guerra-Hernandez6, Works. Distributed
Sizhuo Li1, Maurice Mugabowindekwe1, Sassan Saatchi7, Yuemin Yue8, Zhengchao Chen9, under a Creative
Rasmus Fensholt1 Commons Attribution
NonCommercial
Trees are an integral part in European landscapes, but only forest resources are systematically assessed by na- License 4.0 (CC BY-NC).
tional inventories. The contribution of urban and agricultural trees to national-level carbon stocks remains
largely unknown. Here we produced canopy cover, height and above-ground biomass maps from 3-meter res-
olution nanosatellite imagery across Europe. Our biomass estimates have a systematic bias of 7.6% (overestima-
tion; R = 0.98) compared to national inventories of 30 countries, and our dataset is sufficiently highly resolved
spatially to support the inclusion of tree biomass outside forests, which we quantify to 0.8 petagrams. Although

Downloaded from https://www.science.org at Delft University on September 19, 2023

this represents only 2% of the total tree biomass, large variations between countries are found (10% for UK) and
trees in urban areas contribute substantially to national carbon stocks (8% for the Netherlands). The agreement
with national inventory data, the scalability, and spatial details across landscapes, including trees outside
forests, make our approach attractive for operational implementation to support national carbon stock inven-
tory schemes.

INTRODUCTION example, reports from the United Kingdom have shown that
The quantification of wood resources at the national and local scales woody resources provided by trees outside closed canopy forests
is a prerequisite for sustainable management of timber resources, can exceed those provided by forests (17).
habitats, and carbon stocks (1–5). Local forest management practice Consistent monitoring of both forest and non-forest trees at the
and planning relies to an increasing extent on national maps of national and continental scales remains a challenging endeavor.
forests and forest resources, and the accuracy of national carbon This is because contemporary forest maps derived from medium-
stock estimates may be largely improved by inclusion of such resolution (10 to 30 m) satellite imagery usually miss non-forest
maps at the regional and national levels (6, 7). trees and do not allow for a clear definition and separation of
Most countries focus their inventories and mapping on closed forest and non-forest areas (18). Other maps merge all types of
canopy forests (8–10). Trees outside forests are often associated trees under the variable “tree cover” or “forest cover” (19, 20),
with scattered dryland trees in arid regions where rainfall is not suf- and it often remains unclear which size classes of trees and
ficient for trees to form closed canopies (11). However, in northern shrubs are included and which are excluded. Certain guidelines
countries, a substantial part of the woody resource may be growing provide consistent forest definitions, for example, by the United
in hedgerows, gardens, parks, urban areas, grasslands, and agricul- Nations Food and Agriculture Organization (FAO) (21).
tural lands. These trees outside forests contribute to carbon storage, However, the quality and type of data available at the national
provide resources for local communities, modify the local climate, and continental scales are often not sufficient to apply the guide-
are an important part of habitat networks, affect the hydrological lines consistently in remote sensing studies. Only few countries
cycle, and thus represent an important economic and social value provide a high-quality aerial imagery at the national scale that
(12–14). Many European countries comprise large agricultural and allows for an assessment of all forests and trees (22, 23), but the in-
urban landscapes, and the exclusion of trees outside forests from clusion of canopy height data derived from airborne Light Detec-
systematic carbon stock assessment potentially implies a bias in na- tion and Ranging (LiDAR) campaigns is decisive, to separate
tional inventories and also affects scientific studies and climate large trees from shrubs and bushes and to estimate the wood
models using only closed forest areas as input (15, 16). For volume, dry mass, and carbon stocks. This is particularly the case
for managed forest areas, where small and privately owned lots
1
Department of Geosciences and Natural Resource Management, University of Co-
characterized by different management strategies cause a high
penhagen, Copenhagen, Denmark. 2Laboratoire Evolution et Diversité Biologique, spatial heterogeneity within areas considered as closed canopy
CNRS, UPS, IRD, Université Paul Sabatier, Toulouse, France. 3Department of Com- forest (5). High–spatial resolution tree cover maps are not sufficient
puter Science, University of Copenhagen, Copenhagen, Denmark. 4Laboratoire des for a quantification of wood resources as the tree height is unknown.
Sciences du Climat et de l’Environnement, CEA/CNRS/UVSQ/Université Paris
Saclay, Gif-sur-Yvette, France. 5Department of Geographical Sciences, University Existing global- and continental-scale canopy height (24–26) and
of Maryland, College Park, MD, USA. 6Forest Research Center, School of Agriculture, medium-resolution biomass maps (27, 28) from space-borne
University of Lisbon, Lisbon, Portugal. 7Jet Propulsion Laboratory, California Insti- sensors largely ignore trees outside forests, providing an incomplete
tute of Technology, Pasadena, CA, USA. 8Key Laboratory for Agro-ecological Pro-
cesses in Subtropical Region, Institute of Subtropical Agriculture, Chinese assessment of resources. National canopy height data at a high
Academy of Sciences, Changsha, China. 9Airborne Remote Sensing Center, Aero- spatial resolution (<5 m) is expensive and often only available for
space Information Research Institute, Chinese Academy of Sciences, Beijing, China. parts of the country (29), if at all. It thus remains unknown to
*Corresponding author. Email: sliu@ign.ku.dk (S. Liu); mabr@ign.ku.dk (M.B.)

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 1 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

which extent trees outside forests across the European continent comparing two global canopy height products with airborne
represent a hidden and undervalued resource, as it is the case in LiDAR-derived canopy cover and height, the systematic bias was
the United Kingdom (17). high, particularly for trees outside forests: Here, the bias for
To answer this question and to explore the contribution of non- canopy height and cover was +15.4 and +231.2%, respectively, for
forest trees to tree cover and woody biomass across Europe, here, we a fused Sentinel-2– and spaceborne LiDAR [Global Ecosystem Dy-
generate continental-scale maps of forests, woodlands, and trees namics Investigation (GEDI)]–based map at 10-m resolution (26),
outside forests, including their height at 3-m resolution in 2019, while it was −31.3 and −95.8% for canopy height and cover, respec-
which contains a level of details otherwise only previously accessible tively, for a Landsat- and GEDI-based map at 30-m (Fig. 2D) (24).
from the use of airborne LiDAR surveys. We retrieve these data
from cost-efficient nanosatellites available at a daily basis and at a Quantification of tree cover outside forests across Europe
spatial resolution high enough to map large individual trees (>3-m To refine the forest versus non-forest area definition in relation to
height) using a deep learning approach where a convolutional our dataset, we used the FAO definition (21) to separate forests from
neural network is learned end to end from airborne LiDAR non-forest trees. Here, we aggregated our PlanetScope-based tree
canopy height reference data across Europe. We trained two deep cover and tree height maps into 0.5-ha grid cells. Following the
learning models from PlanetScope optical imagery at 3-m resolu- FAO definition, woody vegetation in 0.5-ha grids where the
tion: (i) a segmentation model predicting binary canopy/no- canopy cover of trees taller than 5-m height exceeded 10% and
canopy (>3-m height) and (ii) a regression model predicting where the land use was not agricultural land or urban according

Downloaded from https://www.science.org at Delft University on September 19, 2023

canopy height. Tree crowns and their shadows are visible in Planet- to (30) was considered as forest. Note that the FAO definition in-
Scope images, representing key features when training deep learn- cludes areas of trees that potentially can grow above 5 m, which
ing methods to retrieve canopy heights. Furthermore, we use cannot be resolved for our map. Conversely, trees in 0.5-ha grids
inventory data to convert tree cover and height into aboveground where the cover of trees taller than 5 m was below 10% or where
biomass maps including both forest and non-forest tree biomass, the land use was urban or agricultural land were considered as
forming the basis for quantifying and comparing forest and non- trees outside forests. Following these definitions, we find that
forest tree resources across European countries. We further apply forests cover 377 million ha, equivalent to 28% of the total land
the model trained in Europe to the boreal and temperate zones of area shown in Fig. 1. The canopy cover of trees outside forests rep-
North America to test the scalability of the approach to different resents 16.1 million ha, which is 1% of the total land area and 4% of
areas evaluated on independent datasets. the total tree cover. The precise percentage of cover due to trees
outside forests depends on a number of factors including land
cover maps and spatial resolution of the product (fig. S4).
RESULTS To analyze the spatial distribution of both forest and trees
A canopy height map at 3-m resolution for Europe outside forests, we aggregated our results at the country scale
We used 0.7 million km2 of canopy height data from airborne (Fig. 3, A and B) and to 1° units (Fig. 3C). Results show that
LiDAR available for Denmark, The Netherlands, Switzerland, almost half of all countries (44%) have a national forest cover
Estonia, Spain, Finland, and Wales (fig. S1B) to train a deep learning between 30 and 45%, and trees outside forests contribute less
model and predict tree height and binary tree cover from Planet- than 5% to the total tree cover for about half (53%) of the countries.
Scope imagery in 2019 at 3-m resolution for all Europe (10 Northern Europe is mostly covered by forest with limited trees cover
million km2; Fig. 1). The training data cover boreal, temperate, outside forests: Finland, Slovenia, and Sweden have the highest per-
and Mediterranean areas (fig. S1A). Model performance was evalu- centage of national forest cover (66, 63, and 61%, respectively) and
ated on 10,000 km2 of randomly sampled airborne LiDAR data (or- the lowest contribution of tree cover outside forests to the total tree
ganized as 1 km–by–1 km tiles), which was not used for model cover (1.8, 2.8, and 1.9%, respectively). On the contrary, the least
training and hyperparameter optimization (Fig. 2 and fig. S2). forested countries (United Kingdom, 12%; The Netherlands, 12%;
At the pixel level, canopy height was predicted with a low and Denmark, 15%) have the highest contribution of tree cover
random error for boreal [root mean square error (RMSE) = 4.33 outside forests to the total tree cover (22.1, 24.5, and 19.5%, respec-
m and relative RMSE (rRMSE) = 21.7%] and temperate zones tively). The highest contribution of trees outside forests at 1° grid
(RMSE = 5.06 m and rRMSE = 27.8%), but the diverse tree structure scale is found in Western Europe along the coast and Southwestern
of Mediterranean regions results in a higher uncertainty (RMSE = Europe (up to 40% in Fig. 3, C and E) where agricultural and resi-
6.44 m and rRMSE = 32.2%), observed as a more severe underesti- dential areas are concentrated and forest coverage is consequently
mation at higher canopies (>30-m height) (Fig. 2A), which can be low. The histogram results show a roughly even distribution of
partially explained by the uneven distribution of the validation data, forest cover percentage, with most of the trees outside forests ac-
with a lower number of points for higher canopies (Fig. 2B). Aggre- counting for less than 10% of the total tree cover (Fig. 3, D and
gating tree cover and height to 1 km–by–1 km grids, we found that E). For forest areas, our results compare well with the national
the overall systematic bias (see Materials and Methods for equation) FAO statistics (fig. S5). Differences can be explained by the slightly
between our PlanetScope predictions and LiDAR data for areas different definitions of forests. For example, the FAO definition in-
outside forests, here defined as trees that we mapped in the cludes areas with shorter trees that can grow taller than 5 m, such as
“urban” and “cropland” classes using a global land cover and land reforestation and afforestation, which were excluded in our map.
use (GLCLU) dataset (30), was −1.6% for canopy height and −2.6% We then quantified tree cover for different land use/cover classes
for canopy cover. For closed canopy forests, here defined as areas (Fig. 3F), by aggregating our data into 0.5-ha grids. For Europe as a
mapped as forest in the GLCLU dataset, the bias was −1.2% for whole, we identify 15.5 million ha of trees located outside the
canopy height and +0.4% for canopy cover (Fig. 2C). When “forest” class of GLCLU, in which half of them are found in

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 2 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Downloaded from https://www.science.org at Delft University on September 19, 2023

Fig. 1. PlanetScope-based canopy height map for Europe at 3-m resolution for 2019. (A) Spatial distribution of average canopy height aggregated to 1-km reso-
lution for visualization purpose. (B) PlanetScope imagery example in Switzerland displayed as near-infrared, red, and green false color composite. (C) Corresponding
canopy height models (CHMs) from airborne LiDAR at 0.5-m resolution. (D) Same area as (B) but showing the PlanetScope-based canopy height prediction at 3 m. The
land base maps in (C) and (D) are from Google Maps satellite imagery (Imagery 2022 CNES/Airbus, Landsat/Copernicus, Maxar Technologies, Map data 2022). The land and
ocean base map in (A) is from www.naturalearthdata.com.

“dense short vegetation” (hereafter referred to as grassland), with a Methods). We aggregated aboveground biomass to the hectare level
median (25 to 75 percentiles) canopy cover of 7.2% (2.1 to 19.2%). (100 m by 100 m) including both forest and non-forest trees (Fig. 4).
Urban areas have on average the highest canopy cover (12.7%; 4.5 to The overall uncertainty is the combined uncertainty due to the
26.9%), with an aggregated tree cover of 5.05 million ha. Tree cover canopy area and height prediction and to the height to biomass con-
of a total of 2.67 million ha is found in cropland with a rather low version and was quantified by comparing our final aboveground
canopy cover percentage per hectare (4.6%; 1.4 to 12.2%). Detailed biomass product (after application of the allometric conversion
statistics are found in table S1. on PlanetScope canopy height and cover data) with field measured
biomass from the Danish (n = 3451) and Spanish (n = 1706) NFI
Quantifying aboveground biomass across Europe data at the plot scale and the country scale (n = 30). At the plot scale,
To quantify the woody resources into aboveground biomass, we the correlation was moderate [correlation coefficient (r) = 0.53 and
used plot-based National Forest Inventory (NFI) data (fig. S13) bias = −23% for Denmark; and r = 0.50 and bias = −25% for Spain]
and airborne LiDAR canopy height model (CHM) data from (fig. S7, C and D), which was expected because of a large uncertainty
Denmark (31) to establish allometric relationships between above- related to both satellite and field data, such as geolocation and image
ground biomass and canopy height, both averaged at the plot level quality errors. These errors were found to be not systematic, as dem-
(n = 11,296). Separate relationships were derived for broadleaf forest onstrated by a comparison of statistics aggregated to the country
(n = 7232; bias = −8.8%), coniferous forest (n = 3768; bias = −6.8%), scale, both from NFI and satellites. Here, the systematic bias was
mixed forest (n = 7536; bias = −9.9%), and plots with sparse tree −10% (underestimation) for Denmark and +5% (overestimation)
cover (n = 409; bias = −5.9%) (fig. S6) using a previously published for Spain, which can be explained by the fact that clear-cut areas
forest-type map from 2018 (32) for the separation (see Materials and and young tree plantations are part of the national statistic but are

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 3 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Downloaded from https://www.science.org at Delft University on September 19, 2023

Fig. 2. Evaluation of canopy height and cover. (A) Residuals at 3 m by 3 m between aerial CHM and our PlanetScope-based predictions, grouped in 5-m height
intervals. (B) Number of evaluation points used in (A). See fig. S1 for biome zones (62). (C and D) Violin plots of canopy height and canopy cover for airborne LiDAR-
based CHMs, PlanetScope predictions (this study), a Sentinel-2–based product from (26), and a Landsat-based product from (24) aggregated to 1 km–by–1 km samples
(that were not used for training and parameter optimization) divided into forest areas (n = 4772) (C) and non-forest areas (n = 5228) (D). Forests and non-forest areas are
defined by the GLCLU dataset. Two-tailed Pearson correlation and bias are calculated between the airborne LiDAR-based CHM samples and the other products (all
aggregated to 1 km by 1 km). See corresponding scatter plots in figs. S2 and S3. The horizontal blue line is the median and the violins show the data distribution.
Pixel-level height distribution of trees outside forests is shown in fig. S12. See Materials and Methods for the bias calculation.

not mapped in our product. The bias was +7.6% over 30 countries 8.5) Mg/ha, and trees in grassland have 0.29 Pg with a median
with a Pearson correlation of 0.98 (Fig. 5D and fig. S9A). The rela- biomass density of 4.6 (2.1 to 11.1) Mg/ha. Overall, 44.6% of the
tively even distribution of the error across countries demonstrates tree biomass outside forests is found in grasslands, 34.0% in
the transferability of the approach beyond Denmark. A previously urban areas, and 16.2% in croplands (Fig. 5C).
published state-of-the-art biomass map from (28) had a larger bias At the country scale, we ranked the top five countries having the
of +17.3% when compared against country statistics. Aggregated to highest contribution of non-forest trees for different GLULCs to the
1 km–by–1 km scale, our map compares reasonably well with exist- national biomass stocks. Ireland was ranked first, with 16.5% of the
ing products (fig. S8) (27, 28, 33), with the advantage of including tree biomass located outside forests, followed by the United
information on non-forest trees (Fig. 4). Kingdom (14.8%), The Netherlands, which has 8.2% of its national
Following the FAO forest definition, European forests have an biomass located in urban areas, and Denmark (9.7%) (Fig. 5C). We
overall aboveground biomass of 35.2 Pg, and trees outside forests found a negative relation between the total tree cover percentage
represent only 0.8 Pg, which is a proportion of 40:1. By further ag- and the biomass contribution by trees outside forest (Table 1).
gregating our biomass maps into forest types and biomes, we found The Netherlands stands out with both the highest tree cover
that broadleaf and coniferous forests in the temperate zone have a outside forest (24.6%) and a high biomass contribution in percent-
biomass density of 140.4 (84.8 to 195.8) Mg/ha and 110.1 (50.4 to age (12.2%). France has the largest total biomass (2388.2 Tg), while
174.6) Mg/ha, respectively (Fig. 5, A and B). Using the global land the tree biomass outside forests is only of 77.2 Tg.
use and land cover class (GLULC) classification, our results reveal
that trees in urban areas have a total biomass of 0.22 Pg with a Beyond Europe
median biomass density of 5.6 (2.5 to 13.2) Mg/ha, trees in crop- To test the scalability of our approach, we applied the tree cover seg-
lands have 0.11 Pg with a median biomass density of 3.8 (1.9 to mentation and canopy height prediction to the boreal and

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 4 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Downloaded from https://www.science.org at Delft University on September 19, 2023

Fig. 3. Quantification of tree cover outside forests across Europe. (A) Total forest area in million hectares and the percentage of surface area at the country scale. (B)
Trees outside forest (TOF) area in million hectares and the percentage of total tree cover at the country scale. (C) Percentage of forest (size of the circle) and TOFs (color of
the circle) cover for 1° by 1° grids. (D) Forest cover percentage histogram at 1° by 1° as shown in (C). (E) TOF cover percentage histogram at 1° by 1° as shown in (C). (F)
Canopy cover of different GLCLU landscape types: The boxplots show the canopy cover percentage per hectare, and the bars denote the total canopy cover area per class.
For the boxplots, the start of the horizontal line represents the minimum value; vertical lines represent first quartile, median, and third quartile values, respectively; and
the end of the horizontal line represents the maximum value.

temperate zones of North America and tested the performance with DISCUSSION
a fully independent LiDAR dataset (1000 km2 at 1-m resolution) Trees outside forests in Europe have always been overlooked, and
that was never seen during training and parameter optimization only the United Kingdom has so far conducted a systematic assess-
(Fig. 6). Errors were found to be comparable to the European test ment of their resources (17). They conclude that, in many areas, the
datasets if aggregated to 1 km–by–1 km samples [coefficient of de- carbon stored in non-forest trees exceeds the forest carbon stocks,
termination (R 2) = 0.71, bias = −15.7%, and rRMSE = 20% for which would imply that current carbon stock assessments include a
height regression; and R 2 = 0.88, bias = +2%, and rRMSE = large bias.
28.6% for tree cover segmentation]. At 3 m by 3 m, the overall Here, we show that, at the continental scale, trees outside forests
RMSE was higher than in Europe (6.5 m). Note that a conversion do not play a critical role for the national aboveground carbon
to biomass would requires region-specific calibrations with local stocks of many Northern European countries. However, this is
field data, which is beyond the scope of this study. because of the dominant role of forest landscapes. We found the
total amount of carbon in non-forest trees (0.37 Pg) is about the

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 5 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Downloaded from https://www.science.org at Delft University on September 19, 2023

Fig. 4. Aboveground biomass estimation using PlanetScope. (A) PlanetScope-based canopy height map at 3-m resolution. (B) Biomass density estimation based on
(A) derived from allometry equation at 30 m by 30 m in kilograms per square meter. (C) Biomass density per hectare aggregated from (B) in megagrams per hectare. The
base map is from Google Maps satellite imagery (Imagery 2022 CNES/Airbus, Landsat/Copernicus, Maxar Technologies, Map data 2022).

same (0.30 Pg) as was previously found in billions of African

dryland trees over the arid and semiarid African Sahel, also with
Table 1. Country statistics of tree cover and biomass. Ten example comparable carbon density values (1.9 to 2.8 Mg/ha versus 1.5
countries sorted by the total tree cover percentage. Full country statistics Mg/ha in semiarid area) (34). This implies that trees outside
are provided in table S1. The class trees outside forests were defined
following the FAO definition.
forests are not only an overlooked resource in dryland area (11),
but they are equally important for European landscapes.
Country Total tree Tree cover Total Tree We also found that there are several countries and regions where
cover (% outside biomass biomass
of forests (% in Tg outside the tree cover outside forests exceeds 20% of the national tree cover
country) of total) forests in Tg (Fig. 3B). In particular, urban areas can have a relatively high tree
(% of total) cover, and their contribution to the national aboveground biomass
Ireland 10.5% 19.7% 36.65 3.24 (8.8%)
can be considerable. In addition, at the subnational scale, many
areas have a large proportion of trees outside forests (Fig. 3C),
United 14.1% 22.1% 264.29 26.29 (9.9%)
Kingdom
and excluding them causes a substantial bias in local carbon
budgets. Moreover, trees outside forests are not only valuable for
The 14.8% 24.6% 47.32 5.79 (12.2%)
Netherlands
their carbon stock, and a systematic identification of their location
could be an integral part of annual monitoring and planning
Denmark 16.9% 19.5% 63.47 4.74 (7.5%)
schemes related to biodiversity, microclimate, habitats, landscape
Ukraine 20.1% 11.4% 1534.2 53.2 (3.5%) values, and hydrological cycles (12, 13).
France 34.6% 9.7% 2388.2 77.2 (3.2%) Our canopy cover, height, and biomass products provide a level
Germany 36.1% 10.0% 2264.3 65.4 (2.9%) of detail that was previously only possible with airborne data from
Italy 38.6% 7.5% 1501.3 37.6 (2.5%) LiDAR, which are costly and rarely available at the national level.
Airborne LiDAR-derived maps of woody resources are not available
Estonia 50.3% 3.0% 259.3 2.7 (1.0%)
at the continental scale, and spaceborne LiDAR from GEDI (35) is
Finland 67.8% 1.7% 1726.2 9.76 (0.6%) only available at a point scale and not over multiple years. Further-
more, GEDI and also the soon-to-be-launched BIOMASS mission

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 6 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Downloaded from https://www.science.org at Delft University on September 19, 2023

Fig. 5. Aboveground tree biomass across tree types, land use/cover and countries. (A) Total tree biomass for different forest and land use/cover types; the color
reflects the biomes. (B) Same as (A) but for biomass density. (C) Top five countries with the highest contribution of non-forest tree biomass to the national biomass stocks,
as well as distribution of biomass among land use/cover classes. (D) Country-level biomass comparison between PlanetScope-based estimations and FAO statistics (5).
Forests are defined by GLULC and divided into broadleaf and coniferous using a previously published forest-type map from 2018 for countries in the European Economic
Area (EEA) and the PROBA-V forest-type map (2019) for countries outside the EEA (Belarus, Moldova, Russia, Ukraine) (32). For more countries, see fig. S9.

do not cover the northern parts of the world (36). Previous studies not free of charge, it is cheaper than aerial imagery, making it fea-
combined sparse GEDI data with Landsat (24) or Sentinel-2 (26) to sible to be acquired annually at the national scales. The daily global
generate global-scale canopy height maps, and, although these coverage facilitates cloud-free mosaics of different periods of the
perform well in forests, they have a large bias in non-forest areas year, and, once a model is trained, it has potential to be applied
(Fig. 2, C and D). Moreover, the major strength of deep learning each year on an operational basis without further need of LiDAR
with 3-m resolution nanosatellite imagery lies in the fact that the data. Our example of North America demonstrates the geographical
model can learn features from the visible tree crown structure, robustness of our models over trained biomes, but canopy height
which is often less clear in Sentinel-2 (10-m) or Landsat (30- estimations in dry and tropical areas require further work. Third,
m) images. LiDAR and aerial imagery differ between countries with regard to
Our maps could be operationally integrated in national carbon spectral, spatial, and temporal resolution, making it difficult to
stock inventory schemes. First, the spatial resolution allows to iden- create biomass maps that are consistent between countries and
tify individual tree crowns and makes it possible to align the satellite years at the continental scale. Although this was not demonstrated
image with field plot data, which can be useful to reduce variances in this study, PlanetScope-based maps have a certain consistency in
when upscaling plot data on, e.g., carbon stocks to the national space and time (37, 38), allowing for large-scale and multiyear
scale. Using different canopy height-to-biomass conversions for biomass assessments. Fourth, integrating trees outside forests in na-
different forest and non-forest types account for the bias caused tional inventories is not only interesting from a carbon stock per-
by the forest understory, which is undetected by classic aerial spective but also allows to quantify a variety of ecosystem services
surveys. Aggregating our biomass maps to the national scales provided by these trees. A previous study used submeter LiDAR and
shows a good alignment with national NFI statistics (bias of 7.7% aerial images to quantify the tree cover outside forests in Denmark,
overestimation for 30 countries) (5), although NFI plot-level data yielding the same proportion as this study (20%) (39), confirming
were only available for Denmark. Second, while the imagery is the reliability of the PlanetScope-based map for the mapping of

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 7 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Downloaded from https://www.science.org at Delft University on September 19, 2023

Fig. 6. Canopy height estimation for the boreal and temperate zones in North America. (A) Canopy height aggregated to 1-km resolution. (B) Scatter plot of canopy
height from aerial CHM and PlanetScope predictions in test areas [red boxes in (A)] using canopy height averaged over 1 km–by–1 km sample sites (n = 1036). (C) Same as
(B) but for canopy cover percentage. (D to F) Zoom-in examples of canopy height predictions in urban (D), cropland (E), and grassland (F) areas at 3-m resolution. The land
base map in (D) to (F) is from Google Maps satellite imagery (Imagery 2022 CNES/Airbus, Landsat/Copernicus, Maxar Technologies, Map data 2022). The land and ocean
base map in (A) is from www.naturalearthdata.com.

non-forest trees. For the United Kingdom, we found the fraction of with higher error values in the Mediterranean zone (Fig. 2A and fig.
urban trees outside forest cover to be 20.2%, which is in general S3), where our map should be used with caution. The inclusion of
agreement with the 16.5% reported in the U.K. NFI (17). Moreover, additional LiDAR data from mixed vegetation types, as well as semi-
the overall forest cover percentage found for the United Kingdom arid and arid regions, will improve future versions of the model. Ac-
(12%) is in line with previous reports (17, 40, 41). knowledging these current shortcomings, recent development of
A drawback of the presented method is the requirement of large state-of-the-art methods from the field of machine learning com-
amounts of airborne LiDAR for training and validation, as well as bined with cost-efficient nanosatellite imagery with a spatial resolu-
NFI data for converting canopy cover and height into biomass, tion below 5 m opens a previously unknown research avenue toward
which now limits the applicability primarily to the Global North. improved monitoring of national tree resources, both in and outside
Although our map has shown robustness across Europe, the gener- forests. This could allow the systematic integration of trees outside
alization to North America yields comparable results when aggre- forests into local, regional, and global carbon budgets, national in-
gated at 1-km resolution, but the performance at the 3-m pixel level ventories, climate models, and carbon credit programs and improve
drops. Nevertheless, fine-tuning on some reference data from a pre- the management of tree resources in countries characterized by pre-
viously unknown region of interest will likely close the performance dominantly non-forest landscapes. While this study concludes that
gap. The conversion of height and cover into biomass can potential- the impact on national statistics of most European countries will not
ly be improved by more region-specific field data, but using allome- be marked due to the dominant role of forests, we also show that the
tric equations from the Spanish NFI data showed that the equations amount of carbon stored in European non-forest trees is compara-
established from the Danish NFI data are robust (see Materials and ble to dryland areas, where trees outside forests are dominating the
Methods). In semiarid regions, small trees form a complex matrix of landscapes (34).
low vegetation, representing more challenging conditions for accu-
rate predictions of the deep learning–based model. We confirm this

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 8 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

MATERIALS AND METHODS Denmark (0.4 m-resolution merged from 2016 to 2021, in total
Summary 43,000 km2), Estonia (1-m resolution from 2011 to 2017, in total
We used CHMs from airborne LiDAR data across Europe to train 45,000 km2), The Netherlands (0.2-m resolution from 2018, in
two deep learning models: a segmentation model predicting binary total 41,000 km2), and Spain (2.5-m resolution from 2008 to
canopy/no-canopy (for trees taller than 3-m height) and a regres- 2015, in total 500,000 km2). Other samples were available for
sion model predicting canopy height, both from PlanetScope patches in Finland (1-m resolution from 2019, in total 50,000
optical imagery at 3-m resolution. The canopy map was used as a km2), Switzerland (0.5-m resolution merged from 2017 to 2020,
mask, and canopy height was only predicted within areas identified in total 20,000 km2), and Wales (1 m-resolution from 2015, in
as tree cover. We included a GEDI- and Sentinel-2–based canopy total 10,000 km2) (fig. S1B). Some data were provided as CHM,
height map from (26) as an additional model input, in both the others as raw cloud points, which were converted to Digital
training and prediction processes, providing a priori knowledge Terrain Model (DTM) and Digital Surface Model (DSM) used to
on the height distribution of forest areas. We observed that this im- calculate CHMs via LAStools (43). CHM data represent the top
proved the underestimation of tall forest trees and that, in case of canopy height in the pixel area but, sometimes, include the height
temporal or spatial mismatches (e.g., clear-cuts and trees outside from buildings and short herbaceous and bush vegetation. Here, we
forests), the prediction followed the PlanetScope image overruling used the building footprints from Microsoft (44) and Open Street
the information from the GEDI/Sentinel-2 map. Here, it is impor- Map (45) to mask buildings. We further visually compared Planet-
tant to mention that the map from (26) severely overestimates trees Scope imagery and the corresponding CHM and decided for a 3-m

Downloaded from https://www.science.org at Delft University on September 19, 2023

in sparsely vegetated areas, which, however, did not affect our result height threshold to separate tree crowns from shrubs, bushes, and
because we applied a binary canopy map at 3 m as mask (fig. S4). We grasses. Woody plants below this threshold were not clearly visible
used Danish NFI data to establish allometric relationships between in the PlanetScope images.
plot-based aboveground biomass and height (>1.3 m) from aerial Training sample generation
CHMs averaged over the field plot. The relationship between To adjust airborne CHM with PlanetScope images, we reprojected
canopy height and biomass was then used to convert canopy the PlanetScope imagery to match the aerial data Universal Trans-
height and cover predicted from PlanetScope imagery into verse Mercator (UTM) projection and resampled the aerial CHM to
biomass at 30-m resolution (same as NFI plot), which was lastly ag- the spatial resolution of PlanetScope data (3 m) by reserving
gregated to the hectare level. We did not use the field measured maximum values (filling small gaps in forests in the aerial CHM,
height but the average vegetation height from the LiDAR data per which are not visible in PlanetScope imagery). We applied a weight-
plot because it comes closest to our predicted height. We also did ed sampling strategy to the aerial CHM data, according to collection
not use PlanetScope predicted height because of alignment uncer- year, mean height, tree cover percentage (area > 3 m), forest type
tainties with the field plots at 3-m resolution. The lack of publicly [coniferous or broadleaf (32)], and dominating land use/cover
available GPS-referenced field plots covering the entire continent classes, to control the training data distribution fitting into the
makes it impossible to train a model that can predict biomass di- model and to minimize possible mismatches between data
rectly from PlanetScope imagery. sources caused by acquisition time differences (fig. S11). The final
dataset contains 100,000 samples, each 328 × 328 pixels correspond-
Data ing to 1 km–by–1 km ground spacing, where 10% was used for val-
PlanetScope imagery idation, 10% as test dataset, and 80% for training the model.
We used PlanetScope imagery from Planet Labs available via a re-
search license. The images are available daily for the entire globe Mapping tree cover and height using deep learning
with approximately 3-m pixel resolution and four multispectral We used the U-Net architecture (46) with an EfficientNet-B4 back-
bands: red, green, blue, and near infrared. We generated yearly 1° bone (47) for both the tree cover segmentation and the tree height
by 1° tiles of mosaics from around 100,000 image scenes for the regression tasks (fig. S10). The modified U-Net has proven to be
study area from 2019. Each mosaic tile consists of roughly 100 able to delineate tree canopies in very high resolution imagery
high-quality cloud-free scenes and was selected within a time (11, 38, 48). EfficientNets have shown to perform well with small
window according to the Moderate Resolution Imaging Spectrora- parameter size by a predesigned network depth, width, and resolu-
diometer (MODIS) phenology product (42). The selected time tion. Because we need to deploy our model to the continental scale
window started 25 days before “senescence” and ends 10 days at 3-m resolution, we chose the EfficientNet-B4 as the trade-offs
before “midgreendown,” which reflects late summertime in between performance and predicting (training) time. A key objec-
Europe, where all woody vegetation still has green leaves and tive of our work was to assess tree resources outside forests as accu-
visible crown structure. For areas with frequent cloud cover, the rately as possible, and the crown size of such trees is often just a few
time window was progressively extended toward earlier dates until pixels in our PlanetScope imagery. In the decoder, we replaced the
the whole tile was fully covered with scenes. To reduce differences simple two convolutional layers in the U-Net by residual blocks
resulting from the mosaicking scenes collected from different (49), which have been demonstrated to improve the model perfor-
sensors and times, we applied a histogram matching to Landsat mance for segmenting small objects (50). Furthermore, spatial and
and Sentinel-2 surface reflectance images from the same date channel squeeze and excitation blocks were integrated in the
range to produce homogeneous mosaics (38). decoder to build attention mechanisms (51). The deep learning
Airborne LiDAR canopy height maps framework is built on an open-source segmentation library (52).
We collected publicly available airborne LiDAR CHM data in Mismatches between aerial CHM and PlanetScope data were
Europe and the United States, with full country coverage for sometimes observed in forested areas, when the images were not ac-
quired from the same date. These inconsistencies cumulatively

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 9 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

contribute to a substantial proportion of loss in the regression, re- relating the total aboveground biomass value per plot with the
sulting in an unstable training progress and large uncertainties in average canopy height derived from airborne LiDAR (for all areas
the tree height estimations. Therefore, we first segmented the tree with woody vegetation > 1.3 m, following the NFI definition) auto-
canopies, which were then used to mask the non-forested areas matically corrects for the bias from undetectable understory. The
for minimizing the regression loss (fig. S10). regression equation in the log-log space can be described as
Tree crown segmentation
We set a 3-m height threshold to the aerial CHM data and converted InðAGBDÞ ¼ α þ β � InðHÞ þ ε ð1Þ
it into a reference binary (tree or background) map. The four-band where AGBD is the aboveground biomass density (in kilograms per
input data from PlanetScope were augmented by random cropping, square meter), H is the average canopy height for woody vegetation
flipping, and brightness distortion and then normalized channel- taller than 1.3 m for the NFI plot area from the aerial CHM, α and β
wise based on training dataset statistics before feeding into the are the regression coefficients, and ɛ is normally distributed error
model. We used the focal Tversky loss (53) as loss function for term with zero mean and standard deviation σ, that is, ɛ ~ N (0,
the segmentation model, and the alpha and beta parameters in σ2). When back-transforming Eq. 1, we corrected for logarithmic
the Tversky index (54) penalizing false positives and false negatives bias using the Baskerville correction (58).
were set to 0.4 and 0.6, respectively. The focal parameter was set to Plot-level biomass data were from the Danish NFI (n = 11,296)
1.2, which incentivizes the model to focus on harder examples collected between 2012 and 2021 (31). Each plot was a circle with a
(when the Tversky index < 0.5), especially for the cases where scat- radius of 15 m, where averaged tree height, estimated aboveground

Downloaded from https://www.science.org at Delft University on September 19, 2023

tered trees or tree lines do not exactly overlap between PlanetScope biomass, and forest fraction were provided. To remove misaligned
and aerial CHM data. We used Adam (55) as optimizer with an plots, we filtered plots where the forest cover difference was more
initial learning rate of 0.0001 and then reduced the learning rate than 30% but kept data regardless of the collection year, assuming
by a factor of 0.1 when the metric has stopped improving for 30 that dynamics are averaged out by the large number of plots. We
epochs. Model training stopped when the learning rate reached 1 derived forest-type–specific relationships using the Copernicus
× 10−7. In our study, each epoch of training took approximately forest map from 2018 (32) to classify the plots in broadleaved
0.5 hour on a GeForce RTX 3090 GPU. The model achieved conver- forest (bf ), coniferous forest (cf ), mixed forest (mf ), and areas of
gence after 180 epochs. Subsequently, the trained model was then sparse tree cover (nf ), which are plots that fall outside of the areas
used to predict a binary tree/non-tree cover map at 3-m for the mapped as forest in the Copernicus forest-type map (fig. S6). We
entire study area. Notably, for each tile (1° by 1°), the prediction acknowledge that the class representing the non-forest trees (nf )
process took only 3 to 5 min. is not ideal, because the Danish NFI does not systematically
Tree height regression measure trees outside forests and other wooded lands. However,
We used the same network architecture and training strategy as for here, we assume that plots being located in areas of sparse tree
the segmentation task, with an inverse sample–frequency weighting cover and having a low biomass are most representative for non-
smooth–L1 loss function to mitigate the underestimation that is forest areas. Because the NFI biomass data include all trees taller
often observed for tall trees (26). The regression loss was only cal- 1.3 m and not only those visible from an aerial perspective and in
culated for areas that overlap between the segmentation tree mask the PlanetScope images (typically >3-m height and 10-m2 crown
and the aerial CHM (>3 m). While this largely eliminated uncer- size), the different relationships include understory or small trees,
tainties from mismatches between the two data sources, the which indirectly correct for the bias of missing small trees in the
model still underestimated taller forests, possibly resulting from PlanetScope-based tree height predictions. This fact makes it im-
unclear training data. For example, coniferous forests vary optically portant to have a separated equation for trees outside forests,
only little in PlanetScope image but can have different heights in the which typically have a lower amount of understory. The equations
CHM, depending on their location. As a consequence, whenever the and the number of NFI plots are as follows
model encounters this situation, it conservatively predicts interme-
diate values. One solution is to encode geo-information into the AGBDbf ¼ exp½ 2:427� � H1:794 � exp½0:1342 =2�
ð2Þ
model (50) or include as input data (26). However, this would ðN bf ¼ 7232Þ
require a more homogeneous distribution of the reference data.
We thus included a GEDI/Sentinel-2–based height map from (26)
as auxiliary data (additional input band) to provide prior knowledge AGBDcf ¼ exp½ 1:652� � H1:547 � exp½0:1522 =2�
to the model about the forested areas, resulting in a total of five ð3Þ
ðN cf ¼ 3768Þ
bands for the tree height regression task. For the regression task,
both the training and prediction times are comparable to tree seg-
mentation. However, the model required a longer time to converge. AGBDmf ¼ exp½ 2:243� � H1:728 � exp½0:1182 =2�
ð4Þ
ðN mf ¼ 7536Þ
Allometry for biomass estimation
Allometric equations provide aboveground biomass estimates from
related tree measurements, such as stem diameter at breast height, AGBDnf ¼ exp½ 2:227� � H1:742 � exp½0:3052 =2�
canopy height, or crown diameter (56). We associated plot-level ð5Þ
biomass derived from species-specific equations and field measured ðN nf ¼ 409Þ
stem diameters with airborne LiDAR CHM (>1.3 m) averaged over We used 20% of the plots randomly selected for each type to
the same plots (56, 57). The plot-based aboveground biomass in- derive error metrics shown in fig. S6. Overall, the equations
cludes trees of all sizes above 1.3 m including understory, and

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 10 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

showed a reasonable performance with a relative low bias (averaged model performance varied across biomes and canopy height
R 2 = 0.62, bias = −8.7%, and RMSE = 5.43 kg/m2). ranges yet was closely related to the data distribution (Fig. 2B). In
Plot-level aboveground biomass in kg was then estimated as general, the model did not show a considerable saturation for taller
canopies, yielding an RMSE of 5.4 m and an mean absolute error of
AGB ¼ ðTCP=100Þ � 900 � AGBD ð6Þ 4.2 m. We found similar metrics as for the segmentation across
where TCP is predicted tree cover percentage aggregated in 30 m– biomes at plot-level evaluation, with R 2 = 0.82, bias = −2.4%, and
by–30 m resolution (rasterized NFI plot size). AGBD is obtained rRMSE = 17.9% for the boreal zone; R 2 = 0.9, bias = −1.6%, and
from Eqs. 2 to 5 according to the dominant forest type within the rRMSE = 14.5% for the temperate zone; and R 2 = 0.45, bias =
plot area. The aboveground biomass at 30-m resolution was lastly −20.1%, and rRMSE = 35.4% for the Mediterranean zone (fig.
aggregated to the hectare level, containing both forest and non- S2D). All slopes below the diagonal line indicated a slight underes-
forest information (Fig. 4). timation, consistent with the pixel-level evaluation (Fig. 2A). As for
trees outside forests, we obtained a bias of −1.6% and rRMSE of
Evaluation 26.4% (fig. S3A). We further evaluated our top height estimation
We randomly selected 10% of the CHM dataset (10,000 1 km–by–1 with field measured plot scale height from Danish NFI data (n =
km samples) as independent test dataset, which was not used for 3451) (fig. S7B), showing an underestimation bias of 9%.
training and to optimize the model parameters during training, North America
neither for the tree cover segmentation nor for the tree height re- To test the generalization of our model, we then predicted tree cover

Downloaded from https://www.science.org at Delft University on September 19, 2023

gression tasks. We evaluated our results at pixel level (3 m by 3 and height for the boreal and temperate zone of North America and
m) using the intersection over union (IoU), which quantifies how compared the results to aerial CHMs that were not used for either
well the predicted boundaries align with the ground truth boundar- training or validation of the model parameters. Here, we acquired
ies, for the segmentation and the RMSE for tree height regression around 1000 km2 of airborne LiDAR CHMs at 1-m resolution from
and aggregated to 1 km–by–1 km grids using the relative systematic (59) and generated 1036 samples with 1 km–by–1 km ground
error (bias %, Eq. 7); positive bias represents an overall overestima- spacing following the same data processing procedure as described
tion, and negative bias represents an overall underestimation (56) in the previous section. These evaluations show that our model can
be deployed to a comparable landscape without a considerable loss
N � �
1 X Ypred Y true of performance (R 2 = 0.71, bias = −15.7%, and rRMSE = 20% for
bias ¼ � 100 ð7Þ
N i¼1 Ytrue height regression; and R 2 = 0.88, bias = +2%, and rRMSE = 28.6%
for tree segmentation) (Fig. 6).
and the rRMSE in percentage described in Eq. 8 Evaluation of the allometry
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi We used NFI data from Spain (fig. S13), mostly located in the region
u N �
u 1 X Ytrue Y pred �2 of Extremadura (768 plots) measured in 2016 and 2017, and the
rRMSE ¼ t � 100 ð8Þ
N i¼1 Ytrue province of Leon in the northwestern part of Spain (1162 plots),
measured in 2019. The Spanish NFI measurements represent an im-
Tree cover segmentation portant remeasurement campaign on which existing NFI plot loca-
Overall, our model achieved an IoU of 0.692 on the test dataset (n = tions have been remeasured with upgraded global navigation
10,000 plots of 1 km by 1 km). Further evaluation revealed a good satellite system equipment to improve the geolocation of plot
alignment between our tree cover map and the tree cover map from centers in existing and new plots (60). Airborne laser scanning
aerial CHM data (using a 3-m height threshold) with overall R 2 data were collected in Leon from October 2019 to May 2020 and
values of 0.92 but with variation across biomes (fig. S2E). Trees in Extremadura from October 2018 to July 2019, which largely
can be segmented more precisely in both boreal (R 2 = 0.97, bias = match with the NFI measurements. In addition, we used 1284
−1.0%, and rRMSE = 13.5%) and temperate zones (R 2 = 0.99, bias = plots that include a diverse mosaic of Mediterranean evergreen
−0.8%, and rRMSE = 18.5%) than in the Mediterranean zone (R 2 = and broadleaf forests (61). We implemented the same procedure
0.39, bias = −8.2%, and rRMSE = 89.8%), which is likely due to two as applied in Denmark, to derive allometric equations using only
reasons: First, a limited amount of training data was usable from the Spanish NFI data
Spain because of the large temporal gap between PlanetScope
imagery and the available aerial CHM data, which is from 2012. AGBDbf ¼ exp½ 3:1218� � H2:072 � exp½0:2892 =2�
ð9Þ
Second, the method of using a constant height threshold to generate ðN bf ¼ 2213Þ
the tree mask is not fully functional in diverse vegetation structures
consisting of trees and shrubs. For areas with trees outside forest,
our predictions compared with aerial CHM data achieve an AGBDcf ¼ exp½ 2:161� � H1:693 � exp½0:2942 =2�
ð10Þ
overall bias = −2.4% and rRMSE = 36.6% (fig. S3 and Fig. 2D). ðN cf ¼ 506Þ
We also calculated a confusion matrix comparing Danish NFI
plots (n = 13,638) classified as either forest or no forest with our
data following the FAO definition, showing a very high agreement AGBDnf ¼ exp½ 1:608� � H1:479 � exp½0:4462 =2�
with an overall accuracy of 0.96 (fig. 7A). ð11Þ
ðN nf ¼ 397Þ
Tree height regression
We evaluated our tree height model by a residual analysis between Using these equations, we estimate 1120.7 Tg biomass for Spain,
reference height from the test dataset and the corresponding pre- which is very close to the estimates using the equations from the
dicted tree height at pixel level (3-m resolution) (Fig. 2A). The Danish NFI data (1139.8 Tg). This consistency gives confidence

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 11 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

in the robustness and ability to generalize our allometric models, 2. G. Ceccherini, G. Duveiller, G. Grassi, G. Lemoine, V. Avitabile, R. Pilli, A. Cescatti, Abrupt
increase in harvested forest area over Europe after 2015. Nature 583, 72–77 (2020).
which were originally developed and validated in Denmark.
3. C. Senf, R. Seidl, Mapping the forest disturbance regimes of Europe. Nat. Sustain. 4,
Biomass estimations 63–70 (2021).
The errors from canopy cover, height, and the allometric equations 4. N. T. Hoang, K. Kanemoto, Mapping the deforestation footprint of nations reveals growing
propagate to the final biomass estimation, so we evaluated the un- threat to tropical forests. Nature Ecology & Evolution 5, 845–853 (2021).
certainty on the biomass maps in three ways. First, we compared the 5. “State of Europe’s forests 2020 report” (Forest Europe, 2020).
predicted biomass at field plot level with estimations from the 6. A. Kangas, R. Astrup, J. Breidenbach, J. Fridman, T. Gobakken, K. T. Korhonen, M. Maltamo,
Danish NFI for 3451 selected plots, resulting in an overall bias of M. Nilsson, T. Nord-Larsen, E. Næsset, H. Olsson, Remote sensing and forest inventories in
Nordic countries–Roadmap for the future. Scand. J. For. Res. 33, 397–412 (2018).
−23% (fig. S7C). Because NFI data were collected in forest areas,
7. S. Magnussen, N.-L. Thomas, R.-N. Torben, Lidar supported estimators of wood volume and
this error is affected by the underestimation of very tall trees, above-ground biomass from the Danish national forest inventory (2012–2016). Remote
where biomass densities can reach 300 Mg/ha (fig. S7C) and possi- Sens. Environ. 211, 146–153 (2018).
ble geolocation uncertainties between the PlanetScope images and 8. J. Breidenbach, A. Granhus, G. Hylen, R. Eriksen, R. Astrup, A century of National Forest
the NFI plot locations recorded by a handheld GPS. Second, we Inventory in Norway–informing past, present, and future decisions. For. Ecosyst. 7,
46 (2020).
compared our biomass estimation with other global biomass prod-
9. E. Tomppo, H. Olsson, G. Ståhl, M. Nilsson, O. Hagner, M. Katila, Combining national forest
ucts (27, 28, 33). The GEDI biomass product (33) does not cover our inventory field plots and remote sensing data for forest databases. Remote Sens. Environ.
entire research area, so we randomly sampled 4000 points (fig. S8), 112, 1982–1999 (2008).
resulting in a bias of 1.8%. The biomass map from (28) shows good 10. J. A. Blackard, M. V. Finco, E. H. Helmer, G. Holden, M. Hoppus, D. M. Jacobs, A. J. Lister,

Downloaded from https://www.science.org at Delft University on September 19, 2023

correlation with ours but with a higher bias (bias = 11.3% and r = G. G. Moisen, M. D. Nelson, R. Riemann, B. Ruefenacht, D. Salajanu, D. Weyermann,
0.83). The biomass map from (27) shows an overall underestimation K. C. Winterberger, T. J. Brandeis, R. Czaplewski, R. E. Mc Roberts, P. L. Patterson, R. P. Tycio,
Mapping U.S. forest biomass using nationwide forest inventory data and moderate reso-
(bias = −37.9% and r = 0.82). Third, we further compared the Plan- lution information. Remote Sens. Environ. 112, 1658–1677 (2008).
etScope-based biomass maps aggregated to the national scale with 11. M. Brandt, C. J. Tucker, A. Kariryaa, K. Rasmussen, C. Abel, J. Small, J. Chave, L. V. Rasmussen,
FAO country statistics (Fig. 5D and fig. S9A) (5): The overall bias P. Hiernaux, A. A. Diouf, L. Kergoat, O. Mertz, C. Igel, F. Gieseke, J. Schöning, S. Li, K. Melocik,
over all 30 countries was 7.6%. Our biomass estimates reached J. Meyer, S. Sinno, E. Romero, E. Glennie, A. Montagu, M. Dendoncker, R. Fensholt, An
similar levels but were slightly higher than the FAO statistics for unexpectedly large count of trees in the West African Sahara and Sahel. Nature 587,
78–82 (2020).
most countries, which could result from trees outside forests that
12. S. Schnell, C. Kleinn, G. Ståhl, Monitoring trees outside forests: A review. Environ. Monit.
are typically not included in national inventory reports. The total Assess. 187, 1–17 (2015).
bias for (28) and (27) at the country scale was calculated as +17.3 13. N. Thomas, P. Baltezar, D. Lagomasino, A. Stovall, Z. Iqbal, L. Fatoyinbo, Trees outside
and +49.9%. forests are an underestimated resource in a country with low forest cover. Sci. Rep. 11,
1–13 (2021).

Environmental and auxiliary data 14. J. Reed, J. van Vianen, S. Foli, J. Clendenning, K. Yang, M. MacDonald, G. Petrokofsky,
C. Padoch, T. Sunderland, Trees for life: The ecosystem service contribution of trees to food
We used a GLCLU map from 2019 at 30 m (30) including forests, production and livelihoods in the tropics. For. Policy Econ. 84, 62–71 (2017).
urban areas, cropland, and grassland to analyze trees outside forests 15. R. J. Zomer, H. Neufeldt, J. Xu, A. Ahrends, D. Bossio, A. Trabucco, M. van Noordwijk,
and trees in different landscapes. We also used global canopy height M. Wang, Global tree cover and biomass carbon on agricultural land: The contribution of
products based on Sentinel-2 from 2020 at 10-m resolution (26) and agroforestry to global and national carbon budgets. Sci. Rep. 6, 29987 (2016).
another product based on Landsat from 2019 at 30-m resolution 16. S. Schnell, D. Altrell, G. Ståhl, C. Kleinn, The contribution of trees outside forests to national
tree biomass and carbon stocks—A comparative study across three continents. Environ.
(24) for comparison. Both products use GEDI spaceborne canopy
Monit. Assess. 187, 1–18 (2015).
height data that were extrapolated with Sentinel-2 and Landsat, re- 17. “Tree cover outside woodland in Great Britain report” (National Forest Inventory, 2017).
spectively. The GEDI/Sentinel-2–based canopy height map was also 18. M. W. Thornton, P. M. Atkinson, D. A. Holland, Sub-pixel mapping of rural land cover
used to assist the height regression task. The Copernicus forest-type objects from fine spatial resolution satellite sensor imagery using super-resolution pixel-
map based on Sentinel-2 data from 2018 at 10-m resolution and a swapping. Int. J. Remote Sens. 27, 473–491 (2006).
PROBA-V–based forest-type map from 2019 at 100 m (32) were 19. D. Zanaga, R. Van De Kerchove, D. Daems, W. De Keersmaecker, C. Brockmann, G. Kirches,
J. Wevers, O. Cartus, M. Santoro, S. Fritz, M. Lesiv, M. Herold, N.-E. Tsendbazar, P. Xu,
used to separate forests into coniferous and broadleaved forest
F. Ramoino, O. Arino, ESA WorldCover 10 m 2021 v200 (2022); https://zenodo.org/
across Europe. We further used the GEDI L4B biomass density record/7254221.
product from 2020 available at 1-km resolution (33) and state-of- 20. M. C. Hansen, P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau,
the-art global biomass density products from 2017 (28) and 2020 S. V. Stehman, S. J. Goetz, T. R. Loveland, A. Kommareddy, A. Egorov, L. Chini, C. O. Justice,
(27) at 100 m for biomass comparisons. We used biome data J. R. G. Townshend, High-resolution global maps of 21st-century forest cover change.
Science 342, 850–853 (2013).
from (62) aggregated into three main biome zones (fig. S2).
21. FAO, FAOSTAT, “Food and agriculture organization of the United Nations, Rome,
Italy” (2018).
22. C. Alexander, P. K. Bøcher, L. Arge, J. C. Svenning, Regional-scale mapping of tree cover,
Supplementary Materials height and main phenological tree types using airborne laser scanning data. Remote Sens.
This PDF file includes: Environ. 147, 156–172 (2014).
Figs. S1 to S13 23. E. Malkoç, M. Rüetschi, C. Ginzler, L. T. Waser, Countrywide mapping of trees outside forests
Table S1 based on remote sensing data in Switzerland. Int. J. Appl. Earth Obs. Geoinf. 100,
References 102336 (2021).
24. P. Potapov, X. Li, A. Hernandez-Serna, A. Tyukavina, M. C. Hansen, A. Kommareddy,
A. Pickens, S. Turubanova, H. Tang, C. E. Silva, J. Armston, R. Dubayah, J. B. Blair, M. Hofton,
REFERENCES AND NOTES Mapping global forest canopy height through integration of GEDI and Landsat data.
1. J.-F. Bastin, Y. Finegold, C. Garcia, D. Mollicone, M. Rezende, D. Routh, C. M. Zohner, Remote Sens. Environ. 253, 112165 (2021).
T. W. Crowther, The global tree restoration potential. Science 365, 76–79 (2019).

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 12 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

25. N. Lang, N. Kalischek, J. Armston, K. Schindler, R. Dubayah, J. D. Wegner, Global canopy 46. O. Ronneberger, P. Fischer, T. Brox, “U-net: Convolutional networks for biomedical image
height regression and uncertainty estimation from GEDI LIDAR waveforms with deep en- segmentation,” in International Conference on Medical Image Computing and Computer-
sembles. Remote Sens. Environ. 268, 112760 (2022). Assisted Intervention (MICCAI) (Springer, 2015).
26. N. Lang, W. Jetz, K. Schindler, J. D. Wegner, A high-resolution canopy height model of the 47. M. Tan, L. Quoc, “EfficientNet: Rethinking model scaling for convolutional neural networks”
Earth. arXiv:2204.08322 (2022); https://doi.org/10.48550/arXiv.2204.08322. in Proceedings of Machine Learning Research (PMLR) (ML Research Press, 2019).
27. L. Xu, S. S. Saatchi, Y. Yang, Y. Yu, J. Pongratz, A. A. Bloom, K. Bowman, J. Worden, J. Liu, 48. M. Mugabowindekwe, M. Brandt, J. Chave, F. Reiner, D. L. Skole, A. Kariryaa, C. Igel,
Y. Yin, G. Domke, R. E. McRoberts, C. Woodall, G.-J. Nabuurs, S. de-Miguel, M. Keller, P. Hiernaux, P. Ciais, O. Mertz, X. Tong, S. Li, G. Rwanyiziri, T. Dushimiyimana, A. Ndoli,
N. Harris, S. Maxwell, D. Schimel, Changes in global terrestrial live biomass over the 21st V. Uwizeyimana, J.-P. B. Lillesø, F. Gieseke, C. J. Tucker, S. Saatchi, R. Fensholt, Nation-wide
century. Sci. Adv. 7, eabe9829 (2021). mapping of tree-level aboveground carbon stocks in Rwanda. Nat. Clim. Chang 13,
28. M. Santoro, O. Cartus, N. Carvalhais, D. M. A. Rozendaal, V. Avitabile, A. Araza, S. de Bruin, 91–97 (2022).
M. Herold, S. Quegan, P. Rodríguez-Veiga, H. Balzter, J. Carreiras, D. Schepaschenko, 49. K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition.
M. Korets, M. Shimada, T. Itoh, Á. Moreno Martínez, J. Cavlovic, R. Cazzolla Gatti, P. da arXiv:1512.03385 (2016); https://doi.org/10.48550/arXiv.1512.03385.
Conceição Bispo, N. Dewnath, N. Labrière, J. Liang, J. Lindsell, E. T. A. Mitchard, A. Morel, 50. W. Sirko, S. Kashubin, M. Ritter, A. Annkah, Yasser Salah Eddine Bouchareb, Y. Dauphin,
A. M. Pacheco Pascagaza, C. M. Ryan, F. Slik, G. Vaglio Laurin, H. Verbeeck, A. Wijaya, D. Keysers, M. Neumann, M. Cisse, J. Quinn, Continental-scale building detection from high
S. Willcock, The global forest above-ground biomass pool for 2010 estimated from high- resolution satellite imagery. arXiv:2107.12283 (2021); https://doi.org/10.48550/arXiv.
resolution satellite observations. Earth Syst. Sci. Data 13, 3927–3950 (2021). 2107.12283.
29. D. J. Luscombe, N. Gatis, K. Anderson, D. Carless, R. E. Brazier, Rapid, repeatable landscape- 51. A. G. Roy, N. Navab, C. Wachinger, “Concurrent spatial and channel ‘squeeze & excitation’in
scale mapping of tree, hedgerow, and woodland habitats (THaW), using airborne LiDAR fully convolutional networks,” in International Conference on Medical Image Computing and
and spaceborne SAR data. Ecol. Evol. 13, e10103 (2023). Computer-Assisted Intervention (Springer, 2018).
30. M. C. Hansen, P. V. Potapov, A. H. Pickens, A. Tyukavina, A. Hernandez-Serna, V. Zalles, 52. Pavel Iakubovskii, et al. Segmentaton Models Pytorch;

Downloaded from https://www.science.org at Delft University on September 19, 2023

S. Turubanova, I. Kommareddy, S. V. Stehman, X. P. Song, A. Kommareddy, Global land use https://github.com/qubvel/segmentation_models.pytorch (2019).
extent and dispersion within natural land cover using Landsat data. Environ. Res. Lett. 17, 53. N. Abraham, N. M. Khan, “A novel focal Tversky loss function with improved attention u-net
034050 (2022). for lesion segmentation,” 2019 IEEE 16th International Symposium on Biomedical Imaging,
31. T. Nord-Larsen, V. K. Johannsen, T. Riis-Nielsen, I. M. Thomsen, B. B. Jørgensen, “Skovstatistik Venice, Italy, 8 to 11 April 2019).
2019: Forest statistics 2019” (Institut for Geovidenskab og Naturforvaltning, Københavns 54. A. Tversky, Features of similarity. Psychol. Rev. 84, 327–352 (1977).
Universitet, 2020). 55. D. P. Kingma, J. Ba, “Adam: A method for stochastic optimization,” Proceedings of the 3rd
32. European Environment Agency, Copernicus Land Monitoring Service (2018). International Conference on Learning Representations (ICLR), San Diego, CA, 7 to 9
33. R. O. Dubayah, J. Armston, S. P. Healey, Z. Yang, P. L. Patterson, S. Saarela, G. Stahl, May 2015).
L. Duncanson, J. R. Kellner, “GEDI L4B Gridded Aboveground Biomass Density, Version 2” 56. T. Jucker, J. Caspersen, J. Chave, C. Antin, N. Barbier, F. Bongers, M. Dalponte, K. Y. van Ewijk,
(ORNL DAAC, 2022). D. I. Forrester, M. Haeni, S. I. Higgins, R. J. Holdaway, Y. Iida, C. Lorimer, P. L. Marshall,
34. C. Tucker, M. Brandt, P. Hiernaux, A. Kariryaa, K. Rasmussen, J. Small, C. Igel, F. Reiner, S. Momo, G. R. Moncrieff, P. Ploton, L. Poorter, K. A. Rahman, M. Schlund, B. Sonké,
K. Melocik, J. Meyer, S. Sinno, E. Romero, E. Glennie, Y. Fitts, A. Morin, J. Pinzon, D. McClain, F. J. Sterck, A. T. Trugman, V. A. Usoltsev, M. C. Vanderwel, P. Waldner, B. M. M. Wedeux,
P. Morin, C. Porter, S. Loeffler, L. Kergoat, B. A. Issoufou, P. Savadogo, J. P. Wigneron, C. Wirth, H. Wöll, M. Woods, W. Xiang, N. E. Zimmermann, D. A. Coomes, Allometric
B. Poulter, P. Ciais, R. Kaufmann, R. Myneni, S. Saatchi, R. Fensholt, Sub-continental-scale equations for integrating remote sensing imagery into forest monitoring programmes.
carbon stocks of individual trees in African drylands. Nature 615, 80–86 (2023). Glob. Chang. Biol. 23, 177–190 (2017).
35. R. Dubayah, J. B. Blair, S. Goetz, L. Fatoyinbo, M. Hansen, S. Healey, M. Hofton, G. Hurtt, 57. N. P. Hanan, L. Prihodko, C. W. Ross, G. Bucini, A. T. Tredennick, “Gridded estimates of woody
J. Kellner, S. Luthcke, J. Armston, H. Tang, L. Duncanson, S. Hancock, P. Jantz, S. Marselis, cover and biomass across Sub-Saharan Africa, 2000–2004” (ORNL DAAC, 2020); https://doi.
P. L. Patterson, W. Qi, C. Silva, The Global Ecosystem Dynamics Investigation: High-reso- org/10.3334/ORNLDAAC/1777.
lution laser ranging of the Earth’s forests and topography. Sci. Remote Sens. 1, 58. G. L. Baskerville, Use of logarithmic regression in the estimation of plant biomass. Can.
100002 (2020). J. For. Res. 2, 49–53 (1972).
36. S. Quegan, T. le Toan, J. Chave, J. Dall, J. F. Exbrayat, D. H. T. Minh, M. Lomas, 59. B. D. Cook, L. A. Corp, R. F. Nelson, E. M. Middleton, D. C. Morton, J. T. McCorkel, J. G. Masek,
M. M. D’Alessandro, P. Paillou, K. Papathanassiou, F. Rocca, S. Saatchi, K. Scipal, H. Shugart, K. J. Ranson, V. Ly, P. M. Montesano, NASA Goddard’s LiDar, hyperspectral and thermal (G-
T. L. Smallman, M. J. Soja, S. Tebaldini, L. Ulander, L. Villard, M. Williams, The European LiHT) airborne imager. Remote Sens. (Basel) 5, 4045–4066 (2013).
Space Agency BIOMASS mission: Measuring forest above-ground biomass from space. 60. A. Pascual, J. Guerra-Hernández, D. N. Cosenza, V. Sandoval-Altelarrea, Using enhanced
Remote Sens. Environ. 227, 44–60 (2019). data co-registration to update Spanish National Forest Inventories (NFI) and to reduce
37. Y. Cheng, A. Vrieling, F. Fava, M. Meroni, M. Marshall, S. Gachoki, Phenology of short training data under LiDAR-assisted inference. Int. J. Remote Sens. 42, 126–147 (2021).
vegetation cycles in a Kenyan rangeland from PlanetScope and Sentinel-2. Remote Sens. 61. A. Pascual, J. Guerra-Hernández, J. Armston, D. M. Minor, L. I. Duncanson, P. B. May,
Environ. 248, 112004 (2020). J. R. Kellner, R. Dubayah, Assessing the performance of NASA’s GEDI L4A footprint
38. F. Reiner, M. Brandt, X. Tong, D. Skole, A. Kariryaa, P. Ciais, A. Davies, P. Hiernaux, J. Chave, aboveground biomass density models using National Forest Inventory and airborne laser
M. Mugabowindekwe, C. Igel, S. Oehmcke, F. Gieseke, S. Li, S. Liu, S. Saatchi, P. Boucher, scanning data in Mediterranean forest ecosystems. For. Ecol. Manage. 538, 120975 (2023).
J. Singh, S. Taugourdeau, M. Dendoncker, X.-P. Song, O. Mertz, C. J. Tucker, R. Fensholt, 62. E. Dinerstein, D. Olson, A. Joshi, C. Vynne, N. D. Burgess, E. Wikramanayake, N. Hahn,
More than one quarter of Africa’s tree cover found outside areas previously classified as S. Palminteri, P. Hedao, R. Noss, M. Hansen, H. Locke, E. C. Ellis, B. Jones, C. V. Barber,
forest. Nat. Commun. 14, 2258 (2022). R. Hayes, C. Kormos, V. Martin, E. Crist, W. Sechrest, L. Price, J. E. M. Baillie, D. Weeden,
39. S. Li, M. Brandt, R. Fensholt, A. Kariryaa, C. Igel, F. Gieseke, T. Nord-Larsen, S. Oehmke, K. Suckling, C. Davis, N. Sizer, R. Moore, D. Thau, T. Birch, P. Potapov, S. Turubanova,
A. Holm-Carlsen, S. Junttila, X. Tong, Alexandre d’Aspremont, P. Ciais, Digital twinning of all A. Tyukavina, N. de Souza, L. Pintea, J. C. Brito, O. A. Llewellyn, A. G. Miller, A. Patzelt,
forest and non-forest trees at national level via deep learning. (2022) 10.21203/rs.3.rs- S. A. Ghazanfar, J. Timberlake, H. Klöser, Y. Shennan-Farpón, R. Kindt, J. P. B. Lillesø, P. van
1661442/v1. Breugel, L. Graudal, M. Voge, K. F. al-Shammari, M. Saleem, An ecoregion-based approach
40. “Woodland cover target: detailed evidence report” (Department for Environment, Food to protecting half the terrestrial realm. Bioscience 67, 534–545 (2017).
and Rural Affairs, 2022). 63. IPCC Guidelines for National Greenhouse Gas Inventories (IPCC, 2006).
41. “The sixth carbon budget: The UK’s path to net zero” (The Climate Change Commit-
tee, 2020). Acknowledgments
42. M. Friedl, J. Gray, D. Sulla-Menashe, “MCD12Q2 MODIS/Terra+ aqua land cover dynamics Funding: This research was supported by the European Research Council (ERC) under the
yearly L3 global 500m SIN Grid V006. NASA EOSDIS land processes DAAC” (2019). European Union’s Horizon 2020 research and innovation programme (grant agreement no.
43. LAStools, Efficient LiDAR processing software (version 141017, academic); 947757 TOFDRY) and a DFF Sapere Aude grant (no. 9064–00049B). S. Liu is funded by China
http://rapidlasso.com/LAStools. Scholarship Council (CSC; grant no. 202104910123). R.F., C.I., and N.L. acknowledge support by
44. Microsoft Building Footprints; https://github.com/microsoft/GlobalMLBuildingFootprints. the Villum Fonden through the project Deep Learning and Remote Sensing for Unlocking
45. OpenStreetMap contributors, Planet dump building footprint database (2015); https:// Global Ecosystem Resource Dynamics (DeReEco; grant no. 34306). C.I. acknowledges support by
planet.openstreetmap.org. the Pioneer Centre for AI, DNRF grant number P1. Y.Y. and M.B. acknowledge support by the
International Partnership Program of Chinese Academy of Sciences (CAS) (092GJHZ2022029GC)
and CAS Interdisciplinary Team (JCTD-2021-16). J.G. acknowledges support by Foundation for

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 13 of 14

S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E

Science and Technology (FCT) through a CEEC contract (CEECIND/02576/2002). Author study: Canopy height, cover, and biomass maps across Europe have been deposited in the
contributions: M.B., R.F., and S. Liu designed the study. M.B. prepared and processed the aerial Zenodo database and are available at https://doi.org/10.5281/zenodo.8154445. Links to
LiDAR CHM data. T.N.-L. provided Danish NFI data. A.P. and J.G.-H. provided Spanish NFI data. download the open LiDAR data can be found under https://publications.jrc.ec.europa.eu/
F.R. and X.T. developed the code for the PlanetScope imagery generation pipeline. J.C. designed repository/bitstream/JRC126223/jrc126223_jrc126223_lidaropensourcedata.pdf. The code for
the biomass allometry, implemented by S. Liu, N.L., C.I., S. Li, M.M. S. Liu wrote the code for the tree segmentation and height regression is built on publicly open-source framework
deep learning framework, supported by F.R., S. Li, and Z.C. Interpretations were done by M.B., Segmentation Models PyTorch, available at https://github.com/qubvel/segmentation_models.
S. Liu, N.L., P.C., J.C., C.I., S.S., Y.Y., R.F., and S. Liu conducted the analyses. S. Liu and M.B. wrote pytorch. Custom code is publicly available at https://doi.org/10.5281/zenodo.8156190.
the first manuscript draft with contributions by all authors. S. Liu designed the figures.
Competing interests: The authors declare that they have no competing interests. Data and Submitted 1 March 2023
materials availability: All data needed to evaluate the conclusions in the paper are present in Accepted 15 August 2023
the paper and/or the Supplementary Materials. PlanetScope imagery was purchased via a Published 15 September 2023
departmental license, and the images cannot be distributed. Derived products produced in this 10.1126/sciadv.adh4097

Downloaded from https://www.science.org at Delft University on September 19, 2023

Liu et al., Sci. Adv. 9, eadh4097 (2023) 15 September 2023 14 of 14

 The overlooked contribution of trees outside forests to tree cover and woody
biomass across Europe
Siyu Liu, Martin Brandt, Thomas Nord-Larsen, Jerome Chave, Florian Reiner, Nico Lang, Xiaoye Tong, Philippe Ciais,
Christian Igel, Adrian Pascual, Juan Guerra-Hernandez, Sizhuo Li, Maurice Mugabowindekwe, Sassan Saatchi, Yuemin
Yue, Zhengchao Chen, and Rasmus Fensholt

Sci. Adv., 9 (37), eadh4097.

DOI: 10.1126/sciadv.adh4097

Downloaded from https://www.science.org at Delft University on September 19, 2023

View the article online
https://www.science.org/doi/10.1126/sciadv.adh4097
Permissions
https://www.science.org/help/reprints-and-permissions

Use of this article is subject to the Terms of service

Science Advances (ISSN ) is published by the American Association for the Advancement of Science. 1200 New York Avenue NW,
Washington, DC 20005. The title Science Advances is a registered trademark of AAAS.
Copyright © 2023 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim
to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).