Remote-sensing support for the Estonian National Forest Inventory, facilitating the construction of maps for forest height, standing-wood volume, and tree species composition

The main forest-forming tree species in Estonia are European aspen (Populus tremula L.), silver birch (Betula pendula Roth), downy birch (B. pubescens Ehrh.), black alder (Alnus glutinosa (L.) Gaertn.), grey alder (A. incana (L.) Moench), Norway spruce (Picea abies (L.) Karst.), Scots pine (Pinus sylvestris L.), and common ash (Fraxinus excelsior L.). Depending on soil fertility and forest management, these species occur either in a variety of mixtures or as pure stands. Norway spruce forms a mid-storey in many stands growing on fertile soils. Forest land (23,308 km²) accounts for more than half (53.6%) of the Estonian terrestrial territory (Forest, 2020). About 47% of the forest land is owned by the state. The rest belongs largely to private owners, with a small share owned by municipalities and other public institutions. According to the forest register, more than 31% of forest land is drained to improve tree growth conditions. On 24.6% of Estonian forest land, use is restricted under a variety of nature-protection regulations (Forest, 2020). Wood is harvested through thinning and through shelter-wood cuttings, and also through the creation of small clearfellings (to a maximum area of 7 ha) with retention trees (Forestman, 2007; Forestact, 2016). More than 50% of the forest land use is classified as retention forestry (Gustafsson et al., 2012). Since the year 2000, the annual regeneration fellings and annual maintenance fellings have been in the respective area ranges of 17.1–38.8 thousand ha and 24.6–42.9 thousand ha (Forest, 2020). The average size of a forest stand is 1.25 ha. Forests are seminatural, with all the above-mentioned tree species regenerating naturally, and with Norway spruce and Scots pine occurring both under cultivation and with natural regeneration. The share of forest land has been increasing since the end of World War II, due to the abandonment of agricultural land. Leading the process of abandonment are those land parcels which are smaller and more remote from roads (Mandel et al., 2019).

The first Estonian-wide NFI with a sample plot layout similar to the current NFI was launched in 1999 (Kohava, 2000). With a rather modest budget and equipment, the NFI was able to give a fairly accurate assessment of the forest area, resources, and cutting volume. The main initial objective of the NFI was the estimation of major characteristics of forest stands. At present, however, the NFI has a wider scope, reporting also such data as the distribution of land by land-use categories and the afforestation and growing stock of non-forest land. The smallest target unit for reporting is the individual county.

Methodologically, the NFI is designed as an annual research effort, which, using optimal methods, must ensure continuous updating of its assembled information, including the forest database. Since 2014, sampling has been conducted with increased frequency, on a network of sample plots covering the entire country (Figure 1). Under this procedure, in each given year approximately 375 clusters (i.e. 20%) from the entire ensemble of sample plots, is measured (Figure 2). This ensures re-measurement of permanent plots once in every 5 years. Point estimates of parameters are calculated with data from the sample plots, as a basis for inferences to the entire population.

Figure 1

Figure 2

The Estonian NFI is designed as a systematic sample without pre-stratification. The sampling grid is established to meet accuracy requirements at the national level. About 5,500 sample plots are measured per year, with the sampling intensity the same throughout the country. The sample (cluster) distribution is based on the national 5×5 km quadrangle grid, determined by the Estonian base-map coordinates system (EPSG:3301).

The observation unit is the individual field plot, determined by its centre coordinates. The method of sampling with partial replacement is used. The sample plot area is subdivided if its area overlaps with different types of forest stands or parcels with different land cover. To enhance the efficiency of the survey, sample plots are concentrated into clusters, defined as 800×800 m squares (Figure 3). Some sample-plot clusters are deemed permanent, others temporary. The radius of a given sample plot depends on the variables selected for assessment, and additionally on variable values (with, e.g., smaller sample plots used for lower or regrowth layers). The radii of the principal sample plots are 10 m and 7 m. For land-use category determination, plots of other radii are taken, although without modification in the scheme for selecting centre coordinates. For the construction of models for stand height and wood volume, we used the sample plots in which a tree layer was present and which were located entirely inside a forest stand.

