Interpreting the uncertainty of model-based and design-based estimation in downscaling estimates from NFI data: a case-study in Extremadura (Spain)

ABSTRACT Remotely sensed data are increasingly used together with National Forest Inventory (NFI) data to improve the spatial precision of forest variable estimates. In this study, we combined data from the 4th Spanish National Forest Inventory (SNFI-4) and from the 2nd nationwide Airborne Laser Scanning (ALS) survey to develop predictive forest inventory variables (total over bark volume (V), basal area (G), and annual increase in total volume (IAVC)) and aboveground biomass (AGB) models for the eight major forest strata in the region of Extremadura that are included in the Spanish Forest Map (SFM). We generated maps at 25 m resolution by applying an area‐based approach (ABA) and 758 sample plots measured with good positional accuracy within the SNFI-4 in Extremadura (Spain). Inventory performance is mainly influenced by spatial scale and vegetation structure. Therefore, in this study, we conducted a comparative analysis of statistical inference methods that can characterize forest inventory variables and AGB uncertainty across multiple spatial scales and types of vegetation structure. Predictions at pixel level were used to produce county, provincial, and regional model-based estimates, which were then compared with design-based estimates at different scales for different types of forest. We developed and tested both methods for forested area (cover, 19,744.15 km2), one province (9126.78 km2), and two counties (1594.42 km2 and 2076.76 km2, respectively) in Extremadura. The resulting relative standard error (SE) for regional level forest type-specific model-based estimates of V, G, IAVC, and AGB ranged from 3.34%–14.46%, 3.22%–12.50%, 4.46%–16.67%, and 3.63%–12.58%, respectively. The performance of the model-based approach, as assessed by the relative SE, was similar to that of the design-based approach at regional and provincial levels. However, the precision of SNFI model-based estimates was higher than that of estimates based on only the plot observations in small areas (e.g. at county level). The standard errors (SE) for model-based inferences were stable across the different scales, while SNFI design-based errors were higher due to the small sample sizes available for small areas. The findings indicate that SNFI-model based maps could be used directly to estimate forest inventory variables and AGB in the major forest strata included in the Spanish Forest Map, leading to potentially large economic savings.


