Incorporating Forest Mapping-Related Uncertainty into the Error Propagation of Wall-to-Wall Biomass Maps: A General Approach for Large and Small Areas

Abstract

The sources of uncertainty in wall-to-wall AGB maps propagate from the tree to pixel, but uncertainty due to forest cover mapping is rarely incorporated into the error propagation process. This study aimed to (1) elaborate an analytical procedure to incorporate forest-mapping-related uncertainty into the error propagation from plot and pixel predictions; (2) develop a stratified estimator with a model-assisted estimator for small and large areas; and (3) estimate the effect of ignoring the mapping uncertainty on the confidence intervals (CIs) for totals. Data consist of a subset of the Brazilian national forest inventory (NFI) database, comprising 75 counties that, once aggregated, served as strata for the stratified estimator. On-ground data were gathered from 152 clusters (plots) and remotely sensed data from Landsat-8 scenes. Four major contributions are highlighted. First, we describe how to incorporate forest-mapping-related uncertainty into the CIs of any forest attribute and spatial resolution. Second, stratified estimators perform better than non-stratified estimators for forest area estimation when the response variable is forest/non-forest. Comparing our stratified estimators, this study indicated greater precision for the stratified estimator than for the regression estimator. Third, using the ratio estimator, we found evidence that the simple field plot information provided by the NFI clusters is sufficient to estimate the proportion forest for large regions as accurately as remote-sensing-based methods, albeit with less precision. Fourth, ignoring forest-mapping-related uncertainty erroneously narrows the CI width as the estimate of proportion forest area decreases. At the small-area level, forest-mapping-related uncertainty led to CIs for total AGB as much as 63% wider in extreme cases. At the large-area level, the CI was 5–7% wider.

1. Introduction

Many countries conduct their national forest inventories (NFIs) following a sampling effort to gain knowledge about their forest resources and to support national forest use and regulation policies. Typically, NFI sampling designs involve systematic distributions of permanent ground plots throughout the country. This system plays a key role by providing information to estimate and map forest attributes at the national level [1], such as aboveground biomass (AGB). Because of the possibility of converting biomass into carbon, information on biomass stock in forests at the national level serves as a basis for supporting the carbon budget, greenhouse gas emissions, and carbon sequestration reports [2,3,4].

Despite of the national emphasis of NFIs, there is an increasing demand for NFI-based estimates of AGB disaggregated into small areas such as counties and administrative regions [5,6]. At these smaller scales, forest information is useful to guide local policies that are part of the national forest policy. Therefore, some countries have already estimated their forest stocks (based on NFI data) by administrative areas [7,8,9,10], thereby playing an essential role in local governmental measures.

In general, NFI plot grids lack the resolution necessary to generate ground-data-based estimates with sufficient precision, requiring refinement through the augmentation of NFI data with auxiliary data [4]. Many techniques and strategies for refining NFI estimates from coarse- to fine-scales have been tested, with most involving remotely sensed data such as satellite imagery and/or airborne LiDAR measurements [3,5,11,12,13,14,15,16]. Research reported by [5] has shown that combining ground and remotely sensed data successfully yields more precise estimates. This is useful both for producing wall-to-wall (large spatial) maps of forest attributes like AGB [17] and for mapping these attributes over smaller areas.

One concern is that mapping forest attributes relies on a long chain of field data measurements that extends predictions from trees to plots and image pixels [18]. Thus, it is crucial to assess and estimate the uncertainties associated with measurements and variable predictions [19] and their propagation to the uncertainty of the mapped forest attributes. While these uncertainties can be estimated, the main role of uncertainty analysis is to provide more robust confidence intervals (CI) for population parameters. The most relevant uncertainties in wall-to-wall biomass maps come basically from six sources [10,18,20]: (1) tree-level measurements (e.g., tree diameter and biomass), (2) tree-level AGB model, (3) plot-level AGB model, (4) pixel-level AGB model, (5) slope correction factor (for fixed-area plots), and (6) forest cover mapping.

It is well-known that these cited uncertainties are propagated to the final estimates [15,18,21]. However, the uncertainty inherent in the modeling process is rarely propagated to uncertainties in the forest area estimates, the latter resulting from image classification. The aim of this study is threefold: (1) elaborate an analytical procedure to incorporate forest-mapping-related uncertainty into the error propagation from plot and pixel predictions; (2) develop a stratified estimator with a model-assisted estimator for small and large areas; and (3) estimate the effect of ignoring the mapping uncertainty on the confidence intervals (CIs) for totals. This study was further motivated by the urgent need for an enhanced procedure for providing CIs for small-area AGB estimates.

2. Materials and Methods

2.1. Study Area and Data Source

The study area covers approximately 55,330 km² or about 28% of the State of Paraná, Brazil, and includes L = 75 counties. This region is part of the Brazilian Atlantic Forest biome, where sub-tropical rainforests cover most of the study area. The elevation ranges from 320 to 1300 m. According to the Köppen classification [22], the region predominantly has a subtropical climate without dry seasons.

The NFI uses a systematic sampling design to allocate clusters across the country, following a regular grid of 20 × 20 km, meaning that each cluster represents a 400 km² cell. Each NFI cluster covers an area of 4000 m² and contains four 1000 m² subunits arranged crosswise. Each subunit is further divided into 10 subplots of 100 m², resulting in a total of 40 subplots per cluster [5]. Data were obtained from a sub-database of the Brazilian NFI, which includes n_c = 152 clusters installed in the study area between February and May of 2013. Figure 1 shows the study area and locations of the clusters.

The sub-database used in this study comprised 10,822 live trees with a diameter at breast height (dbh) ≥ 10 cm. All trees had their dbh and total height measured. The main descriptive statistics of the tree variables are provided in

Table 1. Descriptive statistics of the variables (n = 10,822) for the trees sampled in the NFI sub-dataset used in this study.

Tree Variable	Minimum	Medium	Maximum	CV (%)
Diameter at breast height (cm)	10.0	18.4	118.3	52.0%
Total height (m)	1.4	10.4	28.0	34.9%
Number of trees ha−1 (dbh ≥ 10 cm)	3	237	1,208	91.4%
Aboveground biomass * (kg tree−1)	3.7	151.9	6,476.0	184.7%

* predicted through Equations (25)–(28) of [23].

.

2.2. Analytical Procedure

The forest attribute of interest is forest AGB (Mg ha⁻¹), expressed as the sum of the AGB of woody trees, palms, and ferns. We predicted tree biomass, which was then summed to correspond to predicted plot-level (i.e., cluster-level) AGB per unit area. Next, plot-level AGB was used as the output variable for a pixel-level model. A forest cover map was produced to ultimately map AGB across the study area. This study assumes that the multi-level errors in the modeling and mapping processes are propagated to the final estimate. Figure 2 shows a flowchart of the analytical procedure outlined in this study.

The following subsections describe the analytical procedures for executing satellite image pre-processing, image classification, map accuracy, tree AGB modeling, remote-sensing-based modeling of forest AGB, and uncertainty analysis. The notation in the glossary aids in understading our methods.

2.2.1. Satellite Image Pre-Processing

Satellite image analysis was required to predict AGB at the pixel level in Mg ha⁻¹. Landsat 8 Operational Land Imager (OLI) images were downloaded from the USGS database. Seven scenes were needed to cover the entire study area. The path/row numbers are as follows: 221/077, 221/078, 222/077, 222/078, 223/076, 223/077, and 223/078. The period from December 2013 to February 2014 was chosen for acquiring the satellite images due to the low (0.25%) cloud cover observed in all images and of the alignment with the field data collection period.

The image pre-processing steps necessary for subsequent analyses were as follows: (i) conversion of digital numbers to radiance scale; (ii) conversion of radiance to surface reflectance, the latter performed using the Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) algorithm, which performs the atmospheric correction needed to remove most haze and other atmospheric disturbances [24]; (iii) re-projection to the SIRGAS 2000 reference system; (iv) mosaicking, during which adjustments to image brightness were required; and (v) clipping using the shape of the study area.

2.2.2. Image Classification

