Large scale mapping of forest attributes using heterogeneous sets of airborne laser scanning and National Forest Inventory data

Background: The Norwegian forest resource map (SR16) maps forest attributes by combining national forest inventory (NFI), airborne laser scanning (ALS) and other remotely sensed data. While the ALS data were acquired over a time interval of 10 years using various sensors and settings, the NFI data are continuously collected. Aims of this study were to analyze the effects of stratification on models linking remotely sensed and field data, and assess the accuracy overall and at the ALS project level. Materials and methods: The model dataset consisted of 9203 NFI field plots and data from 367 ALS projects, covering 17 Mha and 2/3 of the productive forest in Norway. Mixed-effects regression models were used to account for differences among ALS projects. Two types of stratification were used to fit models: 1) stratification by the three main tree species groups spruce, pine and deciduous resulted in species-specific models that can utilize a satellite-based species map for improving predictions, and 2) stratification by species and maturity class resulted in stratum-specific models that can be used in forest management inventories where each stand regularly is visually stratified accordingly. Stratified models were compared to general models that were fit without stratifying the data. Results: The species-specific models had relative root-mean-squared errors (RMSEs) of 35%, 34%, 31%, and 12% for volume, aboveground biomass, basal area, and Lorey’s height, respectively. These RMSEs were 2–7 percentage points (pp) smaller than those of general models. When validating using predicted species, RMSEs were 0–4 pp. smaller than those of general models. Models stratified by main species and maturity class further improved RMSEs compared to species-specific models by up to 1.8 pp. Using mixed-effects models over ordinary least squares models resulted in a decrease of RMSE for timber volume of 1.0–3.9 pp., depending on the main tree species. RMSEs for timber volume ranged between 19%–59% among individual ALS projects. Conclusions: The stratification by tree species considerably improved models of forest structural variables. A further stratification by maturity class improved these models only moderately. The accuracy of the models utilized in SR16 were within the range reported from other ALS-based forest inventories, but local variations are apparent.


