REFEWOODLAND:
Tree biomass outside forest areas

Technical report

1 Project overview

The REFEWOODLAND-project was initiated with the intention of complementing existing estimations of above ground tree biomass (AGB) in forest areas with comprehensive information collected in outside forest areas, in order to calibrate a remote sensing data-based model to seamlessly estimate existing AGB on areas outside forests.
With this model a nationwide above ground tree biomass-estimation outside forest areas becomes directly possible, supporting greenhouse gas inventory reporting.

The base remote sensing data used derives from the first national airborne laser-scanning campaign (2018-2025) conducted by Swisstopo and resulting in the swissSURFACE3D (ALS) data. The timing of field recordings was set to be parallel with the existing flight schedule for swissSURFACE3D. For the field measurements, a stratified sampling-design was developed to account for regional differences in biomass-characteristics whilst ensuring a cost-effective sampling procedure. These sample-plots were probed in the following field-campaign for whole Switzerland.

After collecting data for all sample-plots (2019-2023), the above-ground biomass (AGB) was calculated after Chave et al. (2014) and an ALS-based mixed-effects model was set-up and calibrated, using vegetation-volume metrics as explanatory variable, allowing to predict above ground tree biomass-stocks outside forest areas.

2 Reference data collection

2.1 Sampling design

At the beginning of the project, a stratified sampling design was chosen to account for consideration of homogeneous physiological appearance (allometry) of targeted vegetation, as well as providing a cost-effective sampling method rather than a biomass inventory. This sampling approach is documented in Gardi (2018). Price (2023), describes the data derived thresholds for the model application. It is important to point out, that this stratification process does not fulfill the requirements of an inventory design, but is designed, to achieve the best consideration of potential influences in the reference data, for improved calibration of the later created model.

The introduced stratification is based on Swiss Land Use Statistics which are enhanced with additional spatial data in its most recent form. As the most recent form of the Land Use Statistics (Arealstatistik (AREA) 4, 2014-2018) was not yet finalized at the beginning of the project, missing data was completed with the former classification of AREA version 3.

As additional data, a 25 m raster was calculated from the aerial image based vegetation height model (Ginzler and Hobi, 2016) to characterize a vegetation canopy cover (vegetation > 3 m) and the National Forest Inventory (NFI) forest area layer was used to mask our forest areas.

Combined, the AREA data was then filtered by canopy cover > 5 % (NFI vegetation height model (remaining n = 571’783 points)) and a distance to forest areas (NFI) of 30 m (remaining n = 360’334 points), as those regions are either out of interest or thematically covered by the NFI.

Afterwards, a list of parameters that potentially determine growth were selected based on expert input, leading to a stratification-clustering, which is needed for equal sample drawing:

Parameter	Description	Thresholds for defined strata
Landuse	Swiss Greenhouse Gas Inventory Landuse/Landcover Combination Categories (CC) derived from the Swiss Land Use Statistics (AREA 3 & 4)	Settlement (L1) Agriculture (L2) Special crops and shrubs (L3) all others (L4)
Bioregion	Biogeogrpahical regions of Switzerland, defined by FOEN	Jura and Mittelland (B1) Alps (all others) (B2)
Elevation	DHM50 swisstopo	B1: x < 450 m (E1) 450 ≤ x < 570 m (E2) 570 m ≤ x (E3) B2: x < 1100 m (E1) 1100 ≤ x < 1’890 m (E2) 1’890 m ≤ x (E3)
Tree height	Mean per pixel (25 m) tree height from the VHM model NFI (Ginzler and Hobi, 2016) and considering the following boundaries: x < p0.05 (H1) p0.05 ≤ x < p0.95 (H2) p0.95 < x (H3)	H1: x < 6.86 (B1_E1) x < 6.43 (B1_E2) x < 7.08 (B1_E3) x < 6.7 (B2_E1) x < 6.57 (B2_E2) x < 5.01 (B2_E3) H2: 6.86 < x < 15.1 (B1_E1) 6.43 < x < 14.5 (B1_E2) 7.08 < x < 16.3 (B1_E3) 6.7 < x < 15.3 (B2_E1) 6.57 < x < 15.6 (B2_E2) 5.01 < x < 12.7 (B2_E3) H3: x ≥ 15.1 (B1_E1) x ≥ 014.5 (B1_E2) x ≥ 16.3 (B1_E3) x ≥ 15.3 (B2, E1) x ≥ 15.6 (B2_E2) x ≥ 12.7 (B2_E3)
Canopy cover	Mean per pixel (25 m) tree canopy cover from the VHM model NFI (Ginzler and Hobi, 2016) and considering the following boundaries: x < p0.05 (C1) p0.05 ≤ x < p0.95 (C2) p0.95 < x (C3)	C1: x < 0.174 (B1_E1) x < 0.158 (B1_E2) x < 0.162 (B1_E3) x < 0.152 (B2_E1) x < 0.14 (B2_E2) x < 0.11 (B2_E3) C2: 0.174 < x < 0.686 (B1_E1) 0.158 < x < 0.61 (B1_E2) 0.162 < x < 0.636 (B1_E3) 0.152 < x < 0.634 (B2_E1) 0.14 < x < 0.53 (B2_E2) 0.11 < x < 0.39 (B2_E3) C3: x ≥ 0.686 (B1_E1) x ≥ 0.61 (B1_E2) x ≥ 0.636 (B1_E3) x ≥ 0.634 (B2_E1) x ≥ 0.53 (B2_E2) x ≥ 0.39 (B2_E3)

