Estimation of the occurrence, severity, and volume of heartwood rot using airborne laser scanning and optical satellite data

ABSTRACT

Rot in commercial timber reduces the value of the wood substantially and estimating the occurrence, severity, and volume of heartwood rot would be a useful tool in decision-making to minimize economic losses. Remotely sensed data has recently been used for mapping rot on a single-tree level, and although the results have been relatively poor, some potential has been shown. This study applied area-based approaches to predict rot occurrence, rot severity, and rot volume , at an area level. Ground reference data were collected from harvester operations in 2019–2021. Predictor variables were calculated from multi-temporal remotely sensed data together with environmental variables. Response variables from the harvester data and predictor variables from remotely sensed data were aggregated to grid cells and to forest stands. Random Forest models were built for the different combinations of response variables and predictor subsets, and validated with both random- and spatial cross-validation. The results showed that it was not possible to estimate rot occurrence and rot severity with the applied modeling procedure (pR²: 0.00–0.16), without spatially close training data. The better performance of rot volume models (pR²: 0.12–0.37) was mainly due to the correlation between timber volume and rot volume.

KEYWORDS:

• Remote sensing
• multi-temporal remote sensing
• heartwood rot
• Norway spruce

Introduction

Deciding whether to harvest timber from a forest or to allow the trees to continue to grow is one of the most important strategic decisions a forest manager must make. This decision typically relies on information regarding timber quality, volume increment, harvesting costs, and timber prices. It usually involves optimizing the net present value of different harvesting times. The optimal time for harvesting can be influenced by many factors, of which natural disturbances can have a significant effect. By incorporating the risk of disturbances in the decision process, the forest owner can adapt the forest management, especially the decision of when to harvest (Susaeta & Gong, 2019). Without information about these risks, erroneous decisions are made (Möykkynen & Pukkala, 2010), potentially leading to monetary losses for the forest owner.

Due to the effects of global warming, the probability of natural disturbances will likely increase particularly in coniferous forests in the boreal zone (Seidl et al., 2017). Among such disturbances are fungal decay pathogens, which can cause rot in living trees, resulting in reduced growth and decreased wood production. Studies conducted in Norway and Sweden found that 14–25% of Norway spruce (Picea abies (L.) Karst.) trees harvested in final felling were affected by rot (Huse et al., 1994; Noordermeer et al., in press; Stenlid & Wästerlund, 1986; Thor et al., 2005). The national study of Huse et al. (1994), in Norway, showed that most of the rot in Norway spruce, 71%, was caused by the fungus Heterobasidion parviporum Niemelä & Korhonen. The fungi spread through the air by basidiospores that can infect exposed fresh wood tissue on stumps or wounds, or by direct root contact (Garbelotto & Gonthier, 2013). It usually affects the lower and most valuable part of the tree stem, and the occurrence of rot in commercial timber substantially reduces its value. Economic losses in Norway alone are reported to be €18.5 million annually (Noordermeer et al., in press).

Models for the prediction of rot occurrence based on empirical data have been developed in previous studies in Sweden (Thor et al., 2005), Finland (Mattila & Nuutinen, 2007), and Norway (Hylen & Granhus, 2018). These models predict the probability of rot at tree or stand level with explanatory variables describing climatic and edaphic conditions, forest structure, and development stages. The study of Hylen and Granhus (2018) identified diameter at breast height (dbh), age, altitude, growing season temperature sum, and vegetation type as significant explanatory variables for the number of stems with rot observed at breast height. This model is, however, developed for use at national and regional scales, and its performance at an operational management level is currently unknown.

Wall-to-wall remotely sensed information has the potential to enable the mapping and monitoring of rot at an operational level (Lausch et al., 2016). In Norway spruce, rot has been shown to lead to defoliation (Žemaitis & Žemaitė, 2018), and close-range remote sensing has been used to detect differences in the crown condition in Scots pine affected by rot (Pitkänen et al., 2021). Recent research has also shown promising results in using remotely sensed data to detect attacks by bark beetles (Ips spp.) in Norway spruce (Abdullah et al., 2019; Bárta et al., 2021; Huo et al., 2021). These studies have shown that it is possible to detect differences in the spectral signatures of healthy and stressed trees, even in the early stages of an attack using satellite imagery.

Recent studies in Norway have used various sources of remotely sensed information coupled with ground reference data collected from cut-to-length harvesters to model and map rot on single-tree level (Allen et al., 2022a; Allen et al., 2022b; Dalponte et al., 2022a; Dalponte et al., 2022b) and management unit level (Räty et al., 2021). The work on single-tree level identified the classification of rot as a challenging task, with best results in terms of the kappa coefficient of 0.19–0.22 (Dalponte et al., 2022a) and 0.19–0.27 (Allen et al., 2022a). The relatively poor results were explained by the complexity of the problem and the fact that the visual signs of heartwood rot are not detectable to the human eye. Although there are differences in the spectral signatures of healthy and infected trees, there is not enough information to separate the two classes (Dalponte et al., 2022a). Furthermore, the single-tree approach has the drawback of not detecting all trees in the remotely sensed data (Dalponte et al., 2015; Kandare et al., 2017).

Hence, the area-based approach to forest inventory and modeling, whereby the response is modeled on a fixed area, is preferred to avoid systematic errors caused by undetected trees (Coomes et al., 2017). Räty et al. (2021) applied the area-based approach to model the volume of butt rot at the grid cell level and aggregated the predictions to the management unit level. Explanatory information included both environmental and remotely sensed variables. Additionally, information from the harvester was included in the model to analyze this information as a theoretical case. This included variables such as the harvested timber volume, number of harvested stems, quadratic mean dbh of harvested stems of Norway spruce, the difference between the 90th and 10th percentile of the dbh distribution, and the proportion of harvested timber volume of Norway spruce. The study identified the variables describing the maturity of the forest, such as tree height, timber volume, and dbh, as the most important predictor variables.