Introduction
Large-scale forest resource assessments based on National Forest Inventory (NFI) data are essential for supporting national and regional forest planning strategies, especially in the context of climate change (Ene et al. 2013). NFI sampling schemes aim to collect field-based information, from which it may be possible to obtain reliable estimates at different planning scales (Tomppo et al. 2008). At smaller scales, the need for forest information requires downscaling of estimates obtained from NFI data, under the assumption that adaptations are required regarding the data and the inference approaches used in order to produce reliable estimates. Two main methodological approaches can be used for this purpose: (1) a design-based approach in which only field based estimates at the estimation scale are used and (2) a model-based approach in which spatially explicit information is used together with field data in order to increase the efficiency of statistical inference and to enable generation of spatial outcomes (Persson and Ståhl 2020;Lister et al. 2020). These two approaches are very different in terms of the inferential framework and SE estimation. In design-based methodologies, the SE is estimated on the basis of the variability in the study variables in the sampling plots and the sample size in the estimation domain (Kangas 2006). By contrast, in model-based methodologies, SE estimation depends on the performance of a model fitted with data (regarding the variable of interest) estimated from auxiliary variable(s) measured over the study area and the associated variability over the estimation domain (McRoberts 2006). These methodological differences lead to difficulty in interpreting the estimates when working at different scales, even for the same data source (NFI).
Remotely sensed data have been widely used in recent years for mapping and estimating forest inventory variables and biomass (White et al. 2016;Chirici et al. 2020;Saarela et al. 2020). Among remotesensing based techniques, active remote sensing is widely used nowadays due to its capacity to precisely describe the three-dimensional (3D) structure of forests (Wulder et al. 2012;McRoberts, Andersen, and Naesset 2014). The development of forestry applications has mainly focused on the use of airborne laser scanning (ALS) and the implementation of area-based methods (ABA) to estimate forest attributes (McRoberts, Andersen, and Naesset 2014). In ABA, the area under study is discretised in grid cells and sampled with field plots whose locations are chosen on the basis of probability sampling. One common application involves developing models that are applied over the entire forest area of interest (Guerra-Hernández, Tomé, and González-Ferreiro 2016;Guerra-Hernández et al. 2019) and the subsequent use of model-based or model-assisted inference approaches to produce estimates at the desired scale. The increasing availability of remotely sensed data has recently challenged the traditional method of performing forest inventories and has led to interest in these inference techniques (Ståhl et al. 2011;McRoberts, Naesset, and Gobakken 2013;Ståhl et al. 2016;Saarela et al. 2018). Like traditional designbased inference (Saarela et al. 2015), model-based inference provides regional estimates of total and mean values and also large scale mapping of forest characteristics. Methods used to produce extensive maps and to calculate small-area estimates of forest inventory variables are becoming essential in most sophisticated NFIs (Breidenbach and Astrup 2012;Mauro et al. 2017;Breidenbach et al. 2018;Hill, Mandallaz, and Langshausen 2018;Chirici et al. 2020;Vega et al. 2021), as it can be difficult to obtain precise direct estimators of forest attributes within relatively small areas with a few sample plots.
Maps representing the distribution of a variable of interest across a study area are a valuable outcome of forest monitoring systems. Some specific applications require maps with very low uncertainty, which requires great effort from a methodological point of view. However, in most studies, characterization of the uncertainty of mapped or estimated forest inventory variables has depended on either sample-based root mean square error or bias, both of which provide average measures of uncertainty for plot-related predictions and are therefore limited to uncertainty analysis of mean estimates of forest inventory variables for the whole study area (Chen et al. 2016). Thus, it is important to calculate the magnitude of mapping errors, in order to determine the reliability, usage, and interpretation of the map, because map users and producers are interested in the quality of maps (Persson and Ståhl 2020). On the other hand, while error analysis for NFI-based maps has mainly been conducted at stand scale (e.g. (Nilsson et al. 2017;Rahlf et al. 2021)), the errors in NFI-model-based ALS map estimates has not been studied relative to design-based estimates at lower spatial scales. If the stratum-level SE of NFI-model based maps is stable at different scales and forest type, it could easily be used in smaller areas, such as county level, in which the design-based NFI approach will perform less well due to the small sample size.
In this study, the above-mentioned methodologies were used in order to combine NFI and ALS data to predict forest inventory variables in the region of Extremadura, a Mediterranean forest area in Spain (Guerra-Hernández et al. 2019;Pascual et al. 2020;. The Spanish National Forest Inventory (SNFI) requires upgrading of ground positioning protocols to improve data records. In the case of Extremadura, new true positions for samples from the 4th SNFI are enabling the use of remote sensing data for spatial estimation of forest attributes (Guerra-Hernández et al. 2019;Pascual et al. 2020;. To the best of our knowledge, country-level experiences in the Mediterranean forest type have not yet been reported in regard to improving the positional accuracy of SNFI plots and spatially and temporally coincident discrete-return ALS data for a variety of types of forests. Hence, estimation of Mediterranean forest population parameters and rigorous calculation of the associated uncertainties are essential to produce consistent estimates at different scales from SNFI data. Second, the SE associated with design-based inventory and model-based inventory have never been compared in Mediterranean countries at different scales and for different types of forest. Furthermore, this work aimed to extend knowledge about the performance of the SE as a function of the spatial scale and different forest types using NFI-model based maps and estimates from a design-based SNFI approach. The main objectives of the current study were as follows: (i) to construct large scale estimates of forest inventory variables for a large test area (i.e. 27,300 km 2 ) in the region of Extremadura (Spain) by combining SNFI plot data and remotely sensed ALS data as auxiliary information for mapping growing stock volume (V, m 3 ha −1 ), basal area (G, m 2 ha −1 ), annual increase in total volume (IAVC, m 3 ha −1 year −1 ), and aboveground biomass (AGB, Mg ha −1 ) at fine spatial resolution (25 × 25 m); and (ii) to compare the performance of model-based and design-based approaches at regional, provincial and county scales in terms of relative SE for the eight major forest strata included in SNFI plots.

Study area
The study area is located in the region of Extremadura (Spain), which has a total forested area of 19,744.15 km 2 (Figure 1). The Spanish Forest Map of Spain (SFM) (E: 1:25,000) (MAPA 2018) and the sampling design of the Spanish National Forest Inventory (SNFI) were used to cover as much of the whole spectrum of forest structure in the study area as possible (MAGRAMA 2017). The most recent versions of the SFM are based on photointerpretation of aerial photography and digitization of polygons that are classified in relation to the vegetation present in the area. We selected the following eight Mediterranean forest types: Dehesas: sparse old-growth oak forests (Quercus spp.) with low tree density; Encinares: sparse oak forest (Quercus ilex subsp. ballota (Desf.) Samp); Pinaster: even-aged forests of Pinus pinaster subsp. mesogensis Aiton; Quejigares: Mediterranean deciduous oak forests dominated by Quercus faginea and Quercus pyrenaica; Eucalyptus: commercial forests of Eucalyptus spp., dominated by Eucalyptus camaldulensis Dehn., used to produce pulpwood; Alcornocales: dominated by cork oak (Quercus suber L.), from which cork is obtained; Pinea: even-aged forests of Pinus pinea L. managed for cone production; and Riparian forest: composed of a mixed mature forest of riverine species. These eight ecosystems cover an area of 18,335.41 km 2 in Extremadura, i.e. 92.91% of the whole forest area.

NFI data
The design-based inference method used data from 1808 concentric circular plots with 25 m of radius included in the SNFI-4 (Table 1). In the model-based inference approach, 758 of these plots, corresponding to well-georeferenced sample plots were used for model fitting. Field data were collected during 2017 according to the SNFI protocol (Álvarez-González et al. 2014). Forest inventory variables were computed at tree-level using equations from the NFI protocol and were then expanded to units per surface area, by aggregation of the values to plot level. In the case of AGB, the species-specific tree-level allometric equations used by the SNFI were used (Montero, Ruiz-Peinado, and Munoz 2005;Ruiz-Peinado, González, and Del Rio 2012;Ruiz-Peinado, del Rio, and Montero 2011). For more details about the SNFI-4 field data and processing, see (Álvarez-González et al. 2014) and (Alberdi et al. 2017).
The following stand variables were calculated from the tree variable measurements by using tree expansion factors: growing stock volume (V, m 3 ha −1 ), basal area (G, m 2 ha −1 ), annual increase in total volume (IAVC, m 3 ha −1 year −1 ), and aboveground biomass (AGB, Mg ha −1 ). The variables estimated for the 2129 sampled plots used for design-based inference are listed in Table 1. During field measurement of SNFI-4, the data recording uncertainty was addressed by using commercial grade global navigation satellite systems (GNSS) to upgrade positioning information on new samples established in forest areas under expansion and on the set of previously existing SNFI-4 sample plots. The coordinates were measured with a handheld data collection system (TRIMBLE Juno 5B handheld, Trimble Inc. USA) (error range 1-2 m of positioning error after post-processing).

Airborne laser scanning acquisition and processing
Data from two ALS point clouds were used in the study: the first ALS dataset was collected between October 2018 and March 2019 and covers northern Table 1. Summary of ground data collected in the 4 th National Forest Inventory (NFI) for the eight forest types considered. Plot-level estimates are presented for growing stock volume (V, m 3 ha −1 ), stand basal area (G, m 2 ha −1 ), annual increase in total volume (IAVC, m 3 ha −1 year −1 ), and aboveground biomass (AGB, Mg ha −1 ). Extremadura, whereas the second was collected between October 2018 and July 2019 and covers southern Extremadura. The datasets correspond to the second round of countrywide ALS measurements, which are publicly available in Spain through the National Plan for Aerial Orography (hereafter referred to as PNOA). Square ALS blocks of 2 km side, covering the whole region of Extremadura were obtained from the CNIG (Centro Nacional de Información Geográfica) computer server (http://centrodedescargas.cnig.es/ CentroDescargas/index.jsp) to cover the province of Badajoz. The scanning sensors used to collect the ALS data in the study area were a RIEGL LMS-Q1560 for Extremadura-North (EXT-N) and a LEICA ALS80 for Extremadura-South (EXT-S). The nominal laser pulse density varied between 2 points m −2 for EXT-N and 1 point m −2 for EXT-S. The methodology used to filter, process, and generate extensive maps at regional level for each forest strata is summarized in Figure 2. The ALS data were first processed using LAStools software (Isenburg 2020), and the resulting set of 23 ALS statistics (Table 2) were computed for each sample plot. A detailed description of the software parametrization and the workflow involved in processing the ALS point cloud is given by (Pascual et al. 2020). In the final step, wall to wall maps were produced at the SFM forest strata level. For this purpose, the lasclip command was used to clip the normalize point cloud for each forest stratum in the Spanish Forest Map (SFM) (E: 1:25,000) (MAPA 2018). The lascanopy tool was then used to generate the raster metrics at the level of forest strata. Finally, V, G, IAVC, and AGB ALSbased models were applied at 25 m pixel level using the Raster package in R software (Hijmans et al. 2015).

Modeling forest inventory variables
Power-function models were used to establish the relationships between forest inventory variables derived from field measurements and ALS variables. The general expression for the models is as follows: where y i is the objective forest inventory variable for the plot i; X 1;i ; X 2;i : . . . ; �X p;i are the potential explanatory ALS-derived variables belonging to the metrics of height distributions or measurements related to canopy density (Table 2) for plot i; β 0 , β 0 , . . ., β p are the model parameters to be estimated by non-linear regression analysis; and ε i is a zero-mean additive random residual. The first step used in the modeling phase was to select the optimal set of predictor variables to use in estimating the forest inventory variables. The R package "leaps," (Lumley and Miller 2017) was used to select the significant predictors of the regression. Two predictor variables were used to estimate the parameter from the models. Collinearity between regressors was prevented by checking the variance inflation factor (VIF), and regressors with VIF above 10 were disregarded (Belsley, Kuh, and Welsch 2005). The models were fitted using the nls function implemented in the R package BASE (R Core Team 2020). Finally, the model efficiency (Mef, Eq 2), the overall root mean square error (RMSE, Eq. 3), and the relative root mean square error (rRMSE, Eq. 4) were computed to determine the accuracy of ALSderived models for estimating forest inventory variables. The adjusted R 2 was not considered a suitable goodness-of-fit statistic, as the modeling was conducted in non-linear space. We therefore used the Mef statistic, which provides a simple index of performance on a relative scale, where a value of one indicates a "perfect" fit, a zero value reveals that the model is no better than a simple average, and negative values indicate a poor model (Vanclay and Peter Skovsgaard 1997). In addition, leave-one-out cross-validation was performed for each potential regression model using R routines (R Core Team 2020).
RMSE ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where n is the number of plots; yi is the fieldestimated variable of interest for the plots i; � y is the mean value for the field-estimated variable of interest over the n plots; ŷ i is the estimated value of variable of interest derived from the non-linear regression model; and p is the number of parameters in the model.

Model based estimation
For estimation of species-specific V, G, IAVC, and mean AGB, model-based inferential methods that do not require probability sampling were used ( Table 1, No. of plots model-based inference). The available well-georeferenced SNFI samples were considered to have been selected by non-probability, purposive sampling. The model-based unbiased estimator of the species-specific V, G, IAVC, and AGB mean was based on McRoberts (2006): where μ i is the predicted value of the variables of interest V, G, IAVC, or AGB from ALS-based model (Eq. 1) for the ith population unit (pixel, map unit), and N is the total number of population units for each forest type for which ALS-based V, G, IAVC, and AGB models were applied. The corresponding estimator of the standard error (SE) was computed as follows: b SEμ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where z ij = @f X i ; β ð Þ=@β j are the elements of vector Z i , β j are the parameters of the model (Eq. 1), and σ 2β n o is an estimator of the parameter variance- Percentage of all returns above 2.00 m over total all returns covariance matrix as described in McRoberts, Naesset, and Gobakken (2013). The elements σ 2β n o states for the partial derivative of the model in the i position (number of parameter) regarding the j parameter. As Eq. (6) was often computationally intense, 10% of the total number of population units (pixels) were used for the computation through random sample selection. This is consistent with the method used by McRoberts, Naesset, and Gobakken (2013), who demonstrated that using a systematic subsample of the population did not have adverse effects on the variance. In addition, absolute error and relative error were computed as follows: where t is the critical value of the t-student distribution with N-1 degrees of freedom for a 95% confidence level.

Design-based estimation
The SNFI-4 in Extremadura established permanent plots that were systematically distributed at the intersections of a 1 km x 1 km grid in areas identified as forest land in the SNFM. The total samples for design-based inference method at different scales were: 1808 plots for regional level, 802 plots for provincial level, 268 for Siberia county level, and 166 plots for Villuercas county level. The design-based estimator of the species-specific V, G, IAVC, and AGB mean was based on the following: where y i is the field-estimated variable of interest for the plot i, and n is the number of plots. The following is an estimator of the variance of � Y (omitting finite population correction): where The corresponding estimator of the standard error was computed as follows: In addition, absolute and relative errors are computed as indicated in the model-based approach (Eq. 7, Eq. 8). Although the estimators shown above are notunbiased when applied to systematic samples ((Gregoire and Valentine Harry 2008), p 55), they are widely used and are conservative in the sense of producing greater variance in estimating the population mean (Steen et al. 2020). For calculation of design-based estimates, data were included for all plots with centers in the stratum regardless of whether the plots included part of the area outside of the corresponding stratum.

Regional ALS-based models using SNFI plots
The models selected for each variable of interest and forest type are shown in Appendix A, Table A.1, and scatter plots of observed against predicted values for variables with ALS-based models are shown in Figure 3. The performance of the models selected on the basis of density and height metrics differed between the forest type and the variable of interest (Figure 3, Appendix A Table A.1). The Mef values for the resulting regional forest type-specific ALS-based models for predicting V, G, IAVC, and AGB ranged from 0.29-0.89, 0.30-0.80, 0.13-0.74, and 0.27-0.86, respectively (Table A.1, Figure 3). The variation, in terms of rRMSE, for more even-aged Pinea forest to Quercus dominated forest type, ranged from 28.80%-63.12% for V, 29.21%-50.64% for G, 41.03%-94.84% for IAVC, and from 27.22% −51.63% for AGB. In general, regional IAVC-based models were the least accurate models in terms of Mef and rRMSE than V, G, and AGB models at the level of each forest type, except for V in Alcornocales which was slightly worse in terms of Mef.
In terms of Mef, we also found the same pattern in Riparian forest where IAVC-based model was better in terms of Mef than G and AGB models. For the variables modeled using cross-validation, as expected V, G, IAVC, and AGB always yielded the worst goodnessof-fit statistics for all forest types, and the associated conclusions were similar.

Wall-to-wall mapping at regional level
The strata from SFM and the models from Appendix A (Table A.1) were applied to metrics derived from ALS to construct a regular 25 m resolution V, G, IAVC, and AGB map at the regional scale (Extremadura) ( Figure  4). Regarding the V and AGB for the different species, forest production was higher in the central northern counties (Urdes and Gata) and western counties (Siberia and Villuercas) and lower in the interior areas in the province of Badajoz.

Comparison of model-based and design-based inference at regional, provincial and county scale
The results for model-based and design-based estimation of forest inventory variables in terms of relative SE (%) are shown in Figure 5, whereas the overall results in terms of estimate, standard error (SE), absolute error relative (ε), and relative SE (%) -expressed on a per hectare basis by forest type at regional, provincial, and county scale -are shown in Appendix B. Regarding the estimated values, the model-based approach generally yielded lower values than the design-based approach. In terms of errors, the regional model-based results are similar to those obtained for the corresponding designbased estimates from SNFI plots (MITECO 2020) (Appendix B, Table B1). The ALS-based models of IAVC that performed less well for some formations (i.e. Encinares; Mef = 0.16 and Quejigares; 0.26) yielded higher relative values of SE (%) than the design-based methodology at regional, province, and county scale ( Figure 5(b, d)). By contrast, the estimation stabilized, with errors below 10% in homogeneous coniferous forest formations ( Figure 5(c, g)) and based on more accurate ALS-based models for all forest inventory variables. This relative values of SE (%) was also stable for Dehesas formations, where ALS-based models performed less well but the models were fitted with a higher representative sample size than the other Quercus-dominated forests (Encinares, Quejigares, and Alcornocales). The relative SE increased gradually from regional to county scale for the estimates from the design-based approach, as we expected ( Figure 5, Table B3 and B4). This effect was more pronounced for less representative formations across provincial and county scales, especially for Quejigares, Riparian forest ( Figure 5(d, h)), and P. pinea forest in the county of Villuercas (Figure 5g). On the other hand, it is important to note that estimates are not provided by official design-based SNFI aggregated statistics at provincial and county administrative levels (MITECO 2020). The official design-based SNFI provides aggregated statistics at regional level.

Discussion
The main of objectives of this study were (i) to demonstrate that it is possible to produce extensive spatial estimates of the most important forest , Alcornocal (f), P. pinea (g), and Riparian (h) at regional, provincial, and county scales.
inventory variables, even in large Mediterranean forests, by using field data and model-based inference based on predictors obtained using auxiliary ALS variables and (ii) to understand the performance of the model-based and design-based approaches using SNFI and ALS data at different spatial scales.
Regarding the first objective, the performance of new regional forest type-specific ALS-based models for predicting V, G, IAVC, and AGB from SNFI data confirmed the limitation of ALS-based models in differentiating AGB and characterizing vegetation structure under sparse canopy cover conditions as in Dehesas and Encinares (Mef forest inventory models ranging from 0.22-0.29 and 0.16-0.61 in Dehesas and Encinares, respectively). In the case of Quejigares, the performance of ALS-based models, as indicated by Mef (range 0.26-0.54), was also not satisfactory, probably because (i) the data were acquired under leaf-off conditions during the winter in northern Extremadura (Sophie, Donoghue, and Galiatsatos 2020) and (ii) the stand structure of Q. pyrenaica forest is characterized by complex multi-layered, uneven-aged stands. According to the results obtained, the set of models strengthened the idea that the combination of height and vertical canopy structure metrics provides a sufficient and concise quantitative description of a homogeneous vertical structure in Mediterranean areas, as in coniferous (P. pinea and pinaster formations) and Eucalyptus spp forests. In this study, much better goodness-of-fit statistics were obtained than in previous studies carried out with the same type of data in Spain for the main productive forest species (Spanish NFI and laser pulse density of 0.5 points/m 2 ) at the regional level in La Rioja (Pinus sylvestris L.) (Esteban et al. 2020;Fernández-Landa et al. 2018 (Table 3). In general, the performance of the models, as measured by Mef and rRMSE, was slightly better than in previous studies, for three main reasons: (i) more precise geolocation of the SNFI plots in our study, (ii) SNFI coincide temporarily with the ALS data, and (iii) higher density of the ALS data in the present study (1-2 pulses m 2 ) than in Novo-Fernández et al. (2019) and Fernández-Landa et al. (2018) (0.5 points/m 2 ). The Mef associated with ALSbased models (Eucalyptus spp., P. pinaster, and P. pinea) for predicting IAVC, V, and AGB ranged from 0.63 to 0.74, from 0.72 to 0.89, and from 0.76 to 0.86, respectively. The values obtained in the present study were higher than those reported by Novo-Fernández et al. (2019), with Mef for the IAVC, V, and AGB models for the three forest species (E. globulus, P. pinaster, and P. radiata) ranging from 0.51 to 0.69, from 0.71 to 0.83, and from 0.71 to 0.82, respectively. In terms of rRMSE, the values for IAVC, V, and AGB ranged from 41.03 to 46.95%, 28.80% to 46.5%, and from 27.22% to 43.64%, respectively. The accuracy of the ALS-derived forest inventory models was slightly better (except for Eucalyptus spp.) than that obtained by Novo-Fernández et al. (2019)   Regarding the second objective, several studies have reported large area forest inventory variables and AGB estimates using models and remotely sensed data (Chen et al. 2016;Esteban et al. 2020;Chirici et al. 2020). At regional scale, our results in terms of relative SE using model-based inference (ranged from 3.34-14.46%) are in concordance with those reported by (Esteban et al. 2020) who reported relative SE ranging from 2% to 11% using Random forest (RF)-SNFI models to calculate uncertainties in V estimates for the main forest species of La Rioja (Spain). (Esteban et al. 2020) also found model-based stabilities were less for P. halepensis and P. nigra formations for which field plot sample sizes were smaller. Our values of relative SE are slightly higher than (Chirici et al. 2020) who estimated V using NFI field plots and remotely sensed data (Landsat and SAR data) for a large area in Italy. Their SE estimates varied from 2% to 4%, although SE estimates were not reported for specific forest species. It is well known that uncertainty in the estimation of different vegetation strata mapsaffects the determination of the volume stocks, but in practice its effect has been mostly ignored (Margolis et al. 2015;Neigh et al. 2013;Yanjun et al. 2016). Previous studies assessing this source of error have shown very significant accuracy losses for forest stock variables at different scales (Xue et al. 2017;Esteban etal. 2020;Räty et al. 2020;Breidenbach et al. 2021). In our study we have not considered this source of uncertainty, but we assume that it will cause a loss of accuracy for all alternatives studied. It remains for future studies to assess whether the loss of accuracy would be of different importance for the inference methodologies considered in this study.
Several conclusions can be drawn from the second part of this study. First, the errors associated with modelbased estimates by forest types were comparable to those associated with the design-based estimates at regional scale ( Figure 5). At provincial and county levels, the SNFI model-based inference method produced more accurate estimates than the SNFI design-based estimates in small areas. The study findings demonstrated the limitation of SNFI design-based for smaller regions and especially for less representative forest formations area that may contain few plots. Second, the results confirmed that the downscaling potential of model-based approach depends on the statistical strength between field inventory variable and ALS auxiliary variables and also on the structural complexity of the forests (Chen et al. 2016;Vega et al. 2021). The relative SE for the model-based method was stable across different spatial scales, reaching a maximum value of 33.51% for IAVC in Quejigares forest at county level. For the other formations and forest attributes, the estimation stabilized, with errors generally below 25%. Third, model-based estimators benefit from sample sizes for larger areas when used to smaller subsets of the larger areas as long as ALS-based models were fitted with a representative sample size as Dehesas formations which field plot sample sizes were higher than the others Quercus-dominated forests in the Region (McRoberts, Naesset, and Gobakken 2013). The results indicate that greater stability in the sample size forest species formations produces less uncertainty in V, G, IAVC, and AGB estimates from model-based. Table 3. Summary of studies at regional level using countrywide ALS-PNOA coverage and SNFI plots in Spain using a point density of 0.5 points/m 2 . Growing stock volume (V, m 3 ha −1 ), stand basal area (G, m 2 ha −1 ), annual increase in total volume (IAVC, m 3 ha −1 year −1 ), and aboveground biomass (AGB, Mg ha −1 ). In brackets, Independent error assessment using SNFI plots. * not SNFI plots.
In the case of Extremadura, admissible relative SE (%) for sampling design of V (species with vocation for timber production, i.e. Pinaster and Eucalyptus spp) and G (species for protection, recreational, silvopastoral (i.e. Dehesas and Encinares) objectives) are 15% and 30%, respectively. These values from future forest management plans legislation are referred at stand scale (named as "cuarteles") by grouping forest managements units with the same objective. It is important highlight that model-based relative SE (%) for Eucalyptus spp., P. pinaster, and P. pinea was lower than 15% of admissible relative SE (%) of V in Extremadura and lower than 30% of G for Dehesas and Encinares at county level. The results obtained from this work open the door for the use of SNFImodel based maps at stand scale and forest type and to calculate small-area estimates of forest inventory variables. The use of the SNFI-model based map would thus generate considerable economic savings by eliminating or reducing the need for field plot measurement. Based on the estimates of accuracy reported in the tables in Appendix B, the practical response is that the accuracy of the model-based estimates is greater than that of the traditional SNFI design-based estimates at county level.

Conclusions
The potential regression models exhibited no serious lack of fit to the ALS-forest inventory and biomass models, except in regard to characterizing vegetation structure under sparse canopy cover. The comparison demonstrated that design-based SNFI plot errors increased for smaller areas while model-based relative standard errors were stable across the different scales, as long as ALS-based models were appropriate and fitted with a representative sample size. The proposed modelbased approach may be useful in cases where model-assisted inference cannot be applied either because models have been developed from an independent dataset or matching of the samples for model development and application cannot be fully achieved, as in the case of Extremadura where not all SNFI plots have well-georeferenced positions. Finally, the model-based approach could be used operationally to generate stable forest inventory predictions at different scales in the region of Extremadura. 13

Data availability
The data that support the findings of this study are available from the corresponding author, Jurado-Varela A., upon reasonable request. Maps are avaliable to download (http://sitex. gobex.es/SITEX/centrodescargas/viewsubcategoria/56)