<figure><img src="image_1.png" /></figure>

HAL
open science

# Toward a general tropical forest biomass prediction model from very high resolution optical satellite images

Pierre Ploton, Nicolas Barbier, Pierre Couteron, Cécile Antin, Narayanan Ayyappan, Natesan Balachandran, J.-F. Bastin, N. Barathan, G. Chuyong, G. Dauby, et al.

► To cite this version:

Pierre Ploton, Nicolas Barbier, Pierre Couteron, Cécile Antin, Narayanan Ayyappan, et al. Toward a general tropical forest biomass prediction model from very high resolution optical satellite images. Remote Sensing of Environment, 2017, 200, pp.140 - 153. 10.1016/j.rse.2017.08.001. hal-01824235

HAL Id: hal-01824235

https://hal.science/hal-01824235v1

Submitted on 3 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **TITLE.**

2 Toward a general tropical forest biomass prediction model from very high resolution optical

3 satellite images

4

5 **AUTHORS.**

6 P. Ploton¹,², N. Barbier¹, P. Couteron¹, C.M. Antin¹, N. Ayyappan³, N. Balachandran³, N.

7 Barathan³, J.-F. Bastin⁴, G. Chuyong⁵, G. Dauby⁶,⁷, V. Droissart¹,⁸, J.-P. Gastellu-Etchegorry⁹,

8 N.G. Kamdem¹⁰, D. Kenfack¹¹, M. Libalah¹⁰, G.II. Mofack¹⁰, S.T. Momo¹, ¹⁰, S. Pargal¹, P.

9 Petronelli¹², C. Proisy¹,³, M. Réjou-Méchain¹,³, B. Sonké¹⁰, N. Texier¹, ¹⁰, D. Thomas¹³, P.

10 Verley¹, D. Zebaze Dongmo¹⁰, U. Berger¹⁴ and R. Pélissier¹,³

11

12 [1] AMAP Lab, IRD, CIRAD, CNRS, INRA, Montpellier University, Montpellier, France

13 [2] Institut des sciences et industries du vivant et de l'environnement, Montpellier, France

14 [3] French Institute of Pondicherry, Puducherry, India

15 [4] Landscape Ecology and Plant Production Systems Unit, Université libre de Bruxelles,
16 Brussels, Belgium

17 [5] Department of Botany and Plant Physiology, University of Buea, Buea, Cameroon

18 [6] Institut de Recherche pour le Développement, UMR-DIADE, Montpellier, France

19 [7] Evolutionary Biology and Ecology, Faculté des Sciences, Université libre de Bruxelles,
20 Brussels, Belgium
21 [8] Herbarium et Bibliothèque de Botanique africaine, Université libre de Bruxelles, Brussels,

22 Belgium

23 [9] Université Paul Sabatier, CESBIO, Toulouse, France

24 [10] Laboratoire de Botanique systématique et d'Ecologie, Département des Sciences

25 Biologiques, Ecole Normale Supérieure, Université de Yaoundé I, Yaoundé, Cameroon

26 [11] Center for Tropical Forest Science — Forest Global Earth observatory, Smithsonian

27 Tropical Research Institute, Washington, USA

28 [12] Centre de coopération Internationale en Recherche Agronomique pour le

29 Développement, UMR-ECOFOG, Kourou, France

30 [13] Department of Biological Sciences, Washington State University, Vancouver, U.S.A.

31 [14] Technische Universität Dresden, Faculty of Environmental Sciences, Institute of Forest

32 Growth and Forest Computer Sciences, Tharandt, Germany

33

34 Correspondence to:

35 P. Ploton (p.ploton@gmail.com)

36 IRD — UMR AMAP,

37 Botany and Modeling of Architecture of Plants and Vegetations,

38 Bd de la Lironde,

39 TA A51/PS2,

40 34398 Montpellier Cedex 5,

41 France
42

ABSTRACT

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

Very high spatial resolution (VHSR) optical satellite imagery has shown good potential to provide non-saturating proxies of tropical forest aboveground biomass (AGB) from the analysis of canopy texture, for instance through the Fourier Transform Textural Ordination method. Empirical case studies however showed that the relationship between Fourier texture features and forest AGB varies across forest types and regions of the world, limiting model transferability. A better understanding of the biophysical mechanisms on which canopy texture – forest AGB relation relies is a prerequisite to move toward broad scale applications. Here we simulated VHSR optical canopy scenes in identical sun-sensor geometry for 279 1-ha tropical forest inventory plots distributed across the tropics. Our aim was to assess the respective merits and complementarity of two types of texture analysis techniques (i.e. Fourier and lacunarity) on a set of forests with contrasted structure and geographical origin, and develop a general texture-based approach for tropical forest AGB mapping. Across forests, Fourier texture captured a gradient of stands mean crown size reflecting well the progressive changes in stand structure throughout forest aggradation phase (e.g. Pearson's *r* = -0.42 with basal area) while lacunarity texture captured a gradient of canopy openness (, i.e. Pearson's *r* = -0.57 with stand gap fraction). Both types of texture indices were highly complementary for predicting forest AGB at the global level (so-called FL-model). The residual error of the FL-model was structured across sites and could be partially captured with a bioclimatic proxy, further improving the performance of the global model (so-called FLE-model) and reducing site-level biases. The FLE model was tested on a set of real Pleiades images covering a mosaic of high-biomass forests in the Congo basin (mean AGB over 49 field plots: 359±98 Mg.ha⁻¹), leading to a significant relationship (R² = 0.47 on
65 validation data) with reasonable error levels (< 25% rRMSE). The increasing availability of
66 VHSR optical sensors (such as from constellations of small satellite platforms) raises the
67 possibility of routine repeated imaging of the world's tropical forests and suggests that
68 texture-based analyses could become an essential tool in international efforts to monitor
69 carbon emissions from deforestation and forest degradations (REDD+).

70 KEYWORDS

71 Tropical forests, Canopy structure, Forest carbon, REDD, Passive optical imagery, Texture,
72 Fourier transform, Lacunarity

73 1. INTRODUCTION

74 Concerns about the effects of increasing atmospheric carbon dioxide on climate have led to
75 an international initiative aiming at reducing forest-related emissions (Reducing Emissions
76 from Deforestation and forest Degradation, REDD+), notably in the tropics where the bulk of
77 global deforestation occurs (Pan et al., 2011). REDD+ implementation fundamentally relies
78 on our capacity to monitor forest carbon stock and dynamics at multiple spatial scales, from
79 entire countries down to scales at which deforestation and degradation processes occur. In
80 this context, remote sensing naturally becomes an essential tool (Baccini et al., 2012; DeFries
81 et al., 2007; Saatchi et al., 2011). However, remote sensing of forest carbon stocks (often
82 through forest aboveground biomass, hereafter denoted AGB) is challenging in the tropics
83 because most satellite sensors are not sensitive to AGB variation above c. 150 Mg.ha⁻¹, while
84 AGB often exceeds 400 Mg.ha⁻¹ in tropical forests (Slik et al., 2013). This saturation is well
85 documented for passive optical sensors of coarse to intermediate spatial resolution such as
86 the Moderate Resolution Imaging Spectroradiometer (MODIS) or Landsat Thematic Mapper
87 (e.g. Lu, 2006; Lu et al., 2012; Zhao et al., 2016) but also for radar signals (notably L-band SAR, Mermoz et al., 2015), and thus constitutes a strong limit for those data types. In the past years, aircraft-based light detection and ranging systems have become popular for tropical forest AGB mapping, as their data appear to be free of saturation. Small-footprint airborne LiDAR data provide a detailed description of forest three-dimensional (3D) structure from which forest AGB can be estimated with good confidence (e.g. 14% relative error on 1-ha plots in Réjou-Méchain et al., 2015). Unfortunately this information comes at a cost rendering the wall-to-wall and regular coverage of large territories uneconomical (Erdody and Moskal, 2010; Messinger et al., 2016) or is even impossible to collect in some countries because of flight restrictions. An interesting alternative to airborne LiDAR may be found in very high spatial resolution (VHSR) optical images, as they are routinely captured by a variety of satellite platforms and therefore represent an independent, affordable source of data for large-scale monitoring of tropical forests. The improvement brought by the increased spatial resolution over coarser, widely-used optical data (e.g. Landsat), lies on a fundamental change in the analysis: at metric to sub-metric resolutions, individual trees are now discernable in the image, allowing one to exploit the geometric properties of tree crowns (i.e. crown delineation algorithms, e.g. Broadbent et al., 2008; Zhou et al., 2013) or of forest canopies (i.e. texture analysis, e.g. Couteron et al., 2005) that conceptually more directly relate to stand biomass than spectral greenness indices.

106 Texture analysis of panchromatic VHSR canopy images either captured by optical satellite sensors (e.g. Couteron, 2002, 2005) or derived from LiDAR canopy height models (e.g. Frazer et al. 2005, Véga et al. 2015) have shown promising results to retrieve stand and canopy forest attributes, including tropical forest AGB with no sign of saturation (e.g. Bastin et al., 2014; Ploton et al., 2012). Generally speaking, texture analysis (sensu Haralick, 1979) gives us
a measure of the spatial arrangement of grey levels within a canopy image by directly reflecting the contrasts between sunlit and shadowed surfaces, thus providing information on the size and distribution of crowns of canopy trees and inter-crown canopy gaps. Texture is scale-dependent, so that meaningful texture metrics have to be multi-scale by nature or at least measurable at multiple spatial scales. In the Fourier Transform Textural Ordination (FOTO) method (Couteron, 2002; Proisy et al., 2007), texture analysis is carried out within unit-windows large enough (e.g. of c. 1 ha or more) to capture pseudo-periodic repetitions of the patterns of interest, that is local variations in crown and gap sizes. On each unit window, the Fourier transform is applied on the 2D autocovariance function of the grey levels which intrinsically conveys multiscale information. The resulting Fourier spectrum describes the dominating scales of grey levels variation within the unit window. In a second step, Fourier spectra extracted from a large number of unit-windows covering a large forest tract (e.g. hundreds of square kilometers or more), are systematically compared using standard multivariate ordination techniques such as Principal Component Analysis, thus displaying the unit-windows along uncorrelated texture gradients that facilitate their separation and interpretation. Texture indices (i.e. unit-window scores on the main PCA axes) are then used as independent variables in a regression analysis summarizing the statistical relationship between canopy texture and stand structure parameters, for instance AGB, as measured from field plot data (e.g. Couteron et al., 2005; Proisy et al., 2007). The basic idea is that once calibrated, such relationships can be used to predict and map stand structure parameters outside the sampling plots, directly from texture analysis of a canopy image (Fig. 1).

133
<figure><img src="image_3.png" /><figcaption>Figure 1. Typical workflow of texture-based AGB mapping method.</figcaption></figure>