Since the NFI sampling design follows a systematic grid, clusters can be entirely or partially allocated in natural forests or in non-forested areas. Each 100 m² subplot is categorized by the predominant land cover that occupies it. This categorization includes 19 land cover types, although only two were observed across the inventoried area: early successional forest and medium- to late-successional forest. For image classification, these two categories were merged into class ‘forest’, while ‘non-forest’ was assigned to all other land covers. Therefore, the classes ‘forest’ and ‘non-forest’ were our targets. Planted forest was not considered.

Image classification was performed via object-based image analysis (OBIA) to map forest cover (forest/non-forest). OBIA has been widely applied to medium-resolution images (e.g., Landsat) [25,26], providing enhanced land cover classifications, as shown in [27,28,29]. We followed the same segmentation and classification criteria presented by [5]. The salt-and-pepper effect was not corrected because it is minimal in OBIA classifications [30].

2.2.3. Map Accuracy

To assess the accuracy of the forest cover map, we constructed an error matrix and then estimated overall accuracy, user’s accuracy, and producer’s accuracy [31]. Accuracy assessment was conducted at the greater region level, defined as the grand strata, to meet the need for accuracy estimates for these regions in subsequent uncertainty analyses. These greater regions, which are administrative areas within the State of Parana, were chosen due to their relatively uniform levels of agricultural expansion (which is, consequently, related to forest conservation).

Ground reference data were collected from 850 observations, sampled whenever possible from a random point within a 4 km radius around the $n_c$= 152 clusters. Reference data were verified using high-resolution Google Earth imagery as auxiliary data, dating from the same period as the field observations. Although the NFI clusters are systematically distributed on a regular grid, reference data for map accuracy followed a pre-stratified random sampling design.

Table 2. Error matrix expressed in terms of proportion of area and accuracy estimators.

r = Columns: Reference		Non-Forest	Forest	$\mathbf{Total} \left({\mathit{m}}_{\mathit{c}+}\right)$	User’s Accuracy
c = rows: classification	Non-forest	${\displaystyle \raisebox{1ex}{$m_{00}$}\!\left/ \!\raisebox{-1ex}{$m_h$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{01}$}\!\left/ \!\raisebox{-1ex}{$m_h$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{0+}$}\!\left/ \!\raisebox{-1ex}{$m_h$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{00}$}\!\left/ \!\raisebox{-1ex}{$m_{0+}$}\right.}$
	Forest	${\displaystyle \raisebox{1ex}{$m_{10}$}\!\left/ \!\raisebox{-1ex}{$m_h$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{11}$}\!\left/ \!\raisebox{-1ex}{$m_h$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{1+}$}\!\left/ \!\raisebox{-1ex}{$m_h$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{11}$}\!\left/ \!\raisebox{-1ex}{$m_{1+}$}\right.}$
	$\mathrm{Total} \left({m }_{+r}\right)$	$m_{+0}$	$m_{+1}$	$m_h$
	Producer’s accuracy	${\displaystyle \raisebox{1ex}{$m_{00}$}\!\left/ \!\raisebox{-1ex}{${m }_{+0}$}\right.}$	${\displaystyle \raisebox{1ex}{$m_{11}$}\!\left/ \!\raisebox{-1ex}{$m_{+1}$}\right.}$	$OA={\displaystyle \frac{m_{00}+m_{11}}{m_h}}$

OA: overall accuracy. $m_{rc}$: number of sample units (pixels) classified into category r with reference to classification c in a grand stratum h.

presents the error matrix for map accuracy assessment.

Image classification allowed us to estimate (i) the total forest area, $T\widehat{F}A$ (ha), and (ii) the counties’ forest area, ${C\widehat{F}A}_h$ (in ha). Since we have L strata (counties), the total forest area is given by $T\widehat{F}A=\sum _{h=1}^{L }{C\widehat{F}A}_h$. To estimate ${C\widehat{F}A}_h$, we simply clipped the forest/non-forest raster (the final product of the image classification) using vector files representing each county’s boundaries.

2.2.4. Tree AGB Modeling

Dbh and total height were used as predictors in the models of [23] (their Equations (25)–(28)) to predict tree AGB. These authors destructively measured the biomass of 387 trees in forest remnants within the sub-tropical rainforest across southern Brazil, representing the forest typology covering ~90% of our study area. This biomass collection effort is part of a larger project, Tropical Biomass & Carbon [32]. The sum of tree AGB within each cluster corresponded to the cluster AGB.

2.2.5. Remote-Sensing-Based Modeling of Forest AGB

This step involved modeling forest AGB at the pixel level. A variable screening of nine input predictors derived from the OLI sensor was performed prior to modeling. The predictors fall into two categories: (i) surface reflectance of OLI bands (Red, NIR, SWIR1, and SWIR2); and (ii) vegetation indices (

Table 3. Vegetation indices used in the biomass model.

Vegetation Index	Formulation	Reference
DVI	${\rho }_{NIR}-{ \rho }_{Red}$	[35]
EVI	$2.5{\displaystyle \frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+6 {\rho }_{Red}-7.5 {\rho }_{Blue}+1}}$	[36]
NDMI	${\displaystyle \frac{{\rho }_{NIR}-{ \rho }_{SWIR1}}{{\rho }_{NIR}+{\rho }_{SWIR1}}}$	[37]
NDVI	${\displaystyle \frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}}$	[38]
PVI	${\displaystyle \frac{{\rho }_{NIR}-a\times {\rho }_{Red}-b}{\sqrt{1+a2}}}$	[39]
SAVI	$1.5{\displaystyle \frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}+0.5}}$	[40]
SR	${\displaystyle \frac{{\rho }_{Red}}{{\rho }_{NIR}}}$	[41]
TVI	${\displaystyle \frac{120\left({\rho }_{NIR}-{\rho }_{Green}\right)-200({\rho }_{Red}-{ \rho }_{Green})}{2}}$	[42]

ρ: surface reflectance. a = 0.96916, b = 0.084726.

). These categories are commonly combined in biomass modeling [6,29,33]. A known limitation in remote sensing is the saturation of the Normalized Difference Vegetation Index (NDVI) in high-biomass areas, such as dense forests [34]. However, NDVI was used alongside other indices to provide a more comprehensive view of vegetation structure in areas with low biomass.

Landsat scenes were overlaid on NFI cluster locations (see illustration in Figure 3) to align image pixel with ground data. The aim was to extract pixels corresponding to each entirely forested cluster. However, caution is needed when overlaying image and extracting pixel data. First, the dimensions of clusters’ sub-units and image pixels differ: Brazilian NFI cluster sub-units are 20 m wide, meaning they may either fit into a single 30 m pixel or overlap two neighboring pixels, which can cause spectral mixture [6]. Second, Landsat-8 OLI bands have a geodetic accuracy of up to 65 m circular error (about two pixels) at the 90% confidence level [43]. Third, GNSS (Global Navigation Satellite System) devices introduce position errors (in our case, <15 m) that can be larger under canopy cover, so cluster locations in the field may not align precisely with digital maps.

To address these issues, we used AGB per cluster (hereafter ‘plot AGB’) rather than AGB per sub-unit, assuming that an entire cluster better represents mean AGB than its individual sub-units. We also assumed that pixels within a cluster’s 210 m × 210 m window (Figure 3) capture spectral signatures that accurately represent an entirely forested cluster.

Following this procedure, plot AGB in Mg ha⁻¹ was associated with the spectral data of pixels within the cluster’s window. Studies methods, as used by [6,33], are generally regarded as suitable for integrating field clusters (or structures with disjoint sub-plots) with remotely sensed data [44]. The general pixel-based AGB model is expressed in Equation (1):

$${\widehat{B}}_p=\mathrm{f} (\mathit{X}; \widehat{\mathit{\theta }})$$

(1)

${\widehat{B}}_p$: aboveground biomass predicted for a pixel, $\mathit{X}$: surface reflectance of OLI bands and vegetation indices (see Table 3), and $\widehat{\mathit{\theta }}$: model parameter.

The random forest (RF) algorithm was used to (i) construct the pixel-level biomass model, (ii) select the most significant predictors (surface reflectance from OLI and vegetation indices), and (iii) predict pixel-level AGB in Mg ha⁻¹. This nonparametric technique is robust to overfitting and generates multiple decision trees that serve as classifiers, each voting for the most popular class of the output variable [45]. The ‘randomForest’ package [46] in R environment [47] was used to implement the RF approach and then map the AGB (Mg ha⁻¹).