Encouraged by the promising results in using multi-temporal remotely sensed data from nanosatellites (Allen et al., 2022a), the area-based approach used by Räty et al. (2021) was in the present study expanded by incorporating multi-temporal remotely sensed data from other high-resolution optical sensors. The recent results on the mapping of beetle attacks using data from the Sentinel satellites (Abdullah et al., 2019; Bárta et al., 2021; Huo et al., 2021) indicated the importance of the data for the detection of stressed trees. The analysis also included bi-temporal airborne laser scanning (ALS) due to the possible difference in crown density between healthy and infected trees (Pitkänen et al., 2021; Žemaitis & Žemaitė, 2018), and the ability of ALS to measure tree crown density (Popescu, 2007). Bi-temporal ALS can also indicate site index and to an extent soil properties (Noordermeer et al., 2018). Furthermore, rot volume as a response variable was supplemented with proportional responses of rot.

This study aimed to evaluate area-based approaches for estimating rot in Norway spruce trees using remotely sensed data. The specific objectives were to: (i) evaluate the use of multi-temporal remotely sensed data for modeling rot using the area-based approach, (ii) compare three modeling approaches; grid, stand, and grid-to-stand, (iii) assess the effect of using spatially close data, and (iv) evaluate the use of proportional model responses, percentage of trees with heartwood rot (occurrence) and percentage of rot volume to the total timber volume (severity), in addition to rot volume.

Materials and methods

Study area

The study area of 986 square kilometers was located in Etnedal, Nord-Aurdal, Sør-Aurdal, and Nordre Land municipalities in southeastern Norway (60°94”N,09°56”E, Figure 1). The area is characterized by steep terrain with an elevation of the harvested sites ranging from approximately 200 to 980 m above sea level. Norway spruce is the naturally dominant tree species in the area, especially on productive sites (i.e. sites with a high site index). Scots pine (Pinus sylvestris L.) and deciduous species, mainly birch (Betula pubescens Ehrh.), can be found in less productive sites or interspersed in spruce-dominated sites.

Figure 1. Study area and harvesting operations.

Harvester data

Ground reference data were collected from 46 harvesting operations using a single-grip Komatsu 931XC harvester equipped with a 230 H crane and a C144 harvester head, in the period August 2019 – November 2021. A total of 131,375 trees were harvested, 89% of which were Norway spruce, 6% Scots pine, and 5% deciduous, mainly birch. The harvester was fitted with a Septentrio AsteRx-U real-time kinematic Global Navigation Satellite System (GNSS) consisting of two antennas that logged positions and rotations at a rate of one second in the National Marine Electronics Association (NMEA) format. The locations of harvested trees were recorded with a planimetric accuracy of around 1 m using the positions and rotations obtained from the GNSS and the crane extension measured by sensor hardware (Noordermeer et al., 2021). The harvester production report files contained data on geographical coordinates, species and stem profile of harvested trees, and dimensions and assortments of individual logs.

In addition to their regular work related to harvesting and processing trees, the harvester operator visually assessed and recorded the presence or absence of heartwood rot on the butt ends of processed spruce logs. The harvester operator did not, however, assess rot on non-spruce logs.

Commercial height, total tree height, timber volume, and rot volume were computed for each tree. Commercial height was calculated as the height along the stem profile where the tree’s last log was bucked and total tree height was predicted using a taper model. For each log, the log volume was calculated from the stem profile recorded by the harvester. Timber volume was calculated using empirical single-tree volume functions that relied on total tree height and dbh as predictors (Braastad, 1966; Brantseg, 1967; Vestjordet, 1967), and rot volume was calculated in the following procedure: (i) logs for which the subsequent log (i.e. the next log higher up the tree stem) was also classified as having rot, the total volume of the log, calculated from the stem profile, was registered as rot and (ii) logs for which the subsequent log was classified as having no rot, the rot volume was calculated from the stem profile of the lower half of the log length.

A summary of the single-tree harvester data is shown in Table 1.

Table 1. Summary statistics (mean and range) of harvested trees.

Species	n	dbh (mm)	Height (m)	Rot occurrence (%)	Rot severity (%)	Rot volume (m^3)	Timber volume (m^3)
Spruce	116737	197.89 (46–661)	14.96 (2–53)	16.36	4.73 (0–181.68)	0.02 (0–1.61)	0.31 (0.31–6.75)
Pine	8160	269.54 (53–584)	18.81 (5–50)	NA*	NA*	NA*	0.60 (0.01–3.6)
Deciduous	6478	167.34 (50–635)	14.78 (4–53)	NA*	NA*	NA*	0.20 (0.01–4.63)

*Rot was only recorded for spruce.

Aggregation of harvester data

When conducting a forest inventory following the area-based approach, the inventory area is typically tessellated into grid cells of the exact size as the field sample plots. The use of harvester data as ground reference data provides substantial flexibility regarding plot and grid cell sizes, however, research results indicate that sizes between 200 and 400 m² are most appropriate (Maltamo et al., 2019). Accordingly, and conforming to operational practices in Norwegian forest inventories, we used a grid cell size of 250 m² when aggregating the harvester data. Thus, the study area was tessellated into grid cells and clipped by the boundary extent of harvested areas, creating a total of 10,304 cells. Cells with an area <100 m² or with a spruce proportion < 70% were removed from the dataset. Furthermore, the harvester data were aggregated to the stands defined in the forest management plan. The stands were clipped by the boundary extent of harvested areas. A total of 562 stands were assigned harvester data. However, if the cells with aggregated harvester data were in sum less than 0.2 ha at stand level, the cells were discarded. After applying these filtering criteria, 3715 cells and 109 stands were used in the analysis.