All factors combined, result in 217 strata in total considering all parameters with 165 strata being used for the model training. Missing strata were related to canopy cover and vegetation height, which could not be found in the fields. The stratification raster for the field-data collection is provided in the data-package as “strata_comb_code.tif” and the related coding table “strata_LUT_allcomb.csv”.

2.2 Sample selection

2.2.1 Basic principle

To support the economic feasibility of the project, the sample drawing is proceeded randomly, but with higher probability inclusion for positive economic aspects. This is done in a two-phase procedure:

1. Selection of sample-cluster locations
2. selection of related, additional sample points (6-9 per cluster)

This procedure is done to ensure, that each cluster can be handled by the field team in one day, which helps to minimize unnecessary travel and set-up times during field work.

2.2.2 The selection of sample-cluster locations

For the first selection, at least one sample point should be included per stratum and panel, with panel 1 representing Eastern-, Western-, Central- and Southern Switzerland, and panel 2, representing Graubünden, Wallis, Northern Switzerland and Bern.

The distribution of sample plots was designed as random but considering a homogeneous distribution of points. This is achieved by dividing the same number of samples up in both panels, and at least assign one plot per stratum.

Additionally, for each point an inclusion probability P by accessibility, with

\[ P = \frac{1}{t_{ÖV} + 10 \cdot d_{ÖV} + 30 \cdot d_{Str}} \]

\(t_{ÖV}\) : Travel time from Bern, by public transport (Bundesamt für Landestopografie, 2011)
\(d_{ÖV}\) : Distance of the plot to the nearest public transportation station (OSM)
\(d_{Str}\): Distance of the plot to the next street (OSM)

was defined, ranking the plots for sample drawing.

In total, and when applying this method, n₁ = 521 core-plots were selected for further processing.

2.2.3 The selection of additional sample points

For each identified cluster- center, 6-9 additional sample plots were drawn within a radius of 0.2 - 1 km. If no 6 additional sample points were found, the entire cluster was discarded.

This selection results in 466 clusters, with 4644 potential sample plots for recording (Fig. 1).

    Figure 1: Visualization and characterization of all sample points including the respective strata

For the final plot selection and to prepare field recordings, from all sample points, at least 1’500 plots were selected with minimum 750 plots for each panel. For this the points were individually assessed, considering:

• danger assessment of field recording
• the access to the center of the sample plot
• the overall access of the plot, and
• the vegetation cover peculiarities (high cover or vegetation height).

This leaves 350 clusters with an average of 4.2 sample plots.

2.3 Field recordings

2.3.1 Preparation

In advance of the field work, for each plot-cluster a pdf-sheet was created. It includes aerial images, visualizing the location of the plot centers and the respective borders (Gardi, 2019), for navigating, better orientation in the field and to verify the sampling site parameters later in the fields. As the use of GPS solely to locate the plot centers didn’t provide the accuracy needed for the later modeling step, determining the plot center by using the provided aerial image, was defined as standard procedure and delivered faster and more accurate results.