2.3. Uncertainty Analysis

After predicting the AGB at the tree, plot, and pixel levels, uncertainties in AGB predictions were also quantified at these levels. The Monte Carlo method and bootstrap technique were used to assess these sources of uncertainty. Subsequently, the forest-mapping-related uncertainty was analyzed and incorporated into the error propagation from tree to pixel. The R environment [47] was used for the entire uncertainty analysis. Although we acknowledge the importance of all six sources of uncertainty enumerated in the Introduction, data limitations restricted us to assessing only four of them.

2.3.1. Error Propagation from Tree to Plot

As the tree AGB was predicted using the biomass model of [23], the variability of parameters from the models was considered as the main source of uncertainty in the tree AGB predictions. The uncertainty in tree-level AGB predictions, caused by parameter variability, was propagated to the plot level. This propagation was necessary because the plot AGB per unit area served as input for calibrating the pixel-level AGB model (Equation (1)). The error propagation followed these steps:

• Stochastically simulate 5000 values for $\widehat{\beta }$ using the Monte Carlo method, assuming $\widehat{\beta } ~ N\left(\beta , {\sigma }_{\widehat{\beta }}^2\right)$, Note: The variances of $\widehat{\beta }$, ${\mathsf{\sigma }}_{\widehat{\mathsf{\beta }}}^2$’s, were not reported in [23] and were obtained through personal communication.
• Randomly select a sample of $\widehat{\beta }\prime s$ without replacement from the coefficients simulated in the previous step to predict ${\widehat{B}}_{trunk,i}$, ${\widehat{B}}_{branches,i}$, ${\widehat{B}}_{foliage,i}$ and then ${A\widehat{G}B}_{tree,i}$ for all inventoried trees.
• Replicate step (ii) 1000 times.

For every k-th replication generated in steps (i)–(iii), the AGB in the j-th cluster, ${A\widehat{G}B}_j^k$, was predicted as in Equation (2):

$${A\widehat{G}B}_j^k={\displaystyle \raisebox{1ex}{$ \sum _{i=1}^{n_j}{A\widehat{G}B}_{tree,ij}^k$}\!\left/ \!\raisebox{-1ex}{$1000·s$}\right.}$$

(2)

${A\widehat{G}B}_{tree,ij}^k$: AGB of the i-th tree in the j-th cluster for the k-th replication; s = 0.4, i.e., cluster size in ha, and the division by 1000 converts kg to Mg. Therefore, Equation (2) gives AGB in Mg ha⁻¹.

2.3.2. Error Propagation from Plot to Pixel

Uncertainty in the plot-level predictions derived from Equation (2) was propagated to the pixel level. Since Equation (1) is calibrated using the AGB of pixels within a 49 (7 × 7) pixel window (Figure 3), for the k-th replication, the pixel-level AGB predicted through Equation (1) was defined as ${A\widehat{G}B}_p^k={A\widehat{G}B}_j^k$. The uncertainty in the pixel-level predictions was assessed as follows:

• For the j-th forested cluster, randomly select λ pixels from the cluster’s window (Figure 3), and then pair ${A\widehat{G}B}_p^k$ and pixel information (Table 3), where λ = 7, 8, …, 49.
• Replicate step (i) 1000 times (this is $n_{rep})$.
• Using the RF with bootstrap approach, fit Equation (1) using the paired data and predict ${\widehat{B}}_p^k$.
• Let $Q={A\overline{G}B}_p^k={\displaystyle \frac{1}{n_c\times \mathsf{\lambda }}}{\sum }_{j=1}^{n_c}\sum _{p=1}^{\mathsf{\lambda }}{AGB}_p^k$ and $\widehat{Q}={A\overline{\widehat{G}}B}_p^k={\displaystyle \frac{1}{n_c \times \mathsf{\lambda }}}{\sum }_{j=1}^{n_c}\sum _{p=1}^{\mathsf{\lambda }}A{\widehat{G}B}_p^k$. Then, estimate the variance of $(Q-\widehat{Q})$, i.e., $U$.

A formulation for estimating the total variance of ($Q-\widehat{Q}$) is presented in [48], here demonstrated in Equation (6). To calculate this, let

$$\overline{Q}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$n_{rep}$}\right.}{\sum }_{k=1}^{n_{rep}}{\widehat{Q}}_k$$

(3)

be the average of $Q$ over the replications,

$$\overline{U}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$n_{rep}$}\right.}{\sum }_{k=1}^{n_{rep}}{U}_k$$

(4)

be the average of the variances of ($Q-\widehat{Q}$) over the replications, and

$$B={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$(n_{rep}-1)$}\right.}{\sum }_{k=1}^{n_{rep}}{\left({\widehat{Q}}_k-\overline{Q}\right)}^2$$

(5)

be the among-replications variance. The total variance of $(Q-\overline{Q})$ is finally given by

$$T=\overline{U}+\left(1+{\displaystyle \frac{1}{n_{rep}}}\right)·B$$

(6)

Let ${SE}_{AGB}=\sqrt{T}$ be the standard error of AGB in Mg ha⁻¹; then, the general sample-based uncertainty of AGB in Mg ha⁻¹, with error propagated from tree to pixel, is

$$t_v\left(\alpha /2\right)·{SE}_{AGB}$$

(7)

and finally, the relative uncertainty is

$${\delta }_{\overline{Q}}={\displaystyle \frac{t_v\left(\alpha /2\right)·{SE}_{AGB}}{Q}}$$

(8)

where $v=\left(n_{rep}-1\right)\times {\left(1+{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.}\right)}^2$, $r=\left(1+{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$n_{rep}$}\right.}\right){\displaystyle \raisebox{1ex}{$B$}\!\left/ \!\raisebox{-1ex}{$\overline{U}$}\right.}$, and $(1-\alpha )$ is the confidence level.

Due to the funnel-shaped distribution of AGB residuals, which indicates residual heteroscedasticity, Equations (7) and (8) were applied by classes as follows. First, we estimated $Q_M$, i.e., the mean AGB over the $M$ pixels across the study area (we found $Q_M$ = 110.28 Mg ha⁻¹). Second, we estimated $Q_{M,h}$, i.e., the mean AGB over the $M_h$ pixels across the h-th county (with $Q_{M,h}$ ranging from 88.3 to 165.0 Mg ha⁻¹). Based on this range, we defined the following classes: (i) 85–94.9 Mg ha⁻¹, (ii) 95–104.9 Mg ha⁻¹, and so on, up to (viii) 155–165 Mg ha⁻¹, resulting in eight class-specific sampling uncertainties. By substituting $Q$ in Equation (8) with $Q_M$ and $Q_{M,h}$, we obtained the relative sampling uncertainties of AGB for the entire study area and for the h-th county, respectively, as follows:

$${\delta }_{Q_M}={\displaystyle \frac{t_v\left(\alpha /2\right)\times {SE}_{AGB}}{Q_M}}$$

(9)

and

$${\delta }_{Q_{M,h}}={\displaystyle \frac{t_v\left(\alpha /2\right)·{SE}_{AGB}}{Q_{M,h}}}$$

(10)

The estimated CIs are

$$\mathrm{CI} \mathrm{for} \mathrm{AGB} \mathrm{in} \mathrm{Mg} \mathrm{ha}{-}^1 \mathrm{in} \mathrm{the} \mathrm{study} \mathrm{area} {CI}_{Q_{AGB}}=Q_M\left(1\pm {\delta }_{Q_M}\right)$$

(11)

$$\mathrm{CI} \mathrm{for} \mathrm{AGB} \mathrm{in} \mathrm{Mg} \mathrm{ha}{-}^1 \mathrm{in} \mathrm{counties} {CI}_{Q_{AGB,h}}=Q_{M,h}\left(1\pm {\delta }_{Q_{M,h}}\right)$$

(12)

2.3.3. Forest-Mapping-Related Uncertainty