For each grid cell, timber volume and three measures of rot were calculated: rot occurrence, rot severity, and rot volume. Rot occurrence was calculated as the percentage of trees with heartwood rot, rot severity as the percentage of rot volume relative to the total timber volume, and rot volume as the timber volume (m³/ha) with heartwood rot. Formulas used to obtain rot occurrence (eq. 4), rot severity (eq. 3), rot volume (eq. 2), and timber volume (eq. 1) on grid and stand level are shown below where $Rotvolum{e}_j$ and $Timbervolum{e}_j$ are summed rot volume and timber volume for tree $j$,$Presenceofro{t}_j$ is a binary identifier for rot (0 for trees with no rot, 1 for trees with rot) in tree $j$ and $m$ refers to the number of trees in the cell.

(1) $Timbervolume={\sum }_{j=1}^m\left(Timbervolum{e}_j\left(m^{3}h{a}^{-1}\right)\right)$(1)

(2) $Rotvolume={\sum }_{j=1}^m\left(Rotvolum{e}_j\left(m^{3}h{a}^{-1}\right)\right)$(2)

(3) $Rotseverity=\frac{{\sum }_{j=1}^m\left(Rotvolum{e}_j\left(m^{3}h{a}^{-1}\right)\right)}{{\sum }_{j=1}^m\left(Timbervolum{e}_j\left(m^{3}h{a}^{-1}\right)\right)}\cdot 100\mathrm{\%}$(3)

(4) $Rotoccurrence=\frac{{\sum }_{j=1}^m\left(Presenceofro{t}_j\right)}{m}\cdot 100\mathrm{\%}$(4)

To calculate the dominant height (H_dom) for each grid cell, the number of trees to be used in the calculation (N_dom) was determined by dividing the area of the cell in square meters by 100, conforming to the criteria of the 100 largest trees per ha. N_dom was then rounded up to the nearest integer and this number was used to select the n trees with the largest dbh. H_dom was then calculated as the weighted mean of these selected trees, where all but the last tree were assigned a weight of 1, and the last tree was assigned the remainder obtained by subtracting N_dom from the rounded-down nearest integer of N_dom as its weight.

Table 2 shows the statistics of the filtered and aggregated harvester data. The means of properties were quite similar for the different aggregation methods, while the variation was larger for the grid data. The correlation between the rot variables (rot occurrence, rot severity, rot volume) and timber volume was assessed (Figure 2). Rot occurrence and rot severity both had a weak correlation with timber volume. For the grid and stand data, the correlation coefficients were respectively 0.2 and 0.3 for rot occurrence and 0.1 and 0 for rot severity. Rot volume had a moderate correlation with timber volume of 0.5 and 0.6 for the grid and stand data, respectively.

Figure 2. Correlations of rot variables with timber volume.

Table 2. Summary of aggregated harvester data at grid and stand level.

Aggregation	Property	Mean	Standard deviation	Minimum	Maximum
Grid	Area (m^2)	210.2	48.3	100.1	250.0
	n trees (n/ha)	777.5	340.6	40.0	2520.0
	Rot occurrence (%)	16.5	15.3	0.0	100
	Rot severity (%)	6.5	7.1	0.0	50.8
	Rot volume (m^3/ha)	17.0	23.5	0.0	294.3
	Timber volume (m^3/ha)	246.1	151.9	7.7	1404.6
Stand	Area (ha)	1.0	1.1	0.2	7.0
	n trees (m^3/ha)	805.7	246.2	378.7	1672.0
	Rot occurrence (%)	14.3	6.6	3.4	33.4
	Rot severity (%)	6.0	2.9	1.4	15.1
	Rot volume (m^3/ha)	15.1	10.5	2.1	51.8
	Timber volume (m^3/ha)	252.6	124.8	51.5	804.6

Remotely sensed and environmental data

Optical data

Optical satellite images were collected for 2017 and 2018 from Sentinel-2 and the Planet Dove constellation satellites (Planet Team, 2017). Sentinel-2 images have 10 bands with a spatial resolution of 10 m or 20 m while Dove images are of a higher resolution than Sentinel-2 at 3 m, but only contain four bands.

Monthly image composites were generated for Sentinel-2 and Dove. Sentinel-2 image composites were generated in Google Earth Engine, taking the median value of each pixel, considering only shadow-free and cloud-free ones. Analytic Ortho Dove images were collected, corrected to surface reflectance and harmonized to Sentinel-2 to reduce scene-to-scene variability (Kington & Collison, 2022) using Planet’s Data API. Composites for Dove were generated using the R packages terra (Hijmans, 2022) and RStoolbox (Leutner et al., 2022). Images were co-registered before generating composites by taking the median value of each pixel not affected by cloud, haze, shadow, and/or snow. Due to high cloud cover, some months contained cells with missing data. Months in which more than 10% of the cells had missing data were removed. Specifically, this included the months of April, July, September, and October of 2017 and April, September, and October of 2018.

For each of the monthly image composites from Sentinel-2 and Dove, vegetation indices (VIs) were computed. VIs were calculated using the RStoolbox R package (Leutner et al., 2022). In total, 28 VIs were computed from the Sentinel-2 composites and 16 from Dove (Table 3). VIs were computed for Sentinel-2 over seven months and Dove over eight months, resulting in a total of 196 and 128 computed variables, respectively.