Background
Knowledge of status and trends of ecosystem services provided by forests is a key requirement for informed policy and management decisions on various scales, ranging from national to regional and local scales. National forest inventories (NFIs) typically provide this kind of information on the national to regional scale based on a set of sample plots distributed within the country (Tomppo et al. 2010;Vidal et al. 2016). The availability of remotely sensed datasuch as from airborne laser scanning (ALS) or satellite imageryopens up new possibilities for NFIs (McRoberts and Tomppo 2007). By combining the field registrations in the NFI with remotely sensed data, estimates can be improved (McRoberts et al. 2002) which, for example, enables estimates for smaller regions or even local scale (Breidenbach and Astrup 2012;Breidenbach et al. 2021). The principle is to model the relationship between field measured forest attributes and metrics derived from the remotely sensed data. This model is utilized to map information of the forest in the areas between NFI sample locations, which in turn enables estimation for smaller areas than with the use of NFI field plots alone. The use of remotely sensed data does alsoin addition to the ability to provide estimates for smaller areasenable the production of continuous forest resource maps, based on remotely sensed and NFI field data, e.g. (Nord-Larsen and Schumacher 2012;Bohlin et al. 2017;Nilsson et al. 2017).
In 2015, the Norwegian forest resource map (SR16) based on aerial image matching data, was published for the 5 Mha Trøndelag region (Rahlf et al. 2017;Astrup et al. 2019). Following the area-based approach (Naesset 2002), timber volume, biomass, basal area, tree height, site indexand the respective prediction intervalsare provided for maps with a grid cell size of 16 m × 16 m. Stand-like segments were generated using object-based image analysis and segment-level (synthetic) estimates including model-dependent uncertainties are available for the mapped forest attributes (Breidenbach et al. 2016). The map covers the entire forested area in the published region, and provides publicly available forest information at an unprecedented level of detail in the area, compared to the regional statistics previously published based on the NFI field plots. Also in 2015, the Norwegian Mapping Agency started a national ALS campaign that is planned to be finalized in 2021. ALS data from this campaign are used in the current development of SR16 which, as of today, already covers more than 95% of Norway's forest area and is freely available via an online application (NIBIO 2020). SR16 is planned to be updated annually by utilizing the newest NFI data and canopy cover loss maps (Hansen et al. 2013). The latter also opens up for new map layers such as annual volume and biomass loss (NIBIO 2020;Breidenbach et al. 2021) which in Norway is largely driven by harvests.
Stratification is a concept which is often used when predicting forest attributes using ALS-based models in smaller-scale forest management inventories (Maltamo et al. 2021). The relationship between the ALS data and properties of the forest often vary between forest types, and in boreal forest individual strata are typically defined by attributes such as main tree species, site productivity and development class (Naesset 2002(Naesset , 2014. The principle is then to fit separate models for each stratum, and predict using the model corresponding to the stratum which the area of interest belongs to. One aspect of using stratification is that the delineation between strata must be known for the entire area of interest, in order to be able to predict for that area. Whereas stratifications in for example a forest management inventory can be based on attributes derived from manual interpretation of aerial imagery, this will not be possible at a regional or national level. Definition of strata must in those cases therefore be based on information which can more easily be made available for the entire area; for example information derived from remote sensing data. It is of interest to further explore how stratification can be utilised in large scale mapping, and the effect of such a stratification. One aspect that can be a challenge when creating a forest resource map at a national level, is to handle heterogeneous sets of data. In a comparable application ) models were fitted with NFI plots from close-by locations. The ALS data used in the present study consist of hundreds of smaller projects, each with differences in attributes such as acquisition time, sensor model, sensor settings, flying altitude, and point density. The NFI data used as the ground reference are coherent in terms of data acquisition, but are acquired over a period of 5 years, and are also relatively sparsely distributed spatially. The sparse spatial distribution of NFI data in addition to the varying ALS project area sizes lead to ALS projects with only few or no NFI field plots. We therefore used mixed-effects models, in order to account for the diverse nature of the ALS data as well as the spatial distribution of the ground reference data.
The aim of this study was to document and assess the methods used to produce the forest resource map SR16 at national level utilizing heterogeneous sets of data, with a focus on these key issues: Assessing the effectin terms of accuracyof using species-specific models and models stratified by main species and maturity class, by comparing with predictions from general models. Analyzing the effect of using mixed-effects models and quantifying uncertainties at the ALS project level, as well as for the entire study area.

Study area
The 17 Mha study area consists of areas in the southern half of Norway for which ALS data were available ( Fig. 1) and comprises 2/3 of the productive forest area in Norway. The study area is located almost entirely in the boreal climate zone and the forests are mostly dominated by Norway spruce (Picea abies (L.) H. Karst.) or Scots pine (Pinus sylvestris L.). Areas dominated by birch (Betula pubescens Ehrh.) do also occur, typically towards the mountains. More information on the spatial distribution of the main tree species in Norway can be found in Breidenbach et al. (2021).

Field data
The Norwegian NFI's permanent field plots were used as field reference data, with the dataset consisting of all NFI field plots within the study area. Within circular 250 m 2 plots, all trees with a diameter at breast height (DBH) ≥ 5 cm were measured, and DBH and species recorded. A sample of the trees within the plot were selected for an additional height measurement, and the heights of these trees recorded using a Haglöfs Vertex 3 hypsometer. Maturity class was registered on plots in productive forestdefined as having a yearly increment of > 1 m 3 •ha − 1 •year − 1 . The maturity class was assessed according to a classification system ranging from class 2 (young forest) to class 5 (mature forest). Nordic species-specific allometric models were used to predict tree-level timber volume and biomass from which plot-level per-ha values were estimated. Uncertainties due to these allometric models were in the following ignored, as we consider assessment of these outside the scope of this study. More information on the Norwegian NFI can be found in Breidenbach et al. (2020). The following attributes were used in the present study: timber volume (with bark), aboveground biomass, basal area, and Lorey's height.
The species with the largest volume was defined as the main species on the plot.

Forecasting NFI field plot values
The permanent NFI field plots in Norway are remeasured every fifth year, meaning that the age of the newest registrations for a set of plots range from zero to 5 years. In order to have a coherent set of NFI field plots in the modelling dataset, plot attributessuch as volume and basal areawere forecasted to a common field reference date. In the present study this reference date was set to be after the growing season in 2019. The forecasting was carried out by adding an estimated yearly increment to the attributes. This estimated yearly increment was calculated using the observed increment at the plot in the previous five-year period, using measurements from two points in time. The observed increment was determined utilizing the accurate registration of treesincluding local tree idson the NFI field plots. The accurate registration of trees enabled us to select a set consisting of all trees measured at both points in time. Using this set, an annual increment was calculated for each field plot, based on the number of growing seasons between the field measurements. The method was chosen over the use of more general allometric increment models as we believe the locally recorded increment is the best predictor for future growth at a particular location. The method does not take into account ingrowth of trees reaching the 5-cm DBH threshold after the time of the last field registration.
The forecasting method described above is limited to locations where field measurements from two points in time exist, and can therefore only be applied at locations such as NFI field plots.
This procedure ensured that the model dataset, and therefore the model predictions, were valid for the same point in time over the entire mapped area, across multiple ALS projects from multiple years (further described in section "Predictive models"). A summary of the NFI field plot data after the forecasting is given in Table 1.
Field plot data measured in the period 2009-2013 and 2014-2018 were used in this study.

Removal of outliers
NFI field plots with harvest operations after the field work, but before acquisition of ALS data (see section "ALS data") were removed from the dataset ( Table 1). Identification of these plots were based on manual inspection of an initial dataset, and available aerial imagery.

ALS data
Data from the national ALS campaign initiated in 2015 by the Norwegian Mapping Agency were used to produce the current version of SR16. The campaign serves multiple purposes, and is scheduled to be completed in 2021. The main goal of the scanning campaign is to produce a detailed digital terrain model for the entire country.
The acquisition of ALS data in the campaign is organized into projects covering predefined areas. Each project is acquired, processed and delivered by a private contractor to the Norwegian Mapping Agency as one unit. Attributessuch as sensor type and point density are the same for acquisitions within each project area, but vary between projects. For some areas existing ALS data were used, i.e. acquired before 2015. All ALS projects were acquired between 2009 and 2019, some during leaf-off conditions.
The tessellation of the country into project areasas well as the prioritizing of these areaswere based on a range of criteria such as terrain features, municipality borders, existing ALS data and requests from local administrations. The individual projects will be referred to as ALS projects in this paper, and a summary of the characteristics of these ALS projects is given in Table 2. Two products derived from the 367 ALS projects were used: the point cloud, and a 1 m × 1 m resolution digital terrain model. The ALS echo heights in the point cloud were normalized by subtracting the terrain model height, using bi-directional interpolation.

Predictive models
The heterogeneity of the combined ALS datasetin terms of acquisition date and characteristics such as point densityposes a challenge when the ALS data are paired with field observations to create a coherent modelling dataset. The sparse distribution of the NFI field plots on a systematic 3 km × 3 km or 3 km × 9 km grid means that many of the ALS projects will contain relatively few NFI field plots (see Table 2). Fitting models for each ALS project, which is used in many local ALS-based forest management inventories (e.g. Naesset 2004) is then not directly applicable, since the number of field plots contained in each ALS project is too small. In order to utilize a dataset consisting of multiple ALS projects combined with sparsely distributed field plots we therefore used linear mixed-effects regression models, in which the systematic differences between ALS projects are accounted for in the models through random effects (Breidenbach et al. 2008). Mixed-effects models were chosen because some ALS projects did not contain any NFI field plots. The global part of the model can in these cases be used. Model parameters were estimated using the model where y ij is the forecasted field measured attribute at field plot j in ALS project i, k is the number of explanatory variables denoted ×1, ×2, .., ×k, β 0 , β 1 , β 2 , .., β k are the corresponding fixed effects regression parameters, b i is a predicted random effect, n i is the number of field plots in ALS project i and m is the number of ALS projects. σ 2 b is the variance of the random effect and σ 2 ε W i is the residual variance, where σ 2 ε is the mean square residual and W i is a n i × n i matrix. For modelling of volume, basal area and biomass, heteroskedasticity in the residuals was accounted for by letting the diagonal elements of W i contain weights obtained from a variance model. For modelling of Lorey's height all diagonal elements were 1, i.e. W i is an identity matrix (Breidenbach et al. 2016).
For modelling of volume, basal area and biomass we used ALS metric H mean as variable × 1 in Eq. 1. For modelling of Lorey's height we used H 95 as variable × 1. The selection of these most important metrics were based on previous experience. The selection of further variables is described in the next section. The models were fit using the nlme package (Pinheiro et al. 2020) in R (R Core Team 2020).

Model selection
The ALS metrics in the fixed effects terms of the models were chosen by comparing an initial model using a predefined set of metrics to a model with a set of ALS metrics obtained through a stepwise selection: Based on a leave-one-out cross validation (CV), the stepwise models were preferred if the CV revealed an improvement in root mean squared error (RMSE, see section "Accuracy assessment"). We defined a requirement of at least 1 percentage point improvement in RMSE% (see section "Accuracy assessment" for definition) for selecting the stepwise model. Otherwise the initial model with the predefined set of metrics was used. The predefined set of ALS metrics were similar for all modelled attributes, except for Lorey's height (see Table 3). The selection of the predefined metrics was based on previous experience and testing. All metrics are listed in Table 3.

Quantifying the effect of using mixed-effects models
Mixed-effects models were used in order to account for differences among ALS projects. We assessed the potential gain of using mixed-effects models by inspecting the ALS projects containing at least one NFI field plot were included, predictions in ALS projects without any NFI plots were not analyzed in this study variability of the random effects, that is, σ 2 ε in Eq. 1. An estimate for the variance of the random effect, and a 95% confidence interval for this estimate, were obtained using the nlme package (Pinheiro et al. 2020). We quantified the effect by comparing the accuracies of predictions from a mixed-effects model (Eq. 1) with predictions from an ordinary least squares model with the following general form where y j is the forecasted field measured attribute at field plot j, n is the total number of field plots, k is the number of explanatory variables denoted ×1, ×2, .., ×k, β 0 , β 1 , β 2 , .., β k are the regression parameters, and σ 2 ε is the residual variance. The model described in Eq. 2 was fitted to sets of NFI field plots according to main species, and RMSE obtained through a CV. This was compared to RMSE values derived using the same set of field plots and model variables but with mixed-effects models as described in section "Predictive models". The predefined metrics marked in Table 3 were used in both cases.

Time differences between field reference and ALS acquisition date
The acquisition dates for the ALS data ranged from 2009 to 2019 (Table 2), meaning that the time difference between the acquisition of ALS data and the field reference date was up to 10 years. Note that the field reference date denotes the date to which the field data were forecasted, as described in section "Field data". In order to account for the growth between the time of ALS acquisition and the field reference date, information about this time difference was incorporated into the models. The number of growing seasons between the ALS acquisition date and the field reference date is known for the entire dataset since the acquisition date is an integrated part of the ALS data, and the field reference date is one fixed date (Fig. 2). This enabled the time difference information to be included as an additional variable (Time diff ), and to estimate a corresponding regression coefficient in the models which is used for predictions. The expected effect of the Time diff variable in the models would be to obtain a positive coefficient in the model fitting, meaning that the predicted forest attribute values will increase with an increase of Time diff . In other words, when the value of Time diff is large, the particular ALS data represent the forest as it was in the past, and the growth that has since occurred will be taken into account in the predictions through the Time diff variable.
The value of Time diff was in the present study given as the number of growing seasons between the ALS acquisition date and the field reference date (Fig. 2). Growing seasons was used instead of years or other measures of time since it more accurately represents the time period in which the actual growth occurs. For example: the actual growth in two time periods of equal length can differ by up to one full growing season, depending on the start and end of the two periods. Note that in the present study, the value of Time diff was never negative since the field reference date was chosen to be after the acquisition date of all the ALS data used.

Stratification according to main tree species
Dividing the study area into strata and fitting separate models for each stratum has been shown to yield predictions with higher accuracy, and is commonly done when modelling forest attributes using ALS data in forest management inventories (Naesset 2002). The model dataset was stratified according to the main tree species derived from the field measurements: spruce, pine or deciduous trees, and separate models were fit to each of these strata (Fig. 3). Models based on data stratified by tree species will be referred to as species-specific models. For the predictions we used a map of the main species in order to decide which strata a prediction unit (pixel) belonged to. The species map is based on Sentinel-2 and other auxiliary data, and is further described in Breidenbach et al. (2021). As an alternative to using the field observed species, we also tried the predicted tree species for stratifying the model dataset. However, some species misclassifications resulted in large residuals with adverse effects on the models. Results from this test are not further documented here.

Stratum-specific models for use with local forest management inventories
In order to increase the usability of SR16 in local forest management inventories we provide prediction maps based on stratum-specific models. A common practice in ALS-based inventories is to divide the forest into strata, and fit separate models for each stratum (Naesset 2014). In Norway this is typically done through manual delineation and classification of forest stands according to species and maturity class using aerial imagery. To accommodate users of SR16 that already know the main species and maturity class for stands in an area of interestbecause a manual stand-delineation is availablea set of stratum-specific models were fitted. These model parameters were estimated using the same procedure as described above, but the data were split corresponding to six strata in productive forest. The strata were defined by combinations of tree species and maturity class (Fig. 3).
Prediction maps based on these stratum-specific models are available for the full study area, but the use of these data will rely on an additional or existing stratification in the area of interest.

Final models
Several tests were performed before deciding on the final set of explanatory variables. The variable Time diff was Fig. 2 Conceptual overview of the time differences in the data. The field reference date denotes the date to which the field data were forecasted (see section "Field data"). Note that the ALS acquisition date in this figure is one example out of multiple acquisition dates used in the study Fig. 3 Overview of the division of the model dataset according to main tree species strata (left), and strata defined by main tree species and maturity class (right). The number of plots is given for each subset. The stratum-specific datasets consist of plots from productive forest only, because the maturity class is registered in the field only in productive forest (see section "Field data") inspected and removed if the parameter estimate was negative. Because the field reference date was set to be after the acquisition time for all ALS data, a decrease in the predicted attribute value with increasing Time diff is not logical. We further removed the variable Terrain slope if its p-value was > 0.05.
Final models for prediction were fitted to each of three subsets of the model data, split according to main species, as described in section "Time differences between field reference and ALS acquisition date". The application of the appropriate model was based on a map of predicted species, described in Breidenbach et al. (2021).
The level of uncertainty for the predicted attribute values in the final forest resource map is important information for the user. We therefore calculated prediction intervals at the 16 m × 16 m pixel level as described by Breidenbach et al. (2016). The identity was used in the calculation of the prediction intervals for Lorey's height because the residual variance was homogeneous (no heteroscedasticity).

Accuracy assessment
Predictions were compared to field reference values at the plot level by computing the RMSE and mean difference (MD) as where y ij and ŷ ij is the field measured and predicted attribute at field plot j in ALS project i, n i is the number of field plots in ALS project i, n is the total number of field plots and m is the number of ALS projects. Predicted attributes were based on a CV. Relative RMSE and MD were calculated as a percentage of the mean field measured value, and denoted RMSE% and MD%, respectively. The accuracy was assessed overall, and for subsets of plots according to a stratum. RMSE and MD were also calculated for individual ALS projects containing at least 30 NFI field plots.

Variable selection
The comparison of predictions from models with a predefined set of ALS metrics and from models with ALS metrics derived through a stepwise selection procedure resulted only in small differences. This resulted in the selection and use of models with the predefined set of ALS metrics in all cases (the metrics are listed in Table 3).

The effect of using mixed-effects models
The main effect of using a mixed-effects model is to have a model which will account for differences between the ALS projects, and have the ability to adapt to project-specific conditions. The variability of the random effects shows that there are differences between the ALS projects that are not reflected in the fixed-effects terms, and the comparison of RMSEs show that mixed-effects models performed better than ordinary least squares models ( Table 4). The gain of using a mixed-effects model varied between the modelled attributes with the decrease in RMSE when using a mixed-effects modelin percentage of the mean field measured valueranging from 0.3% for basal area in pine forest, to 6.3% for Lorey's height in deciduous forest (Table 4).

Species-specific models
The RMSEs of the species-specific models were 44 m 3 •ha − 1 (35%) for volume, 1.5 m (12%) for Lorey's height, 5.5 m 2 •ha − 1 (31%) for basal area and 26 t•ha − 1 (34%) for aboveground biomass (Table 5 and Fig. 4). The MD% was low for all six combinations of main tree species and attributes, and were all in the range between 0 and 0.3%.
The use of species-specific modelscompared to using one general modelresulted in a reduction in RMSE% of 7.1, 2.5, 2.2 and 5.1 percentage points for volume, Lorey's height, basal area and aboveground biomass, respectively (Table 5).
In SR16, the models are ultimately used to produce predictions for a set of 16 m × 16 m pixels covering the entire study area. The species-specific models require information about the main species for each pixel, and we conducted a CV in which the main species for the leftout plot was the predicted main species for the 16 m × 16 m pixel containing the plot center. The species map used here is based on Sentinel-2 and other auxiliary data which is further described in Breidenbach et al. (2021). This procedure resulted in predictions which incorporate the errors introduced by using predicted species, in contrast to the use of the main species determined from field measurements. A slight increase in RMSE% was observed when using the predicted species, with a 3.4 percentage points increase overall for volume (Table 5). For Lorey's height, basal area and biomass the overall increase in RMSE% when using the predicted species in the validation were 2.1, 1.4 and 3.1 percentage points, respectively. For pine dominated forest, results in Table  5 imply that a general model is better than the speciesspecific models. Contrary to this, deciduous forest profits most from stratifying the dataset by tree species.

Stratum-specific models for local forest management inventory applications
Model parameters were estimated using data stratified by main species and maturity class (Fig. 3). Using the cross-validated predictions, the RMSE% and MD% were calculated for each stratum. In order to assess the effect of fitting separate models for each stratum we also calculatedfor the plots within each stratum -RMSE% and MD% with predictions from the species-specific models. The results show that the stratum-specific models in most cases yielded lower RMSE% than the species-specific models applied within the corresponding stratum ( Table 6). The RMSE% of the 6 × 4 combinations of strata and response variables listed in Table 6, ranged between 8.1% -39.0% with an average of 25.1%. The MD% for the stratum-specific models ranged from − 0.1 -0.6% with a mean of 0.2%. The MD% were for all 6 × 4 combinations of strata and attributes closer to zero with the stratum-specific models than with the general models (Table 7). The average decrease in MD% due to the usage of stratum-specific models was found to be 1.2 percentage points.

Accuracy at the ALS project level
In addition to the overall and stratum-wise accuracy, we assessed the accuracy at the individual ALS project level, using the species-specific model predictions and observed values for the NFI field plots within each ALS project. The per-project RMSEs and MDs for the 115 ALS projects with more than 30 NFI field plots are summarised in Table 8. Note that predicted values used in the calculations are derived using the species-specific models, and the predicted main species (see section "Stratification according to main tree species"). On average the RMSE and MD at the individual ALS project level corresponded with the observed overall accuracy. For the 115 ALS projects, the mean per-project RMSE% were 34.7%, 14.0%, 29.9% and 34.4% for volume, Lorey's height, basal area and biomass, respectively (Table 8). RMSE% values within individual ALS projects ranged up to 59.0%, 17.6%, 46.6% and 54.2% for volume, Lorey's height, basal area and biomass, respectively. Correspondingly, the highest MD% values at the individual ALS project level were 15.9%, 3.7%, 24.3% and 19% for volume, Lorey's height, basal area and biomass (Table 8). Plots of observed vs predicted volume for the individual ALS projects with the largest RMSE and MD are shown in Fig. 5.

Discussion
Large-scale forest resource mapping using 3D remotelysensed data is an active endeavor in several countries that offers new opportunities for forest management and research. The plot-level RMSEs in this study were comparable to RMSEs found in other similar studies: Nilsson et al. (2017) describe a nationwide mapping of forest in Sweden using ALS and NFI field plot data, and report RMSEs after a CV in three of 397 ALS projects: 31-44 Table 4 Effect of using mixed-effects models at the large scale: 1) Variability of the random effects, shown as the estimated standard deviation of the random effect (σ b ) and a 95% confidence interval for this estimate in parenthesis. 2) Differences in RMSE between using a mixed-effects model, and an ordinary least squares model, ΔRMSE = RMSE (ordinary least squares model) -RMSE (mixed model). The ΔRMSE is also given as a percentage of the mean field measured value, in parenthesis. More details in section "Quantifying the effect of using mixed-effects models"   . 4 Observed vs. predicted volume, Lorey's height, basal area and aboveground biomass for NFI field plots with spruce, pine and deciduous trees as the main species. Predictions are based on a CV using the species-specific models m 3 •ha − 1 for volume, 1.2-1.9 m for Lorey's height and 4.9-5.9 m 2 •ha − 1 for basal area. The mean per-project RMSEs in the present study were all within the range reported by Nilsson et al. (2017). It could be noted that the results from the Swedish study were from only three of 397 ALS projects, and that the overall RMSE for all ALS projects might be different. This comparison only gives a cursory impression because RMSEs are affected by the properties of the studied population. Differences between ALS projects were handled differently in the Swedish approach, namely by using local regression models. Nilsson et al. (2017) fitted separate models using data from geographical sub-regions. Predictions were then made for each of the 397 ALS project areas separately, using a local model. Data from the local ALS project as well as the closest similar projects were used. The Swedish ALS data were collected in a more systematic fashion than in Norway, with rectangular and regularly sized project areas, and were therefore more suitable for using an approach with local models. The mixed-model approach used by us does, however, simplifies modelling, since separate models are not needed for each individual ALS project.
The mixed-effects model approach is intended to incorporate in the models differences between ALS projects, and the results show that the mixed-effects models performed better than ordinary least squares models. The improved performance was however most pronounced for volume and biomass, and less for Lorey's height and basal area.
Nord-Larsen and Schumacher (2012) report in a similar study using ALS and NFI data from Denmark an overall RMSE of 80.4 m 3 •ha − 1 for volume and 7.4 m 2 •ha − 1 for basal area. This corresponds to 44.8% and 47.7% of the reported mean field measured values for volume and basal area, respectively. The RMSEs reported in the Danish study areboth in absolute and relative termsslightly higher than the RMSEs for volume and basal area we found. Several differences between the data in the two studies may be the underlying cause for the observed difference in accuracy. The mean ALS point density was, for example, lower in the study by Nord-Larsen and Schumacher (2012).
Using separate models fitted to stratified sets of data has been shown to improve the prediction accuracy (e.g. Naesset and Gobakken 2008), and it is commonly used in forest management inventories in the Nordic countries (Maltamo et al. 2020). Stratification did also improve the predictions in our case, with the use of separate models according to main tree species resulting in a lower overall RMSE, even when including the errors introduced by using remote sensing based predictions of main species. For pixels in which the predicted main species is wrong, this approach might result in larger errors for the modelled attributes, compared to using one general model. One related and positive side-effect is,  however, that any improvements in the species predictionsfor example by using multitemporal instead of bi-temporal satellite data (Persson et al. 2018;Breidenbach et al. 2021)will also lead to improvements in the modelled attributes.
The results further show that stratum-specific models using species and maturity class improve the accuracy. The improvements in terms of RMSE over using species-specific models were however modest, with slightly larger improvements in MD than in RMSE.
The RMSE and MD at the individual ALS project level was on average similar to the overall values, but larger values of RMSE and MD were observed for some individual ALS projects. Factors such as acquisition specifications, flying date and species predictions contribute to the prediction accuracy and uncertainty at the project level. While our approach shares many aspects of the methodology used in local forest management inventories, application at a national level has to rely on automated processing to a greater degree. Due to the large number of field plots, only limited manual inspection of the dataset was carried out. Comparison of predicted versus observed attribute values at the individual ALS project level (e.g. Fig. 5) indicate that further manual inspection of individual observations in the modelling dataset might reveal additional outliers. Identification and removal of such outliers could improve the prediction at the ALS project level.
Systematic errors are of concern when applying models in one ALS project which have been developed with data from other ALS projects (Naesset 2009;Noordermeer et al. 2019;Tompalski et al. 2019). The average and per-project MD% for the projects with 30 or more NFI field plots show that the systematic errors overall were small, but with systematic errors at the project level of up to e.g. 16% for volume. Measures could be taken to increase the accuracy locally, for example by combining the NFI field plots with a small set of local field plots .
The current study was conducted in boreal forest with a structure and species composition which differ from forests in temperate and tropical regions. Modeling of forest attributes using ALS data has however also been demonstrated in other climate zones such as in temperate (e.g. van Leeuwen and Nieuwenhuis 2010) and tropical forest (e.g. Asner et al. 2012;Ioki et al. 2014;Hansen et al. 2015;Mauya et al. 2015). Building on the experience from such studies, an application of the mapping described in the current study should be possible also for these forest types. The availability of representative field measurements is however a requirement to ensure consistent and reliable maps. Because field and ALS data are relatively costly to acquire, it is likely the main obstacle for development of similar forest resource maps in many countries. An alternative to ALS data would be to use available satellite imagery, as described by Reese et al. (2003). This would however influence the accuracy and types of attributes possible to include in the map. While further research is needed, there are promising results from predictions in temperate and tropical forest using data derived from Sentinel2 imagery (e.g. Gonçalves et al. 2019;Ahmadi et al. 2020).
The combination of NFI and remotely sensed data is useful for mapping forest attributes and the improvement of estimates on various scales. It is currently analysed whether further variables, for example those related to stand age (Schumacher et al. 2020) and the protective functions of forests, can be included in the forest resource map SR16. Because of their fine spatial resolution and relatively high accuracy, forest resource maps offer new insights and decision support for example with respect to harvest monitoring ) and road network planning. Table 8 Summary of RMSE and MD calculated using NFI field plots within individual ALS projects. Predictions from a CV, using the species-specific models and predicted species as described in section "Stratification according to main tree species". The 115 ALS projects with 30 or more NFI field plots were included