To date, the FOTO method has been applied in a number of tropical forest case studies to derive plot-level AGB estimates (e.g. Ploton et al., 2012; Proisy et al., 2007; Singh et al., 2014). Study areas generally covered a few hundreds of square kilometers, corresponding to one or two VHSR images. At such a scale the texture gradient provides good relationships with forest stand structure parameters, with a relative error on 1-ha plot AGB estimations
generally below 20%. It has indeed been demonstrated from simulations that in closed-
canopy forests, FOTO-texture indices characterize the size distribution of canopy crowns
(Barbier et al., 2012), which are allometrically linked to tree size and AGB (Jucker et al., 2016;
Ploton et al., 2016). However, it is increasingly evident that the texture-structure
relationship is forest type- and/or site-dependent (Bastin et al., 2014), thus limiting broad-
scale applications. As a consequence, mapping forest AGB on large mixed-forest landscapes
currently requires a stratification into homogeneous texture-AGB strata (e.g. Bastin et al.,
2014). It is worth mentioning that similar limitations are encountered when mapping forest
AGB with airborne LiDAR data, requiring complementary information layers to allow multi-
site comparisons (e.g. Asner et al., 2011; Vincent et al., 2014, 2012). Much like the universal
airborne LiDAR approach (Asner et al., 2011), generalizing the texture-structure relationship
across contrasted forest types from different regions of the world would be a major step
toward an operational method to monitor forest structure and AGB from VHSR optical data.
This paper is a first attempt in this direction, which tries to overcome two major limitations
of the current FOTO approach.

The first limitation is that in a VHSR canopy image, tree crowns and gaps are visible to the
naked-eye from local brightness variations, but texture patterns captured by FOTO cannot
discriminate large tree crowns from canopy gaps as they both appear as aggregates of bright
or dark pixels of similar size in the image. Coarse FOTO-texture thus indistinctly results from
either closed canopies comprising large tree crowns or from open canopies where
aggregates of tree crowns alternate with canopy openings. Thus, a same FOTO-texture may
represent high or low AGB values, respectively (see for instance, Pargal et al., 2017). Our
hypothesis is that environmental conditions and type of gap dynamics influence canopy
structure and heterogeneity, and in turn makes texture-structure relationships vary across
sites. To overcome this drawback, we sought to complement FOTO texture with additional texture features characterizing other aspects of the canopy structure. We focused in particular on lacunarity analysis and derived texture features, which provides a quantitative measure of image 'gapiness' and proved to correlate well with canopy cover and gap volume on simulated forest stands (Frazer et al., 2005).

A second major constraint limiting the empirical exploration of the link between canopy texture and the underlying forest stand structure over large regions covered by several VHSR images is of instrumental nature. Forest texture indices are indeed sensitive to sensor optical properties and spatial resolution (Ploton et al., 2012; Proisy et al., 2007), but also to sun-sensor geometry (Barbier et al., 2011; Barbier and Couteron, 2015). Disentangling the effects of forest structure and instrumental configuration on canopy texture requires having a sufficiently large set of VHSR images acquired in controlled sun-sensor configurations and covering an extensive network of ground-truth field plots of varying stand structure and canopy aspect. An efficient workaround is to adopt a modelling approach coupling the simulation of spatially explicit 3D forest stand models (referred to as "mockups") and the generation of virtual canopy images by applying a radiative transfer model onto the 3D mockups (Barbier et al., 2012). By controlling parameters of the radiative transfer model, canopy images in identical acquisition configurations can be generated for a variety of 3D stand mockups. Building forest mockups from simple tree shapes and stand inventory data (e.g. tree diameter distributions, local tree size allometries) may appear as a coarse representation of reality, but it is sufficient to capture how the 3D arrangement of stems and leaves interacts with important processes of the forest dynamics (Stark et al., 2015; Taubert et al., 2015). The simulation approach proved a valuable tool to bridge the gap between ground observations and remote sensing data (e.g. Frazer et al. 2005, Schneider et al., 2014),
and ultimately help us understand how signal information can be reliably interpreted in terms of biophysical parameters.

In this study, we used a simulation approach based on extensive field inventory data to test whether a generalized canopy texture model could provide consistent estimates of AGB across structurally contrasted forest types. We first generated virtual canopy scenes in homogeneous acquisition configurations from 279 1-ha field plots collected in several tropical forest sites on three continents. These virtual canopy scenes allowed us to address how canopy texture features in contrasted forest dynamics and environmental conditions vary with the underlying stand structure, notably AGB. More specifically, we investigated from these simulated images the respective merits of the FOTO method and the lacunarity analysis in predicting AGB both locally (within sites) and globally (across sites), and explored whether combining the two methods improved AGB prediction models. Finally, we tested the texture-structure inversion model developed from simulations on three real Pleiades images acquired in comparable configurations over 49 1-ha field plots from a complex mosaic of forest types in Central Africa.

## 2. MATERIAL AND METHODS

### 2.1 Forest inventory data

As a basis of both empirical estimations of forest plot AGB and simulations of 3D stand mockups for generating virtual canopy images, we used inventory data from a set of 279 1-ha forest sample plots in tropical Africa (157 ha), India (37 ha), and French Guiana (85 ha). Forest inventory data contained tree diameter at breast height (D) for all trees with $D \ge 10$ cm, for a total of 142,791 trees. Tree taxonomic identification was available at the species
level for 78% of the trees, at the genus level for 6%, at the family level for 5% and 11% of the trees were left unidentified. We used the Dryad Global Wood Density Database (Chave et al., 2009; Zanne et al., 2009) and the World Agroforestry Center's Wood Density Database (ICRAF, 2007) to attribute to each individual tree the mean wood density value of its species.

For those known only to the genus or family level, the average wood density at that taxonomic level was used (Chave et al., 2006). Unidentified individuals were attributed the average wood density of the 1-ha plot. Besides *D*, other tree dimensions were recorded on a subset of trees per plot, namely total tree height (*H*, n=23,237), trunk height (*Ht*, n=6,502) defined as the height to the lowest main living branch (as in Ploton et al., 2016, see also interactive discussion on this paper for more information) and tree crown diameter (*Cd*, n = 4,438). The distribution of field inventory data among sampling sites is provided in Table A1.

We established allometric models for tree *H*, *Ht* and *Cd* in order to predict the dimensions of unmeasured trees. In each plot, we used allometric models established at the plot-, site- or regional-levels (i.e. Africa, India and French Guiana) depending on data availability.

Following Feldpausch et al. (2012), we used a three-parameter Weibull function to predict tree *H* from *D*: $H = a(1 - \exp(-bD^c))$. Trunk height was modeled as a power function of *H*: $Ht = a * H^b$, fitted in logarithmic units and accounting for Baskerville correction (Baskerville, 1972). The same model and transformation were applied to fit the crown diameter, *Cd*, to *D*, at the exception that a 2-segments model was used when a significant breaking point was found in the relationship (following the procedure described in Antin et al., 2013). Finally, we used the pantropical model of Chave et al. (2014), including *D*, *H* and wood density to compute tree AGB estimates. In order to avoid increasing technical complexity, errors on AGB estimates from field data were neglected in this paper, though
246 they would have to be accounted for and propagated through the up-scaling process in
247 operational applications (e.g. Pargal et al., 2017; Réjou-Méchain et al., 2017).

248 Overall, the sample plots spanned wide gradients of tree density (139 to 848 trees.ha⁻¹),
249 stand basal area (10.3 to 57.7 m².ha⁻¹) and dominant height (17.8 to 50.5 m), and thus
250 captured a large array of forest age and 3D organizations.

251 **2.2 Generation of 3D forest mockups**

252 In order to simulate, using a radiative transfer model, realistic virtual canopy images in
253 homogenous acquisition configurations, we constructed 3D representations of the 279 1-ha
254 sample plots using local tree allometries and tree location data (either exact coordinates or
255 by sub-quadrats of 20 m sides). The 3D modelling procedure of a sample plot can be
256 summarized as follows: (1) we built a simplified representation of each tree in the plot with
257 the crown modelled as an ellipsoid of diameter Cd and depth Hc (i.e. H – Ht, in m) using field
258 measurements if available or local allometric models otherwise; (2) if tree coordinates within
259 the plot were known (i.e. for 135 of the 279 1-ha plots), we generated a single plot mockup
260 by placing trees at their actual position; if tree positions were only known at the quadrat-
261 level (i.e. for the remaining 144 plots): (i) trees were sorted from the tallest to the smallest;
262 (ii) starting from the tallest tree downward, trees were tentatively located in their respective
263 quadrats by generating random X and Y coordinates within an uniform distribution, and; (iii)
264 retained if their crown volume did not overlap crowns of already placed trees by more than
265 25% of tree crown volume. If the algorithm failed to place a tree after 50 iterations (in 5.3%
266 of the cases), its crown diameter was reduced by randomly sampling a value in the
267 distribution of the crown diameter allometric model residuals. For those 144 plots, we
generated 3 mockups per plot so as to capture some variability in canopy texture emerging
from the random component of the mockup generation algorithm.

In order to help the interpretation of the canopy texture analyses introduced further in this
paper, we also extracted from the 3D stand mockups simple parameters of stand structure
and heterogeneity, such as mean and maximum tree sizes, stand gap fraction (defined as the
cumulated area of pixels < 2 m in height) and canopy roughness (defined as within-plot
variation of canopy height, Table 2).

## 2.3 Simulation of canopy images

We simulated very-high spatial resolution optical images of the 567 forest mockups using
the Discrete Anisotropic Radiative Transfer model (DART, version 5.5.3, Gastellu-Etchegorry
et al., 2015). DART is able to handle the forest mockups as a 3D array of cells (voxels) of
different sizes. We used 1-m³ voxels so to match with the metric resolution at which FOTO
analyses are typically performed. Voxels contain a mixture of air and vegetation elements,
either as solid facets or as turbid medium with a predefined distribution of leaf volume
density. As leaf and bark optical properties of most tree species in our plots were not
available, we used the optical properties of a single African canopy tree species (Terminalia
superba Engl. & Diels) as input parameters for all the trees in DART. Tree trunks were
represented as solid cylinders and crowns as turbid ellipsoids with a spherical distribution of
leaf angles. Volume density within tree crown cells was adjusted to represent a realistic
range of Leaf Area Index at the plot scale: (i) we computed the cumulated crowns volume of
each plot; (ii) calibrated the vector of plot cumulated crowns volume to plot LAI by assigning
a LAI of 4 and 8 to its 5th and 95th percentiles, respectively, and; (iii) used this relationship to
derive the LAI of each plot and associated crown cell densities. Plot-scale LAI varied from 2.2
291 in the most open-canopy, low-density plots to 9.9 in the most highly packed plot (with
292 mean±sd of 6.3±1.4), in good agreement with tropical forest LAI estimations available in the
293 literature (Asner and Martin, 2008; Vincent et al., 2017). To represent a vegetated
294 understory in open canopy forests, we additionally attributed to the lowest 1-m layer above
295 the ground the optical properties of a pioneer tree species (*Musanga cecropioides* R.Br. ex
296 Tedlie) with a spherical distribution of leaf angles and leaf area density of 0.3 m².m⁻³.