Table 3. Overview of vegetation indices computed for Sentinel-2 (S2) and Planet Dove (Dove).

Vegetation index	Bands	Satellite	Reference
CLG	Red edge 3, green	S2	(A. A. Gitelson et al., 2003)
CLRE	Red edge 3, red edge 1	S2	(A. A. Gitelson et al., 2003)
CTVI	Red, nir	S2 and Dove	(Perry & Lautenschlager, 1984)
DVI	Red, nir	S2 and Dove	(Richardson & Wiegand, 1977)
GEMI	Red, nir	S2 and Dove	(Pinty & Verstraete, 1992)
GNDVI	Green, nir	S2 and Dove	(A. A. Gitelson et al., 1996)
KNDVI	Red, nir	S2 and Dove	(Camps-Valls et al., 2021)
MCARI	Green, red, red edge 1	S2	(Nagler et al., 2000)
MNDWI	Green, swir 2	S2	(Xu, 2006)
MSAVI	Red, nir	S2 and Dove	(Qi et al., 1994)
MSAVI2	Red, nir	S2 and Dove	(Qi et al., 1994)
MTCI	Red, red edge 1, red edge 2	S2	(Dash & Curran, 2004)
NBRI	Nir, swir 3	S2	(Lozano-Garcia et al., 1991)
NDREI1	Red edge 2, red edge 1	S2	(A. Gitelson & Merzlyak, 1994)
NDREI2	Red edge 3, red edge 1	S2	(Barnes et al., 2000)
NDVI	Red, nir	S2 and Dove	(Rouse et al., 1974)
NDWI	Green, nir	S2 and Dove	(McFeeters, 1996)
NDWI2	Nir, swir 2	S2	(Gao, 1996)
NRVI	Red, nir	S2 and Dove	(Baret & Guyot, 1991)
REIP	Red, red edge 1, red edge 2, red edge 3	S2	(Guyot & Baret, 1988)
RVI	Red, nir	S2 and Dove	(Birth & McVey, 1968)
SATVI	Red, swir 2, swir 3	S2	(Marsett et al., 2006)
SAVI	Red, nir	S2 and Dove	(Huete, 1988)
SLAVI	Red, nir, swir 2	S2	(Lymburner et al., 2000)
SR	Red, nir	S2 and Dove	(Birth & McVey, 1968)
TTVI	Red, nir	S2 and Dove	(Thiam, 1998)
TVI	Red, nir	S2 and Dove	(Deering, 1975)
WDVI	Red, nir	S2 and Dove	(Richardson & Wiegand, 1977)

Airborne lasscanning data

Two sets of ALS data were prepared for the study area. The first set (T1) was collected in 2007 with a point density of 0.7 points per m², while the second set (T2) consisted of data from several acquisitions covering different parts of the study area in the years 2013–2017. The point density in these acquisitions was either two or five points per m².

ALS data were processed with the lidR package (Roussel & Auty, 2022) in R. Point cloud normalization was performed using a triangulated irregular network (TIN) created from echoes classified as ground. For each echo, the height was normalized by subtracting the interpolated height of the TIN.

For each grid cell, canopy height variables including percentiles at 10% intervals (H10, H20, … , H90) were calculated from laser echoes above a threshold of 2 m above the ground. In addition, the mean height (Hmean), standard deviation (Hsd), coefficient of variation (Hcv), kurtosis (Hkurt), and skewness (Hskewness) were calculated for echoes above the threshold of 2 m. Canopy density variables were calculated by first dividing the range between the 95% height percentile and the 2 m threshold into ten vertical layers of equal height. The proportion of echoes above each layer to the total number of echoes was calculated, resulting in ten canopy density metrics (D0, D1, … , D9). ALS metrics were calculated for T1 and T2, as well as the differences between corresponding T1 and T2 metrics. Since the collection of ALS data for T2 was spread over several years (2013–2017) these differences were calculated as the relative change per year. The ALS metrics are hereafter referred to as ALS variables.

Environmental data

Forest management plans are typically revised every 10–15 years for the majority of the productive forest area in Norway. The plans provide stand-level information and are usually produced using a combination of remotely sensed information and field surveys. We collected two variables obtained from manual photo interpretation, i.e. site index and age, from the forest management plan covering the study area. Some grid cells intersected multiple stands, in which case the values of the stand with the largest overlap were selected.

A water table model referred to as the depth to water (DTW) is available at a national level (Hoffmann et al., 2022). The mean value in centimeters from the soil surface to the water table was extracted for each grid cell. Furthermore, ALS data from T1 were used to create a digital terrain model (DTM) with one-meter resolution using a TIN. The R package Terra (Hijmans, 2022) was then used to calculate slope, aspect, terrain ruggedness index (TRI), topographic position index (TPI), roughness, and flow direction from the DTM. In addition, the height above sea level was extracted from the DTM, and mean values were extracted for each cell.

Modeling and validation

The complete dataset contained all VIs computed from Sentinel-2 and Dove in the years 2017 and 2018 (324 variables), ALS variables from T2 and the ALS difference variables (48 variables), environmental variables (8 variables), and variables collected from the forest management plan (2 variables). The variables were divided into four different subsets of predictor variables: all variables (all), only ALS variables (als), only VIs variables (optical), and only environmental variables (env).