2.3.2 Field measurements protocol and data collection

For each plot, all vegetation types were recorded in a concentric sample circle design, covering two radii: 12.62 m (a05) and 21.85 m (a15). The larger sample circle is taken as sampling plot with higher vegetation-inclusion probability, in case the inner circle contains less than 6 trees above 3 m height and a dbh ≥ 5 cm.

The detailed instructions for the field work of which workflow to follow is documented in the “field measurements protocol” (Gardi, 2019), including a list of parameters and the use of the data-collection app, which was setup specifically for this project, based on open data kit (ODK)- services. This setup of the ODK-service allows to directly submit the filled forms from personalized mobile handhelds to a server and to collect all plot-recordings from different mobile devices (8 field workers, working also in parallel between 2019-2023) in one database.

2.3.3 Field recordings timeline

Within the REFEWOODLAND-project, timed field recordings were conducted (Fig. 2), synchronized with the parallel LiDAR recordings conducted by Swisstopo as basis to provide the national swissSURFACE3D LiDAR dataset (Swisstopo, 2021). This allows to calibrate a LiDAR-based model, based on the field data collected and with low influences of temporal effects as pruning, removal of trees but also vegetation growth.

For each field worker, a selection of plots were loaded in sets on their ODK mobile service, which were timed with the ALS flight schedule and was from then on visible in their personalized field computers and ready for recording.

In the field data collection process, an overall number of 1738 plots were being recorded between 2019-03-30 and 2023-10-06.

Figure 2: Timeline of field-recordings for reference data, including systematical errors (cf. chapter 2.4) and plots, where the access was denied despite arrival.

2.3.4 Timing of the ALS and field recordings

These field recordings were timed after the ALS flight schedule and thus the las-recording date (Fig. 3) and resulted in an average early recording of 157 days.

    Figure 3: Timing overview of field-data recordings (delay of field-recordings: negative values [d], delay of las-recordings: positive values [d]))

The distribution of time-delay-differences (Fig. 4) show both issues faced during data collection: A late start of the project, with ALS recordings already in progress, as well as unknown delays of the ALS recording-schedule and data collection issues due to heavy snow-cover (Swisstopo, 2024).

Figure 4: Frequency distribution of timing success of field and LiDAR recordings (delay of field-recordings: negative values [d], delay of las-recordings: positive values [d]), of all usable field recordings (access granted)

2.4 Data preparation - cleaning for systematic errors

The first data-cleaning procedure was conducted continuously when recordings were submitted to the database. This leaves the opportunity to discuss with the field team about peculiarities of the plot or the recording procedure which could influence the data quality.
In this step, changes in the initial situation (absence of trees), localization problems (non-accessible, missing references or alternatively missing GPS in the field) or errors in the las-files were adressed.

The second data-cleaning was done after the LiDAR data was matched with the field data, and a plausibility control could be made (input: rfwdlnd_XXX.csv –> cleaned output: XXX_data.csv). In case of height measurements, data was marked as not plausible, if the height differences between field data and LiDAR data exceeded 40 % related to the lower nominator, without potential vegetation influencing effects possible. Reasons for deviations were related to personal measuring errors (invisible crown top) or failures of the device, observed in high air-humidity situations, among others.

All data cleaning procedures (Fig. 5) were conducted and documented in the “3_field_data_cleaning.R” - R script, provided also in the data package for documentation and to enable retrospective adjustments.

In total and in course of this procedure, 112 plot-recordings were removed from the data. All removed plots are shown in the following interactive map, commenting the reason of removal (see also removed_plots.csv).

    Figure 5: Removed plots and reason for removal after data-cleaning

2.5 Field data storage and basic characteristics

The final collected and cleaned data is stored in three data tables, which are linked to a generative ID (Cluster_ID, Plot_ID, Tree_ID and Shoot_ID) as the primary identifier:

• plot_data.csv : Containing plot-related data, like submission date
• tree_data.csv : All tree-related data, as species, tree-height and diameter
• shoots_data.csv : In case, shoots were identified in the tree-data, the individual shoots are characterized in this table.