In this phase, we developed a procedure to incorporate the forest-cover-mapping-related uncertainty into the error propagation from tree to pixel. The uncertainty in forest cover predictions arises from the supervised image classification, which inherently carries a degree of uncertainty (e.g., due to calibration samples and sample size). This mapping-related uncertainty is estimated using the variance of a model-assisted estimator, specifically a regression estimator [49,50]. A regression estimator typically includes a correction for systematic sampling errors added to the mean population unit estimate [49]. In this study, the regression estimator was calibrated with values derived from the forest/non-forest classification, and the correction term was taken from the classification error matrix. The variance of this estimator was subsequently used to estimate CIs for the mean proportion forest, as described in [51,52,53].

For simplicity, we divided this process into two sub-topics: large-area approach, where the proportion forest and CIs are estimated using simple random sampling, and the small-area approach, where general estimators are adapted for use at the stratum (county) level. The model-assisted estimator provides the uncertainty of the estimate for the mean proportion forest.

• Large-area approach

According to [51], the synthetic estimator of the mean proportion forest $\left({\overline{\widehat{y}}}_{syn}\right)$ is the mean of the forest/non-forest classification for all pixels, where the mean takes values between 0 and 1. This estimator also represents the relationship between the $T\widehat{F}A$ and the total area, as follows:

$${\overline{\widehat{y}}}_{syn}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$M$}\right.}\sum _{p=1}^M{\widehat{y}}_p={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$A$}\right.}T\widehat{F}A{\widehat{y}}_p=\left\{\begin{array}{c}1,p-\mathrm{t}\mathrm{h} \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l} \mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{d} \mathrm{a}\mathrm{s} \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\\ 0,p-\mathrm{t}\mathrm{h} \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l} \mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{d} \mathrm{a}\mathrm{s} \mathrm{n}\mathrm{o}\mathrm{n}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\end{array}\right.$$

(13)

${\widehat{y}}_p$: predicted class for the p-th pixel.

The mean systematic classification error in ${\overline{\widehat{y}}}_{syn}$ is expressed by the following:

$$B\widehat{ia}s\left({\overline{\widehat{y}}}_{syn}\right)={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}{\sum }_{p=1}^m{\epsilon }_p={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}{\sum }_{p=1}^m\left({\widehat{y}}_p-y_p\right)$$

(14)

$y_p$: reference class for the p-th pixel.

The difference between ${\overline{\widehat{y}}}_{syn}$ and $B\widehat{ia}s\left({\overline{\widehat{y}}}_{syn}\right)$ is the regression estimator [51] for proportion forest (${\overline{\widehat{y}}}_{reg}$), as follows:

$${\overline{\widehat{y}}}_{reg}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$M$}\right.}{\sum }_{p=1}^M{\widehat{y}}_p-{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}{\sum }_{p=1}^m\left({\widehat{y}}_p-y_p\right)$$

(15)

The variance of ${\overline{\widehat{y}}}_{reg}$ is

$$V\widehat{a}r\left({\overline{\widehat{y}}}_{reg}\right)={\displaystyle \frac{(1-f)}{m(m-1)}}{\sum }_{p=1}^m{({\epsilon }_p-{\overline{\epsilon }}_p)}^2$$

(16)

where ${\overline{\epsilon }}_p=B\widehat{ia}s\left({\overline{\widehat{y}}}_{syn}\right)$; $f=m/M$.

• Small-area approach

If Equation (15) is an unbiased estimator at the large area level, it may not be unbiased at the small area. In this context, a small area refers to a county, where stratified sampling [54] is proposed to address this issue. We treated counties as strata and assigned weights (${\widehat{W}}_h$) to each h-th county. Additionally, studies such as [51,55,56] have demonstrated that stratified estimators are more precise than, for example, simple expansion estimators, as follows:

$${\widehat{W}}_h={\displaystyle \frac{{\sum }_{p=1}^{M_h}{\widehat{y}}_{p,h}}{\sum _{p=1}^M{\widehat{y}}_p}}={\displaystyle \frac{{C\widehat{F}A}_h}{T\widehat{F}A}} {\widehat{y}}_{p,h}=\left\{\begin{array}{c}\begin{array}{c}1,p-\mathrm{t}\mathrm{h} \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l} \mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{d} \mathrm{a}\mathrm{s} \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t} \mathrm{i}\mathrm{n} \mathrm{t}\mathrm{h}\mathrm{e} h-\mathrm{t}\mathrm{h} \mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{y}\\ 0,p-\mathrm{t}\mathrm{h} \mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l} \mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{d} \mathrm{a}\mathrm{s} \mathrm{n}\mathrm{o}\mathrm{n}-\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t} \mathrm{i}\mathrm{n} \mathrm{t}\mathrm{h}\mathrm{e} h-\mathrm{t}\mathrm{h} \mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{y}\end{array}\end{array}\right.$$

(17)

In stratified-sampling designs, ${\overline{\widehat{y}}}_{syn}$ must be replaced by the stratified proportion forest, ${\overline{\widehat{y}}}_{st}$ (st for stratified). The ${\overline{\widehat{y}}}_{st}$ is an unbiased estimator of the proportion forest and is generally not the same as ${\overline{\widehat{y}}}_{syn}$ [54]. It is given by

$$\begin{gathered} {\overline{\widehat{y}}}_{st}={\sum }_{h=1}^L\left[{\widehat{W}}_h·{\overline{\widehat{y}}}_h\right] \\ \mathrm{being} {\overline{\widehat{y}}}_h={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$M_h$}\right.}{\sum }_{p=1}^M{\widehat{y}}_{p,h}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$A_h$}\right.}{C\widehat{F}A}_h \end{gathered}$$

(18)

${\overline{\widehat{y}}}_h$: estimate of the proportion forest of the h-th county.

The estimate of the stratified bias, ${Bias(\overline{\widehat{y}}}_{st})$, can be obtained as in Equation (19):

$$\begin{gathered} {\overline{\widehat{y}}}_{st}={\sum }_{h=1}^L\left[{\widehat{W}}_h·{\overline{\widehat{y}}}_h\right] \\ \mathrm{being} {Bias(\overline{\widehat{y}}}_{h^{*}})={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m_{h^{*}}$}\right.}{\sum }_{p=1}^{m_{h^{*}}}{\epsilon }_{p,h^{*}}={\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m_{h^{*}}$}\right.}{\sum }_{p=1}^{m_{h^{*}}}\left({\widehat{y}}_{p,h^{*}}-y_{p,h^{*}}\right) \end{gathered}$$

(19)

${Bias(\overline{\widehat{y}}}_{h^{*}})$: estimate of the bias of the h^*-th grand stratum.

Note that ${Bias(\overline{\widehat{y}}}_{h^{\mathrm{*}}})$ is distinguished by (*), indicating that it was estimated per grand stratum. This distinction is made in cases where $m_h$ is small (≤20) or unavailable for a given stratum, in which [54] suggests combining two or more strata to obtain estimates. Since we had, on average, approximately 11.3 (=850/75) sample units per stratum (i.e., per county), we aggregated counties into grand strata.

The regression estimator for proportion forest $\left({\overline{\widehat{y}}}_{reg,st}\right)$ is as follows:

$$\begin{gathered} {\overline{\widehat{y}}}_{reg,st}={\overline{\widehat{y}}}_{st}-B\widehat{ia}s\left({\overline{\widehat{y}}}_{st}\right) \\ ={\sum }_{h=1}^L\left[{\widehat{W}}_h·{\overline{\widehat{y}}}_h\right]-\sum _{h=1}^L\left[{\widehat{W}}_h·B\widehat{ia}s\left({\overline{\widehat{y}}}_{h^{*}}\right)\right] \\ ={\sum }_{h=1}^L\left\{{\widehat{W}}_h·\left[{\overline{\widehat{y}}}_h-B\widehat{ia}s\left({\overline{\widehat{y}}}_{h^{*}}\right)\right]\right\} \end{gathered}$$

(20)

The within-grand-stratum variance of the estimate ${\overline{\widehat{y}}}_h$ is

$$V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)={\displaystyle \frac{(1-f_{h^{*}})}{m_{h^{*}}(m_{h^{*}}-1)}}{\sum }_{p=1}^{m_{h^{*}}}{({\epsilon }_{p,h^{*}}-{\overline{\epsilon }}_{p,h^{*}})}^2$$