Figure 3

The remote-sensing support system of the Estonian NFI is based on data from the European Union Copernicus programme, from the NASA/USGS Landsat programme, and from the airborne photography and laser scanning programme of the Estonian Land Board. The medium spatial resolution (10–30 m) multispectral images from Sentinel-2 MSI (ESA, 2015) and Land-sat-8 OLI (USGS, 2019) sensors and SAR images from Sentinel-1 (ESA, 2020) are used for the prediction of tree species composition and for the detection of changes. The data are available from the regional data centre ESTHub (2016).

Aerial photography and airborne laser scanning are conducted by the Estonian Land Board under a repetition schema that provides, for any given location, either summer or springtime data in each second year, and additionally produces measurements from each similar growth season after every four years (Figure 4) (Maa-amet, 2020). The point density of the archived ALS data ranges from 0.15 m⁻² to 2 m⁻². The ALS pulse footprint diameter at canopy level is about 0.5 m, with a scanning angle that does not exceed 30° from nadir. Data are distributed according to a 1 km² map sheet system (Maa-amet, 2019). Before the year 2017, the Estonian Land Board used the Leica ALS50-II airborne scanning system. From 2017 onward, measurements have instead been taken with the Riegl VQ-1560i. In the Estonian NFI, orthophotos (in the RGB + NIR bands) are currently used only for visual interpretation, in the estimation of land cover type during the preparation of fieldwork. In the remote-sensing module offered in this paper for the Estonian NFI, ALS data are used for the prediction of forest height and wood volume. Ancillary data sources for our offered module are a 1:10,000 base map and soil map of Estonia, the FMI database with its stand-level forest-management inventory data, and a digital terrain module (DTM) provided by the Estonian Land Board.

Figure 4