2.6 Biomass estimation from terrestrial data

For the model, the biomass as response variable has to be calculated from the field recordings as reference. A previously conducted project, investigating tree-biomass estimation for urban areas in Switzerland (Gardi et al., 2016; Mathys et al., 2019), showed that the biomass estimation following the function presented by Chave et al. (2014), shows good results for Switzerland. It is applied in the following form:

\[ agb = \left(0.0673 \cdot (\rho \cdot d^2 \cdot h)\right)^{0.976} \]

with ρ as the wood density of the respective tree species, d the diameter of the tree and h the tree height. The values of d and h were determined in the field, while wood density values of individual tree species were taken from the wood-density database provided by Zanne et al. (2009), which was further completed for missing values and documented in the respective comments column of the data collection (wood-densities.csv). In case a species was missing, the density of the respective wood-genus was taken for calculation.

For 2203 trees, which are composed of several shoots, the equation was applied separately for each shoot individually and the biomass was summed up to represent the whole tree.

2.7 LiDAR Indices

To create an explanatory variable, LiDAR indices were derived from the swissSurface3D point-cloud dataset (Swisstopo, 2021). This data contains a pre-classification, providing the used classes “vegetation” and “buildings” for filtering purposes.

To derive the metrics used for the modeling step, a filter on vegetation data was set and voxel metrics in a resolution of 0.5 m were calculated. This focus on the vegetation data and the voxel-based metrics ensure, that the density of returns in the LiDAR data remain consistent for all strata and is not influenced by e.g. water-areas where reduced returns can be expected.

With this, the mean height of the vegetation top-voxels (tox.avg) were calculated and the number of voxels were assigned to the sample circle with the smaller radius (a05) (all.050.veg.tox).

All calculated LiDAR indices are documented in the respective “las_metrics.csv” table.

2.8 Calculating the green volume

Based on this voxel-metrics, for each plot, a green volume-value is calculated with:

\[ \text{gv} = \text{tox.avg} \cdot \left( \frac{\text{all.050.veg.tox}}{2000} \right) \]

as product of the mean vegetation height of the top voxels identified, and the vegetation cover share per a05 sample-circle (500 m²), containing a maximum of 2000 voxels possible. This green volume parameter is further used as the explanatory variable for the modeling step.

2.9 Combined data for training

The results of the two prior working-steps, needed for the model training step and including the calculation of additional LiDAR indices for testing, are documented in the following .csv tables:

• plot_results.csv
• shoot_results.csv
• tree_results.csv

2.10 Data filter to address border- effects

To train the model, the data was finally selected according to the technical training and prediction scenarios, specifically the consideration of inclusion and exclusion of vegetation at the spatial prediction borders (Fig. 6) .

For each detected tree in the field, a probability is given, that part of its biomass is virtually cut-off the sample area as it lies outside its 500%nbsp;m² area (a05). As the reference biomass is calculated based on the center of the tree, and assumes that all branches are related to the LiDAR indices derived, and as the ratio between cut-off vegetation-points and reaching in vegetation-points is not balanced, as in all cases a tree is detected inside, not outside, a potential overestimation of the biomass modeling could occur. This effect can be amplified in the prediction process, when branches from adjoining trees outside the prediction area are detected inside the prediction area.
As the model already references a higher amount of biomass, to a lower count of LiDAR returns and predicts additional biomass also for branches, this can lead to excessive prediction values.

To minimize this effect, all sample areas where no related vegetation (e.g. trunk center) was recorded in the field measurements, but whose vegetation was detected from the LiDAR data, are also retained in the training data for the model calibration in order to take the double counting effect into account.

Figure 6: Visualization of border effect in the training (left) and prediction (right) process considering different scenarios: (a) self compensating double-border effect situation, (b) overestimation of the predicted biomass (border effect) due to a reduced green volume detected in the LiDAR data during the training process and (c) training data adjustment to compensate for systematic prediction errors (border effect)

2.11 Plot reports

For all data, which was used for the modeling step, plot reports were generated for 1566 plots, to facilitate plausibility control and the visual assessment of the model performance of the REFEWOODLAND biomass model. These plot-reports contain both, field-recording elements and LAS-metrics derived from the matching area for comparison and serve as a platform to include collected pictures of the sample-plots and peculiarities in the plot-documentation (Starke, 2025a).