(21)

where $f_{h^{*}}=m_{h^{*}}/M_{h^{*}}$ represents the sampling fraction that can be ignored when is ≤0.02.

Like the estimator ${Bias(\overline{\widehat{y}}}_{h^{*}})$, the estimator of the variance of ${\overline{\widehat{y}}}_h$, $V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)$, was also estimated per grand stratum rather than per stratum. As long as $m_h$ is large in all strata, ${Bias(\overline{\widehat{y}}}_h)$ and $V\widehat{a}r\left({\overline{\widehat{y}}}_h\right)$ should be estimated per stratum.

The variance of ${\overline{\widehat{y}}}_{reg,st}$ in the stratified sampling for proportion forest is

$$V\widehat{a}r\left({\overline{\widehat{y}}}_{reg,st}\right)={\sum }_{h=1}^L\left[{\widehat{W}}_h^2 ·V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)\right]$$

(22)

Note that Equation (22) can be calculated even with h and h^* in the summation, with the peculiarity that all strata within a grand stratum will share the same variance. The CI for the proportion forest by stratum (county) and study area are, respectively,

$$\begin{gathered} \hspace{1em}\hspace{1em}\hspace{1em}{CI}_{{\widehat{y}}_h}={\overline{\widehat{y}}}_{reg,h}\pm t_{m_{h^{*}}}(\alpha /2)·\sqrt{V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\overline{\widehat{y}}}_{reg,h}\left(1\pm {\displaystyle \frac{t_{m_{h^{*}}}(\alpha /2)·\sqrt{V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)}}{{\overline{\widehat{y}}}_{reg,h}}}\right) \\ ={\overline{\widehat{y}}}_{reg,h}\left(1\pm {\delta }_{{\overline{\widehat{y}}}_h}\right) \end{gathered}$$

(23)

and

$$\begin{gathered} \hspace{1em}\hspace{1em}\hspace{1em}{CI}_{\widehat{y}}={\overline{\widehat{y}}}_{reg,st}\pm t_m(\alpha /2)·\sqrt{V\widehat{a}r\left({\overline{\widehat{y}}}_{reg,st}\right)} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\overline{\widehat{y}}}_{reg,st}\left(1\pm {\displaystyle \frac{t_m(\alpha /2)·\sqrt{V\widehat{a}r\left({\overline{\widehat{y}}}_{reg,st}\right)}}{{\overline{\widehat{y}}}_{reg,st}}}\right) \\ ={\overline{\widehat{y}}}_{reg,st}\left(1\pm {\delta }_{{\overline{\widehat{y}}}_{st}}\right) \end{gathered}$$

(24)

where (1 − α) is the confidence level; ${\overline{\widehat{y}}}_{reg,h}={\overline{\widehat{y}}}_h-B\widehat{ia}s\left({\overline{\widehat{y}}}_h\right)$.

The CIs for the $CFA$ (in ha) and for the $TFA$ (in ha) can be directly estimated by multiplying the lower and upper limits of the Cis by the total area of a county ($A_h$) and the total study area ($A$), respectively, as follows:

$${CI}_{CFA}=A_h\left[{\overline{\widehat{y}}}_{reg,h}\left(1\pm {\delta }_{{\overline{\widehat{y}}}_h}\right)\right]$$

(25)

$${CI}_{TFA}=A\left[{\overline{\widehat{y}}}_{reg,st}\left(1\pm {\delta }_{{\overline{\widehat{y}}}_{st}}\right)\right]$$

(26)

2.3.4. Uncertainty Analysis with Error Propagated from Pixel to Forest Mapping

To estimate the CI for the total AGB in Mg, considering both the AGB CIs and the forest area CIs, a multiplicative approach is used. This approach is known as “uncertainty of a product” as described by [57]. The CIs for the total AGB are derived by multiplying the CIs for the mean AGB in Mg ha⁻¹ (from Equations (11)–(12)) with the CIs for the forest area in hectares (from Equations (25) and (26)).

LL of the total AGB in counties is

$${CI}_{AGB,h}^{-}=A_h\left[Q_{M,h}·{\overline{\widehat{y}}}_{reg,h}\left(1-{\delta }_{Q_{M,h}}·{\delta }_{{\overline{\widehat{y}}}_h}\right)\right]$$

(27a)

UL of the total AGB in counties is

$${CI}_{AGB,h}^{+}=A_h\left[Q_{M,h}·{\overline{\widehat{y}}}_{reg,h}\left(1+{\delta }_{Q_{M,h}}·{\delta }_{{\overline{\widehat{y}}}_h}\right)\right]$$

(27b)

LL of the total AGB in the study area is

$${CI}_{AGB}^{-}=A\left[Q_M·{\overline{\widehat{y}}}_{reg,st}\left(1-{\delta }_{Q_M}·{\delta }_{{\overline{\widehat{y}}}_{st}}\right)\right]$$

(28a)

and UL of the total AGB in the study area is

$${CI}_{AGB}^{+}=A\left[Q_M·{\overline{\widehat{y}}}_{reg,st}\left(1+{\delta }_{Q_M}·{\delta }_{{\overline{\widehat{y}}}_{st}}\right)\right]$$

(28b)

LL: lower limit; UL: upper limit.

The final analysis aimed to quantify the impact of ignoring the forest-mapping-related uncertainty. To achieve this, we compared the CIs estimated using the method that incorporates forest-mapping-related uncertainty (as in Equations (27) and (28)) with those estimated assuming no classification uncertainty (as in Equations (11) and (12)). The difference between these two approaches for estimating the CIs was calculated using Equations (29) and (30):

$${Diff}_{LL}\left(\%\right)=100 {\displaystyle \frac{y_{LL1}-y_{LL0}}{y_{LL1}}}\hspace{1em}\mathrm{Difference} \mathrm{for} \mathrm{the} \mathrm{lower} \mathrm{limit}$$

(29)

$${Diff}_{UL}\left(\%\right)=100 {\displaystyle \frac{y_{UL1}-y_{UL0}}{y_{UL1}}}\hspace{1em}\mathrm{Difference} \mathrm{for} \mathrm{the} \mathrm{upper} \mathrm{limit}$$

(30)

$y_{LL1}$ and $y_{UL1}$ represent the lower and upper limits, respectively, that account for the mapping-related uncertainty. $y_{LL0}$ and $y_{UL0}$ represent the lower and upper limits, respectively, that exclude the mapping-related uncertainty.

2.4. Variance Estimators

A variety of estimators of the mean proportion forest and the variance of the mean were compared to the regression estimator, as outlined below.

2.4.1. Stratified Random Sampling for Proportions (SRSP) Estimator

The study of [55] shows how stratified random sampling for proportions (SRSP) can be implemented to estimate the mean and CIs for the proportion forest. These authors’ Equation (9) is equivalent to Equation (5.52) in [54] and provides the estimate for the mean proportion forest across the entire population. This model was used to estimate ${\overline{\widehat{y}}}_{h^{*}}$, i.e., mean proportion forest per grand stratum. The stratified proportion forest, ${\overline{\widehat{y}}}_{st}$, was estimated as in Equation (18). It is important to note that Equation (10) in [55] gives the SRSP estimator for the variance of the proportion forest, which corresponds to Equation (5.57) in [54] when considering the sample variance estimate. The SRSP estimator requires three inputs: (1) $m_{rc}$ (their $n_{ik}$), the sample count at cell (r, c) in the error matrix; (2) $m_{c.}$ (their $n_{i.}$), the sum of sample counts of the map class c in the error matrix; and (3) $W_c$ (their $W_i$), the proportion of area mapped as class c. One should note that $W_c$ differs from ${\widehat{W}}_h$ in our Equation (17). The first difference is that index c refers to the map class, while index h refers to the stratum. Second, $W_i$ indicates the relative area of a class among the mapped classes, whereas ${\widehat{W}}_h$ denotes the proportion of forest area in a stratum relative to the $T\widehat{F}A$. Since our study involves a forest/non-forest map, index c refers to these two classes. Therefore, Equation (10) in [55] provides a weighted average of the variance across the forest/non-forest classes. To obtain a stratified variance suitable for our study, Equation (10) in [55] was estimated per grand stratum, providing within-grand-stratum variances for ${\overline{\widehat{y}}}_h$, $V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)$. The variance of ${\overline{\widehat{y}}}_{st}$, $V\widehat{a}r\left({\overline{\widehat{y}}}_{st}\right)$, was estimated according to Equation (5.7) in [54], where ${\widehat{W}}_h$ was used. The CI for the proportion forest was estimated as in [55], and CIs for the final estimates of AGB (in Mg) were obtained similarly to those in Equations (27) and (28).