The first two variables to be predicted are basal-area-weighted forest height (Lorey's mean height) H and wood volume M. The prediction of H and M is based on ALS data. The data processing system is designed in such a way that the user can develop and apply any models that take as their inputs ALS point-cloud metrics in combination with ancillary variables. However, for Estonian forests a simple linear regression model can be used with sufficient predictive power (Lang et al., 2012) for forest height: (1)HALS=a+b⋅HPx+ɛ.{H_{ALS}} = a + b \cdot {H_{Px}} + \epsilon . Here a and b are model parameters, H_Px, is one of the upper percentiles of the point-cloud height distribution, and ɛ is the model residual error. In the calculation of point-cloud height distribution percentiles, points near the ground are excluded. Wood volume M can be predicted with a model that is based on a relationship involving basal area G, forest height H, and stand form factor F (M = G·H·F): (2)MALS=(a⋅HP80b+c⋅HP25)⋅Kd+ɛ,{M_{ALS}} = \left( {a \cdot H_{P80}^b + c \cdot {H_{P25}}} \right) \cdot {K^d} + \epsilon , where a, b, c, and d are the model parameters; H_P80 is the 80th percentile; H_P25 is the lower quartile of the point-cloud height distribution; and K = C_above/C is a proxy for canopy cover calculated using the ALS pulse return count C_above above a height threshold (1–2 m from the ground) and total count of returns C per target pixel or sample plot (Lang et al., 2012). In calculating the height percentiles of point clouds, returns from near the ground are excluded, and K is calculated using all returns recorded for each emitted pulse. The FUSION toolbox (McGaughey, 2018) is used for ALS data processing in this study, with ALS metrics calculated for 10 m target pixels. Point clouds for sample plots were extracted using a circle with a radius of 15 m.

The parameters for height and wood volume models are estimated using NFI sample-plot field measurement data for each particular flight campaign, with each campaign covering roughly one-quarter of Estonia's area (Figure 4). No correction is currently made for springtime phenology and the consequent increased gap fraction in the canopy, and tree species are not currently distinguished for H and M models. The current species-independent modelling also ignores differences in the interaction of laser pulse with tree canopy which may arise from variations in scanner settings and from variations in the summertime structure of the canopy. The species-independent model was selected because it was considered that the number of field measurements would become insufficient if the set of field measurements were partitioned into subsets according to dominant species. In the current study, the wood volume corresponds to the upper canopy layer.

Prediction of tree species composition for target units with a size of 10–30 m is a challenging task, since most forests in Estonia are of mixed composition, since the spectral signatures of broadleaf deciduous trees are quite similar to each other, and since variables that are related to stand age and structure have a substantial influence on forest spectral signature (Nilson & Peterson, 1994) as calculated from Land-sat-8 OLI and Sentinel-2 MSI images. The proposed solution for the Estonian NFI is based on machine learning, using a large number of samples from the FMI database and data from a 1:10,000 soil map. Lang et al. (2018) found the random forest algorithm implementation in the GRASS GIS to be suitable for the construction of tree-species composition maps.

The machine-learning procedure starts with the selection of a potential training set of stands from the FMI database. The next step after the cloud masking of image data is the removal of outliers, based on relationships of spectral radiance with forest age and wood volume in red, near infrared, and shortwave infrared bands (Lang et al., 2016). Additional pairwise comparison of images is applied to remove stands with detectable disturbances that may have occurred during the period covered by the selection of satellite images. Experience has shown that this winnowing procedure yields a final set of about 100,000 samples, appropriate for machine learning for the whole of Estonia. The target area is divided into 8 overlapping subareas, to account for regional characteristics of forests in Estonia and to exclude training samples which have been measured over a long distance (Lang et al., 2018). Predictions for each subarea are made using either single images or combinations of images that cover more than half of the subarea, with a maximum of 3 images allowed in any one combination.

The model fitting on empirical data starts with the selection of informative features. The next step is model hyper-parameter optimization. Finally, a prediction map is constructed. Although the variable being predicted is the dominant species code, additional useful information is obtained from the entire vector of class probabilities (Lang et al., 2018) for each target pixel. Finally, the predicted class probabilities (where each class corresponds either to one of the possible dominant forest species or to a classification as non-forest) are averaged for each pixel, and dominant species information is extracted according to the estimated proportions.

For the construction of a tree species composition map, we used Landsat-8 OLI and Sentinel-2 images from April 2018 to September 2018. The images were resampled to 25 m spatial resolution and converted to the EPSG:3301 coordinate system. The methods for image processing and data analysis resemble those recently described by Lang et al. (2018).

The data processing system is implemented with existing software, but with the addition of a calling script written in Python. Each task (e.g. the prediction of forest height, of wood volume, or of tree species composition; or the detection of changes; or the appraisal of predictions) is structured as a sequence of specific steps (e.g. a preparation phase, the extraction of features from remote sensing data, the fitting of a model, the application of the selected model, and validation). The user edits the script template of a specific step, inserts proper parameter values, and executes the script. Each processing execution is assigned a unique ID, which is then used in a file-storage folder structure. The scripts produce an execution history log as a directed acyclic graph (DAG). Version control, the generated data-folder structure, and the generated DAG jointly ensure traceability and repeatability of the data-processing steps. In considering alternatives to this scheme, we evaluated some formal workflow tools (Airflow, Luigi, Jenkins), but found that they failed to add sufficient value. The principal software tools used in our scheme are QGIS, GRASS GIS, Python, GDAL, ESA SNAP, the Orfeo Toolbox (OTB), LAStools, FUSION, R, RStudio, and Git.

During execution, all intermediate results are stored in the local computer. In addition, the intermediate and final results are copied to a backup file server. Intermediate results are stored for 2 years (this being the longest interval between successive measurement updates). Final results are stored for some longer period, as decided by the system administrator. The hardware configuration comprises an 8-core workstation CPU, 24 GB of RAM, and 6 TB of workstation disk storage, with additionally 2.3 TB of backup storage. All hardware is virtualized.

For characterization of forest height and wood volume models, the mean residual squared error (3)RSE=∑(Yˆ-Y)2NRSE = \sqrt {{{\sum {{\left( {\hat Y - Y} \right)}^2}} \over N}} and mean residual error (4)MRE=∑(Yˆ-Y)NMRE = {{\sum \left( {\hat Y - Y} \right)} \over N} were calculated, with N the number of observations, Ŷ the predicted value for each sample plot in the sum, and Y the corresponding measured value.

For the assessment of forest height and wood volume predictions, we used a set of forest stands from the overlap area of the NE and NW blocks, and additionally a set of forest stands from the overlap area of the SE and SW blocks (Figure 4). We applied the following selection criteria: the stand polygon area was required to lie in the range of 1.2–8.0 ha; the count of 10 m pixels within the stand polygon was required to be > 100; the number of pixels within the stand polygon without value (i.e. with NoData label) was required to be < 10 in the corresponding comparison maps; and the proportion of evergreen tree species in the stand was required to be either ≥ 75% or ≤ 25%. The last of these filters (the disjunctive criterion regarding the proportion of evergreens) was applied to select contrasting sets of stands according to tree species. On the overlap area of NW and NE blocks were 1,157 evergreen and 605 deciduous stands. On the overlap area of SW and SE blocks were 644 evergreen and 1,275 deciduous stands. For each stand, the mean value of pixels was calculated from forest-height and wood-volume prediction maps.

The tree species composition was validated on NFI sample plots, because we used FMI data for the prediction model, and the NFI and FMI datasets are independent. Two selection criteria were used: it was required that the sample plot be described as forest in the NFI, and the probability of the plot being non-forest was required to be < 25% according to predicted pixel values in the tree species map. This pair of filters helped to remove disturbed areas and validation points in which the spectral signature is mainly influenced by objects other than forest trees. We analyzed the prediction for dominant tree species and the prediction for the proportion of evergreen coniferous trees in the species composition.

The NFI is based on contact measurements of trees on a set of small sample plots. For spatial units smaller than the individual county, feature variables from remote-sensing data can be used to construct maps of forest-inventory variables. As in the case of many other countries, Estonia has a national programme of springtime aerial photography and ALS for topographic mapping, with special flights additionally performed at higher altitudes in the summertime for forest-inventory purposes.

To construct maps of forest height based on the 3-dimensional ALS point clouds and NFI sample-plot measurements, it suffices to apply a linear model that takes as argument an upper percentile of the point-cloud height distribution. Upon taking a model independent of tree species and comparing forest height predictions based on springtime against predictions based on summer ALS data, we found a slight bias in the predicted forest height, dependent on the proportion of evergreen coniferous species in the subject forests. The difference proved more pronounced in deciduous forests, where greater heights were predicted when summer ALS data were used.

Lidar pulse returns are rarely recorded from tree stems during airborne measurements. The prediction of standing-wood volume therefore requires a model that accounts for tree height and forest density (defined as the number of trees per unit area). With detection of single trees not reliable from sparse point clouds with about 1 point per m², canopy cover can be taken as a proxy for forest density. However, the interaction of laser pulse with forest canopy, including the degree of pulse penetration to the ground, is substantially dependent on the amount of foliage, in addition to suffering dependence on flight parameters and scanner settings. There are two options for incorporating canopy cover information into the standing-wood volume model. The first option is to use a variable that predicts canopy cover. The second option is to calculate forest height metrics from point clouds which include near-to-ground points. The former option has been adopted for Estonian forests, while jurisdictions adopting the second option include Finland (Kotivuori et al., 2016). Arumäe and Lang (2018) analyzed ALS-based canopy cover (CC_ALS) estimation errors using digital hemispherical images and found the random variability of CC_ALS to be about 10–12% for sample plots 30 m in radius. In view of the random variability in CC_ALS, the combined effect of location errors and small sample-plot sizes, and the failure of tree stems to exert direct influence on the sparse ALS point clouds, the fitted for wood volume prediction are reliable for practical use for Estonian forests. The relative residual error of the models remained in the range of 30–38% of the mean measured value. A further decrease of random errors in predicted values, although possible, is difficult to achieve at pixel level (Arumäe & Lang, 2016).

A bias similar to the result encountered with evergreen coniferous-dominated forest height prediction was found in the standing-wood volume prediction, with the common model yielding an smaller by about 40 m³ ha⁻¹ for summer than for springtime ALS data. On the other hand, a systematically greater, with a discrepancy of the same order of magnitude, was found for forests dominated by deciduous trees for summer as compared against springtime ALS data. The problem could be addressed by dividing the empirical data into subsets according to the dominant tree species and fitting model parameters for each subset. However, in this approach the number of observations decreases and it becomes necessary to know the tree species distribution with high accuracy before the model can be applied. The number of observations for model fitting can be increased by including sample-plot measurements from several years. If information about forest disturbances is reliable, then forest height or wood-volume data for the sample plots measured a few years before the collection of ALS data can be updated with a forest growth model.

Tree species composition prediction was tested in this study for 25 m target pixels. Although the uncertainty in species was substantial at the pixel level, the proportion of aggregated evergreen coniferous species in the composition was predicted reliably, and therefore could in future versions of the system be used for correcting phenology effects.

The NFI is intended to provide timely forest data for the whole country and to furnish published maps of the main variables, namely forest height, standing-wood volume, tree species composition, forest age, basal area, relative density, and location of disturbances. However, even for a small country, such as Estonia it is not economically feasible to carry out ALS measurements every year for the whole area, making it necessary in the case of many Estonian locations to update the data by other means. With a geographically comprehensive prediction of tree species composition, it becomes possible to use forest growth modelling to update forest height and wood-volume predictions for those areas lacking current-year ALS measurements. The current Estonian Land Board flight schedule makes it necessary to apply the modelling over a period of up to two years. For the determination of forest height, an algebraic difference model is used to simulate the growth of Estonian forests (Kiviste, 1997): (5)HALS,updated=f(A1,H1,A,OHOR,PE,TEKE).{H_{ALS,updated}} = f\left( {{A_1},{H_1},A,OHOR,PE,TEKE} \right). Here H_ALS,updated, the updated height for a forest of some specified age A, is calculated from the known value of height H₁ at age A₁, the known soil organogenic layer thickness OHOR (in cm), and the known stand origin TEKE (given as natural or cultivated) for a known dominant species PE. A model similar to (5) is available for wood volume. For areas without current FMI data, forest age can be calculated from a site fertility-index model used in Estonia, based, in turn, on forest height and stand age (Forestinv, 2017).

The remaining two variables to be predicted for forest management planning are stand basal area G and relative density T. For forest inventory, the relative density is calculated as the ratio of stand basal area to a reference (a normal stand in the given place). It is, in principle, possible to predict G and T from ALS data. However, we propose to use the relationship M = G·H·F where wood volume is the value predicted from ALS data, and where the stand form factor height H·F is predicted with a forest-height-based model taken from forest inventory regulations (Forestinv, 2017). For each target pixel, the forest height is then predicted from ALS data with model (1).

As an essential part of our remote-sensing-based prediction system for forest variables, we designed a module for the estimation of prediction errors. This error-budget module is based on the k-NN method and on the assumption that NFI sample-plot field measurement data provide ‘true’ values for the forest variables. According to the proposed method, the prediction error of a forest variable for a given pixel is calculated from known errors of similar pixels (found by the k-NN algorithm) in the NFI data. The validation of prediction models is still based on simple validation datasets. Nevertheless, for practical use, the described method makes it possible to construct uncertainty maps of predicted values for the entire target area.

In this paper, we have presented a fully functional remote-sensing support module for the Estonian NFI, supporting the construction of maps of forest height, standing-wood volume, and tree species composition. Our study has, however, revealed the need for correction of phenology effects in the ALS data. This can be achieved by increasing the accuracy of tree species composition estimates, using longer time series of multispectral satellite data and dense time series of Sentinel-1 SAR data as proposed by Dostálová et al. (2018).

When the system becomes fully operational, it will be possible to determine also forest age at high accuracy, by using time series of multi-temporal satellite images for the detection of forest regeneration fellings (Peterson et al., 2004; Liira et al., 2006). Weak disturbances, such as normal thinnings, can also be detected from bi-temporal multispectral satellite image pairs (Uiga et al., 2003), as well as from sparse point clouds obtained through repeated ALS measurements (Arumäe, et al., 2020). With the combination of these various refinements in technique, Estonian forestry is entering an era in which data for forests at every point of the country are updated yearly and can be used for sustainable forest-management planning, on the basis of 10–30 m spatial units.