2.12 R-scripts for data-processing

All tasks of data preparation, modeling and evaluation steps were conducted with the statistics software R (R Core Team, 2023), for spatial prediction and to create the presented rasterdata the software ArcGIS Pro was used.
All relevant scripts developed, are described in the “Metadata REFEWOODLAND R-Scripts” document (Starke, 2025b).

3 Results and data provided

3.1 The mixed-effects model

We used above ground tree biomass (field inventory data) and green volume (ALS data) per sample plot to fit a linear mixed-effects model using the R function “lmer” from the package “lme4” (Bates et al., 2015). Response (above ground biomass) and explanatory variable (green volume) were square root transformed to obtain approximately normally distributed residuals. To account for the hierarchical data structure and spatially related sample plots, strata was included as a random effect term. Thus, the mixed-effects model formula is given by:

\[ \sqrt{\text{agb}} \sim (\beta_0 + \beta_1 \cdot \sqrt{gv_{ij}}) + (\gamma_j + \text{strata}_j) + \varepsilon_{ij} \]

where \(agb\) is the above ground biomass in t/plot (500 m², a05), \((\beta_0 + \beta_1 \cdot \sqrt{gv_{ij}})\) is the fixed effect of green volume, \((\gamma_j + \text{strata}_j)\) is the random effect of strata (random intercepts), and \(ε_{ij}\) is the error term.

As the number of observations per stratum is in some cases strongly limited, random slopes are not included in the model to avoid the over-fitting of random effects.

3.2 Model description

The presented model is characterized with a conditional \(R^{2} = 0.644\) (fixed and random effects) and a marginal \(R^{2} = 0.628\) (variance of fixed effects only) (Lüdecke et al., 2021), which highlights the improvement when considering random effects and supports the stratification success. The normalized visualization of the terrestrial AGB values with the model predictions show a potential underestimation of biomass on higher AGB plots, but with low frequency of observation (Fig. 7).

Figure 7: Normalized visualization of above ground biomass field measurements (AGB) and predicted AGB values by the REFEWOODLAND model

To validate this model, a cross-validation method was applied without considering random effects (leave-one-out, 10-fold). The results show high error percentages, especially for low biomass proportions (Fig. 8). Partly, this is to be seen as a result of the consideration of the border effect error in combination with marginal vegetation volumes (negative values mean overestimation of the biomass), as discussed in the section 2.10. One in this way affected sample point, is shown exemplary in Fig. 9.
Overall and with a median of 18.4 %, the predicted AGB remains underestimated in the majority of observations, especially predicting larger AGB values. Overall this leads to an average underestimation of the predicted AGB values of 13.26 % (cross-validated results).

Figure 8: Excerpt of the observed cross-validated errors of the above ground biomass prediction without random effects (100 % error boundaries)

Figure 9: Illustration of the canopy height model of the highest observed error (plot R-064-6). An AGB field estimate of 20 kg / plot (inner circle) is superimposed by strong edge effects of green volume, resulting in a prediction of 3300 kg / plot

3.3 Error prediction

Next to the AGB prediction value, and to estimate a confidence interval of the conducted prediction, we used a bootstrapping method when creating the final raster. This method supports the estimation of confidence intervals for mixed effects model by incorporating both the variance of the fixed and random effects (Bates et al., 2015).

As result, a 95% confidence interval for the prediction of the biomass in each 25x25 m pixel (excluding forests, buildings, water bodies) was determined by parametric bootstrapping (1000 iterations). For this, the function “bootMer” from the package “lme4” (Bates et al., 2015; cf. model formulation in S1) over a dataframe with a value for each 25 m pixel was applied.

3.4 Spatial model application

The modelling and bootstrapping results were applied spatially to the 25 m resolution raster data (corresponds closely to the diameter of the circular 500 m² (12.62 m radius) sample plots of the field inventory) of green volume derived from swissSURFACE3D, whereby green volume is calculated as described in section 2.8, but 2000 is replaced by 2500 to obtain pixel-based predictions. Green volume is only calculated for pixels with average vegetation height greater than 3 m.