Three approaches to modeling were used for the aggregated data. In the grid and stand approach, predictor and response variables (rot occurrence, rot severity, rot volume, and timber volume) aggregated to the respective levels were used for training and validating models. The last approach was a grid-to-stand approach where the models were first trained at grid cell-level. CV predictions and observed values at grid cell-level were then averaged to the respective stands weighted by the area of the grid cells to create estimates at stand-level. These stand-level values were then used to validate the model.

A model was created for each combination of the four response variables and the four subsets of predictor variables, each model was fitted with and without outlier removal for grid data, totaling 32 models for grid-level and 16 models for stand-level. All models’ performance was assessed with random cross-validation (CV) and spatial CV. An overview of the data sources and modeling steps is shown in Figure 3.

Figure 3. Schematic diagram of the data sources, modeling, and validation steps.

Regression method

The random forest algorithm has become a popular and widely used algorithm in both classification and regression problems (Biau & Scornet, 2016). The algorithm presented by Breiman (2001) can handle datasets with a large number of variables, non-linearity, and at the same time having few parameters that require tuning. Several random forest implementations are available, of which we used one of the R package ranger (Wright & Ziegler, 2017) for computational efficiency.

Preliminary modeling using the R package tuneRanger (Probst et al., 2019) showed that hyperparameter tuning had little effect on model performance. Therefore, the hyperparameters were fixed to their default values: Number of trees to 500, minimum node size to five, and mtry to the largest integer value, not greater than or equal to the square root of the number of predictor variables. Impurity importance was used for estimating variable importance measures. Variable importance measures were averaged across the folds used in the validation and then scaled according to the greatest value on a scale from 0 to 1.

Outlier removal

To improve model robustness, outlier removal was attempted to remove grid cells where harvester data and ALS data did not match. Outlier removal was only applied to the grid dataset as the inaccuracy of harvester data would have little effect on the stand dataset. A linear model with H_dom as the response and the H90 ALS variable from T2 as the predictor was fitted to the harvester grid data. Observations with a Cook’s distance higher than 0.5 of the mean Cook’s distance were flagged as outliers. These observations were removed from the training data but not from the validation data. Modeling was performed with and without outlier removal.

Validation

Two validation methods were used for assessing model performance. The first method was random k-fold CV, by which observations were randomly assigned to folds. The second validation method was spatial clustering CV, by which spatially close samples were assigned to the same folds. Clusters based on the distances between observations were created with k-means clustering. The k in the k-fold CV was set to 10 for all datasets. For the spatial CV, the average distance between folds was 13.0 km for the grid approach and 11.4 km for the stand approach, with minimum and maximum distances of 1.4 km and 28.0 km for grid, and 0.5 km and 24.1 km for stand.

Machine learning models applied on spatially dependent datasets have been shown to overestimate performance when random CV is used compared to spatial CV (Meyer et al., 2018, 2019). By including both random and spatial CV we can assess the effects of spatially close training data.

The results were evaluated using root mean squared error (RMSE), relative root mean squared error (RMSE%), pseudo coefficient of determination (pR²) calculated as the squared correlation between observed and predicted values, and mean difference (MD). Additionally, a fuzzy set evaluation, as proposed by Gopal and Woodcock (1994), was used on a scale from five to one to complement the preceding metrics. The fuzzy set score from five to one gives a percentage-wise distribution of the deviation between the predicted value and the observed value, corresponding to ≤ 10%, ≤20%, ≤30%, ≤40%, and > 40%, respectively. This can be helpful in operational forest management by assessing the model performances in terms of deviation from the true value within a range, eg., 0–20%. For each response variable, a null-model was created. Null-models for rot severity and rot occurrence used the average value of the response as a constant estimate, while rot volume was estimated from true timber volume and average rot severity.

Formulas for RMSE, RMSE%, pR², and MD are shown in formulas 5–8 where $y_i$ and $\widehat{y_i}$ is the observed and predicted values of the response variable for cell $i$, $\underset{\hspace{0.17em}}{y}$ is the observed mean of the response variable and $n$ is the number of cells.

(5) $RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{n}}$(5)

(6) $RMSE\mathrm{\%}=\frac{RMSE}{\bar{y}}\ast 100$(6)

(7) $p{R}^2=cor{\left(y,\hat{y}\right)}^2$(7)

(8) $MD=\frac{{\sum }_{i=1}^n\left(\widehat{y_i}-y_i\right)}{n}$(8)

Results

Performance in terms of pR² are shown in Figure 4 for all models. Random CV led to higher performances compared to the spatial CV for all approaches. The drop in performance, from random- to spatial CV, was reduced using ALS variables as predictors for timber volume- and rot volume models.

Figure 4. Performance of models.

The effect of outlier removal depended on the validation method. Results using the random CV, showed reduced performance for every model evaluated. Applying the spatial CV, the effect of outlier removal was small. However, it improved the results for some rot volume and rot severity models.

Performance metrics for the best-performing models for each response and modeling approach are shown in Table 4. Of the rot responses, the rot volume models achieved the highest pR² values for all approaches with 0.12, 0.14, and 0.37 and RMSE-values of 22.4, 9.9, and 11.4 using spatial CV for grid-, stand- and grid-to-stand approach, respectively. The RMSEs for the grid approach were clearly larger than the other two, while the difference between stand and grid-to-stand was smaller. The performance of rot occurrence and rot severity models was low for all approaches. The pR² of the rot severity models for all approaches ranged from 0.00 to 0.01 with spatial CV. Rot occurrence models performed better for the grid and grid-to-stand approach with pR² of 0.06 and 0.16, and RMSE-values of 15.0 and 7.7, respectively, while the stand approach had pR² of 0.03 and RMSE of 6.6.