2.4.2. Synthetic Estimator

The difference between the synthetic and regression estimators is that, in the first, the proportion forest per grand stratum, ${\overline{\widehat{y}}}_h$, and the stratified proportion forest, ${\overline{\widehat{y}}}_{st}$, were estimated in their synthetic form shown in Equation (18). For both methods (synthetic and regression), estimators of variance of ${\overline{\widehat{y}}}_h$, $V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)$, and variance of ${\overline{\widehat{y}}}_{st}$, $V\widehat{a}r\left({\overline{\widehat{y}}}_{st}\right)$, were estimated as in Equations (21) and (22), respectively. Although $V\widehat{a}r\left({\overline{\widehat{y}}}_{st}\right)$ is the same for both approaches, the CI for proportion forest differs between them, due to the replacement of ${\overline{\widehat{y}}}_{reg,st}$ by ${\overline{\widehat{y}}}_{st}$. Thus, the magnitude of $B\widehat{ia}s\left({\overline{\widehat{y}}}_{st}\right)$ is the unique factor that can differ the CIs estimated from the variance of the synthetic and regression estimators. The idea is to estimate such magnitude. The CI for proportion forest was estimated as in Equation (24), substituting ${\overline{\widehat{y}}}_{reg,st}$ by ${\overline{\widehat{y}}}_{st}$. CIs for the final estimates of AGB (in Mg) were obtained similarly as in Equations (27) and (28).

2.4.3. Ratio Estimator

Ratio estimator in stratified random sampling (see pp. 164–166 of [54]) was the final approach we compared. This method is sample-based and does not rely on results from the error matrix, as in [56], or from image classification, as in [55]. In this case, the NFI clusters serve as the sample units. The area sampled in a grand stratum, in ha, is determined by $x_{h^{*}}=0.4\times n_{c,h^{*}}$, and the mean among clusters is ${\overline{x}}_{h^{*}}=0.4$, where $n_{c,h^{*}}$ is the number of 4000 m² clusters installed in the h^*-th grand stratum. The total forest area sampled in a grand stratum, in ha, is $y_{h^{*}}=0.4\times n_{c,h^{*}}\times {\overline{p}}_{c,h}$, and the mean among clusters is ${\overline{y}}_{h^{*}}=0.4\times {\overline{p}}_{c,h}$, where ${\overline{p}}_c$ is average proportion forest over the clusters in a grand stratum. The distribution of the clusters across grand strata is as follows: grand stratum I, $n_{c,h^{*}=I}=34$; grand stratum II, $n_{c,h^{*}=II}=71$; and grand stratum III, $n_{c,h^{*}=III}=47$.

For both simple random and systematic samplings, the ratio estimator $\widehat{R}={\displaystyle \raisebox{1ex}{$\overline{y}$}\!\left/ \!\raisebox{-1ex}{$\overline{x}$}\right.}$ is valid for the grand stratum h^* [54], which results in ${\widehat{R}}_{h^{*}}={\displaystyle \raisebox{1ex}{${\overline{y}}_{h^{*}}$}\!\left/ \!\raisebox{-1ex}{${\overline{x}}_{h^{*}}$}\right.}$. The same authors also present an estimate derived from a single combined ratio, ${\widehat{R}}_c$ (c for combined), for stratified sampling. Equation (31) serves as an unbiased estimator of the proportion forest for the large area.

$${\widehat{R}}_c=\sum _{h=1}^L{\displaystyle \frac{\raisebox{1ex}{${N_{h^{*}}\times \overline{y}}_{h^{*}}$}\!\left/ \!\raisebox{-1ex}{$N$}\right.}{\raisebox{1ex}{$N_{h^{*}}\times {\overline{x}}_{h^{*}}$}\!\left/ \!\raisebox{-1ex}{$N$}\right.}}=\sum _{h=1}^L{W}_{h^{*}}{\displaystyle \frac{{\overline{y}}_{h^{*}}}{{\overline{x}}_{h^{*}}}}$$

(31)

where the grand-stratum population sizes are as follows: grand stratum I, $N_{h^{*}=I}=\mathrm{3,075,216}$; grand stratum II, $N_{h^{*}=II}=\mathrm{6,560,065}$; and grand stratum III, $N_{h^{*}=III}=\mathrm{4,197,368}$.

To gain an initial understanding of ${\widehat{R}}_{h^{*}}$ and ${\widehat{R}}_c$, the data comprise 19 entirely forested clusters, 77 partially forested clusters, and 56 non-forest clusters. One remark is that because our ratio estimator is based on plot sampling, $W_{h^{*}}$ takes a different form compared to the other methods. Specifically,$W_{h^{*}}=N_{h^{*}}/N$ is more familiar than the form shown in Equation (17), though both terms are statistically correct and ensures that ${\sum _{h^{*}}W}_{h^{*}}=1$.

To mitigate statistical complexities, we assumed that sample units (i.e., clusters) are randomly located within the grand stratum; as such, the within-grand-stratum variance of ${\widehat{R}}_{h^{*}}$ (Equation (23) of [54]) is an approximation of $V\widehat{a}r \left({\widehat{R}}_{h^{*}}\right)$, which may lead to an overestimate of the variance.

$$V\widehat{a}r\left({\widehat{R}}_{h^{*}}\right)\approx {\displaystyle \frac{1-f_{h^{*}}}{n_{c,h^{*}}\times {\overline{x}}_{h^{*}}^2}}\left[{\displaystyle \frac{{\sum }_{j=1}^{n_{c,h^{*}}}{\left(y_{h^{*}}-{\widehat{R}}_{h^{*}}{x}_{h^{*}}\right)}^2}{n_{c,h^{*}}-1}}\right]$$

(32)

where $f_{h^{*}}$ is the sampling fraction ($n_{h^{*}}/N_{h^{*}}$) and can be ignored when $f_{h^{*}}$ ≤ 0.02.

For stratified random sampling, the variance of ${\widehat{R}}_c$ can be approximated analogously to applications in [54].

$$V\widehat{a}r\left({\widehat{R}}_c\right)\approx \sum _{h=1}^L\left[{\widehat{W}}_{h^{*}}^2 ·V\widehat{a}r\left({\widehat{R}}_{h^{*}}\right)\right]$$

(33)

Based on the method we developed for the ratio estimators, ${\widehat{R}}_{h^{*}}$ and ${\widehat{R}}_c$ are essentially proportion forest estimators. Therefore, we will hereafter refer to both ratios as ${\overline{\widehat{y}}}_h$ and ${\overline{\widehat{y}}}_{st}$, respectively. Similarly, the variance of ${\widehat{R}}_{h^{*}}$, $V\widehat{a}r\left({\widehat{R}}_{h^{*}}\right)$, and the variance of ${\widehat{R}}_c$, $V\widehat{a}r\left({\widehat{R}}_c\right)$, will be referred to as $V\widehat{a}r\left({\overline{\widehat{y}}}_{h^{*}}\right)$ and $V\widehat{a}r\left({\overline{\widehat{y}}}_{st}\right)$, respectively. The CI for the proportion forest was obtained as in Equation (24), substituting $m$ by $n_{{c,h}^{*}}$ and ${\overline{\widehat{y}}}_{reg.st}$ by ${\overline{\widehat{y}}}_{st}$. CIs for the total AGB (in Mg) were estimated similarly to Equations (27) and (28).

3. Results

3.1. Forest Cover Map

Overall, user’s and producer’s accuracies are shown in

Table 4. Error matrix of the image classification for the grand strata.