297 We used DART in ray tracing, reflectance mode with the sun and the atmosphere as the only
298 radiation sources. We simulated top of atmosphere (TOA) panchromatic reflectance images
299 (see Gastellu-Etchegorry et al., 2015 for further details) in the 0.35 µm - 1.1 µm domain that
300 we decomposed into 75 spectral bands of 0.01 µm each. The radiative transfer calculations
301 were conducted to represent viewing and solar illumination elevation angles set to 70° and
302 56.5°, respectively, and a solar-viewing azimuth difference set to -92°, to match the
303 acquisition parameters of our real satellite images (see section 2.4). For each DART image
304 pixel, the reflectance channel of any VHSR sensor can then be obtained from the convolution
305 of reflectance values provided in the 75 bands by the sensor's modulation transfer function.
306 Computation time was about 40 minutes per image with about 12 GB of RAM on 64-bit
307 Windows with an Intel Core i7 GHz processor.

## 2.4 Real satellite images

309 To test on real VHSR optical images the texture-structure inversion models developed with
310 the simulated canopy images, we used a set of three Pleiades images acquired over a
311 network of 49 1-ha sample plots in the transition area between evergreen and semi-
312 deciduous forests of Eastern Cameroon, in central Africa. Images were taken approximately
313 two years after the establishment of the sample plots. Although the three images were in
backward configuration (i.e. with sun behind the sensor), textural bias between images induced by image acquisition configurations is expected to be low thanks to fairly homogeneous acquisition angles (Table 1; Barbier and Couteron, 2015). Following previous applications, texture analysis was based on panchromatic band (450–830 nm) images, resampled from 0.7 m to 1 m per pixel through mean pixel aggregation.

Table 1. Date and acquisition parameters of Pleiades panchromatic satellite images over Eastern Cameroun, central Africa

<table><thead><tr><th colspan="2">Date</th><th colspan="2">Sensor</th><th colspan="2">Sun</th></tr><tr><th>month</th><th>year</th><th>azimuth</th><th>elevation</th><th>azimuth</th><th>elevation</th></tr></thead><tbody><tr><td>12</td><td>2014</td><td>227.6</td><td>64.5</td><td>145.3</td><td>57.0</td></tr><tr><td>12</td><td>2014</td><td>248.5</td><td>79.7</td><td>146.4</td><td>55.3</td></tr><tr><td>1</td><td>2015</td><td>229.6</td><td>65.7</td><td>138.3</td><td>57.3</td></tr></tbody></table>

321
322
323
324

#1 catalog ID: DS_PHR1A_201412210944201_FR1_PX_E013N03_0909_01764

#2 catalog ID: DS_PHR1A_201412210943258_FR1_PX_E013N05_0606_01395

#3 catalog ID: DS_PHR1A_201501160944349_FR1_PX_E013N03_0702_00899

325

## 2.5 Extraction of texture features from real or simulated canopy images

We performed texture analysis of panchromatic VHSR canopy images using both FOTO (Couteron, 2002) and lacunarity (Frazer et al., 2005) analyses as they are expected to provide complementary information on stand structure attributes. FOTO method has been extensively described elsewhere (e.g. Couteron, 2002; Couteron et al., 2005) and we only give a brief outline of the procedure hereafter. The first step consisted in dividing a panchromatic satellite image into square 1-ha unit canopy windows (with 1-m spatial resolution thus a spacing of $S = 100$ in both X and Y directions), a size which proved suitable from previous studies to capture several repetitions of the largest tree crowns in closed forest stands (Couteron, 2002; Proisy et al., 2007). On each unit-window, the two-dimensional Fast Fourier Transform (fft2) function was applied to transpose the spectral
radiance from the spatial domain to the frequency domain, using sine-cosine functions at
integer frequencies (i.e. wavenumbers, 1, 2, ..., S/2) along the X and Y directions of the plane.
The squared amplitude of the fft2 yields a 2D periodogram, which represents an
apportionment of the variance in spectral radiance among spatial frequency bins in all
possible planar directions within the unit canopy window. The orientation information was
here neglected by averaging the periodogram across directions, which provided the one-
dimensional so-called radial- or r-spectrum. Knowing image spatial resolution (1 m),
wavenumbers can equivalently be expressed in wavelengths (λ, in m) so that a r-spectrum
gives the frequency distribution of the number of times a pattern of size λ repeats itself in a
1-ha unit canopy window. A characteristic of the method is that sampling according to
harmonic Fourier frequencies (stemming from fft2) results in denser spectrum values as λ
decreases, so that intermediate scales of patterns relevant for canopy characterization may
be poorly sampled. To overcome this we used the modified algorithm presented in Barbier
and Couteron (2015), in which, after signal centering, unit-window size is doubled in each
direction using zero-padding (Cressie, 1993), which allows increasing intermediate λ
sampling. The resulting list of λ was pruned to conserve a sampling of λ values as regular as
possible, thus decreasing information redundancy at small wavelengths. Finally, a set of 25 λ
values was retained from 2 up to 99 m. All r-spectra were assembled into a single table, F,
with the individual 1-ha unit-windows as rows and Fourier spatial frequencies as columns.
Unit-windows with r-spectra dominated by high Fourier frequencies (low λ) are expected to
display fine grain canopy textures due to the succession of small tree crowns, while coarse
grain canopy textures correspond to r-spectra dominated by low Fourier frequencies
expected to reflect the repetition of large canopy trees (see e.g. Barbier et al., 2012).
359 Lacunarity was formally defined by Mandelbrot (1983) as the deviation of a fractal pattern
360 from translational invariance, a concept further expanded to non-fractal patterns (Allain and
361 Cloitre, 1991). It basically measures how patterns fill space, those having more or larger gaps
362 having generally a higher lacunarity or gapiness. In the present paper, we computed
363 lacunarity following Frazer et al. (2005) to which the reader should refer for more details.
364 Lacunarity was computed in the same 1-ha unit canopy windows as Fourier r-spectra. On
365 each unit-window, a moving square box of size s (in pixels) is glided by one pixel at a time
366 along all rows and columns. For each box of size s, the mass, M, is computed as the sum of
367 all pixel spectral radiance (previously normalized per unit-window), so that the probability
368 distribution function Q(M,s) represents the frequency distribution of M divided by the
369 number of boxes at size s. The lacunarity statistic is then computed as Δ(s) =
370 Z<sub>Q</sub><sup>(2)</sup>(s)/(Z<sub>Q</sub><sup>(1)</sup>(s))<sup>2</sup>, where Z<sub>Q</sub><sup>(1)</sup>(s) and Z<sub>Q</sub><sup>(2)</sup>(s) are the mean and mean of squared values of
371 Q(M, s), respectively. In our study we used 10 discrete box sizes, s, from 1 to 99 m, and
372 normalized all lacunarity values per unit canopy windows by dividing them by Δ(1). We then
373 denote as a lacunarity spectrum the series of normalized lacunarity statistics computed at
374 each box size for a given unit-window. All normalized lacunarity spectra for a given satellite
375 image were assembled into a single table, L, with the unit-windows as rows and the box sizes
376 as columns. The decreasing pattern of lacunarity with increasing box size reflects the rate of
377 increase in canopy heterogeneity from the finer to the coarser scale patterns, which in the
378 context of forest canopies is related to the degree of inter-crowns canopy openness (see
379 Frazer et al., 2005).

380 **2.6 Statistical analyses**