Table 4. Performance metrics for best-performing models.

								Fuzzy set evaluation
Response	Predictors	CV	RMSE	RMSE%	pR^2	MD	Outlier removal	5	4	3	2	1
Grid
Rot occurrence	all	Random	13.59	82.58	0.22	−0.58	F	0.09	0.20	0.28	0.72	0.63
Rot occurrence	als	Spatial	14.96	90.93	0.06	−0.86	F	0.09	0.18	0.26	0.74	0.65
Rot severity	all	Random	6.65	102.96	0.11	−0.33	F	0.07	0.15	0.23	0.77	0.69
Rot severity	all	Spatial	7.26	112.41	0.01	−1.12	F	0.08	0.15	0.22	0.78	0.71
Rot volume	all	Random	20.09	118.90	0.27	−0.99	F	0.07	0.15	0.22	0.78	0.71
Rot volume	als	Spatial	22.38	132.45	0.12	−0.93	T	0.06	0.12	0.18	0.82	0.75
Timber volume	all	Random	84.20	34.22	0.70	−1.66	F	0.30	0.55	0.72	0.28	0.19
Timber volume	als	Spatial	93.19	37.87	0.62	−0.11	F	0.26	0.48	0.65	0.35	0.24
Stand
Rot occurrence	all	Random	6.01	42.06	0.16	−0.14	F	0.13	0.37	0.54	0.46	0.37
Rot occurrence	als	Spatial	6.55	45.86	0.03	0.02	F	0.17	0.31	0.47	0.53	0.38
Rot severity	all	Random	2.67	44.52	0.17	−0.14	F	0.20	0.34	0.51	0.49	0.39
Rot severity	env	Spatial	3.07	51.10	0.00	−0.19	F	0.17	0.29	0.41	0.59	0.46
Rot volume	all	Random	9.23	61.02	0.23	−0.51	F	0.11	0.24	0.35	0.65	0.50
Rot volume	als	Spatial	9.86	65.22	0.14	−0.19	F	0.13	0.25	0.34	0.66	0.59
Timber volume	als	Random	62.51	24.74	0.75	1.54	F	0.39	0.65	0.83	0.17	0.07
Timber volume	als	Spatial	73.19	28.97	0.66	6.75	F	0.34	0.60	0.82	0.18	0.11
Grid-to-stand
Rot occurrence	all	Random	5.27	34.45	0.66	−1.28	F	0.34	0.52	0.65	0.35	0.25
Rot occurrence	als	Spatial	7.66	50.04	0.16	−1.91	F	0.16	0.31	0.46	0.54	0.37
Rot severity	env	Random	2.30	37.84	0.64	−0.36	F	0.31	0.50	0.64	0.36	0.28
Rot severity	optical	Spatial	4.08	67.00	0.01	−1.00	T	0.19	0.29	0.41	0.59	0.48
Rot volume	all	Random	7.43	46.70	0.70	−2.15	F	0.20	0.43	0.51	0.49	0.41
Rot volume	als	Spatial	11.37	71.44	0.37	−3.28	F	0.14	0.30	0.42	0.58	0.50
Timber volume	all	Random	40.40	15.94	0.90	−2.38	F	0.61	0.84	0.94	0.06	0.05
Timber volume	als	Spatial	56.83	22.42	0.80	1.09	F	0.41	0.65	0.85	0.15	0.04

The predicted versus observed values for spatial CV of the best-performing models are shown in Figure 5. Timber volume models performed well, although they underpredicted observations with large volumes. The predictions for rot models were centered around the mean, underestimating small values and overestimating large ones.

Figure 5. Observed versus predicted values of best-performing models for spatial CV.

Variable importance plots are shown in Figure 6. The best-performing models for rot occurrence and rot volume models are dependent on ALS variables. The variable Hsd and height percentiles are assigned high importance. The importance of ALS variables is lower for rot severity models and is dependent on the approach. The grid approach assigned high importance to environmental and optical variables, while the stand approach only selected environmental variables.

Figure 6. Variable importance of the best models for spatial CV.

The timber volume models showed substantially better performance than any of the rot models, with pR²-values of 0.62, 0.66, and 0.80 for grid-, stand- and grid-to-stand approaches, respectively, using spatial CV. The drop in performance between validation methods was small compared to any of the rot measures when ALS variables were included as predictors. Height percentiles and Hmean from the ALS variables scored high on variable importance for models.

The fuzzy set evaluation scores show that, using spatial CV, 65% of the timber volume predictions were within ± 30% of the true value on grid level and 85% at stand level using the grid-to-stand aggregation approach. For rot occurrence, rot severity, and rot volume the comparable values were 26%, 22%, and 18% at grid level, respectively, and 47%, 41%, and 34% at the stand level.

Comparing the results from the null-models in Table 5 with the best regression models (Table 4), the RMSEs were reduced by all regression models using random CV, with the exception of rot volume models. For rot occurrence and rot severity the RMSEs were reduced by 11% and 5% respectively on grid level and 8% and 9% on stand level. The difference in RMSEs was small on grid level (0%) but 14% lower on stand level for rot volume. For spatial CV, the RMSE were reduced with 2% for rot occurrence at grid cell level, while for rot severity and rot volume the RMSEs increased with 3% and 11%. At the stand level, there was no difference for rot occurrence while the RMSEs of rot severity and rot volume increased with 5% and 22% respectively.

Table 5. Null-model results. Average values of response used for prediction.

					Fuzzy set evaluation