j = Columns: Reference		Non-Forest	Forest	$\mathbf{Total} \left({\mathit{m}}_{\mathit{i}+}\right)$	User’s Accuracy	${\mathit{W}}_{\mathit{i}}$
Grand stratum * ‘Mid-west’
i = rows: classification	Non-forest	0.85	0.01	0.86	0.99	0.867
	Forest	0.03	0.11	0.14	0.78	0.133
	Total	0.88	0.12	1		1.000
	Producer’s accuracy	0.96	0.94		$OA=0.96$
Grand stratum * ‘Mid-south’
i = rows: classification	Non-forest	0.59	0.01	0.60	0.98	0.689
	Forest	0.10	0.30	0.40	0.75	0.311
	Total	0.69	0.31	1		1.000
	Producer’s accuracy	0.86	0.97		$OA=0.89$
Grand stratum * ‘Southeast’
i = rows: classification	Non-forest	0.39	0.06	0.45	0.87	0.466
	Forest	0.05	0.50	0.55	0.91	0.534
	Total	0.44	0.56	1		1.000
	Producer’s accuracy	0.89	0.89		$OA=0.89$

* Grand stratum can be visualized in Figure 1 as a greater region.

. This table shows that the forest/non-forest classification achieved an overall accuracy of 0.90. The forest cover map resulting from the image classification is displayed in Figure 5. The estimated total forest area ($T\widehat{F}A$) is 1,762,569 ha.

3.2. Mapping and Confidence Interval of Biomass

The (RF-based) pixel-level AGB model (Equation (1)) produced AGB predictions with an RMSE of 12.26 Mg ha⁻¹ (=11.87%). Figure 4 shows the relationship-observed vs. -predicted AGBs. The output data used to fit Equation (1) are AGB (Mg ha⁻¹) with errors propagated from trees to plot and pixel levels. Among the input variables in the RF training, EVI (Table 3) was the most important.

The gap shown in Figure 4 is due to the absence of clusters with observed AGB values between 170 and 280 Mg ha⁻¹, similar to the AGB variations found in [3] for southern Brazil. The fitted RF-based model was then used to predict pixel-level AGB. The funnel-shaped relationship between observed vs. predicted AGBs (Figure 4) indicates residual heteroscedasticity, a common issue in tree variables such as volume and biomass [58,59]. To address potential model residual issues, we estimated the AGB sampling errors (Mg ha⁻¹) by classes of estimated mean AGB. These sampling errors ranged from 19.6 to 27.6 Mg ha⁻¹ and are presented in

Table 5. Standard errors of estimates of AGB by class of mean AGB.

Class of Mean AGB	Standard Error of AGB (Mg ha−1)
85–94.9 Mg ha−1	20.23
95–114.9 Mg ha−1	23.03
105–114.9 Mg ha−1	27.62
115–124.9 Mg ha−1	20.56
125–134.9 Mg ha−1	21.48
135–144.9 Mg ha−1	23.78
145–154.9 Mg ha−1	22.74
155–165.0 Mg ha−1	19.69

Standard errors originated from $\sqrt{T}$ (Equation (6)).

. Figure 5 displays the spatial distribution and AGB stocks across the study area.

The mean AGB predicted at pixel level allowed us to determine the class of mean AGB (Table 5) to which each county and the entire study area belong. We then used the class-specific standard errors presented in Table 5 to estimate the CIs for the total AGB in each county and the study area. These CIs were estimated with the forest-mapping-related uncertainty incorporated into the error propagation from trees to plot and pixel levels.

Table 6. Estimated forest area and confidence intervals for the total AGB (in Mg) stocked in the 10 most biomass-stocked counties and in the study area.

Ranking of Counties by AGB Stock	$\mathit{C}\widehat{\mathit{F}}\mathit{A}$ (ha)	${\overline{\widehat{\mathit{y}}}}_{\mathit{r}\mathit{e}\mathit{g},\mathit{h}}$	Confidence Interval (95%)
			Lower Limit	Total	Upper Limit
1st-Prudentópolis	89,818	0.414	6,559,415	9,841,894	13,644,139
2nd-Guarapuava	109,535	0.257	5,433,926	8,548,769	12,286,599
3rd-Cruz Machado	82,189	0.570	5,830,320	8,220,944	10,897,854
4th-São Mateus do Sul	57,061	0.438	5,845,709	7,857,820	10,138,573
5 th-Bituruna	77,361	0.651	5,537,339	7,728,680	10,155,447
6th-General Carneiro	77,424	0.739	5,467,624	7,385,800	9,486,184
7th-Coronel Domingos Soares	77,099	0.409	4,572,440	6,753,636	9,240,165
8th-Pinhão	77,406	0.300	4,182,012	6,117,693	8,380,487
9th-Teixeira Soares	32,879	0.376	3,867,473	5,124,560	6,554,368
10th-São João do Triunfo	33,305	0.475	3,949,145	5,107,761	6,404,205
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
Study area	$T\widehat{F}A$=	${\overline{\widehat{y}}}_{reg,st}$ = 0.357	152,007,465	218,128,166	291,927,114
	1,762,569

$C\widehat{F}A$: county’s forest area. $T\widehat{F}A$: total forest area. ${\overline{\widehat{y}}}_{reg,h}$: estimate of proportion forest in the h-th county. ${\overline{\widehat{y}}}_{reg,st}$: estimate of proportion forest in the large area.

shows these CIs for the 10 counties with the highest biomass stocks, as well as for the entire study area, alongside estimates of $C\widehat{F}A\prime s$, $T\widehat{F}A$, ${\overline{\widehat{y}}}_h$, and ${\overline{\widehat{y}}}_{reg,st}$. Table S1 (Supplementary file) gives estimates for the remaining counties.

An important finding of this study is the effect of ignoring the uncertainty due to the forest cover classification. The differences in the lower and upper limits of the total AGB, with and without adding this uncertainty, were graphed (Figure 6) as a function of the counties’ estimated proportion forest. At the county level, the percentage differences ranged from ±5% to ±63%. The x-axis shows the counties’ estimated proportion forest, ranging from 4% to 75%. The differences at the large-area level are summarized at the bottom of Table 8.

3.3. Comparing Variance Estimators

Table 7. Estimates of mean, bias, and standard error of proportion forest calculated through different estimators.

Origin of the Estimator	Grand Stratum a	${\widehat{\mathit{W}}}_{\mathit{h}}$	${\overline{\widehat{\mathit{y}}}}_{\mathit{h}}$	$\mathit{B}\widehat{\mathit{i}\mathit{a}}\mathit{s}\left({\overline{\widehat{\mathit{y}}}}_{\mathit{h}}\right)$c	$\sqrt{\mathit{V}\widehat{\mathit{a}}\mathit{r}\left({\overline{\widehat{\mathit{y}}}}_{\mathit{h}}\right)}$
Regression estimator	Mid-west	0.093	0.147	0.025	0.012
	Mid-south	0.463	0.342	0.091	0.018
	Southeast	0.443	0.511	−0.008	0.022
	(Large area)	(1.00)	0.399 b	0.041	0.013
SRSP estimator	Mid-west	0.093	0.187	-	0.010
	Mid-south	0.463	0.279	-	0.013
	Southeast	0.443	0.436	-	0.021
	(Large area)	(1.00)	0.340 b	-	0.011
Synthetic estimator	Mid-west	0.093	0.147	-	0.012
	Mid-south	0.463	0.342	-	0.018
	Southeast	0.443	0.511	-	0.022
	(Large area)	(1.00)	0.399 b	-	0.013
Ratio estimator	Mid-west	0.222	0.099	-	0.031
	Mid-south	0.474	0.399	-	0.045
	Southeast	0.303	0.483	-	0.054
	(Large area)	(1.00)	0.358 b	-	0.028

^a: Grand stratum can be visualized in Figure 1 as a greater region. ^b: estimate of ${\overline{\widehat{y}}}_{st}$. ^c: bias term to subtract from the estimate of the proportion forest, ${\overline{\widehat{y}}}_h$.