381 2.6.1 Canopy texture indices
In order to systematically compare unit-windows textural properties (as characterized by Fourier and lacunarity spectra), both the **F** and **L** tables were submitted independently to Principal Component Analysis (PCA) as in Couteron et al. (2005) and Frazer et al. (2005). The **F** and **L** tables were column-wise normalized prior to perform PCAs, which allows uncovering texture gradients even when the absolute variation between unit-windows is small. This led to two independent sets of texture gradients, synthetizing FOTO-texture indices on the one hand (that is, unit-windows scores along the axes of the PCA on the **F** table, or F-PCA) and Lacunarity-texture indices on the other hand (i.e. unit-windows scores along the axes of L-PCA). We then investigated the agreement between these two analyses using a co-inertia analysis (COIA; Dolédec and Chessel, 1994, Dray et al. 2003). COIA is a simple and robust method for the simultaneous analysis of two data tables matching by observations (here, unit canopy windows. The method maximizes the square covariance between the projected coordinates of the unit-windows on the PCA axes of **F** and **L**. It thus simultaneously maximizes the variance in **F** and **L** separately (the two PCAs) and the correlation between the two sets of principal components. The texture gradients captured by FOTO and by the lacunarity analysis can thus be projected on common orthogonal components. The unit-window scores on these components are canopy texture indices used in further analyses as image-derived predictors of plot-level AGB.

### 2.6.2. Modelling global and local texture – AGB relationships

We used the Random Forest (RF) algorithm as implemented in the package "randomForest" version 4.6.12 of R statistical software (R Core Team, 2016) to develop global (i.e. across sites) forest AGB regression models based on canopy texture indices. RF is an ensemble of learning methods for classification and regression based on decision trees (Breiman, 2001)
405 increasingly used in remote-sensing studies, notably for carbon mapping (Baccini et al., 2012; Mascaro et al., 2014; Vieilledent et al., 2016). The aim is to overcome over-fitting
406
407 problems that can occur using individual decision trees. The algorithm builds a large number
408 of decision trees that are trained on random samples of both the training set (i.e. the
409 observations) and the predictor variables, and outputs the average of individual trees
410
411 predictions. An interesting feature of such machine learning techniques is that it allows
411
integrating a large number of variables of potentially different statistical distributions
412 (Marvin et al., 2016) and does not require specifying an explicit model form.

413 We used a similar (RF) regression framework, but with forest "site" taken as a random
414 variable in regression models, to assess the predictive power of canopy texture indices on
415 plot AGB at the local level (i.e. within sites). The mixed-effects RF model (MERF, Hajjem et
416 al., 2014) is defined as a classical linear mixed-effects model (Wu and Zhang, 2006), except
417 that a more general and unknown function, f(Xᵢ), where i denotes a forest site, is used to
418 estimate the within-site fixed effect using RF in place of the usual linear predictor, the
419 random part of the model being still assumed linear. The model is fitted with an expectation-
420 maximization (EM) algorithm (see Hajjem et al., 2014 for details and R script). For simplicity,
421 across site variations in texture – AGB relationships were accounted for with site-specific
422 random intercepts only.

423 2.6.3 AGB models validation

424 We evaluated regression model errors using a five-fold cross-validation procedure
425 independent of the internal validation scheme of RF. The dataset was randomly split into 5
426 folds of approximately equal size and the AGB of forest plots in each fold was predicted with
427 the regression model calibrated on the four other folds. This procedure was repeated 10
428 times so to compute mean and standard deviation of goodness-of-fit statistics. For each
429 model, we report the amount of variation explained as computed from calibration (R²c) and
430 validation (R²v) data. Using validation data only, we further provide the root-mean-square
431 error (RMSE) and relative RMSE (rRMSE in %), along with the mean signed deviation (MSD),
432 used to compare site-level prediction bias in RF and MERF models (Xu et al., 2016).

433 **3. RESULTS**

434 **3.1 Texture analysis of virtual canopy images**

435 We applied both the FOTO method and the lacunarity analysis on all unit canopy windows
436 simulated from the 279 1-ha plots of tropical forest inventory. The first factorial plane of
437 FOTO-PCA explained 51.1% (37.2% and 13.9% for axes 1 and 2, respectively) of the total
438 variability among Fourier r-spectra (F-table) and produced a typical correlation circle with
439 spatial patterns sorted from short to long wavelengths (or from high to low frequency/
440 wavenumber) (Fig. 2A). Wavelengths (λ) of about 5 to 10 m (visually relating to small crown
441 sizes) were found on the positive side of the first PCA axis (F-PCA1), intermediate λ (15 to 20
442 m) on the positive side of F-PCA2 and λ of 25 to 35 m (largest crown sizes) on the negative
443 side of F-PCA1. The ordination gradient on the first FOTO-PCA plane thus corresponded to a
444 fineness-coarseness canopy texture gradient, with aperiodic canopies found close to the
445 origin and heterogeneous ones (i.e. mixing small and large crown patches) found on the
446 lower-left part of the plane. The general decrease of canopy grain size along F-PCA1
447 negatively correlated with stand maximum tree size (Pearson's r ≈ -0.7 with DBH<sub>m</sub> and H<sub>m</sub> in
448 Table 2) and canopy roughness (r = -0.65 with H<sub>v10</sub>), which coincided with a moderate
449 decrease of parameters reflecting stand aggradation (r = -0.48 and -0.42 with DBH<sub>mean</sub> and G, respectively).
450

451 The first Lacunarity-PCA plane explained 98.1% (81.2% and 17% for axes 1 and 2, respectively) of the variance among normalized lacunarity spectra (L-table). The first PCA 452 axis (L-PCA1) sorted unit-windows according to their degree of gapiness (Fig. 2B) as shown 453 by a negative correlation with stand gap fraction (r=-0.57 with GF in Table 2), which 454 coincided with a positive correlation with stand basal area (r=0.6) and mean tree diameter 455 (r=0.48). L-PCA1 thus sorted unit-windows from fairly open stands (featuring frequent 456 patches of small trees, less than 2 m in height) to fully-stocked stands having a canopy made 457 of larger trees. L-PCA2 captured variations in the rate of decrease in lacunarity with box size, 458 which reflected a gradient in canopy grain and heterogeneity. Fairly homogeneous canopy 459 (smooth canopy grain) characterized by a low rate of decrease in lacunarity with box size 460 were found on the negative side of L-PCA2, while heterogeneous canopy were found on the 461 positive side of L-PCA2. As a matter of fact, L-PCA2 was positively correlated with the canopy 462 roughness index, H<sub>v10</sub> (r=0.61) and with maximum tree size (r=0.76 and 0.66 with H<sub>m</sub> and 463 DBH<sub>m</sub>, respectively). F-PCA1 and L-PCA2 displayed a very high correlation of 0.93 (axes 464 orientations are arbitrary and not informative in PCA) and fundamentally depicted the same 465 coarseness–fineness gradient of canopy aspect. On the other hand, correlations between F- 466 PCA2 and both Lacunarity-based PCA axes were weak, suggesting complementarity. 467

468

469

470
471

472

473

474

475

476

477

478

479

480

481

482

483

484

<figure><img src="image_4.png" /></figure>

485

486

487

488

489

490

Figure 2. Canopy texture ordinations based on (A) the FOTO method and (B) the lacunarity analysis. In both cases, scatter plots of PCA scores along the first two principal axes are shown (left panels), with 3 example sites highlighted with particular symbols (Paracou, Uppangala and Yellapur). Top-right panels represent correlation circles with wavelength, λ (A) or box size, s (B) in meter. Bottom-right panels represent histograms of eigenvalues in % of total variance.
491

492

493

494

495

Table 2. Correlation between stand structure parameters extracted from three-dimensional mockups and canopy window scores on the texture ordination axes based on the FOTO method (F-PCA1 and F-PCA2) and the lacunarity analysis (L-PCA1 and L-PCA2). Probability value of Pearson correlation test are provided between brackets and coded following standard notation (*** P ≤ 0.001, ** P ≤ 0.01, * P ≤ 0.05, ns = non-significant).

<table><thead><tr><th></th><th>Stand structure parameters</th><th>F-PCA1</th><th>F-PCA2</th><th>L-PCA1</th><th>L-PCA2</th></tr></thead><tbody><tr><td>GF</td><td>Gap fraction (cumulated area of pixels &lt; 2 m in height)</td><td>-0.04 (ns)</td><td>0.15 (***)</td><td>-0.57 (***)</td><td>0.07 (*)</td></tr><tr><td>H<sub>v10</sub></td><td>Standard deviation of maximum tree height among 10 m quadrats</td><td>-0.65 (***)</td><td>-0.24 (***)</td><td>-0.24 (***)</td><td>0.61 (***)</td></tr><tr><td>DBH<sub>m</sub></td><td>Stand maximum tree DBH</td><td>-0.69 (***)</td><td>-0.21 (***)</td><td>0.23 (***)</td><td>0.66 (***)</td></tr><tr><td>H<sub>m</sub></td><td>Stand maximum tree height</td><td>-0.71 (***)</td><td>0.1 (**)</td><td>-0.12 (**)</td><td>0.76 (***)</td></tr><tr><td>Cd<sub>m</sub></td><td>Stand maximum crown diameter</td><td>-0.56 (***)</td><td>-0.2 (***)</td><td>0.28 (***)</td><td>0.53 (***)</td></tr><tr><td>N</td><td>Stand tree density</td><td>0.2 (***)</td><td>0.1 (**)</td><td>-0.01 (ns)</td><td>-0.12 (**)</td></tr><tr><td>G</td><td>Stand basal area</td><td>-0.42 (***)</td><td>0.06 (ns)</td><td>0.6 (***)</td><td>0.44 (***)</td></tr><tr><td>DBH<sub>mean</sub></td><td>Quadratic mean of tree DBH</td><td>-0.48 (***)</td><td>0 (ns)</td><td>0.48 (***)</td><td>0.43 (***)</td></tr><tr><td>CD<sub>mean</sub></td><td>Quadratic mean of tree crown diameters</td><td>-0.29 (***)</td><td>-0.08 (*)</td><td>0.45 (***)</td><td>0.22 (***)</td></tr></tbody></table>

496

497
Combining the F and L tables in a common co-inertia analysis (Fig. 3) confirmed a slight but significant co-structure, with a RV coefficient (equivalent to the R² between two variables) of 0.27 (P < 0.001). The first COIA plane captured 99% of the co-inertia, which was mostly one-dimensional (61.9% on COIA-1). The first COIA plane (i.e. COIA-1 and 2) also captured 91.5% and 98.6% of the cumulated inertia projected on the first F-PCA and L-PCA planes, respectively. Figure 3 illustrates that while the co-inertia analysis preserves the structure of each separate PCA, i.e. the canopy texture gradients described by the FOTO (Fig. 3A) and lacunarity analyses (Fig. 3D), some disagreements between the two methods appeared at both the site and unit-window level. On the first COIA axis (COIA-1) for instance, Yellapur and Paracou sites are very close to each other and far from Uppangala in terms of mean FOTO-texture properties (large empty circles in Fig. 3 C), while Yellapur is intermediate between Paracou and Uppangala in terms of mean lacunarity-texture properties (large full circles, Fig. 3 C). In Paracou, canopy texture described by FOTO and lacunarity provided very similar information for all the unit-windows (short arrows in Fig. 3F), while for Yellapur and Uppangala, some unit-windows that exhibited very different FOTO textures (i.e. found at the opposite on COIA-1) showed a convergence towards similar lacunarity-derived textures (long convergent arrows in Fig. 3F). This demonstrates that FOTO and Lacunarity texture gradients, although globally correlated, can combine differently at the scale of unit-windows, and thus bring relevant information for discriminating individual canopy features.

517

518

519

520
<figure><img src="image_5.png" /><figcaption>Figure 3. Co-inertia analysis. Position on the first co-inertia plane of the FOTO r-spectra wavelengths, λ (A) and the lacunarity box size, s (D). Components of the F-PCA (B) and the L-PCA (E) projected onto the co-inertia axes. Ordination of windows from the 3 example sites (Paracou, Uppangala and Yellapur) on COIA-1 (C) with large empty and full circles representing the average site-level score for FOTO and Lacunarity features, respectively. Normed scores of 10 randomly sampled canopy windows from the 3 example sites on the first co-inertia plane (F), with each arrow linking a unit-window position for FOTO and Lacunarity characteristics, respectively.</figcaption></figure>

## 3.2 Canopy texture - AGB models

We built global stand AGB prediction models based on the two first axes of either F-PCA (F-model) or L-PCA (L-model) derived from the virtual canopy images (Fig. 4-A, C). Both models provided quite low goodness of fit (Rc² = 0.20 and 0.42 for F- and L-model, respectively). Model R² computed from independent validation data were even lower (R²v = 0.09 and 0.31), indicating that the internal validation scheme of RF tends to overestimate model
performance. The two models showed a general bias with respect to the 1:1 line, with a tendency to overestimate low biomass levels and underestimate high biomass levels, as well as more or less pronounced positive or negative biases occurring at the site-level (up to MSD = -32% for the L-model at Paracou). Using a mixed-effect RF algorithm (MERF) with a random site effect largely improved goodness-of-fit statistics (R²v = 0.31 and 0.61 for F- and L-model, respectively, Fig 4-B, D), indicating that a substantial gain in AGB prediction can be achieved when accounting for local (i.e. site-level) variation in texture – AGB relationships.

Combining both FOTO and Lacunarity texture features from COIA in a single global model (FL-model in Fig. 4-E) increased the prediction power on stand AGB (R²v = 0.49 with rRMSE = 32.3%). This global model also improved the prediction accuracy at Paracou (MSD = -18.1%, Fig. 4-E), but site-level MSDs remained substantial. It is noteworthy that site-level MSDs seemed to decrease with the range of biomass encompassed across plots in a site: while plots in Yellapur (MSD = 10.3%) and Paracou (MSD = -18.1%) are restricted to low and high biomass levels, respectively, in Uppangala (MSD = -5.3 %) plots have been sampled along a biomass gradient spanning from c. 150 up to > 600 Mg.ha⁻¹. This indicates that the underlying texture-AGB relationship varies across sites as a function of sampling, but it may also vary with respect to local site characteristics. For instance, while tree density in Paracou correlated well with FOTO and lacunarity scores on COIA-1 (r = 0.74 and 0.63, respectively), it did not in Yellapur (r < 0.18). Similar observations can be made for maximum tree slenderness (i.e. Hm/Dm) which correlated well with scores in Uppangala (r = 0.91 and 0.88, respectively) but not so clearly in Yellapur (r = 0.42 and 0.53, respectively). Again, accounting for across site variation in the texture – structure model by including a random site effect largely improved model performance (Fig. 4-F). The MERF FL-model outperformed all other models, with R²v = 0.69 and rRMSE = 23.5%.
559

<figure><img src="image_6.png" /></figure>

560

561

562

563

L-model w/o site effect

<figure><img src="image_7.png" /></figure>

564

565

566

567

568

FL-model w/o site effect

<figure><img src="image_8.png" /></figure>

569

570

571

572

F-model w. site effect

<figure><img src="image_9.png" /></figure>

L-model w. site effect

<figure><img src="image_10.png" /></figure>

FL-model w. site effect

<figure><img src="image_11.png" /></figure>

573

<caption>Figure 4. Aboveground biomass prediction models based on FOTO texture (F-model),</caption>

574

Lacunarity texture (L-model) and the two sources of texture information (FL-model), either

575

without site effect (RF algorithm, panels A, C and E) or with site effect (MERF algorithm,

576

panels B, D and F). Texture features were extracted from virtual canopy scenes. Goodness of

577

fit statistics are defined in the Methods.
Although the mixed modeling strategy proved powerful at sampled sites (i.e. where sample plots are available to account for a random site effect), it is not relevant to predict plot AGB in an unsampled new site. We therefore introduced in the FL-model a bioclimatic proxy, E, proposed by Chave et al. (2014), which is a compound variable based on water deficit, temperature and precipitation seasonality available as a global gridded layer at 2.5 arc sec resolution (see http://chave.ups-tlse.fr/pantropical_allometry.htm). E can be seen as expressing a bioclimatic stress on forest growth and was demonstrated to capture pantropical variation in forest AGB through its effect on tree height-diameter allometries (Chave et al., 2014), a source of variation than can hardly be captured by canopy texture.

Figure 5-A shows indeed that E largely captured the random site effect of the FL-model, with a downward adjustment of plot AGB (i.e. negative random intercept) at sites with high bioclimatic stress (i.e. when E is positive and high) and vice versa. Incorporating E along with FOTO and lacunarity texture features (FLE-model, Fig. 5-B) provided very similar results to the MERF FL-model (R² v= 0.69, with rRMSE = 24.9%) while maintaining local prediction errors (MSD) of the 3 example sites below 10%.

593

594

595

596

597

598

599
<figure><img src="image_12.png" /><figcaption>Figure 5. (A) Bioclimatic proxy E vs random site intercepts of FL-model (in Fig. 4-E); (B) Random forest AGB prediction model based on FOTO texture, lacunarity texture and E (FLE-model). Texture features were extracted from virtual canopy scenes. Goodness of fit statistics are defined in the Methods.</figcaption></figure>

## 3.3 Application to real satellite images

As the virtual canopy scenes are necessarily imperfect representations of canopy images, we also tested the texture – AGB models on a set of three real VHSR optical images acquired in homogeneous acquisition configurations over our sampling plots in central Africa. The F- and L-models did not, or weakly discriminate AGB variations ($R^2$ = 0.08 and $R^2$ = 0.1, respectively) among the 49 1-ha forest inventory plots. Combining F and L texture features (FL-model) led to a significant improvement of model goodness-of-fit ($R^2v$ = 0.42 with rRMSE = 25.1%), which further improved by the inclusion of E as an additional predictor (FLE model, $R^2v$ = 0.46 with rRMSE = 23.2%, Fig. 6). Though a limited sample of canopy images in homogeneous acquisition configurations were available to us, these results from central
620 Africa confirm the power of the FLE model to predict regional variation in plot AGB from real
canopy images too.
621

622

623

624

625

626

627

<figure><img src="image_13.png" /></figure>

628 Figure 6. Aboveground biomass prediction models based on FOTO texture, Lacunarity
629 texture and the bioclimatic stress proxy E (FLE-model). Texture features were extracted from
630 three Pleiades covering 49 1-ha forest inventory plots in central Africa. Goodness of fit
631 statistics are defined in Methods section.
632

632

633 **4. DISCUSSION**

634 Over the past decade, about two dozen studies successfully used canopy texture analysis
635 applied on VHSR optical or LiDAR canopy height models data to uncover spatial gradients of
636 forest structure and AGB, including in high-biomass tropical forests (e.g. Couteron et al.,
637 2005; Frazer et al., 2005; Malhi and Román-Cuesta, 2008; Proisy et al., 2007). These studies
638 were however limited to relatively small geographical areas and often retrieved structure
639 gradients within a single, homogeneous forest type and image. The few attempts made in
mosaics of heterogeneous forest patches (e.g. Bastin et al., 2014; Pargal et al., 2017; Singh et al., 2014) have shown that the link between canopy texture and stand structure depends on forest type, hindering broad-scale applications of texture-based methods. Using 279 1-ha plots distributed among different forest types across the tropics, we evaluated whether a generalized biomass prediction model based on canopy texture indices could provide consistent predictions at both local and global scales. The results presented here based on virtual canopy images show that it is worth complementing FOTO with lacunarity texture indices capturing canopy features related to canopy openness. Introducing a bioclimatic stress variable that captured regional variation in tree height-diameter allometry also substantially improved the accuracy and precision of forest AGB retrievals among forest sites. A practical application of the method to a mosaic of forest types in the Congo basin showed that forest AGB inferences could be made with reasonable precision (i.e. ≤ 25% of error) up to 600 Mg.ha⁻¹, without saturation.

## 4.1. Contrasted canopy texture - stand AGB relationships across sites

Global models based only on FOTO or lacunarity texture indices (F- and L-model, respectively) significantly explained AGB variations but presented a relatively high uncertainty (rRMSE = 38-45%), consistent with what other studies have reported for heterogeneous forest landscapes (e.g. Bastin et al. 2014). To circumvent the problem and improve AGB inferences, previous empirical studies typically used a forest type stratification step prior to the calibration of within-class texture-AGB models, suggesting that the relationship is dependent on forest type (Bastin et al., 2014; Singh et al., 2014). Our results indeed showed that forest plots from different geographical sites were often clustered along
the texture gradients and that texture indices did not always correlate with the same stand structure parameters at different sites, generating systematic biases in global models.

The relationship between FOTO-texture indices and forest AGB is expected to hold well for periodic canopy patterns with homogeneous grain size (Proisy et al., 2007). In this ideal case, Fourier r-spectra peak around the scale of the mean crown size of canopy trees (Barbier et al., 2010). Because crown size is allometrically related to tree AGB (Jucker et al., 2016; Ploton et al., 2016) and the AGB of the 20 biggest trees in 1-ha stands capture on average c. 85% of the variation across whole stands AGB (Bastin et al., 2015), FOTO texture accurately predicts variation in forest stand AGB. This is for instance the case along the successional development pathway of mangrove tree cohorts (Proisy et al., 2007) but can also apply to mixed-forests where the coarseness-fineness canopy texture gradient reflects differences in the stature of mature trees (e.g. due to contrasting soil fertility, Couteron et al., 2005) or along forest biomass gradients. At the Paracou study site for instance, a good linear correlation (Pearson's r = -0.88) existed between tree density (N) and the mean quadratic diameter (DBH<sub>mean</sub>) suggesting that sample plots align along a self-thinning trajectory (see also Vincent et al., 2012). In this case, our simulation procedure based on simple allometric relationships thus produced fully-stocked mockups along with virtual canopy images often displaying a periodic aspect (see Fig. 2) and for which a local F-model (i.e. based on FOTO texture only) provided very accurate AGB predictions (rRMSE = c. 10%). However, when the forest enters the gap phase dynamics (Withmore, 1975), stand structure becomes heterogeneous in both the horizontal and vertical dimensions, leading to higher canopy roughness, more frequent canopy gaps alternating with clusters of large trees (e.g. Franklin et al., 2002; Guariguata and Ostertag, 2001; Spies, 1998). In Uppangala for instance, where N and DBH<sub>mean</sub> did not correlate over the sampling plots (Pearson's r = -0.09), simulated
canopy images were highly heterogeneous (see Fig. 2), and a local L-model (i.e. based on lacunarity texture only) produced more accurate AGB predictions than a local F-model (rRMSE< 10% vs. c. 20%, respectively). Our study showed that lacunarity is an efficient method to reveal canopy openness and heterogeneity gradients, which confirms previous findings (e.g. Frazer et al., 2005; Malhi and Román-Cuesta, 2008) and stresses the importance of those canopy structural descriptors for predicting local variation of forest AGB at some sites (e.g. Uppangla). It is therefore easy to understand that, as hypothesized, complementing FOTO texture with lacunarity texture provided a more comprehensive description of the potential range of canopy structures and thus largely improved our multi-site prediction model (FL-model) at both the local and global levels.

At the global scale however, variation in canopy height (Fayad et al., 2016; Saatchi et al., 2011) and tree slenderness (Feldpausch et al., 2012), which proved critical for accurate AGB estimations (Chave et al., 2014), cannot be directly accounted for by a 2D analysis of canopy texture, whatever the method. Along a local forest successional gradient, canopy height variations follow increment in tree size visible from above through the mean crown sizes. But, the maximal height reachable by dominant trees is known to generally reveal the growth potential of a forest in relation to soil and regional bioclimatic constraints, in particular with regard to water stress (Chave et al., 2014). Introducing the E bioclimatic variable, which combines water deficit with temperature and precipitation seasonality (see Chave et al., 2014) improved goodness of fit of our global regression model (FLE-model) and significantly reduced local prediction errors. More studies are however needed to further test the robustness of the bioclimatic E variable as a proxy for forest potential height to be used in large-scale assessments of forest AGB from 2D optical data.
On the set of real images from the Congo basin, the generalized texture-based model was able to detect spatial variations in AGB within a tropical forest mosaic which may be seen as a difficult case since being characterized by both high AGB levels (i.e. average AGB over 49 field plots of 359±98 Mg.ha⁻¹) and forest types showing important variations in stand 3D structure, from closed-canopy mixed-species mature stands to open-canopy *Marantaceae* forests. We indeed found that the multi-site FLE-model was the most robust with uncertainty of about 70 Mg.ha⁻¹ at 1-ha scale (RMSE) corresponding to a relative error of less than 25% (rRMSE). This result confirms and expands the results previously obtained with a variety of satellite sensors (Quickbird, GeoEye, IKONOS, SPOT-5), aerial images and even images freely available from the Google Earth engine (Bastin et al., 2014; Meng et al., 2016; Ploton et al., 2012; Proisy et al., 2007; Singh et al., 2014, 2015) in temperate or tropical forests. This represents an important step toward using VHSR optical imagery for broad-scale assessments of forest AGB, as it eliminates the need of a preliminary forest classification (see for instance Pargal et al., 2017). Recent progress also opened the perspective to perform texture analysis on inter-calibrated satellite images from various sensors and/or in various configurations, provided they partly overlap (Barbier and Couteron, 2015).

## 4.2. On 3D stand mockups and virtual canopy images for model calibration

As canopy texture is sensitive to variation in shadowing pattern created by the various combinations of sun-sensor-scene angles (Barbier et al. 2011), we used in this paper a simulation approach to assemble a large array of above-plot VHRS canopy images in homogeneous acquisition configurations, which remains challenging from real satellite sensors. The approach based on 3D stand mockups derived from field data and allometric
relationships, on which a radiative transfer simulation is applied, follows previously published studies on image texture simulations, notably Barbier et al. (2012, 2010), Barbier and Couteron (2015) and Proisy et al. (2016). At the scale of 1-ha forest plots, the size distribution of large objects (individual and aggregated tree crowns) and their spatial arrangements (e.g. inter-crown gaps and associated shadows, variations in tree density and size between subplots, etc.) determine the coarsest patterns visible on canopy scenes and are crucial determinants of spatial variations of the apparent reflectance. Using a unique, simplified tree shape representation (i.e. cylindrical trunks and non-plastic, ellipsoid crowns) allows generating sufficiently realistic brightness variation patterns at those coarse scales for the interpretation of canopy texture gradients. For instance, Proisy et al. (2016) showed that this approach allows producing texture r-spectra that have similar frequency peaks than the real image counterparts over a wide range of mangrove successional stages. In the present study, a particular care was taken to develop local tree size allometries and integrate field information on the spatial positions of trees within stands, so as to mimic the coarse-scale heterogeneity observed in real field plots. However, our stand simulation procedure could still be improved by accounting for crown shapes and plasticity (like for instance in Boudon and Le Moguédec, 2007) or for between-species variations in trees inner-crown properties, either geometric (e.g. foliage clumping and porosity, leaf angle distribution) or optical (Schneider et al., 2014), that may influence grayscale variations at smaller spatial scales. Considering these different aspects is a long term research effort, and will likely benefit from technological advances in LiDAR scanning and unmanned aerial vehicles (e.g. Morton and Cook, 2016). We suspect lacunarity analysis to be more sensitive than FOTO to these geometric and optical parameters, making lacunarity a complement to FOTO appearing more relevant when analyzing simulated images than when analyzing real images.
An invaluable advantage of virtual canopy images is that they can be simulated in homogeneous acquisition configurations on the basis of field data from different sites and regions. Texture indeed quantifies pattern characteristics related to the contrasts between sunlit and shadowed surfaces and is thus highly sensitive to the sun to sensor and scene angles (Barbier et al. 2011). Even if inter-calibration procedures can help, particular configurations (e.g. backward scattering modes, near hotspot directions) are intrinsically detrimental to canopy texture analysis and fatal to large-scale applications over a mosaic of several images (Barbier and Couteron, 2015). Besides, programming image acquisition in controlled configurations, though possible, remains difficult, sometimes costly and often very long, particularly over tropical zones where optical images frequently suffer atmospheric pollutions related to dense cloud cover or persistent nebulosity. Another pitfall of the detection of texture gradients through ordination is that the principal axes extracted from the data are highly dependent on available samples and particular image sets may not reflect the entire possible gradients of forest structure and canopy texture variations. In our example sites for instance, it clearly appears that when the texture gradient is short, as it was the case in Paracou, the texture-structure relationship (as measured by R²) was not as strong as it was in areas covering a more extended gradient (as in Uppangala). Simulated images thus allowed us to place all plots and sites along a large gradient of canopy texture in order to stabilize the PCAs axes and thus to increase robustness of the texture-structure relationship, making it more universal than when fitted within sites. We therefore call for the collaborative construction of a pantropical database of 3D forest mockups, that could be incremented from plot data and local allometries in order to document all the particular situations encountered in tropical forests. These mockups could then be processed with DART in homogeneous acquisition conditions to produce arrays of canopy images of
sufficient generality. The relevant simulated spectral bands could finally be assembled
according to the specific characteristics of any of the VHSR sensors from which it could
appear realistic to gather a very large array of real world, homogeneous images as to
perform global scale canopy texture analysis. This would be a way to consolidate the
reference gradients in canopy texture that emerged from our simulation results as to
benchmark real observations, i.e. canopy texture analyses performed from real canopy
images. This perspective is a straightforward extension of the approach we introduced in the
present paper, and that would lead us towards a robust operational texture-structure
relationship from which accurate assessments of tropical forests AGB and carbon stocks
could be produced.

790

791 Acknowledgments

Pierre Ploton was supported by an Erasmus Mundus PhD grant from the 2013–2016 Forest,
Nature and Society (FONASO) doctoral program. Satellite images were acquired through the
Forest project funded by the European Institute of Technology-Climate Knowledge &
Innovation Community (grant N° PIN0040_2015-3.1-044_P032_01-01). Forest inventory data
were collected with the support of the CoForTips project as part of the ERA-Net BiodivERsA
2011-2012 European joint call (ANR-12-EBID-0002), the IRD project PPR FTH-AC
“Changements globaux, biodiversité et santé en zone forestière d’Afrique Centrale”, the
IFPCAR project “Controlling for Uncertainty in Assessment of Forest Aboveground Biomass in
the Western Ghats of India” (grant N° 4509-1), Eramet, the World Bank, WWF, the African
Development Bank and the Center for Tropical Forest Science – Forest Global Earth
Observatory (CTFS-ForestGEO) of the Smithsonian Tropical Research Institute. We are
803 thankfull to the Pallisco company for granting access to their concession, as well as local and
804 national authorities in Cameroon, India and DRC. Numerical experiments benefited from
805 high-performance computing resources of the Montpellier Bioinformatics Biodiversity
806 platform and the SouthGreen bioinformatics platform.

807

808 Reference

809 Allain, C., Cloitre, M., 1991. Characterizing the lacunarity of random and deterministic fractal
810 sets. Phys. Rev. A 44, 3552.

811 Antin, C., Pélissier, R., Vincent, G., Couteron, P., 2013. Crown allometries are less responsive
812 than stem allometry to tree size and habitat variations in an Indian monsoon forest.
813 Trees 27, 1485-1495. doi:10.1007/s00468-013-0896-7

814 Asner, G.P., Martin, R.E., 2008. Spectral and chemical analysis of tropical forests: Scaling
815 from leaf to canopy levels. Remote Sens. Environ. 112, 3958-3970.

816 Asner, G.P., Mascaro, J., Muller-Landau, H.C., Vieilledent, G., Vaudry, R., Rasamoelina, M.,
817 Hall, J.S., Breugel, M. van, 2011. A universal airborne LiDAR approach for tropical
818 forest carbon mapping. Oecologia 168, 1147-1160. doi:10.1007/s00442-011-2165-z

819 Baccini, A., Goetz, S.J., Walker, W.S., Laporte, N.T., Sun, M., Sulla-Menashe, D., Hackler, J.,
820 Beck, P.S.A., Dubayah, R., Friedl, M.A., 2012. Estimated carbon dioxide emissions
821 from tropical deforestation improved by carbon-density maps. Nat. Clim. Change 2,
822 182-185.

823 Barbier, N., Couteron, P., 2015. Attenuating the bidirectional texture variation of satellite
824 images of tropical forest canopies. Remote Sens. Environ. 171, 245-260.
825 doi:10.1016/j.rse.2015.10.007
826 Barbier, N., Couteron, P., Gastelly-Etchegorry, J.-P., Proisy, C., 2012. Linking canopy images
to forest structural parameters: potential of a modeling framework. Ann. For. Sci. 69,
827
828 305-311. doi:10.1007/s13595-011-0116-9

829 Barbier, N., Couteron, P., Proisy, C., Malhi, Y., Gastellu-Etchegorry, J.-P., 2010. The variation
830 of apparent crown size and canopy heterogeneity across lowland Amazonian forests.
831 Glob. Ecol. Biogeogr. 19, 72-84.

832 Barbier, N., Proisy, C., Véga, C., Sabatier, D., Couteron, P., 2011. Bidirectional texture
833 function of high resolution optical images of tropical forest: An approach using LiDAR
834 hillshade simulations. Remote Sens. Environ. 115, 167-179.

835 Baskerville, G.L., 1972. Use of Logarithmic Regression in the Estimation of Plant Biomass.
836 Can. J. For. Res. 2, 49-53. doi:10.1139/x72-009

837 Bastin, J.-F., Barbier, N., Couteron, P., Adams, B., Shapiro, A., Bogaert, J., De Cannière, C.,
838 2014. Aboveground biomass mapping of African forest mosaics using canopy texture
839 analysis: toward a regional approach. Ecol. Appl. 24, 1984-2001.

840 Bastin, J.-F., Barbier, N., Réjou-Méchain, M., Fayolle, A., Gourlet-Fleury, S., Maniatis, D., de
841 Haulleville, T., Baya, F., Beeckman, H., Beina, D., 2015. Seeing Central African forests
842 through their largest trees. Sci. Rep. 5. doi:doi:10.1038/srep13156

843 Boudon, F., Le Moguédec, G.L., 2007. Déformation asymétrique de houppiers pour la
génération de représentations paysagères réalistes. Rev. Electron. Francoph. Inform.
844
845 Graph. 1.

846 Breiman, L., 2001. Random forests. Mach. Learn. 45, 5-32.

847 Broadbent, E.N., Asner, G.P., Peña-Claros, M., Palace, M., Soriano, M., 2008. Spatial
848 partitioning of biomass and diversity in a lowland Bolivian forest: Linking field and
849 remote sensing measurements. For. Ecol. Manag. 255, 2602-2616.

850 doi:10.1016/j.foreco.2008.01.044

851 Chave, J., Coomes, D., Jansen, S., Lewis, S.L., Swenson, N.G., Zanne, A.E., 2009. Towards a

852 worldwide wood economics spectrum. Ecol. Lett. 12, 351-366. doi:10.1111/j.1461-

853 0248.2009.01285.x

854 Chave, J., Muller-Landau, H.C., Baker, T.R., Easdale, T.A., Steege, H. ter, Webb, C.O., 2006.

855 Regional and phylogenetic variation of wood density across 2456 neotropical tree

856 species. Ecol. Appl. 16, 2356-2367. doi:10.1890/1051-

857 0761(2006)016[2356:RAPVOW]2.0.CO;2

858 Chave, J., Réjou-Méchain, M., Búrquez, A., Chidumayo, E., Colgan, M.S., Delitti, W.B.C.,

859 Duque, A., Eid, T., Fearnside, P.M., Goodman, R.C., Henry, M., Martínez-Yrízar, A.,

860 Mugasha, W.A., Muller-Landau, H.C., Mencuccini, M., Nelson, B.W., Ngomanda, A.,

861 Nogueira, E.M., Ortiz-Malavassi, E., Pélissier, R., Ploton, P., Ryan, C.M., Saldarriaga,

862 J.G., Vieilledent, G., 2014. Improved allometric models to estimate the aboveground

863 biomass of tropical trees. Glob. Change Biol. 20, 3177-3190. doi:10.1111/gcb.12629

864 Chuyong, G.B., Condit, R., Kenfack, D., Losos, E.C., Moses, S.N., Songwe, N.C., Thomas, D.W.,

865 2004. Korup forest dynamics plot, Cameroon. Trop. For. Divers. Dynamism 506-516.

866 Couteron, P., 2002. Quantifying change in patterned semi-arid vegetation by Fourier analysis

867 of digitized aerial photographs. Int. J. Remote Sens. 23, 3407-3425.

868 Couteron, P., Pelissier, R., Nicolini, E.A., Paget, D., 2005. Predicting tropical forest stand

869 structure parameters from Fourier transform of very high-resolution remotely sensed

870 canopy images. J. Appl. Ecol. 42, 1121-1128.

871 Cressie, N.A., 1993. Statistics for spatial data: Wiley Series in Probability and Mathematical

872 Statistics. Find This Article Online.
873 DeFries, R., Achard, F., Brown, S., Herold, M., Murdiyarso, D., Schlamadinger, B., de Souza, C., 2007. Earth observations for estimating greenhouse gas emissions from deforestation in developing countries. Environ. Sci. Policy 10, 385-394.

876 Dolédec, S., Chessel, D., 1994. Co-inertia analysis: an alternative method for studying species-environment relationships. Freshw. Biol. 31, 277-294.

878 Erdody, T.L., Moskal, L.M., 2010. Fusion of LiDAR and imagery for estimating forest canopy fuels. Remote Sens. Environ. 114, 725-737.

880 Fayad, I., Baghdadi, N., Guitet, S., Bailly, J.-S., Hérault, B., Gond, V., El Hajj, M., Minh, D.H.T., 2016. Aboveground biomass mapping in French Guiana by combining remote sensing, forest inventories and environmental data. Int. J. Appl. Earth Obs. Geoinformation 52, 502-514.

884 Feldpausch, T.R., Lloyd, J., Lewis, S.L., Brienen, R.J., Gloor, M., Monteagudo Mendoza, A., Lopez-Gonzalez, G., Banin, L., Abu Salim, K., Affum-Baffoe, K., 2012. Tree height integrated into pantropical forest biomass estimates. Biogeosciences 3381-3403.

887 Franklin, J.F., Spies, T.A., Van Pelt, R., Carey, A.B., Thornburgh, D.A., Berg, D.R., Lindenmayer, D.B., Harmon, M.E., Keeton, W.S., Shaw, D.C., 2002. Disturbances and structural development of natural forest ecosystems with silvicultural implications, using Douglas-fir forests as an example. For. Ecol. Manag. 155, 399-423.

891 Frazer, G.W., Wulder, M.A., Niemann, K.O., 2005. Simulation and quantification of the fine-scale spatial pattern and heterogeneity of forest canopy structure: A lacunarity-based method designed for analysis of continuous canopy heights. For. Ecol. Manag. 214, 65-90.

895 Gastellu-Etchegorry, J.-P., Yin, T., Lauret, N., Cajgfinger, T., Gregoire, T., Grau, E., Feret, J.-B., Lopes, M., Guilleux, J., Dedieu, G., 2015. Discrete Anisotropic Radiative Transfer
897 (DART 5) for modeling airborne and satellite spectroradiometer and LIDAR
898 acquisitions of natural and urban landscapes. Remote Sens. 7, 1667-1701.

899 Guariguata, M.R., Ostertag, R., 2001. Neotropical secondary forest succession: changes in
900 structural and functional characteristics. For. Ecol. Manag. 148, 185-206.

901 Hajjem, A., Bellavance, F., Larocque, D., 2014. Mixed-effects random forest for clustered
902 data. J. Stat. Comput. Simul. 84, 1313-1328.

903 Haralick, R.M., 1979. Statistical and structural approaches to texture. Proc. IEEE 67, 786-804.

904 ICRAF, 2007. Wood density database. World Agrofor. Cent. Nairobi Kenya
905 http://db.worldagroforestry.org//wd.

906 Jeyakumar, S., Ayyappan, N., Muthuramkumar, S., Rajarathinam, K., in press. Impacts of
907 selective logging on diversity, species composition and biomass of residual lowland
908 dipterocarp forest in central Western Ghats, India. Trop. Ecol.

909 Jucker, T., Caspersen, J., Chave, J., Antin, C., Barbier, N., Bongers, F., Dalponte, M., van Ewijk,
910 K.Y., Forrester, D.I., Haeni, M., Higgins, S.I., Holdaway, R.J., Iida, Y., Lorimer, C.,
911 Marshall, P.L., Momo, S., Moncrieff, G.R., Ploton, P., Poorter, L., Rahman, K.A.,
912 Schlund, M., Sonké, B., Sterck, F.J., Trugman, A.T., Usoltsev, V.A., Vanderwel, M.C.,
913 Waldner, P., Wedeux, B.M.M., Wirth, C., Wöll, H., Woods, M., Xiang, W.,
914 Zimmermann, N.E., Coomes, D.A., 2016. Allometric equations for integrating remote
915 sensing imagery into forest monitoring programmes. Glob. Change Biol.
916 doi:10.1111/gcb.13388

917 Lu, D., 2006. The potential and challenge of remote sensing-based biomass estimation. Int. J.
918 Remote Sens. 27, 1297-1328.

919 Lu, D., Chen, Q., Wang, G., Moran, E., Batistella, M., Zhang, M., Vaglio Laurin, G., Saah, D.,
920 2012. Aboveground Forest Biomass Estimation with Landsat and LiDAR Data and
921 Uncertainty Analysis of the Estimates. Int. J. For. Res. 2012, 1–16.

922 doi:10.1155/2012/436537

923 Malhi, Y., Román-Cuesta, R.M., 2008. Analysis of lacunarity and scales of spatial

924 homogeneity in IKONOS images of Amazonian tropical forest canopies. Remote Sens.

925 Environ. 112, 2074–2087.

926 Mandelbrot, B.B., 1983. The fractal geometry of nature. Macmillan.

927 Marvin, D.C., Koh, L.P., Lynam, A.J., Wich, S., Davies, A.B., Krishnamurthy, R., Stokes, E.,

928 Starkey, R., Asner, G.P., 2016. Integrating technologies for scalable ecology and

929 conservation. Glob. Ecol. Conserv. 7, 262–275. doi:10.1016/j.gecco.2016.07.002

930 Mascaro, J., Asner, G.P., Knapp, D.E., Kennedy-Bowdoin, T., Martin, R.E., Anderson, C.,

931 Higgins, M., Chadwick, K.D., 2014. A Tale of Two “Forests”: Random Forest Machine

932 Learning Aids Tropical Forest Carbon Mapping. PLoS ONE 9, e85993.

933 doi:10.1371/journal.pone.0085993

934 Meng, S., Pang, Y., Zhang, Z., Jia, W., Li, Z., 2016. Mapping Aboveground Biomass using

935 Texture Indices from Aerial Photos in a Temperate Forest of Northeastern China.

936 Remote Sens. 8, 230. doi:10.3390/rs8030230

937 Mermoz, S., Réjou-Méchain, M., Villard, L., Le Toan, T., Rossi, V., Gourlet-Fleury, S., 2015.

938 Decrease of L-band SAR backscatter with biomass of dense forests. Remote Sens.

939 Environ. 159, 307–317. doi:10.1016/j.rse.2014.12.019

940 Messinger, M., Asner, G., Silman, M., 2016. Rapid Assessments of Amazon Forest Structure

941 and Biomass Using Small Unmanned Aerial Systems. Remote Sens. 8, 615.

942 doi:10.3390/rs8080615

943 Morton, D.C., Cook, B.D., 2016. Amazon forest structure generates diurnal and seasonal

944 variability in light utilization. Biogeosciences 13, 2195.
945 Pan, Y., Birdsey, R.A., Fang, J., Houghton, R., Kauppi, P.E., Kurz, W.A., Phillips, O.L.,
946 Shvidenko, A., Lewis, S.L., Canadell, J.G., 2011. A large and persistent carbon sink in
947 the world's forests. Science 333, 988-993.

948 Pargal, S., Fararoda, R., Rajashekar, G., Balachandran, N., Réjou-Méchain, M., Barbier, N.,
949 Jha, C.S., Pélissier, R., Dadhwal, V.K., Couteron, P., 2017. Inverting Aboveground
950 Biomass-Canopy Texture Relationships in a Landscape of Forest Mosaic in the
951 Western Ghats of India Using Very High Resolution Cartosat Imagery. Remote Sens. 9,
952 228.

953 Ploton, P., Barbier, N., Takoudjou Momo, S., Réjou-Méchain, M., Boyemba Bosela, F.,
954 Chuyong, G., Dauby, G., Droissart, V., Fayolle, A., Goodman, R.C., Henry, M.,
955 Kamdem, N.G., Mukirania, J.K., Kenfack, D., Libalah, M., Ngomanda, A., Rossi, V.,
956 Sonké, B., Texier, N., Thomas, D., Zebaze, D., Couteron, P., Berger, U., Pélissier, R.,
957 2016. Closing a gap in tropical forest biomass estimation: taking crown mass variation
958 into account in pantropical allometries. Biogeosciences 13, 1571-1585.
959 doi:10.5194/bg-13-1571-2016

960 Ploton, P., Pélissier, R., Proisy, C., Flavenot, T., Barbier, N., Rai, S.N., Couteron, P., 2012.
961 Assessing aboveground tropical forest biomass using Google Earth canopy images.
962 Ecol. Appl. 22, 993-1003. doi:10.1890/11-1606.1

963 Proisy, C., Couteron, P., Fromard, F., 2007. Predicting and mapping mangrove biomass from
964 canopy grain analysis using Fourier-based textural ordination of IKONOS images.
965 Remote Sens. Environ. 109, 379-392.

966 Proisy, C., Féret, J.-B., Lauret, N., Gastellu-Etchegorry, J.-P., 2016. Mangrove forest dynamics
967 using very high spatial resolution optical remote sensing., in: Remote Sensing of Land
968 Surfaces: Urban and Coastal Area. N.N. Baghdadi & M. Zribi, Paris, pp. 274-300.
969 R Core Team, 2016. R: A language and environment for statistical computing. R Foundation
970 for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

971 Réjou-Méchain, M., Tanguy, A., Piponiot, C., Chave, J., Hérault, B., 2017. biomass: an r
972 package for estimating above-ground biomass and its uncertainty in tropical forests.
973 Methods Ecol. Evol.

974 Réjou-Méchain, M., Tymen, B., Blanc, L., Fauset, S., Feldpausch, T.R., Monteagudo, A.,
975 Phillips, O.L., Richard, H., Chave, J., 2015. Using repeated small-footprint LiDAR
976 acquisitions to infer spatial and temporal variations of a high-biomass Neotropical
977 forest. Remote Sens. Environ. 169, 93–101.

978 Rutishauser, E., Herault, B., Petronelli, P., Sist, P., 2015. Tree height reduction after selective
979 logging in a tropical forest. Biotropica.

980 Saatchi, S.S., Harris, N.L., Brown, S., Lefsky, M., Mitchard, E.T., Salas, W., Zutta, B.R.,
981 Buermann, W., Lewis, S.L., Hagen, S., 2011. Benchmark map of forest carbon stocks in
tropical regions across three continents. Proc. Natl. Acad. Sci. 108, 9899–9904.
982

983 Schneider, F.D., Leiterer, R., Morsdorf, F., Gastellu-Etchegorry, J.-P., Lauret, N., Pfeifer, N.,
984 Schaepman, M.E., 2014. Simulating imaging spectrometer data: 3D forest modeling
985 based on LiDAR and in situ data. Remote Sens. Environ. 152, 235–250.
986 doi:10.1016/j.rse.2014.06.015

987 Singh, M., Evans, D., Friess, D.A., Tan, B.S., Nin, C.S., 2015. Mapping Above-Ground Biomass
988 in a Tropical Forest in Cambodia Using Canopy Textures Derived from Google Earth.
989 Remote Sens. 7, 5057–5076. doi:10.3390/rs70505057

990 Singh, M., Malhi, Y., Bhagwat, S., 2014. Biomass estimation of mixed forest landscape using a
991 Fourier transform texture-based approach on very-high-resolution optical satellite
992 imagery. Int. J. Remote Sens. 35, 3331–3349. doi:10.1080/01431161.2014.903441
993 Slik, J.W., Paoli, G., McGuire, K., Amaral, I., Barroso, J., Bastian, M., Blanc, L., Bongers, F.,
994 Boundja, P., Clark, C., 2013. Large trees drive forest aboveground biomass variation
995 in moist lowland forests across the tropics. Glob. Ecol. Biogeogr. 22, 1261–1271.

996 Spies, T.A., 1998. Forest structure: a key to the ecosystem. Northwest Sci. 72, 34–36.

997 Stark, S.C., Enquist, B.J., Saleska, S.R., Leitold, V., Schietti, J., Longo, M., Alves, L.F., Camargo,
998 P.B., Oliveira, R.C., 2015. Linking canopy leaf area and light environments with tree
999 size distributions to explain Amazon forest demography. Ecol. Lett. n/a-n/a.
1000 doi:10.1111/ele.12440

1001 Taubert, F., Jahn, M.W., Dobner, H.-J., Wiegand, T., Huth, A., 2015. The structure of tropical
1002 forests and sphere packings. Proc. Natl. Acad. Sci. 112, 15125–15129.

1003 Tucker, L.R., 1958. An inter-battery method of factor analysis. Psychometrika 23, 111–136.

1004 Vieilledent, G., Gardi, O., Grinand, C., Burren, C., Andriamanjato, M., Camara, C., Gardner,
1005 C.J., Glass, L., Rasolohery, A., Rakoto Ratsimba, H., Gond, V., Rakotoarijaona, J.-R.,
1006 2016. Bioclimatic envelope models predict a decrease in tropical forest carbon stocks
1007 with climate change in Madagascar. J. Ecol. 104, 703–715. doi:10.1111/1365-
1008 2745.12548

1009 Vincent, G., Antin, C., Laurans, M., Heurtebize, J., Durrieu, S., Lavalley, C., Dauzat, J., 2017.
1010 Mapping plant area index of tropical evergreen forest by airborne laser scanning. A
1011 cross-validation study using LAI2200 optical sensor. Remote Sens. Environ. 198, 254–
1012 266.

1013 Vincent, G., Sabatier, D., Blanc, L., Chave, J., Weissenbacher, E., Pélissier, R., Fonty, E.,
1014 Molino, J.-F., Couteron, P., 2012. Accuracy of small footprint airborne LiDAR in its
1015 predictions of tropical moist forest stand structure. Remote Sens. Environ. 125, 23–
1016 33.
1017 Vincent, G., Sabatier, D., Rutishauser, E., 2014. Revisiting a universal airborne light detection and ranging approach for tropical forest carbon mapping: scaling-up from tree to stand to landscape. Oecologia 175, 439-443. doi:10.1007/s00442-014-2913-y

1020 Withmore, T.C., 1975. Tropical rain forest of the far east. Claredon Prees Oxf. Univ. Prees Lond.

1022 Wu, H., Zhang, J.-T., 2006. Nonparametric regression methods for longitudinal data analysis: mixed-effects modeling approaches. John Wiley & Sons.

1024 Xu, L., Saatchi, S.S., Yang, Y., Yu, Y., White, L., 2016. Performance of non-parametric algorithms for spatial mapping of tropical forest structure. Carbon Balance Manag. 11. doi:10.1186/s13021-016-0062-9

1027 Zanne, A.E., Lopez-Gonzalez, G., Coomes, D.A., Ilic, J., Jansen, S., Lewis, S.L., Miller, R.B., Swenson, N.G., Wiemann, M.C., Chave, J., 2009. Data from: towards a worldwide wood economics spectrum. Dryad Digital Repository.

1030 Zhao, P., Lu, D., Wang, G., Wu, C., Huang, Y., Yu, S., 2016. Examining Spectral Reflectance Saturation in Landsat Imagery and Corresponding Solutions to Improve Forest Aboveground Biomass Estimation. Remote Sens. 8, 469. doi:10.3390/rs8060469

1033 Zhou, J., Proisy, C., Descombes, X., Le Maire, G., Nouvellon, Y., Stape, J.-L., Viennois, G., Zerubia, J., Couteron, P., 2013. Mapping local density of young Eucalyptus plantations by individual tree detection in high spatial resolution satellite images. For. Ecol. Manag. 301, 129-141.
APPENDIX

1037

1038

1039 Table A1. Distribution of forest inventory data among sampling sites. Tree dimensions collected in field plots include the diameter at breast

1040 height (D), tree height (H), trunk height (Ht) and crown diameter (Cd).

<table><thead><tr><th rowspan="2">Source</th><th rowspan="2">Country</th><th rowspan="2">Site</th><th rowspan="2">Plot (ha)</th><th colspan="4">Tree dimensions</th><th colspan="4">Taxonomic identification (in %)</th></tr><tr><th>D</th><th>H</th><th>Ht</th><th>Cd</th><th>Species</th><th>Genus</th><th>Family</th><th>unidentified</th></tr></thead><tbody><tr><td>Bastin et al. (2014)</td><td>DRC</td><td>Malebo</td><td>28</td><td>10092</td><td>10092</td><td>0</td><td>0</td><td>96,5</td><td>3,5</td><td>0</td><td>0</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Deng-Deng</td><td>15</td><td>7847</td><td>343</td><td>342</td><td>42</td><td>62,1</td><td>17</td><td>17</td><td>3,9</td></tr><tr><td>Rutishauser et al. (2015) & Vincent et al. (2012)</td><td>French Guyana</td><td>Paracou</td><td>85</td><td>53685</td><td>2191</td><td>1667</td><td>1740</td><td>85,3</td><td>0</td><td>0</td><td>14,7</td></tr><tr><td>CTFS - Chuyong et al. (2004)</td><td>Cameroon</td><td>Korup</td><td>50</td><td>24590</td><td>0</td><td>0</td><td>0</td><td>95,8</td><td>4</td><td>0,2</td><td>0</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Djomedjo</td><td>10</td><td>4959</td><td>505</td><td>504</td><td>186</td><td>65,4</td><td>14,4</td><td>17,7</td><td>2,5</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Lomié</td><td>8</td><td>3571</td><td>539</td><td>128</td><td>0</td><td>75,5</td><td>7,7</td><td>11,1</td><td>5,7</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Mabounié</td><td>12</td><td>4475</td><td>592</td><td>525</td><td>113</td><td>53,4</td><td>19,2</td><td>24,4</td><td>3</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Mindourou I</td><td>11</td><td>5163</td><td>908</td><td>559</td><td>306</td><td>61,6</td><td>16,8</td><td>15,1</td><td>6,5</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Mindourou II</td><td>12</td><td>5598</td><td>920</td><td>667</td><td>233</td><td>56,4</td><td>26,2</td><td>11,7</td><td>5,7</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Ngoila</td><td>6</td><td>2611</td><td>126</td><td>161</td><td>0</td><td>69,6</td><td>18,5</td><td>6,1</td><td>5,8</td></tr><tr><td>Unpublished (IRD-AMAP/ENS)</td><td>Cameroon</td><td>Nemeyong</td><td>5</td><td>2824</td><td>260</td><td>227</td><td>90</td><td>65,6</td><td>17,2</td><td>14,6</td><td>2,7</td></tr><tr><td>Ploton et al. (2012)</td><td>India</td><td>Uppangala</td><td>15</td><td>8220</td><td>1729</td><td>1722</td><td>1728</td><td>0</td><td>0</td><td>0</td><td>100<sup>1</sup></td></tr><tr><td>Pargal et al. (2017)</td><td>India</td><td>Yellapur</td><td>22</td><td>9156</td><td>5032</td><td>0</td><td>0</td><td>97,8</td><td>0,4</td><td>0</td><td>1,8</td></tr></tbody></table>

1041

1042 ¹ AGB values at Uppangala were computed using the average wood density of four 1-ha plots in the area (from Jeyakumar et al., *in press*)
1043

List of figure captions

1044

Figure 1. Typical workflow of texture-based AGB mapping method.

1045

1046

1047

1048

1049

1050

Figure 2. Canopy texture ordinations based on (A) the FOTO method and (B) the lacunarity analysis. In both cases, scatter plots of PCA scores along the first two principal axes are shown (left panels), with 3 example sites highlighted with particular symbols (Paracou, Uppangala and Yellapur). Top-right panels represent correlation circles with wavelength, $\lambda$ (A) or box size, $s$ (B) in meter. Bottom-right panels represent histograms of eigenvalues in % of total variance.

1051

1052

1053

1054

1055

1056

1057

1058

Figure 3. Co-inertia analysis. Position on the first co-inertia plane of the FOTO r-spectra wavelengths, $\lambda$ (A) and the lacunarity box size, $s$ (D). Components of the F-PCA (B) and the L-PCA (E) projected onto the co-inertia axes. Ordination of windows from the 3 example sites (Paracou, Uppangala and Yellapur) on COIA-1 (C) with large empty and full circles representing the average site-level score for FOTO and Lacunarity features, respectively. Normed scores of 10 randomly sampled canopy windows from the 3 example sites on the first co-inertia plane (F), with each arrow linking a unit-window position for FOTO and Lacunarity characteristics, respectively.

1059

1060

1061

1062

1063

Figure 4. Aboveground biomass prediction models based on FOTO texture (F-model), Lacunarity texture (L-model) and the two sources of texture information (FL-model), either without site effect (RF algorithm, panels A, C and E) or with site effect (MERF algorithm, panels B, D and F). Texture features were extracted from virtual canopy scenes. Goodness of fit statistics are defined in the Methods.
1064 Figure 5. (A) Bioclimatic proxy E vs random site intercepts of FL-model (in Fig. 4-E) ; (B)
1065 Random forest AGB prediction model based on FOTO texture, lacunarity texture and E (FLE-
1066 model). Texture features were extracted from virtual canopy scenes. Goodness of fit
1067 statistics are defined in the Methods.

1068 Figure 6. Aboveground biomass prediction models based on FOTO texture, Lacunarity
1069 texture and the bioclimatic stress proxy E (FLE-model). Texture features were extracted from
1070 three Pleiades covering 49 1-ha forest inventory plots in central Africa. Goodness of fit
1071 statistics are defined in Methods section.

1072

1073