The spatial prediction rasterizes the dataframe resulting from the “bootMer” bootstrapping to a 3-band raster (Fig. 10) where band 1 is the median prediction, band 2 the 0.025^th quantile and band 3 the 0.975^th quantile, giving a 95 % prediction interval for the predicted (median) value. Any values predicted to be below 0 (including in the 0.025^th quantile values) were replaced with 0 since it is not logical to have negative tee biomass in the field.

The NFI forest mask, which applies to the raster data labeled as “masked”, was buffered by 30 m to ensure that any forest trees which occur at the forest edge, were not included in the estimates of TOF AGB. Without this masking buffer, TOF AGB would be inflated in areas neighboring forests, since the 2017 version of the forest mask tends to underestimate forest area. A 30 m buffer was chosen to ensure a full pixel width from the forest edge is covered.

File name	Spatial resolution	Units	Value no green volume	Value forest mask	Usage
agb_T_pixel_CH.tif	25 m	t/pixel	NA	prediction	Model output (internal)
agb_T_pixel_CH_masked.tif	25 m	t/pixel	NA	NA	Model evaluation
agb_T_pixel_CH_masked_forestNA.tif	25 m	t/pixel	0	NA	Forest mask assessment
agb_T_pixel_CH_masked_noNAs_GHGI.tif	25 m	t/pixel	0	0	GHGI
agb_T_ha_CH_masked_forestNA_100m.tif	100 m	t/ha	0	NA	Public data, sharing

Figure 10: Screenshot of the provided 3-band raster-layer (agb_T_pixel_CH_masked_forestNA.tif) showing all values for a sample raster cell with a resolution of 25 m (background: SWISSIMAGE)

An additional 1 band raster (Fig. 11) is also produced which gives the prediction interval for each pixel = the 0.975^th percentile value minus the 0.025^th percentile value. Units of the prediction are tonnes of above ground tree biomass per 25 m resolution pixel (625 m²). Forest areas are masked out (to NoData) using the NFI forest mask (Waser et al., 2015), version 2017, with a 30 m buffer.

File name	Spatial resolution	Units	Value no green volume	Value forest mask	Usage
PI_range_T_pixel_CH_masked.tif	25 m	t/pixel	NA	NA	Model evaluation
PI_range_T_pixel_CH_noNAs_GHGI.tif	25 m	t/pixel	0	0	GHGI

Figure 11: Screenshot of the provided 1-band raster-layer (PI_range_T_pixel_CH_masked.tif) showing the forecast interval for a sample raster cell with a resolution of 25 m (background: SWISSIMAGE)

3.5 Model-based evaluation

For the following evaluation steps, the predicted biomass was spatially summarized by the first three levels of the strata (land use (AREA 4, 25.11.2021), bioregion and elevation) to average values per hectare of the given strata (Fig. 12).

The model-based evaluation contains all input data, with a green volume value > 0. Effects of missing vegetation and thus no prediction error are ruled-out by this evaluation and highlights the direct performance of the model and characterizes the vegetation detected for the respective stratum observed.

Figure 12: AGB prediction statistics for raster cells with AGB detected

3.6 Land-use based evaluation

For the AREA land-use based evaluation (AREA 4, 25.11.2021), areas without green volume detected are included in the evaluation to seamlessly represent the area cover (Fig. 13). For this, these areas are given an AGB value of 0 (as opposed to NoData). To obtain a measure of prediction variance per strata, the average value for the 0.025^th percentile and the 0.975^th percentile were also calculated for each strata with 0 values, after below 0 predictions were converted to 0.

Figure 13: Land-use and area-related AGB prediction statistics related to the first three levels of the strata (land use, production region and altitude)

The values shown in Fig. 14 apply to the first level of stratification.

Figure 14: Land-use and area-related AGB prediction statistics related to the first level of the strata

4 Discussion

4.1 Field data

4.1.1 Sampling design