presents all estimates of the mean, bias term, and variance of proportion forest estimated using the tested estimators. In this table, ${\overline{\widehat{y}}}_h$ represents the synthetic form, and the regression estimator is derived by subtracting $B\widehat{ia}s\left({\overline{\widehat{y}}}_h\right)$ from ${\overline{\widehat{y}}}_h$. Three main observations can be made from Table 7. First, the estimates of proportion forest for the large area showed minimal variation across the different estimators. Second, the ratio estimator had a variance that was twice as large as the others, although it accurately estimated mean proportion forest. Third, related to the ratio estimator, which is plot sample-based, ${\widehat{W}}_h$ was significantly overestimated in the mid-west, the most deforested grand stratum (see $W_i$ in Table 4), in comparison to the other three estimators.

The estimates presented in Table 7 were used to estimate the CI for proportion forest, allowing us to estimate the CI for forest cover as in Equations (12) and (13). After also estimating the CIs for Q of AGB (Equations (14) and (15)), the forest-mapping-related uncertainty was incorporated (Equations (16)–(19)) into the model-related uncertainties represented by the CIs for Q of AGB.

Table 8. CIs for AGB (in Mg) stocked in the study area calculated through different estimators and differences with and without adding the forest-mapping-related uncertainty.

Origin of the Estimator	Confidence Interval (95%)		Range
	Lower Limit	Upper Limit
Adding forest-mapping-related uncertainty (A)
Regression estimator	152,007,465	291,927,114	139,919,649
SRSP estimator	146,227,162	277,613,198	131,386,036
Synthetic estimator	170,618,677	322,975,900	152,357,223
Ratio estimator	138,387,641	315,158,106	176,770,465
Without adding forest-mapping-related uncertainty (B)
-	182,108,175	303,808,153	121,699,978
Difference a between (A) and (B)
Regression estimator	−19.8%	−4.1%	-

^a: difference calculated as in Equations (22) and (23).

provides the final lower and upper limits of AGB stored across the large area, with results provided for the four estimators. The effect of ignoring the mapping-related uncertainty is also provided at the bottom of Table 8.

4. Discussion

4.1. Sensitivity, Generalization, and Application of Estimators

The first issue with the regression estimator is its mathematical sensitivity to the distribution and sufficiency of the sample [49] used to estimate components of the error matrix (Table 2). Sample units must effectively capture the mean accuracy achieved in the forest/non-forest classification for the model-assisted estimator of proportion forest $\left({\overline{\widehat{y}}}_{reg,st}\right)$ to remain unbiased and for its variance $\left(V\widehat{a}r\left({\overline{\widehat{y}}}_{reg,st}\right)\right)$ to be minimized [49], thereby producing more reliable CIs for AGB. This is especially crucial at the small area level, where the sample size is limited. In such cases, a limited reference data sample can lead to unreliable estimates for the variance of ${\overline{\widehat{y}}}_h$ and the bias term $\left(Bias\left({\overline{\widehat{y}}}_{h^{*}}\right)\right)$. This can result in misleading conclusions about the CI for strata with small sample sizes. As suggested, small strata should be aggregated, as reported in [54], which also mentions that little reduction is gained beyond six strata, although using more is reasonable for geographic subdivisions such as counties or administrative regions. A second issue is that the regression estimator remains unaffected by the sensor’s spatial resolution. Research by [6] observed that the spatial resolution of the employed sensor is not a significant source of uncertainty in forest maps at small scales, like counties. They reached this conclusion by comparing fine (30) and coarse (1000 m) spatial resolutions and found similar accuracy in county-level AGB estimates.

While our analytical procedure assigns counties as small areas, the framework outlined in this paper can be broadly applied to assess uncertainty for any small area and for various forest attributes. An ideal application of our approach would be when estimates are needed at geographic subdivision levels, such as counties or administrative regions. Users may opt to apply our full procedure for error propagation in the modeling process and incorporate forest-mapping uncertainty using their own method. In addition, airborne laser scanning (ALS) data are frequently used to map forest attributes across both large and small areas [21,60,61,62,63,64,65]. In such cases, ALS data can be combined with optical images (e.g., satellite images) or used independently to model and map forest attributes. The procedure for incorporating the forest-mapping-related uncertainty into the error propagation when biomass from ALS data are predicted can be applied accordingly.

4.2. Limitation and Comparison with Existent Methods

Two primary limitations are noted. First, we assumed that the map units (pixels) are randomly selected for accuracy assessment across and within small areas. However, remote-sensing researchers often systematically select sample units for accuracy purposes. In such cases, $V\widehat{a}r\left({\overline{\widehat{y}}}_h\right)$ in Equation (10) becomes a biased estimator of the variance of ${\overline{\widehat{y}}}_h$ [55]; therefore, it should be adjusted accordingly [54]. Second, the estimators for the mean and variance of the mean were developed for proportional allocation within stratified sampling [54]. In practice, these estimators perform optimally when the strata, as administrative regions, states, or counties, are known in advance. Knowing the size of each stratum allows users to proportionally allocate reference data for map accuracy assessment. On the other hand, users who classify maps into categories and then stratify their collected reference data into those categories should use a post-stratified estimator, as outlined in [51].

While comparing the variance of the regression estimator with that of the other three estimators used to assess forest-mapping-related uncertainty, we noted that overall precision is numerically reflected in the CI width shown in Table 8. The CI width increases as a consequence of rising variance in ${\overline{\widehat{y}}}_{reg,st}$ and $B\widehat{ia}s\left({\overline{\widehat{y}}}_{st}\right)$. Basically, this means that the more accurate and precise the classification, the smaller the variance and, consequently, the narrower the CI width, as shown in Table 8. In this regard, Table 8 reveals that all estimators produced wider ranges compared to the method that excluded the forest-mapping-related uncertainty, as expected when an additional source of uncertainty is incorporated into the error propagation. Among the estimators, the regression and SRSP produced the closest CIs and the narrowest ranges, with the SRSP estimator offering the highest precision (smallest variance).

One potential reason for the superior performance of the SRSP estimator is its use of forest/non-forest map classes as strata, effectively employing map classes as weights ($W_i$) in the estimator. This advantage has also been reported in [56]. Other relevant findings emerged from the comparison between estimates with and without incorporating the forest-mapping-related uncertainty. For the regression estimator, the lower limit of the CI was shifted leftward as expected; however, the upper limit also moved leftward. This shift occurs because the CI limits adjust in response to the bias introduced by image classification, which is accounted for in Equation (20). When forest cover is overestimated (positive bias), the model-assisted estimator of the variance causes both CI limits to move leftward relative to the scenario without forest-mapping-related uncertainty, and vice versa. This bias correction is achieved through a bias adjustment term, which shifts the CI position accordingly. Comparing the regression and synthetic estimators in Table 8 highlights this effect, as the synthetic estimator lacks this bias adjustment. Thus, the model-assisted estimator offers an advantage by producing CIs that are robust in terms of both precision and accuracy, while other estimators typically maintain only precision.

The synthetic and regression estimators exhibited equal variances due to their similarity in calculation. However, the magnitude of the estimated bias term (Table 7) was substantial enough to reduce the CI width (Table 8) obtained from the regression estimator by approximately 13 million Mg compared to the synthetic form. Additionally, the estimated bias term shifted the CI limits roughly 20 million Mg to the right relative to the synthetic estimator. This resulted in a lower total population estimate using the regression estimator. Lastly, the ratio estimator closely approximated the mean proportion forest estimated by the regression estimator. However, it did not produce variance estimates with the same accuracy as the other methods. Considering that the ratio estimator is plot-sample-based and does not rely on remotely sensed data or image processing, we could even conclude: “Not bad!”

5. Conclusions

Ignoring the forest-mapping-related uncertainty erroneously narrows the CI width as the estimate of proportion forest area decreases. At the small-area level, incorporating this uncertainty can widen the CI for total biomass by up to 63% in extreme cases, while at the large-area level, the CI expands by 5–7%. One recommendation is therefore to consider mapping-related uncertainty particularly in studies that involve image classification of forest/non-forest with low accuracy and in research conducted in areas with highly fragmented forests. Studies with low classification accuracy are prone to produce misleading conclusions about the CI for forest population parameters, leading to biased CIs influenced by the discrepancy between the user’s and producer’s accuracy, and by the overall classification accuracy. In areas with fragmented forests, the smaller the forest remnants, the more unreliable the CI for total AGB becomes. This issue is especially critical when forest cover is below 50%.