Approach	Response	RMSE	RMSE%	pR2	5	4	3	2	1
Grid	Rot occurrence	15.31	93.05	NA*	0.08	0.17	0.24	0.76	0.66
	Rot severity	7.06	109.31	NA*	0.07	0.15	0.22	0.78	0.70
	Rot volume	20.12	119.09	0.28	0.07	0.15	0.22	0.78	0.70
Stand	Rot occurrence	6.54	45.76	NA*	0.16	0.39	0.50	0.50	0.38
	Rot severity	2.93	48.80	NA*	0.18	0.32	0.46	0.54	0.44
	Rot volume	8.07	53.33	0.41	0.18	0.32	0.46	0.54	0.44

*Constant estimate was used for rot occurrence and rot severity.

Discussion

By using an area-based approach combined with multi-temporal data from high-resolution optical sensors, we were not able to build a robust model for the estimation of rot. With spatial CV the relative responses performed poorly with pR² of rot severity ranging from 0.00 to 0.01 and rot occurrence ranging from 0.03 to 0.16. Rot volume models performed better than the relative responses with pR² ranging from 0.12 to 0.37.

The poor performance of the relative measures of rot may have been caused by a lack of predictive power of the used predictor variables. The classification of rot at single-tree level has been shown to be challenging when using multi-temporal imagery (Dalponte et al., 2022b) and hyperspectral imagery (Allen et al., 2022a; Allen et al., 2022b; Dalponte et al., 2022a). Žemaitis and Žemaitė (2018) found significant differences in crown defoliation between trees with and without rot damage, but noted that almost all of the morphological indicators were non-specific. Thus, despite the existence of detectable differences between trees with and without rot, other stressors may hinder our ability to detect these differences using remote sensing. The effect of such stressors may be augmented in the present study, which was conducted at the spatial scale of grid cells and stands as opposed to smaller individual trees in the mentioned studies.

Although rot volume was the best-performing response of rot, its correlation with timber volume may cause it to be misleading. In practice, estimating rot volume involves two parts: estimating timber volume and estimating rot severity (i.e. the extent to which the timber is damaged by rot). The poor performance of the rot severity models indicated that our rot volume models were mainly driven by volume estimation rather than the estimation of rot. In this study, both rot volume and rot occurrence models shared the high importance of height percentiles with timber volume models, further highlighting the link with timber volume. Contrary to Räty et al. (2021) however, our rot volume models did not outperform the null-models for rot volume.

Interestingly, ALS variables such as Hsd and Hkurt were assigned high importance for both rot occurrence, rot severity, and rot volume, but not for the timber volume models. Räty et al. (2021) also found high importance of Hvar for their rot volume models. Hsd can be associated with stand density, stand variation, or used for identifying thinning, although the low performance of spatial CV makes it likely to be a result of spatial overfit. Rot severity had the weakest correlation with timber volume and was assigned high importance to environmental, optical, and non-volume-related ALS variables.

Compared to rot occurrence or rot severity, rot volume is less useful for making decisions on forest management due to the correlation to timber volume. If the models cannot differentiate between stands with similar volumes, they also cannot assist in prioritizing. Comparing the null-models and the rot volume models shows that using the true volume and the average rot severity result in a pR² of 0.41 compared to a rot volume model which resulted in a pR² of 0.37. We therefore suggest that future studies should, in addition to rot volume, include rot occurrence or rot severity, as both are suitable alternatives that are useful in decision-making.

Dalponte et al. (2022b) found that classification accuracy increased when only individual tree crowns with an area greater than 9 m² were included instead of all individual tree crowns. This was in part explained by the more representative spectral information, avoiding a mixed spectral signature from multiple trees within a satellite pixel. Aggregating harvester data to larger units such as grid cells or stands could result in a loss of representativeness of the data, especially for larger units. Since the effect of rot on the vigor of a tree is quite low at the individual tree level, aggregating data to large units may obscure any response present at the tree level. When predicting rot at the stand level, factors that explain the risk of rot may be more important than detecting the presence of rot. This would explain the shift from predominantly optical variables to exclusively environmental variables when predicting rot severity using the grid- and stand approach, respectively. The stronger correlation with timber volume may obscure this effect for rot occurrence.

Spatial overfitting was a significant issue for all models estimating rot, as indicated by the drop in performance with spatial CV. Similarly, Räty et al. (2021) reported a drop in performance from a pR² of 0.60 to 0.37 for rot volume models when using spatial CV. They also found high variable importance measures for the coordinates of stands. However, we observed a substantial improvement in model performance for random CV of all rot responses when using the grid-to-stand approach, and the performance was better than the stand approach. This improvement could be attributed to the use of spatially close training data, which could potentially be implemented at the operational management level, as local data from harvesters with accurate positioning systems are becoming more widely available.

To our knowledge, the only comparative study is that of Räty et al. (2021), which used an approach similar to our stand approach to predict rot volume. They reported values of RMSE, RMSE%, and pR² of 15.6, 65.4, and 0.37, respectively, for their spatial CV. The results were similar to ours, and their increased RMSE and pR² may partially be explained by the larger spread in rot volume and rot severity. Based on their variable importance plots, it is unlikely that their model would be able to describe rot severity better, as important variables are related to the volume or position of stands. The study of Räty et al. (2021) contains observations covering a broader geographical area of Norway than the present study. The high importance of stand position might indicate that their model can describe differences in rot between regions but not on a local scale.

Although the study area in this study is quite large, it is limited in the range of environmental variables compared to the study of Hylen and Granhus (2018) who identified growing season temperature sum, altitude, and age as important predictor variables. In the present study, data was only collected from mature and harvested stands, and the models might not apply to younger or older forests.