In the first step of the project, data on vegetation was collected to describe the biomass and carbon stock outside forest areas. This data takes into account potential variations in vegetation type that can be found across predefined strata in Switzerland. In the modeling step, this differentiation proved to be an improvement in model fit, although it should be noted that some strata are characterized by a small number of observations, which contain a high sampling error and thus lead to a tendency for the model to be overfitted. However, the method was used to ensure a cost-effective distribution of sample areas and was not intended as a standalone inventory method. Therefore, the regionally observed sampling error was not considered problematic for the further procedure. Consequently, the stratification and sampling procedure were retained in their intended form, but further measures were implemented in the modeling step by limiting the power of the associated random effects.

4.1.2 Field data quality

To improve the field data quality and quantity, it must be mentioned, that the field- measuring error varies in an unpredictable way. During the field-data collection steps it was frequently observed, that the height measurements which were made with laser measuring devices were not correct at foggy air-conditions. Other problems were identified in finding the centre of the sample plot accurately. Solutions for further data acquisition steps could thus be the integration of mobile laser-scanning devices in the data collection step, spatially-referenced with high accurate GPS modules. This would assure the precise localization of the plot-center and offer the opportunity to integrate high point density evaluation steps in the biomass estimation for more precise training data.

4.2 Model

4.2.1 Biomass references

The derived above ground biomass values were based on Chave et al. (2014), which had proven itself in prior research for trees in urban environments (Mathys et al., 2019). However, it must be noted, that its application for other land-use systems remains unclear what raises the uncertainty of results.

This approach is also highly dependent of the vegetation height and the dbh of the trees observed. Especially in high dbh-classes, allometric functions show a larger variability, than expected for e.g. forest areas. Reasons could be found with long-term modifications of crowns (pruning), or strongly influenced changing growth-conditions for trees, especially considering light conditions (dense vegetation scenarios (shoots) and influencing buildings). For the model development and due to the missing individualization of trees for the model application, adjustments could be considered with adjustments of the vegetation height and crown cover, which was partially addressed in the stratification process, but still leaves the error for human-based influences and the separation of individual trees and thus allometric variability, especially in closed crown cover situations.

4.2.2 Explanatory variables, resilience and sensitivity

As an explanatory variable, a variable should be developed that is as independent as possible from the point density of the collected data without compromising the accuracy of the prediction. This has two advantages: resistance to disturbances during collection, such as water surfaces, and applicability of the model to data from different sources, as used for time analysis.

The first potential advantage was tested using a variable from an earlier phase based on the number of LiDAR returns (counters), which resulted in high errors along riverbanks and other water-influenced areas. This effect was completely avoided with the modified approach and therefore does not need to be considered further.

The second potential advantage of sensor resilience lies in the possible interaction with the sensitivity of the variables and thus of the model. This sensitivity to change is directly related to the identification of vegetation per voxel in the vertical direction (voxel stack) and to the identification of the highest vegetation points in the scanning process (top voxel). While the vegetation area required to calculate the vegetation cover has a high detection probability, as the scan is also predominantly vertical, the detection of the highest vegetation points and thus the voxels may be subject to greater uncertainty.

In order to detect changes, the vegetation within the next voxel height range must be detected, which could easily be overlooked due to its appearance as branches. However, since this feature remains the same during the initial detection and during a possible follow-up measurement, deviations in the detection rate of the outer top voxels could be caused by higher point densities in the secondary LiDAR recording. Since there is also a certain probability of detection in older ALS measurements, this effect appears to be negligible, but could be given further consideration in time series analysis.

4.2.3 Spatial prediction quality and model application

Even though individual tree biomass measurements lie within the ranges of the model predictions (Gardi et al., 2016) and the model in general shows a higher performance compared to existing methods (Price et al., 2017), prediction results must be handled with care.

In this case, the predicted value of a single raster cell contains a model-related error, building up for multiple reasons, such as diverse vegetation appearances outside forest areas, pruning effects or also false ground detection in the LiDAR data. Due to the large data behind the model, handling these influences, the model is especially useful for estimates on a larger scale. For the model use, it is recommended to evaluate larger areas of interest only, which thus guarantees the best model performance and the best mean biomass estimate possible.

4.3 Conclusions