The accuracies of the volume models, with values of RMSE of 38% and 22% at grid- and stand level, respectively, were found to be similar or somewhat smaller in comparison to earlier studies using ALS and harvester data for estimating timber volumes. For instance, Hauglin et al. (2018) reported a relative RMSE of 20% and 32% at high and medium productivity sites, respectively, using a grid cell size of 400 m². In comparison, Maltamo et al. (2019) reported a relative RMSE of 33% at a grid cell size of 200 m² and only 9% at stand level.

The results presented in this study may be partially affected by errors in the harvester data, as discussed by Hauglin et al. (2018). The mentioned study noted that errors in tree positions could result in co-registration errors between harvester and ALS data. This impact may be expected to be greater at grid cell level due to the smaller area and number of trees compared to the stand level. Trees that were not harvested and therefore not registered in the data, such as retention trees, could lead to an underestimation of the actual volume. Additionally, some trees may be too large for the harvester head to process, and therefore manually felled and thus missing in the harvester data. However, the small effect of the outlier removal indicated that errors in the harvester data had little impact on the models. Moreover, the point density of the ALS acquisitions varied between two and five echoes per m², and the time of acquisition ranged from two to eight years before harvest. Previous studies have shown that the use of different sensors and pulse densities could significantly affect the height, density, and intensity variables calculated from ALS data, however, the impact on the precision forest attributes, such as volume and tree height has been shown to be minor (Næsset, 2009; Ørka et al., 2010). While Noordermeer et al. (2022) found no significant issues when using data from multiple ALS acquisitions to model timber volume, it remains unclear how potential sensor effects could impact the utility of ALS data for modeling rot properties. However, using multiple ALS acquisitions could add additional noise to the variables used for modeling rot properties, and a more uniform dataset may therefore provide better results. Furthermore, an eight-year time difference could have significant effects due to changes in both timber volume and rot.

Constellations of nanosatellites, such as Planet Dove provide imagery of high spatial and temporal resolutions but suffers in geometric and radiometric lower quality compared to Sentinel-2 (Frazier & Hemingway, 2021). Radiometric inconsistencies between images, geometric errors, and unusable data missed by the provided mask, e.g. undetected clouds, will all introduce noise in the data. The applied harmonization tool, the co-registration of images, and the aggregation of pixels to a higher level (grid or stand) should reduce noise introduced from these sources. Lower image quality, wide overlapping bands of the Dove Classic sensor (Collison et al., 2021), and the lack of bands can explain the weaker performance of Dove imagery compared to Sentinel-2. Using newer generations of Planet’s Dove satellites and improved pre-processing methods of imagery could be applied in future studies to improve the performance.

Forest management activities can have large effects on the infection risk and spread of H. parviporum (Garbelotto & Gonthier, 2013; Pukkala et al., 2005), and the fungus can remain active in infected stumps for decades (Greig & Pratt, 1976). Thus, records of management history, or techniques to derive previous activities, could be a viable route for future studies to improve the performance of rot models.

Conclusions

In this study, we evaluated area-based methods for estimating rot in Norway spruce trees using remotely sensed data. Models for rot volume, rot severity, rot occurrence, and timber volume were developed using stem rot data collected by a cut-to-length harvester and multi-temporal remotely sensed data from the Sentinel-2 and Planet Dove constellation satellites and ALS. The models’ performance for timber volume was consistent with previous studies. However, for rot responses, the models showed relatively poor performance, with RMSEs of 91%, 112%, and 119% at the grid cell level, for rot occurrence, rot severity, and rot volume, respectively. At the stand level, the best results using spatial CV were 46%, 51%, and 65% RMSE. The only comparative study known to us by Räty et al. (2021) reported an RMSE of 65% for rot volume using a similar approach.

The variable importance analysis showed that ALS was the most important source of information for both rot occurrence and volume. For rot severity at the grid cell level, the models also included optical data from the two satellite sources and environmental information. At the stand level, the models were dominated by environmental variables.

Including spatially close data in training improved the performance of models. The results showed reduced accuracy when predicting new areas and improved when observations in nearby stands or within stands were allowed in the training data.

Compared to estimates based solely on the average values of rot, the models only slightly improved the estimates, and only for rot volume. This effect can be explained by the positive correlation between timber volume and rot volume, and the fact that ALS data has proven to be a valuable data source for timber volume estimation. For responses less correlated to timber volume, such as rot severity and occurrence, the results indicated that the models were less useful.

Author contributions

Conceptualization: E.H., J.W., M.D., T.G., H.O.; Methodology: J.W., M.D., T.G., H.O.; Formal analysis: E.H., J.W.; Investigation: E.H., J.W.; Data Curation: J.W., M.D., T.G., L.N., H.O.; Writing – original draft: E.H., J.W.; Writing – review & editing: E.H., J.W., M.D., T.G., L.N., H.O.; Visualization: J.W.; Supervision: T.G., H.O.; Project administration: T.G., H.O.; Funding acquisition: T.G., H.O.

Acknowledgments

We would like to thank Valdres Skog AS for collecting and sharing the harvester data and Planet Labs Inc. for providing the images used in this study under the “Education and Research Program”. Furthermore, we are grateful for the constructive and valuable comments from three anonymous reviewers on previous versions of this study.

Data availability statement

The data that support the findings of this study are owned by Valdres Skog AS. Restrictions apply to the availability of these data, which were used under license for this study. Data are available upon reasonable request from the corresponding author, [H.O.Ø.], with the permission of Valdres Skog AS.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported by the Research Council of Norway under the project PRECISION (project no. 281140).