The model successfully achieved its goal of predicting above-ground biomass outside forest areas and improves existing methods with new data usable. The model performed well for its intended purpose to support the LULUCF AGB reporting for TOF, but it has limitations that are based on the green-volume-based method applied, substituting the terrestrial dbh-value with a LiDAR-based green-volume variable in the prediction process, as sensor-robust AGB prediction solution. In further development steps, thus the sensitivity of this model approach to representing changes over time and based on different data sources is of particular interest, to highlight the additional value of the model in AGB-change observation, based on multiple input data.

5 References

Bates, D., Mächler, M., Bolker, B., Walker, S., 2015. Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67. https://doi.org/10.18637/jss.v067.i01

Bundesamt für Landestopografie, 2011. Verkehrsmodellierung VM-UVEK (ARE), INFOPLAN-ARE, swisstopo.

Chave, J., Réjou-Méchain, M., Búrquez, A., Chidumayo, E., Colgan, M.S., Delitti, W.B.C., Duque, A., Eid, T., Fearnside, P.M., Goodman, R.C., Henry, M., Martínez-Yrízar, A., Mugasha, W.A., Muller-Landau, H.C., Mencuccini, M., Nelson, B.W., Ngomanda, A., Nogueira, E.M., Ortiz-Malavassi, E., Pélissier, R., Ploton, P., Ryan, C.M., Saldarriaga, J.G., Vieilledent, G., 2014. Improved allometric models to estimate the aboveground biomass of tropical trees. Global change biology 20, 3177–3190. https://doi.org/10.1111/gcb.12629

Gardi, O., 2019. Anleitung zur aufnahme der felddaten und qualitätssicherung V2. REFEWOODLAND data-package 16p.

Gardi, O., 2018. REFEWOODLAND: Erhebung von referenzdaten für baumbiomasse ausserhalb des waldes. REFEWOODLAND Presentation 44p.

Gardi, O., Schaller, G., Neuner, M., Mack, S., 2016. Ermittlung der kohlenstoffspeicherung von bäumen im siedlungsgebiet am beispiel der stadt bern. Schweizerische Zeitschrift fur Forstwesen 167, 90–97. https://doi.org/10.3188/szf.2016.0090

Ginzler, C., Hobi, M.L., 2016. Das aktuelle vegetationshöhenmodell der schweiz: Spezifische anwendungen im waldbereich. Schweizerische Zeitschrift fur Forstwesen 167, 128–135. https://doi.org/10.3188/szf.2016.0128

Lüdecke, D., Ben-Shachar, M.S., Patil, I., Waggoner, P., Makowski, D., 2021. performance: An R package for assessment, comparison and testing of statistical models. Journal of Open Source Software 6, 3139. https://doi.org/10.21105/joss.03139

Mathys, L., Gardi, O., Kükenbrink, D., Morsdorf, F., 2019. Non-forest tree biomass references. Federal Office for the Environment (FOEN), Bern, Switzerland.

Price, B., 2023. REFEWOODLAND: Spatial data for prediction of AGB for trees outside forest (TOF), metadata. REFEWOODLAND data-package 3p.

Price, B., Gomez, A., Mathys, L., Gardi, O., Schellenberger, A., Ginzler, C., Thürig, E., 2017. Tree biomass in the swiss landscape: Nationwide modelling for improved accounting for forest and non-forest trees. Environ Monit Assess 189, 106. https://doi.org/10.1007/s10661-017-5816-7

R Core Team, 2023. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Starke, M., 2025a. REFEWOODLAND - plot reports REFEWOODLAND data-package.

Starke, M., 2025b. REFEWOODLAND - description of r-scripts created for the project REFEWOODLAND data-package.

Swisstopo, 2024. Mail from swisstopo about the current state of LiDAR data collection, 24.08.2024 1p.

Swisstopo, 2021. Geodata of switzerland: Comprehensive overview (Technical Report). Federal Office of Topography swisstopo, Wabern, Switzerland.

Waser, L., Fischer, C., Wang, Z., Ginzler, C., 2015. Wall-to-wall forest mapping based on digital surface models from image-based point clouds and a NFI forest definition. Forests 6, 4510–4528. https://doi.org/10.3390/f6124386

Zanne, A.E., Lopez-Gonzalez, G., Coomes, D.A., Ilic, J., Jansen, S., Lewis, S.L., Miller, R.B., Swenson, N.G., Wiemann, M.C., Chave, J., 2009. Global wood density database.