Open AccessArticle
Combining Camera Relascope-Measured Field Plots and Multi-Seasonal Landsat 8 Imagery for Enhancing the Forest Inventory of Boreal Forests in Central Russia
by
and
Georgy Rybakov
1,*,
Jussi Peuhkurinen
2, Petri Latva-Käyrä
2, Maria Villikka
2, Sanna Sirparanta
2, Alexander Kolesnikov
2, Virpi Junttila
3 and Tuomo Kauranne
31
RusFor Consult Oy Ab, Raseborg, FI-10600 Ekenäs, Finland
2
Arbonaut Ltd., Kaislakatu 2, FI-80130 Joensuu, Finland
3
LUT School of Engineering Science, Lappeenranta University of Technology, FI-53851 Lappeenranta, Finland
*
Author to whom correspondence should be addressed.
Remote Sens.2018,10(11), 1796;https://doi.org/10.3390/rs10111796
Submission received: 23 October 2018 /Revised: 8 November 2018 /Accepted: 8 November 2018 /Published: 13 November 2018
(This article belongs to the SectionForest Remote Sensing)
Abstract
:The study considers a forest inventory for the mean volume, basal area, and coniferous/deciduous mapping of a large territory in central Siberia (Russia), employing a camera relascope at arbitrary sized sample plots and medium resolution satellite imagery Landsat 8 from the leaf-on and leaf-off seasons. The research bases are on field plots and satellite data that are acquired for the real operational forest inventory, performed for industrial purposes during summer–fall 2015. Sparse Bayesian regression was used to estimate linear regression models between field-measured variables and features derived from satellite data. Coniferous/deciduous mapping was done, applying maximum likelihood classification. The study reported the root mean square error for the mean volume and basal area under 25% for both the plot level and compartment level. The overall accuracy of the forest-type classification in coniferous, mixed coniferous/deciduous, and deciduous classes was 71.6%. The features of Landsat 8 images from both seasons were selected in almost every model, indicating that the use of satellite imagery from different seasons improved the estimation accuracy. It has been shown that the combination of camera relascope-based field data and medium-resolution satellite imagery gives accurate enough results that compare well with previous studies in that field, and provide fast and solid data about forests of large areas for efficient investment decision making.
1. Introduction
Assessment and monitoring of forests’ growing stock is important for sustainable use of forest resources, planning investments in forest sector and mitigating effects of climate change. Forests’ growing stock is systematically monitored with different kinds of forest inventory methods. For model based forest inventory using remotely sensed data, there are different sources of data that could be applied, such as satellite imagery, aerial photography, airborne laser scanning (ALS) or a combination of these. The choice of the data source used is application-dependent. It depends on the required accuracy and the level of detail in the resulting data and outcomes. The methods based on airborne data are rather complex in terms of practical arrangements, since they require a mission-specific data acquisition campaign. Due to their relatively high accuracy and detailed data deliverables, airborne data are convenient for the purpose of operational, detailed forest management planning. A forest inventory based on satellite imagery is economical for large area mapping, as in the general case, medium-resolution satellite imagery such as Landsat 8 Operational Land Imager (OLI) is freely available. However, a forest inventory based on medium resolution satellite images has lower accuracy than one that is based on airborne measurements [1], but its accuracy can still be adequate for general forest growing stock assessment at a low spatial resolution. For a concise but comprehensive overview of the relative merits of different remote sensing modalities and estimation methods, please see the review [2] and references therein.
According to the National Forest inventory guidelines of Russia [3], the period of repetitive forest data observations is set to 10–15 years. For some of Russia’s forested areas, this recommendation is not always followed, due to remoteness, lack of transport infrastructure, high field work costs, and the low economic value of the forests, considering the overall forest industry development level in such regions. In practice, this means that some of the existing official forest data are extremely outdated, being acquired, e.g., in the 1980s–90s of the 20th century, and they do not consider, for example, recent forest fire events and growth increments. At the very same time, there are a number of forest industry development initiatives that currently require information about the state of the forests, in order to provide efficient investment decisions both on the level of private investors and the government. In most cases, forest industry projects in Russia consider large-size forest concessions, ranging from hundreds of thousands to millions of hectares of forestland. Due to the financial structuring of such projects, it is often required that forest assessment results for an entire area under consideration are delivered in a short time frame, with accurate and objective results. From that point of view, one possible approach is the use of satellite imagery-based forest inventory.
Satellite imagery-based forest inventory is indeed already a rather old technique; the fundamentals of the approach were introduced already in the 90s; see, e.g., [4,5]. However, despite the fact that satellite imagery-based forest inventory techniques have been well elaborated, one of the most critical phases of the estimation process, field data acquisition, has not been changing much over time. Field data measurements, on the one hand, influence the quality of the forest inventory results, but on the other hand, a field campaign carries a large portion of the forest inventory campaign’s costs. Typically, employing conventional tools such as calipers and a hypsometer/clinometer in field data acquisition has been rather labor intensive. Furthermore, the approach using conventional tools, in some specific cases, suffers from low objectivity and demands field data cross-validation by revisiting a certain number of sample plots, in order to independently confirm field data quality. Modern-type calipers and hypsometers/clinometers are indeed rather sophisticated tools nowadays; see e.g., [6], possessing Global Positioning System (GPS) readings that improve the objectivity of the measurements. Also, the accuracy and precision of measurements is expected to be high when trained mensurationists are using the tools [7]. Even though such tools are advanced equipment, and even though they can provide high accuracy and precision, it is known from their practical application that they remain complex and time consuming, as well as requiring well-trained and experienced field crews. In that specific context, there is a need to employ a field sampling tools that enable faster field campaigns, maintaining a good level of acquired data accuracy, and providing a high level of data objectivity in harsh conditions. Objectivity requirement of the data acquisition can be achieved by automated non-contact field data collection methods. These include terrestrial laser scanning (TLS) [8,9,10] and close-range photogrammetry [11,12]. These methods, however, require purpose-built equipment and computationally heavy post-processing, and their performance is related to the complexity of the forest structure, and their measurement accuracy in field tests has not always been convincing. Mobile phone technology offers an alternative to automated field data collection, with a more generally available platform. Mobile devices and applications have been successfully tested in participatory projects for collecting data for forest fuel estimation [13], in forest monitoring [14], and in creating forest cover maps [15]. Besides their use as a data collection device designed for professional and non-professional users, mobile phones can be used as measuring devices. The camera relascope [16] is an alternative method to automate and streamline field data collection, and to allow non-professionals to participate in data acquisition campaigns.
The idea of using data collected using mobile applications as training data in satellite-based mapping of forest related variables has been introduced in [13,15,16]. In this study, we push the idea forward by testing a commercially available camera relascope tool, Trestima (www.trestima.com/w/en/) [17], in collecting field sample data at arbitrary sized sample plots, and we use the collected field data as training data in Landsat 8 OLI based forest inventory in a large area in central Siberia (Russia). Trestima differs from another camera relascope tool, Relasphone, described in [16] in its level of automatization. While in Relasphone, the user must count the measured trees and define the species for each tree, Trestima detects the species and delivers the basic metrics of a forest stand as basal area, number of stems, and diameter distribution, based on automatic image interpretation.
The way in which this article uses Landsat 8 OLI data is relatively simple. We apply satellite images from two seasons, leaf-on (summer) and leaf-off (spring/autumn) seasons. By this, we aim at a better separation of tree species. However, we would expect further improvements of the results to be obtainable from modern methods that use a time series of satellite images (see e.g., [18]), or even the difference of heights derived from different space-born instruments, such as the method proposed in [19].
This study considers the forest inventory for mean volume, basal area, and coniferous/deciduous mapping of a large area in central Siberia, approaching 2.6 million hectares (ha), employing modern camera relascope-based field measurements at arbitrary sized sample plots and medium-resolution freely sourced satellite imagery from Landsat 8 OLI. The applied camera relascope tool was Trestima, which is a commercial and operationally used forest inventory system used mainly in stand-level assessments. The assessment accuracy is evaluated earlier in Finnish conditions for stand-level mapping [20]. However, the data collected with the tool has not been used for training data in satellite-based inventory.
The objective of the study was to validate the applicability of camera relascope measurements in field plots in the environment of the large-scale forest inventory, and to test the suitability of these field data in medium-resolution satellite-based forest inventory. Also, the study validates the use of Landsat 8 OLI images from different seasons, with the aim to improve estimation accuracy.
2. Materials and Methods
This section and its subsections describe the materials used, and the methods divided in the characterization of the study area, acquisition and processing of input data, data analysis and modeling, and finally the description of the validation methods.Figure 1 presents a workflow diagram of the whole process.
2.1. Study Area
The forest inventory took place in the summer–fall season of 2015. The study area, later the area of interest (AOI), is located in central Siberia (Krasnoyarsk region, Russia); seeFigure 2.
The dimensions of the AOI are approximately 450 km along and 350 km across. In total, the AOI comprises 2,691,316 ha. The AOI is featured by lack of road infrastructure, and mainly represents intact forests (taiga). The AOI consists of forest units, the majority of which has a rectangular form of so-called forest kvartals scattered in either large or small groups and/or lying all along over the territory (see [21] for the definition of a kvartal). Geographically, the AOI is located in the proximity of the Yenisei river on its west and east banks.
Forests on both banks of the river were represented by the coniferous species of Scots pine (Pinus sylvestris L.), Norway spruce (Picea abies L.), Fir (Abies Mill.), Larch (Larix sibirica Lebed.), Siberian pine (Pinus sibirica Du Tour), and by the deciduous species of Birch (Betula L.) and Aspen (Populus tremula L.), all being of different ages at the different parts of the AOI. The most southern part of the AOI differs from the rest of the area. It has been under exploitative forest use, and deciduous species are more dominant. At the time of the study, actual dominance and distribution of species over the area were not well known. Some parts of the AOI have undergone recent severe forest fires. The environment on the west bank is characterized by a flat terrain of lowlands and widespread swamps. The environment on the east bank of the river is rather mountainous, with elevations of up to 700–900 m above sea level, and characterized by frequent large tributaries of the Yenisei river that penetrate deep into the territory over many kilometers.
2.2. Remote Sensing Data
The remote sensing data employed in the study were cloud-free Landsat 8 OLI imagery of both leaf-on (summer) and leaf-off (autumn) seasons for the entire AOI. Processing level-1 images were obtained from the U.S. Geological Survey, Earth explorer service [22]. The image search was first conducted for images from 2015 with less than 10% cloud cover; however, the search criteria were later adjusted to include images with higher cloud cover percentages, and the search was expanded to also include the years 2013 and 2014. Leaf-off images were searched from both autumn and spring seasons. Mountainous areas was problematic because of snow cover: suitable leaf-off images were not detected for that area. The downloaded imagery was first analyzed visually using LandsatLook images as aids. This was done to identify the most suitable images, possible mosaicking the order of images, and checking that the images covered the entire area, and also after masking out the clouds. The total number of images used was 18 and 9 for leaf-on and leaf-off mosaics, respectively. Acquisition dates for leaf-on images were 1 July 2014, 12 July 2014, 4 August 2014, 13 August 2014, 6 June 2015, 15 June 2015, 22 June 2015, 25 August 2015, and 1 September 2015. Acquisition dates for leaf-off images were 29 September 2013, 28 September 2014, and 30 September 2014. For further processing, the original 12 bit images were resampled to 32 bit images, because the ArboliDAR tools [23] used for feature calculation required 32 bit raster data as an input.
The mosaicking of images was done in ArcGIS 10.0 software. Due to variation in satellite image acquisition conditions, the same ground object on two overlapping images can result in different spectral values [24]. Because of this, radiometrically uniform mosaics using multi-temporal scenes should be created prior to employing the satellite imagery into image analysis. To overcome radiometric differences between the multi-temporal scenes, relative normalization by pre-processing of the input images was performed. Relative radiometric normalization is based on the assumption that the difference between two satellite images from different dates but from the same area can be explained by a linear relationship [25]. To resolve this problem, we used the variance-preserving mosaic (VPM) algorithm that minimizes the overall error of the normalization, and aims to preserve the average variance of input scenes [26]. Normalization was done using the VPM algorithm implemented into the ArboLIDAR toolbox, developed by Arbonaut Ltd., Finland [23]. The normalization procedure utilizes overlapping portions of the images, and prior to normalization, clouded and shadowed areas as well as image edges with black or clearly abnormal pixel values were delineated. Portions of the image that do not contain relevant information and that might affect the normalization negatively, need to be masked out; i.e., pixel values need to be set as No Data. Polygon shapefiles were used as masks. The AOI was buffered by 50 km. Areas further than 50 km from AOI were set as No Data. The black edges of each image needed to be set to No Data. In addition to the black edges, some pixels on the edge of the scenes had abnormal values. Worldwide Reference System (WRS-2) scene boundaries downloaded from U.S. Geological Survey Earth explorer service [22] were used to remove areas on the edges of composite images. Other masks were created mostly with manual editing. Images were analyzed visually to detect clouds and shadows, or other areas with inferior quality, and then these areas were delineated by hand. For some of the Landsat scenes, LandsatLook quality images were further processed to mask out very scattered small clouds and their shadows. If there were clouds outside the AOI, they needed to be masked out as well, so that they would not affect the normalization result negatively. Several images adjacent to each other had no visible image borders. This was the case, for example, for when the images were from the same date. Based on visual analysis, it was decided that several of the masked composite images could be mosaicked without normalization. These mosaics and other individual composite images were then normalized and mosaicked into one image mosaic covering the entire AOI (Figure 3).
2.3. Sample Design Plan
Initially, the sample plot design was an L-shaped cluster format with four plots in each cluster. The distance between center coordinates of plots was 200 m. The AOI defined by the rectangular planning compartment, kvartal boundaries, which were preliminarily digitized in GIS shape file format, was divided into two areas—on the west bank and on the east bank of Yenisei, respectively. The forested area was defined by delineating out non-forest areas based on Hansen tree cover data [27], as well as rivers and roads, with a 100 m buffer. The accessible forest was defined by using roads and rivers (at the field campaign the rivers were also used as transport lines) that were preliminarily evaluated in terms of their usability. The maximum distance from a usable road to a plot was 2000 m, and no rivers could be crossed on the way to the plot. Cluster locations were randomly sampled in accessible forest area. Each initial sample plot was of Landsat 8 OLI -pixel size of 30 × 30 m.
Next, the initial sample plots were compared over the forest area of AOI variation. The aim was to have field data representing the variation of the whole AOI, in order to obtain good estimation accuracy and precision with a low number of field plots; see e.g., [28,29,30]. The representativeness of the initial sample was analyzed by comparing and plotting green, red and NIR (near infrared) band values and NDVI (normalized difference vegetation index) image values from the initial sample plots, and from the whole forested area of the AOI. The values for the whole forest area of AOI were from a large random sample (n = 15,355). The materials used for the analysis were cloud-free image mosaics produced from the best available summer season Landsat 8 OLI images that were obtained before August 2015. Image parts with cloud or cloud shadow were removed and supplemented with data from other images.
Arbitrary sized forest patches around the plot center coordinates were then drawn manually. A patch was not allowed to be in a road or a river and the land use/land cover class, tree species composition and volume were to be visually assessed as being homogeneous inside the patch. When it comes to the volume, the assessment was done in the sense that, provided the land use/land cover class, tree species composition is the same, and then if the patch of forest seemed visually homogeneous, it was assumed that the forest is also similar in volume. Thus, in manual segmentation, a plot could be moved to a more homogeneous location, or removed if it could not be moved. The manual segmentation was performed based on the visual interpretation of Landsat 8 OLI imagery, with the support of satellite imagery available on Google Earth. On average, an arbitrary-sized plot area was set to be approximately 0.5–1 ha (6–12 Landsat 8 OLI pixels), which could vary depending on the homogeneity.
The size of the initial sample was 848 plots, with 424 plots on each side of the river Yenisei. Some of the plots from the initial sample were removed because they were not accessible after reviewing them, based on the information from the field team. After the removal of the inaccessible plots, the revised sample consisted of a total of 792 plots. All in all, 392 and 400 plots were left at the east and west banks of the river Yenisei, respectively. The representativeness of the revised sample against Landsat 8 OLI satellite image spectral values is presented in inTable 1.
The analysis showed that on the east bank of the river, the revised sample well represented the variation of the whole inventory area. There were some extreme values (e.g., in green and red bands) that were not covered in the sample. It was assumed that these extremes did not represent the forests, but some other land use classes. On the west bank of the river, the sample covered the variation well; however, the means were shifted towards high NIR and NDVI values. Also, the variation of the inventory area differed significantly, compared to the areas on the east bank of the river. The most southern planning compartments (kvartals) were well-represented in the sample, and the rest of the kvartals on the west bank of the river Yenisei were underrepresented. The reason for this was that the accessible area on the west bank was quite limited.
Finally, during manual segmentation, the number of sample plots was reduced to 385 and 399 at the west and east banks of the river Yenisei, respectively, due to some inaccessible plots found in visual interpretation. An average area of a forest patch was 0.8 ha, with a minimum area of 0.4 ha and a maximum area of 1.8 ha. The representativeness of the final sample was investigated by calculating NIR and NDVI values as the mean of pixel values in the patch. This slightly averaged the values compared to the pixel-based method used in analyzing the revised sample. However, the averaging effect was minor, since the patches were segmented in homogeneous areas. The representativeness of the final sample was checked again against the random sample. Scatter plots inFigure 4 show how NIR and NDVI values of the final sample compared with the random sample presenting the variation in AOI.
2.4. Field Plot Data Acquisition
The fieldwork took place in the period between late August and early October 2015. The Trestima camera relascope system was used to measure the field plots. The Trestima camera relascope is a new generation of terrestrial forest inventory tools, which is based on using a mobile application. As input data, Trestima uses terrestrial imagery of two types: basal area images and individual tree images, captured with the mobile application. Trestima automatically detects species and delivers such basic metrics of a forest stand such as basal area, the number of stems, and diameter distribution, which is defined as the number of trees in each diameter class, as well as detecting diameter and height, when trees are measured individually [31,32]. Because there was no assured information available on the accuracy of the automatic species detection in the study area, each image at each sample plot was reviewed manually, and if the tree species was wrong, it was changed.
The field work was performed using the mobile phones Sony Xperia Z3 (Android) and Nokia Lumia 720 (Microsoft). Prior to field measurements, each sample plot was covered with a systematic net of points (30 × 30 m). The net was visible on each mobile device. The points in the mobile devices served as navigation orienteers for the field crew. The assignment was to approach to as near as possible at each such point, taking into consideration the local environment and the mobile device’s GPS accuracy, and to capture two sample images of basal area in approximately perpendicular directions inward a sample plot (Figure 5). On a typical plot, there were seven points, which resulted in 14 basal area images (7 × 2 images). Simultaneously, at each sample plot, the field crew selected a minimum of two model trees per each tree species, and measured the individual trees’ diameter and height (seeFigure 6). This means that the number of model tree images was dependent on how many tree species were present on the plot. The average value of basal area from basal area images was used as the field measured basal area of the plot. Model trees were used for estimating the height–diameter curve. All in all, 16,516 sample images of either type, i.e., the basal area image or the model tree image, were taken from a total of 599 field plots.
By default, Trestima calculates the sample plot volume automatically, based on the basal area sample images from the plot, as well as the height and diameter measurements of subjectively selected trees per tree species. The volumes are calculated using volume models per tree species, configured in Trestima [31]. However, as the default Trestima volume models were not tailored to this project, the estimation of the forest volume per each sample plot was done using Trestima diameter distribution from each plot, which defined the number of trees on each diameter class by species. Heights were estimated with height–diameter models (h(d)-models) generated from the model tree data, and tree level allometric volume functions specific to the forest region of the AOI were used to calculate diameter class level volumes. The height per tree species was estimated for each diameter class by using h(d)-models estimated for the southern part, and the rest (western, south-western, and eastern) separately. The models were created in the statistical program R [33] and in the package nlme [34] using nonlinear mixed-effects modeling. The mixed-effects models were applied for modeling the relationship between tree height and tree diameter, because the variables were clustered and thus spatially correlated, see e.g., [35]. The plot was used as the random variable in the modeling, assuming that different plots have somewhat different growth conditions, and thus plot affects the h(d)-model curve. The h(d)-models were created separately for each major tree species of the following areas: in the southern part pine, spruce, birch, and aspen each had their own models, and in the western, southwestern, and eastern parts, cedar and other (deciduous) tree species each had their own, in addition to those mentioned before. Some minor tree species data were combined to suitable major tree species groups: in the western, southwestern, and eastern parts larch and fir were merged with spruce, in the southern part, cedar and other species with pine. All h(d)-models were made with the Näslund function [36] after testing several options.
Finally, the volume of the growing stock per each sample plot was estimated by summing up the volumes of the diameter classes per tree species, which were calculated based on tree-specific allometric volume models [37], taking as the input data the estimated quantity of the trees of each respective diameter class.
2.5. Land Use Classification
Maximum Likelihood Classification tool in ArcGIS 10.0 toolbox Spatial Analyst was used to classify a Landsat 8 OLI mosaic with supervised maximum likelihood method. The area was classified into nine land use classes (Conifer forest, Mixed conifer forest, Deciduous forest, Swamp, Grass, Harvested, Regrowth, Open land, Water and Burnt area). Training areas consisted of field plots measured with Trestima, and additional training areas established based on visual interpretation by creating several polygons for each land use class. RapidEye satellite imagery and satellite imagery available on Google Earth were used as ancillary information for establishing training areas. An independent random sample of 207 visually interpreted plots was acquired for validation.
2.6. Estimation of Forest Volume and Basal Area
All in all, data from a total of 599 plots were delivered from the field campaign. Of the 599 plots, 544 were used for further modeling. Some plots were left outside modeling, due to poor satellite image quality on the plot region, or due to problems in measurements. The highest possible model quality was assured by not using plots with some quality issues that might distort the models. The forest volume modeling was done separately for the four major parts of the AOI: South, West, Southwest, and East parts (Figure 7). The division was based on satellite imagery and forest characteristics, i.e., the modeling area was the area with homogenous satellite data and similar forest types (lowlands, mountainous).
The area of each respective part is estimated as follows:
- South—ca. 50,000 ha.
- West—ca. 1,000,000 ha;
- Southwest—ca. 175,000 ha;
- East—ca. 1,100,000 ha.
At the estimation, different field plot combinations were used in the different areas. By default, both leaf-on and leaf-off Landsat 8 OLI images were used in models, except for the mountainous area of the east bank, where different image compositions were applied because good quality leaf-off images were not available for the whole area. Also, because of mountains, the East area had large shadowed areas, and therefore, separate shadow area models were estimated. Field data statistics for the modeling areas are presented inTable 2. The used field plot and Landsat 8 OLI combinations were:
- South area (both leaf-on and leaf-off imagery):
- -
- 228 plots (all plots from the Southern part);
- West area (both leaf-on and leaf-off imagery):
- -
- 119 plots (all plots from the Western part);
- Southwest area (both leaf-on and leaf-off imagery):
- -
- 217 plots (119 of plots from the Western and 76 of plots from the Southern part);
- East area:
- -
- Leaf-on and leaf-off imagery (model 1): 160 plots;
- -
- Leaf-on imagery only (model 2): 160 plots;
- -
- Shadowed leaf-on and leaf-off imagery (model 3): 43 plots;
- -
- Shadowed leaf-on imagery (model 4): 43 plots.
Models for volume and basal area estimation were built using field plot data and satellite image features calculated for sample plots using Sparse Bayesian Regression (SBR) method [38,39]. SBR is a linear regression method, which automatically selects the needed covariates from the set of candidate covariates to estimate a forest variable. By default, in SBR, a specific set of covariates is employed in order to estimate each forest variable.
The modeling and estimation were done for totals and for different forest classes (coniferous, mixed-coniferous, and deciduous) separately, meaning that deciduous and coniferous forest classes had different models. In the case of forest class-based modeling, all models except East area models 3 and 4, had different models for each forest class (coniferous, mixed-coniferous, and deciduous). East area models 3 and 4 were mainly coniferous, and therefore, forest class-specific models were not possible to estimate. The modeling was executed in the statistical program “R”, and in ArcMap 10.0 with the ArboLIDAR toolbox developed by Arbonaut Ltd., Finland [23]. Results were calculated for a grid with a cell size of 0.81 ha (90 × 90 m), which is 3 × 3 Landsat pixels, and corresponded close to the average size of the field sample plot.
The list of satellite image features derived from remote sensing data used as candidate covariates in SBR and their description is presented inAppendix A. Features such as the mean and standard deviation of band values were calculated from Landsat 8 OLI bands 2–6. Features based on NDVI images were calculated from both leaf-on and leaf-off images. NDVI images were used for calculating 14 Haralick textural features to capture the textural variation of images [40].
2.7. Accuracy Assessment
The accuracy assessment of volume and basal area estimations was conducted with plot-level leave-one-plot-out validation, and additionally at the compartment level for each respective modeling area (seeFigure 7).
There were no independent reference data for compartment-level estimates. In the distribution of forest variables, there is very often a significant spatial correlation [41]. Although the field plots in each cluster were not necessarily inside the same stand, we can expect that forest variables tend to be more similar when they are spatially closer together. We merged the field plots in clusters, and compartment-level accuracies were calculated by comparing the estimated and field measurements values at the cluster level. The clusters had on average, an area of around 3 ha, and this emulates the size of an average compartment.
An accuracy assessment was performed using the root mean square error (RMSE), which shows the probability of an estimate to deviate from the observed (measured) value (Equation (1)):
and relative RMSE (Equation (2)):
where is an estimated value for the variabley in sample ploti, is a measured value for the variabley in sample ploti, is an average of the measured variablesy, and is the number of field sample plots.
Simultaneously, bias—the difference of the predicted mean and the measured mean, which indicates a systematic error in the estimation, was calculated using Equation (3):
and relative bias with Equation (4):
The accuracy assessment of land use classification was done using separate random sample plot data withN randomly sampled plots. From the original sample, the plots with sufficient reference imagery were used in the accuracy assessment process only. Each sample plot was given a land use class from the land use classification map, based on the 30 × 30 m map cell that the sample point fell in. The correct class was then given to each point visually. Accuracies were calculated for all nine classes found from the project area, for forest classes only, and for forest–non-forest classification.
Overall accuracy (OA) and Cohen’s Kappa coefficient (Kappa) were calculated from the error matrix to express the accuracy of the land use classification (Equation (5)):
and (Equation (6)):
whereNcorrect is the number of correct observations,N is the total number of observations,p0 is the relative number of correct observations, andpe is the relative number of correct observations that would be expected by the relative frequency of each class. Kappa value can range from −1 to 1. A Kappa value of 1 implies complete agreement, and Kappa values of less than 0.4 imply a less than moderate agreement.
3. Results
Plot-level validation results show low and relatively similar RMSE and bias values for total value estimates in all modeling areas (Table 3). Also, forest class-based estimates show similar accuracy levels in terms of relative RMSE (Table 4,Table 5 andTable 6). However, absolute RMSE values indicate that the estimation error for the deciduous forest class is smaller in the South area than in other modeling areas. Based on field plot data, the South area is dominated by the deciduous forest class, whereas in other areas, coniferous and mixed coniferous classes represent the majority of forests.
Satellite features from leaf-on images were selected in all models, except the South area, with the mixed coniferous model (Appendix A). Features from leaf-off images were selected in 11/13 models for which Landsat 8 OLI images from both seasons were available. The features were not selected in the West area coniferous and southwest mixed coniferous models.
Coefficients of determination (R2) for timber volume vary between 0.55–0.13, 0.65–0.012, and 0.53–0.22 for coniferous, deciduous, and mixed coniferous classes, respectively. The highest R2 values were in the South area, with the highest number of field plots. The lowest values were in the mountainous East area. The scatter plot inFigure 8 presents the measured field plot values against the estimated values in forest class-based estimation.
On the compartment level, the validation demonstrated smaller RMSE values than on the plot level (Table 7).
The maximum likelihood classification of different land use classes resulted in an overall accuracy of 66.2% and a Kappa coefficient of 0.55, which indicates moderate agreement with ground truth classes (Table 8). Classification results between forest and non-forest classes were slightly better, yielding an overall accuracy of 90.3% and a Kappa coefficient 0.66, which indicates a good agreement with ground truth classes. Results between three forest classes were 71.6% for overall accuracy and 0.56 for the Kappa coefficient, indicating moderate agreement with the ground truth.
4. Discussion
Satellite imagery-based forest inventories have been studied extensively and used for forest inventories all over the world for several decades. The camera relascope field measurement technique in turn has already been in use for some years as a forest inventory tool, e.g., at pre-harvesting inventory, measuring at specific forest compartments, and at other small- or medium-sized forest inventories [42]. The concept of using a camera relascope for measuring field reference plots in satellite-based inventory was introduced few years ago (also [16]). This research is unique in the sense that the data was acquired in a real operational large-scale forest inventory (ca. 2.6 million ha) conducted for industrial purposes. Also, it is the first study where camera relascope measurements based on automatic image interpretation are used in combination with satellite-based forest inventory.
Nevertheless, when it comes to similar studies in boreal forests, the results of such studies concerned with satellite imagery-based forest inventories over large areas have shown quite large variations in prediction accuracy, depending on the mean compartment size, study area, field reference data, and satellite data used. For example, [43] achieved an RMSE of 36% at the compartment level when using national forest inventory field reference plots and medium-resolution satellite data with rather large average compartment sizes (19 ha) in Sweden. Ref [44] reported a RMSE of 48% for the compartment-level total mean volume by using Landsat TM images with compartment-level field inventory data as the reference data in Finland. Ref [45] used Landsat ETM+ and field sample plots in Alberta, Canada to estimate above-ground biomass and volume. The plot level RMSE in the modeling data was 40% when the volume was estimated using forest structural parameters (height and crown closure estimates), and 59% when direct estimation was used. The reported R2 values, 0.71 and 0.30 for estimation through forest structural parameters and direct estimation, respectively, were significantly high [45]. Unfortunately, compartment-level RMSEs were not reported. As [46] points out in a discussion paper, the relative RMSE tends to decrease, when the reference area size increases. According to [46], the RMSE of the total volume in satellite-based forest inventories has been between 40% and 60% at typical Finnish compartment levels (average compartment size 1–2 ha). More recent studies applying Sentinel-2 satellite images in the estimation of total volume have reported promising accuracies, although in small study areas. Research carried out in Norway in forest areas dominated by coniferous forests reported a relative RMSE of 30.1% and a R2 of 0.48 when using field plots measured using traditional field methods [47]. In the same study, training plots based on analyzing point cloud data from an unmanned aerial vehicle resulted in a relative RMSE of 37.5% and a R2 of 0.30. The study was carried out in a small research area of 7330 ha, and field data comprised 33 plots of size 500 m2. Similar results were reported in research done in Eastern Poland from a study area of 7800 ha dominated by Scots pine [1]. Studies reported a relative RMSE for a volume of 35.1% and a R2 of 0.24. The number of field plots were 94, and the size of the plot was 400 m2. In that respect, the results of the present study, with RMSE being reported for both the plot level and the compartment level (with mean size of approximately 3 ha) in all major parts of the AOI (South, West and East) under 25%, have shown that they compare well to previous satellite-based forest inventories in coniferous forests of the boreal zone.
It is important to have a good spatial match between the field reference data and the corresponding satellite features. Ref [45] presented a somewhat similar method to our study in improving the match between the field reference data and satellite image features. They used automatic segmentation to remove the spectrally uncharacteristic pixels from the 3 by 3 pixel window surrounding the field plot. In their case, the field plot size was only about half of the image pixel size. However, the segmentation improved the correlations between the image features, and the field measurements consistently compared with the full 3 by 3-pixel window. We used manual segmentation to delineate the plot and camera relascope to measure the sample plots. The camera relascope approach made it possible to measure larger reference plots, and not only one sample plot inside the segment, compared to manual field measurement methods with the given budget. However, it would have been interesting to compare the results between the manually measured smaller samples plot and camera relascope plots with the use of the segmentation method. This was unfortunately not possible, due to the budget limitations in our field campaign. Ref [20] reported RMSEs of 19.7 to 29.7% for the plot level (plot size 32 by 32 m) basal area when comparing Trestima measurements with accurate manual field measurements. The most accurate results were achieved when the images were taken toward plot center. In [20], the increase from four images per plot to eight images per plot improved the results only marginally. The mean diameter and height were measured more accurately: 5.2–11.6 and 10.0–13.6 for mean diameter and height, respectively, depending on the tree species [20]. We measured larger plots in different conditions so that the results from [20] are only directive for the measurement accuracy for our study.
Our results strengthen the findings of earlier studies concerning the suitability of camera relascope data in satellite-based forest inventory. Ref [16] reported a RMSE of 29.7% for the basal area when comparing basal area measured with a Relasphone application and accurate field measurements in a small study area in Finland. The plot size in their data was between 0.2 and 1.0 ha, and on average, two Relasphone sample plots were measured in each plot. In the same study, the RMSE of the basal area in a test carried out in Durango, Mexico, was 17.2% and 18.0% for two observers. The reference plot size was 0.25 ha. Furthermore, the field plot data measured in Durango with a Relasphone was used as training data in a Landsat 8-based mean volume model. The standard error of the model estimate was 34% of the mean measured stem volumes [16]. The reported RMSEs in [16] cannot be used as a reference for the accuracy of our study; the Relasphone application uses a different sampling technique, and field plot size used was smaller. Despite the challenges in comparing the results of different studies, applying a camera relascope in measuring field training plots for satellite-based forest inventory seems to produce estimation accuracies that are comparable with methods applying manual field measurement methods for the same purpose, i.e., estimating the volume in boreal forests.
Applying multi-seasonal satellite data proved to be useful. Satellite features from Landsat 8 OLI images of both seasons were selected for almost every model. This indicates that multi-seasonal satellite images have more information than the imagery from one season in the context of the forest inventory. Thus, when applicable, it is recommended to search for multi-seasonal imagery, and to test case-by-case if leaf-on, leaf-off or a combination of these produced the best results.
The presented concept and methodology was very promising for being used in large areas of satellite-based forest inventory and mapping projects. Future research and testing were still needed to validate the methodology against other field measurement techniques, and to fine tune the optimal plot type and image sample design. Also, the applicability in more challenging environments with more heterogenous forest structures and a wider variety of species should be tested.
5. Conclusions
The results obtained in this research effort demonstrated that it was possible to attain useful and accurate forest estimation results, and to produce the corresponding forest maps for planning the use of forest resources, applying a combination of camera relascope field measurements and medium-resolution satellite imagery in inaccessible northern boreal forests. Using the camera relascope in measuring the training data allowed for measuring substantially larger field plots, compared to manual field measuring methods. The achieved compartment level accuracy (average compartment size 3 ha) was 11.8–24.9% and 14.0–35.0%, in terms of RMSE, relative to the measured mean value for the basal area and mean volume, respectively. These can be considered as good results compared to satellite-based forest inventory studies in general. The method needs further testing and validation in finding the optimal field plot type, including the size of the field plot and the number of mobile phone images per plot. Also, manual field measurements should be carried out in the same study area to compare the accuracy and efficiency of different field measurement methods more carefully.
Author Contributions
Conceptualization, G.R. and J.P.; methodology, G.R., J.P., P.L.-K., S.S., V.J., and T.K.; software, A.K.; validation, G.R., J.P., and P.L.-K.; formal analysis, J.P., P.L.-K., S.S., A.K., V.J., and T.K.; investigation, G.R., J.P., P.L.-K., M.V., S.S., A.K., and V.J.; resources, G.R., J.P., and T.K.; data curation, G.R., P.L.-K., and M.V.; writing—original draft preparation, G.R., J.P., P.L.-K., M.V., S.S., A.K., V.J., and T.K.; writing—review and editing, G.R., J.P., V.J., and T.K.; visualization, P.L.-K. J.P., and M.V.; supervision, J.P. and T.K.; project administration, G.R., J. P., and T.K.; funding acquisition, T.K.
Funding
This research received no external funding.
Conflicts of Interest
The authors declare no conflict of interest.
Appendix A
Table A1. List of satellite image features used for respective areas and forest types. D = deciduous, M = mixed coniferous, C = coniferous, S = shadow.
Table A1. List of satellite image features used for respective areas and forest types. D = deciduous, M = mixed coniferous, C = coniferous, S = shadow.
| Feature | South | West | Southwest | East Leaf-On | East Leaf-On and Leaf-Of | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| D | M | C | D | M | C | D | M | C | D | M | C | S | D | M | C | S | |
| Lon_B2_MEAN | x | x | x | x | |||||||||||||
| Lon_B2_SD | x | x | x | ||||||||||||||
| Lon_B3_MEAN | x | x | x | x | x | x | x | x | |||||||||
| Lon_B3_SD | x | x | x | x | |||||||||||||
| Lon_B4_MEAN | x | x | |||||||||||||||
| Lon_B4_SD | x | x | |||||||||||||||
| Lon_B5_MEAN | x | x | x | x | x | ||||||||||||
| Lon_B6_MEAN | x | ||||||||||||||||
| Lon_B6_SD | x | ||||||||||||||||
| Loff_B2_SD | x | ||||||||||||||||
| Loff_B3_MEAN | x | x | x | x | x | x | |||||||||||
| Loff_B3_SD | x | ||||||||||||||||
| Loff_B4_MEAN | x | ||||||||||||||||
| Loff_B5_MEAN | x | x | x | x | x | x | |||||||||||
| Loff_B5_SD | x | x | x | ||||||||||||||
| LOff_B6_MEAN | x | x | |||||||||||||||
| Loff_B6_SD | x | x | |||||||||||||||
| NDVI_On_HR_01 | x | x | |||||||||||||||
| NDVI_On_HR_03 | x | ||||||||||||||||
| NDVI_On_HR_04 | x | ||||||||||||||||
| NDVI_On_HR_05 | x | x | x | ||||||||||||||
| NDVI_On_HR_07 | x | ||||||||||||||||
| NDVI_On_HR_09 | x | ||||||||||||||||
| NDVI_On_HR_10 | x | ||||||||||||||||
| NDVI_On_HR_11 | |||||||||||||||||
| NDVI_On_L8_SD | x | ||||||||||||||||
| NDVI_Off_HR_01 | x | ||||||||||||||||
| NDVI_Off_HR_03 | x | x | |||||||||||||||
| NDVI_Off_HR_05 | x | ||||||||||||||||
| NDVI_Off_HR_06 | x | ||||||||||||||||
| NDVI_Off_HR_10 | x | ||||||||||||||||
| NDVI_Off_HR_12 | x | ||||||||||||||||
| NDVI_Off_HR_13 | x | ||||||||||||||||
| NDVI_Off_HR_14 | x | ||||||||||||||||
| NDVI_OffL8_MEAN | x | ||||||||||||||||
| Loff_B4_SD^0.4 | x | ||||||||||||||||
| LOff_B5_SD^0.2 | x | ||||||||||||||||
| LOn_B2_MEAN/Lon_B3_MEAN | x | ||||||||||||||||
Table A2. Description of the satellite features.
| ## | Feature | Description |
|---|---|---|
| 1 | Lon_BX_MEAN | Mean value of band x from Leaf-on Landsat 8 OLI image |
| 2 | Lon_BX_SD | Standard deviation of band x from leaf-on Landsat 8 OLI image |
| 3 | Loff_BX_MEAN | Mean value of band x from Leaf-off Landsat 8 OLI image |
| 4 | Loff_BX_SD | Standard deviation of band x from leaf-off Landsat 8 OLI image |
| 5 | NDVI_On/Off_HR_01 | Angular second moment from leaf-on/leaf-off NDVI image |
| 6 | NDVI_On/Off_HR_02 | Contrast from leaf-on/leaf-off NDVI image |
| 7 | NDVI_On/Off_HR_03 | Correlation from leaf-on/leaf-off NDVI image |
| 8 | NDVI_On/Off_HR_04 | Sum of squares from leaf-on/leaf-off NDVI image |
| 9 | NDVI_On/Off_HR_05 | Inverse difference moment from leaf-on/leaf-off NDVI image |
| 10 | NDVI_On/Off_HR_06 | Sum average from leaf-on/leaf-off NDVI image |
| 11 | NDVI_On/Off_HR_07 | Sum variance from leaf-on/leaf-off NDVI image |
| 12 | NDVI_On/Off_HR_08 | Sum entropy from leaf-on/leaf-off NDVI image |
| 13 | NDVI_On/Off_HR_09 | Entropy from leaf-on/leaf-off NDVI image |
| 14 | NDVI_On/Off_HR_10 | Difference variance from leaf-on/leaf-off NDVI image |
| 15 | NDVI_On/Off_HR_11 | Difference Entropy from leaf-on/leaf-off NDVI image |
| 16 | NDVI_On/Off_HR_12 | Information Measure of correlation |
| 17 | NDVI_On/Off_HR_13 | Information Measure of correlation |
| 18 | NDVI_On/Off_HR_14 | Maximal correlation Coefficient |
| 19 | NDVI_On/OffL8_MEAN | Mean of Leaf-on/Leaf-off NDVI indexes |
| 20 | NDVI_On/OffL8_SD | Standard deviation of Leaf-on/Leaf-off NDVI indexes |
References
- Hawryło, P.; Wezyk, P. Predicting growing stock volume of Scots pine stands using Sentinel-2 satellite imagery and airborne image-derived point clouds.Forests2018,9, 274. [Google Scholar] [CrossRef]
- Koch, B. Remote Sensing supporting national forest inventories. InKnowledge Reference for National Forest Assessments; Food and Agriculture Organization of the United Nations: Rome, Italy, 2015; pp. 77–92. ISBN 978-92-5-308832-4. [Google Scholar]
- Federal Forestry Agency of the Russian Federation. Order No 207, Dated 6 June 2011. On Adoption of Methodological Recommendations for Conduction of Forest State Inventories.Rossiiskaya Gazeta. 2011. Federal Issue 5500(174). Available online:https://rg.ru/2011/08/10/leshoz-dok.html (accessed on 10 April 2017).
- Tomppo, E. Satellite Image-Based National Forest Inventory of Finland.Photogramm. J. Finl.1990,12, 115–120. [Google Scholar]
- Tomppo, E.; Katila, M. Satellite Image-Based National Forest Inventory of Finland—First Large Area Results. In Proceedings of the European Collaborative Programme Workshop on Remote Sensing for Forestry Applications, Copenhagen, Denmark, 13–15 November 1991; pp. 63–67. [Google Scholar]
- Haglöf Company Group. Available online:http://www.haglofcg.com/index.php/en/products/instruments (accessed on 11 May 2017).
- Luoma, V.; Saarinen, N.; Wulder, M.A.; White, J.C.; Vastaranta, M.; Holopainen, M.; Hyyppä, J. Assessing precision in conventional field measurements of individual trees.Forests2017,8, 38. [Google Scholar] [CrossRef]
- Liang, X.; Kankare, V.; Hyyppä, J.; Wang, Y.; Kukko, A.; Haggrén, H.; Yu, X.; Kaartinen, H.; Jaakkola, A.; Guan, F.; et al. Terrestrial laser scanning in forest inventories.ISPRS J. Photogramm. Remote Sens.2016,115, 63–77. [Google Scholar] [CrossRef] [Green Version]
- Liang, X.; Hyyppä, J.; Kaartinen, H.; Lehtomäki, M.; Pyörälä, J.; Pfeifer, N.; Holopainen, M.; Brolly, G.; Francesco, P.; Hackenberg, J.; et al. International benchmarking of terrestrial laser scanning approaches for forest inventories.ISPRS J. Photogramm. Remote Sens.2018,144, 137–179. [Google Scholar] [CrossRef]
- Kvochkin, D.; Ustyugov, V. Mobile device of the automated forest inventory.J. Ind. Pollut. Control2017,33, 976–980. [Google Scholar]
- Mikita, T.; Janata, P.; Surový, P. Forest stand inventory based on combined aerial and terrestrial close-range photogrammetry.Forests2016,7, 165. [Google Scholar] [CrossRef]
- Mokroš, M.; Liang, X.; Surový, P.; Valent, P.; Čerňava, J.; Chudý, F.; Tunák, D.; Saloň, S.; Merganič, J. Evaluation of close-range photogrammetry image collection methods for estimating tree diameters.ISPRS Int. J. Geo-Inf.2018,7, 93. [Google Scholar] [CrossRef]
- Ferster, C.J.; Coops, N.C. Integrating volunteered smartphone data with multispectral remote sensing to estimate forest fuels.Int. J. Digit. Earth2016,9, 171–196. [Google Scholar] [CrossRef]
- Pratihast, A.; Herold, M.; Avitabile, V.; De Bruin, S.; Bartholomeus, H.; Souza, C., Jr.; Ribbe, L. Mobile devices for community-based REDD+ monitoring: A case study for central Vietnam.Sensors2013,13, 21–38. [Google Scholar] [CrossRef] [PubMed]
- Lesiv, M.; Moltchanova, E.; Shchepashchenko, D.; See, L.; Shvidenko, A.; Comber, A.; Fritz, S. Comparison of Data Fusion Methods Using Crowdsourced Data in Creating a Hybrid Forest Cover Map.Remote Sens.2016,8, 261. [Google Scholar] [CrossRef]
- Molinier, M.; López-Sánchez, C.A.; Toivanen, T.; Korpela, I.; Corral-Rivas, J.J.; Tergujeff, R.; Häme, T. Relasphone—Mobile and participative in situ forest biomass measurements supporting satellite image mapping.Remote Sens.2016,8, 869. [Google Scholar] [CrossRef]
- Trestima Forest Inventory System. Available online:https://www.trestima.com/w/en/ (accessed on 31 October 2018).
- White, J.C.; Wulder, M.A.; Hobart, G.W.; Luther, J.E.; Hermosilla, T.; Griffiths, P.; Coops, N.C.; Hall, R.J.; Hostert, P.; Dyk, A.; Guindon, L. Pixel-based image composition for large-area dense time series applications and science.Can. J. Remote Sens.2014,40, 192–212. [Google Scholar] [CrossRef]
- Sadeghi, Y.; St-Onge, B.; Leblon, B.; Prieur, J.-F.; Simard, M. Mapping boreal forest biomass from SRTM and TanDEM-X based on canopy height model and Landsat spectral indices.Int. J. Appl. Earth Obs. Geoinf.2018,68, 202–213. [Google Scholar] [CrossRef]
- Vastaranta, M.; Latorre, E.G.; Luoma, V.; Saarinen, N.; Holopainen, M.; Hyyppä, J. Evaluation of a Smartphone App for Forest Sample Plot Measurements.Forests2015,6, 1179–1194. [Google Scholar] [CrossRef] [Green Version]
- Karvinen, S.; Välkky, E.; Gerasimov, Y.; Dobrovolsky, A.Northewest Russian Forest Sector in a Nutshel; Finnish Forest Research Institute: Joensuu, Finland, 2011; 144p, ISBN 978-951-40-2343-9. [Google Scholar]
- U.S. Geological Survey. Earth Explorer Service. Available online:https://earthexplorer.usgs.gov (accessed on 17 August 2017).
- Arbonaut Products. Available online:http://www.arbonaut.com/en/products (accessed on 1 November 2017).
- Yuan, D.; Elvidge, C.D. Comparison of relative radiometric normalization techniques.ISPRS J. Photogramm. Remote Sens.1996,51, 117–126. [Google Scholar] [CrossRef]
- Vicente-Serrano, S.M.; Prez-Cabello, F.; Lasanta, T. Assessment of radiometric correction techniques in analyzing vegetation variability and change using time series of Landsat images.Remote Sens. Environ.2008,112, 3916–3934. [Google Scholar] [CrossRef]
- Eivazi, A.; Kolesnikov, A.; Junttila, V.; Kauranne, T. Variance-preserving mosaicing of multiple satellite images for forest parameter estimation: Radiometric normalization.ISPRS J. Photogramm. Remote Sens.2015,105, 120–127. [Google Scholar] [CrossRef]
- Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R.; et al. High-Resolution Global Maps of 21st-Century Forest Cover Change.Science2013,342, 850–853. [Google Scholar] [CrossRef] [PubMed]
- Demaerschalk, J.P.; Kozak, A. Suggestions and criteria for more effective regression sampling.Can. J. For. Res.1974,4, 341–348. [Google Scholar] [CrossRef]
- Hawbaker, T.J.; Keuler, N.S.; Lesak, A.A.; Gobakken, T.; Contrucci, K.; Radeloff, V.C. Improved estimates of forest vegetation structure and biomass with LiDAR-optimized sampling design.Geophys. Res.2009,114. [Google Scholar] [CrossRef]
- Maltamo, M.; Bollandsås, O.M.; Næsset, E.; Gobakken, T.; Packalén, P. Different plot selection strategies for field training data in ALS-assisted forest inventory.Forestry2011,84, 23–31. [Google Scholar] [CrossRef]
- Rouvinen, T. Kuvia metsästä.Metsätieteen Aikakauskirja2014,2, 119–122. [Google Scholar] [CrossRef]
- Rouvinen, T. Trestima—Digital Photographs for Forest Inventory.Sibirskij Lesnoj Žurnal (Siberian J. For. Sci.)2014,1, 69–76. [Google Scholar]
- R Core Team.R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018. [Google Scholar]
- Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D.nlme: Linear and Nonlinear Mixed Effects Models; R Package Version 3.1-137; Core Team: Vienna, Austria, 2018. [Google Scholar]
- Mehtätalo, L. A longitudinal height-diameter model for Norway spruce in Finland.Can. J. For. Res.2004,34, 131–140. [Google Scholar] [CrossRef]
- Näslund, M. Skogsförsöksanstaltens gallringsförsök i tallskog.Meddelanden från Statens skogsförsöksanstalt1936,29. Available online:https://pub.epsilon.slu.se/10159/1/medd_statens_skogsforskningsanst_029_01.pdf (accessed on 1 November 2017).
- Shevelev, S.L.; Kuzmichev, V.V.; Pavlov, N.V.; Smolyaninov, A.S.Лесoтаксациoнный справoчник для южнo-таежных лесoв Средней Сибири (Forest Inventory Handbook for the South Taiga Forests of Central Sibiria); M. VNIILM: Pushkino, Russia, 2002; 166p. [Google Scholar]
- Junttila, V.; Maltamo, M.; Kauranne, T. Sparse Bayesian estimation of forest stand characteristics from airborne laser scanning.For. Sci.2008,54, 543–552. [Google Scholar]
- Tipping, M.E. Sparse Bayesian learning and the relevance vector machine.J. Mach. Learn. Res.2001,1, 211–244. [Google Scholar]
- Haralick, R.M.; Shanmugam, K.; Dinstein, I. Textural features for image classification.IEEE Trans. Syst. Man Cybern.1973,SMC-3, 610–621. [Google Scholar] [CrossRef]
- McRoberts, R.E.; Tomppo, E.O. Remote sensing support for national forest inventories.Remote Sens. Environ.2007,110, 412–419. [Google Scholar] [CrossRef]
- Metsä Group. Metsä Group’s Owner-Members Can Update Their Growing Stock Information by Taking Photos with Their Mobile Phones. Available online:http://www.metsagroup.com/en/media/Pages/Case-Updating-growing-stock-information-by-taking-photos-with-mobile-phones.aspx (accessed on 14 October 2016).
- Holmgren, J.; Joyce, S.; Nilsson, M.; Olsson, H. Estimating stem volume and basal area in forest compartments by combining satellite image data with field data.Scand. J. For. Res.2000,15, 103–111. [Google Scholar] [CrossRef]
- Mäkelä, H.; Pekkarinen, A. Estimation of forest stand volumes by Landsat TM imagery and stand-level field-inventory data.For. Ecol. Manag.2004,196, 245–255. [Google Scholar]
- Hall, R.J.; Skakun, R.S.; Arsenault, E.J.; Case, B.S. Modeling forest stand structure attributes using Landsat ETM+ data: Application to mapping of aboveground biomass and stand volume.For. Ecol. Manag.2006,225, 378–390. [Google Scholar] [CrossRef]
- Päivinen, R.; Anttila, P. How reliable is a satellite forest inventory?Silva Fenn.2001,35, 125–127. [Google Scholar] [CrossRef]
- Puliti, S.; Saarela, S.; Gobakken, T.; Ståhl, G.; Næsset, E. Combining UAV and Sentinel-2 auxiliary data for forest growing stock volume estimation through hierarchical model-based inference.Remote Sens. Environ.2018,204, 485–497. [Google Scholar] [CrossRef]
Figure 1. Workflow diagram of the process. V = mean volume, G = basal area.
Figure 2. Location and delineation of the study area (green area). Rivers (blue) and roads (grey) are illustrated in the image at higher scale. Background map data source and copyright DeLorme/ESRI.
Figure 2. Location and delineation of the study area (green area). Rivers (blue) and roads (grey) are illustrated in the image at higher scale. Background map data source and copyright DeLorme/ESRI.
Figure 3. False color image presenting a leaf-on mosaic covering the area of interest (area delineated in yellow).
Figure 3. False color image presenting a leaf-on mosaic covering the area of interest (area delineated in yellow).
Figure 4. Scatter plots of NDVI and NIR values of random samples representing the variation of the inventory areas (black crosses) and the final sample (gray triangles) in the west and east banks of the river.
Figure 4. Scatter plots of NDVI and NIR values of random samples representing the variation of the inventory areas (black crosses) and the final sample (gray triangles) in the west and east banks of the river.
Figure 5. Image capture from the Trestima web service interface showing an example of a field sample plot shape and the samples’ visualization (green arrows correspond to basal area samples, and red arrows correspond to model tree diameter and height samples). The example plot includes several tree species, and therefore the number of model tree samples (red arrows) is large.
Figure 5. Image capture from the Trestima web service interface showing an example of a field sample plot shape and the samples’ visualization (green arrows correspond to basal area samples, and red arrows correspond to model tree diameter and height samples). The example plot includes several tree species, and therefore the number of model tree samples (red arrows) is large.
Figure 6. Example of model tree diameter and height measurements (Norway spruce).
Figure 7. Modeling areas and locations of measured field sample plots (red), and locations of not measured plots of the final sample (black). Note that the symbols of the measured sample plots cover the modeling area South at the image scale.
Figure 7. Modeling areas and locations of measured field sample plots (red), and locations of not measured plots of the final sample (black). Note that the symbols of the measured sample plots cover the modeling area South at the image scale.
Figure 8. Scatter plot of field-measured and estimated mean volume at plot level.
Table 1. Mean, standard deviation (std), and minimum (min) and maximum (max) of Landsat 8 Operational Land Imager (OLI) bands and normalized difference vegetation index (NDVI) values from random sample presenting whole AOI and from revised field sample locations (n = 792) for the east and the west bank of the river.
Table 1. Mean, standard deviation (std), and minimum (min) and maximum (max) of Landsat 8 Operational Land Imager (OLI) bands and normalized difference vegetation index (NDVI) values from random sample presenting whole AOI and from revised field sample locations (n = 792) for the east and the west bank of the river.
| Satellite Imagery Values | East Bank of Yenisei | West Bank of Yenisei | ||
|---|---|---|---|---|
| Inventory Area | Initial Sample | Inventory Area | Initial Sample | |
| Green mean | 7380 | 7245 | 7530 | 7470 |
| Green std | 382 | 279 | 316 | 284 |
| Green min | 6802 | 6887 | 6895 | 7025 |
| Green max | 22,388 | 9071 | 9865 | 9606 |
| Red mean | 6509 | 6361 | 6564 | 6443 |
| Red std | 438 | 340 | 324 | 275 |
| Red min | 6041 | 6092 | 6096 | 6148 |
| Red max | 23,393 | 9013 | 8982 | 8654 |
| NIR mean | 12,848 | 12,894 | 14,677 | 17,074 |
| NIR std | 1580 | 1761 | 2782 | 4348 |
| NIR min | 5855 | 7071 | 8590 | 9970 |
| NIR max | 23,684 | 18,600 | 26,962 | 31,418 |
| NDVI mean | 0.324 | 0.334 | 0.373 | 0.432 |
| NDVI std | 0.055 | 0.062 | 0.078 | 0.108 |
| NDVI min | −0.188 | 0.039 | 0.061 | 0.192 |
| NDVI max | 0.558 | 0.503 | 0.620 | 0.654 |
Table 2. Field data statistics for modeling areas. MV = mean value, Std = standard deviation, n = number of plots. D = mean diameter at breast height, H = mean tree height, V = mean volume, G = basal area, Vc = mean volume of the coniferous forest class, Vd = mean volume of the deciduous forest class, Vm = mean volume of the mixed coniferous forest class.
Table 2. Field data statistics for modeling areas. MV = mean value, Std = standard deviation, n = number of plots. D = mean diameter at breast height, H = mean tree height, V = mean volume, G = basal area, Vc = mean volume of the coniferous forest class, Vd = mean volume of the deciduous forest class, Vm = mean volume of the mixed coniferous forest class.
| Modeling Area | D, cm | H, m | V, m3 ha−1 | G, m2 ha−1 | Vc, m3 ha−1 | Vd, m3 ha−1 | Vm, m3 ha−1 | |||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MV | Std | MV | Std | MV | Std | MV | Std | MV | Std | MV | Std | MV | Std | |
| South | 21.1 | 5.6 | 18.8 | 3.3 | 169.9 | 69.9 | 20.2 | 7.8 | 249.9 | 64.0 | 140.1 | 54.6 | 192.6 | 69.3 |
| West | 21.2 | 4.3 | 17.2 | 2.9 | 195.5 | 50.3 | 24.1 | 5.2 | 187.9 | 53.1 | 169.3 | 65.6 | 202.0 | 46.9 |
| Southwest | 21.1 | 5.0 | 17.6 | 3.5 | 186.7 | 61.6 | 22.8 | 6.7 | 202.8 | 62.7 | 141.9 | 51.0 | 201.6 | 54.6 |
| East | 19.8 | 3.7 | 16.1 | 2.7 | 223.5 | 44.3 | 26.4 | 4.7 | 231.4 | 32.2 | 210.3 | 41.3 | 226.8 | 47.6 |
| East, shadow1 | 21.1 | 5.0 | 15.9 | 3.0 | 205.0 | 41.7 | NA | NA | NA | NA | NA | NA | NA | NA |
1 Forest class specific models were not estimated.
Table 3. Summary of plot-level validation results (total values). RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
Table 3. Summary of plot-level validation results (total values). RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
| Modeling Area/Model | n | Basal Area, m2 ha−1 | Volume, m3 ha−1 | ||||||
|---|---|---|---|---|---|---|---|---|---|
| RMSE | RMSE% | Bias | Bias% | RMSE | RMSE% | Bias | Bias% | ||
| South | 210 | 4.1 | 20.3 | −0.3 | −1.3 | 39.2 | 23.6 | −1.2 | −0.7 |
| West | 110 | 4.3 | 18.0 | 0.1 | 0.3 | 40.3 | 20.6 | 2.4 | 1.2 |
| Southwest | 176 | 4.3 | 19.0 | 0.6 | 2.5 | 42.6 | 22.8 | 6.1 | 3.3 |
| East, model 11 | 155 | 4.6 | 17.2 | 0.9 | 3.2 | 42.2 | 18.9 | 10.7 | 4.8 |
| East, model 22 | 155 | 4.6 | 17.5 | 0.7 | 2.5 | 43.5 | 19.5 | 8.1 | 3.6 |
| East, model 33 | 16 | 3.8 | 16.0 | 1.4 | 4.6 | 34.5 | 16.8 | 3.7 | 1.8 |
| East, model 44 | 16 | 4.3 | 18.2 | −0.3 | −1.1 | 40.0 | 19.5 | −0.2 | −0.1 |
1 Model with leaf-on and leaf-off images.2 Model with leaf-on images.3 Model with leaf-on and leaf-off images in shadow area.4 Model with leaf-on images in shadow area.
Table 4. Summary of plot-level validation results in the coniferous forest class. RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
Table 4. Summary of plot-level validation results in the coniferous forest class. RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
| Modeling Area/Model | n | Basal Area, m2 ha−1 | Volume, m3 ha−1 | ||||||
|---|---|---|---|---|---|---|---|---|---|
| RMSE | RMSE% | Bias | Bias% | RMSE | RMSE% | Bias | Bias% | ||
| South | 21 | 5.1 | 17.3 | −0.8 | −2.7 | 48.2 | 19.3 | −4.7 | −1.9 |
| West | 39 | 4.6 | 19.5 | −0.7 | −2.9 | 43.0 | 22.9 | −3.3 | −1.8 |
| Southwest | 47 | 4.5 | 18.2 | −0.8 | −3.2 | 41.7 | 20.6 | −3.0 | −1.5 |
| East, model 11 | 26 | 3.2 | 12.1 | −0.0 | −0.1 | 32.8 | 14.2 | −1.8 | −0.8 |
| East, model 22 | 26 | 3.1 | 11.8 | 0.0 | 0.1 | 32.4 | 14.0 | −2.4 | −1.0 |
1 Model with leaf-on and leaf-off images.2 Model with leaf-on images.
Table 5. Summary of plot-level validation results in the mixed coniferous forest class. RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
Table 5. Summary of plot-level validation results in the mixed coniferous forest class. RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
| Modeling Area/Model | n | Basal Area, m2 ha−1 | Volume, m3 ha−1 | ||||||
|---|---|---|---|---|---|---|---|---|---|
| RMSE | RMSE% | Bias | Bias% | RMSE | RMSE% | Bias | Bias% | ||
| South | 59 | 4.9 | 20.6 | −0.7 | −2.8 | 48.6 | 25.2 | −4.6 | −2.4 |
| West | 66 | 4.1 | 16.5 | 0.4 | 1.8 | 36.8 | 18.2 | 5.3 | 2.6 |
| Southwest | 84 | 4.6 | 18.9 | 1.3 | 5.2 | 46.7 | 23.2 | 12.2 | 6.0 |
| East, model 11 | 91 | 4.5 | 17.0 | 0.9 | 3.5 | 42.8 | 18.8 | 13.4 | 5.9 |
| East, model 22 | 91 | 4.6 | 17.5 | 0.6 | 2.2 | 42.6 | 18.8 | 8.7 | 3.8 |
1 Model with leaf-on and leaf-off images.2 Model with leaf-on images.
Table 6. Summary of plot-level validation results in the deciduous forest class. RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
Table 6. Summary of plot-level validation results in the deciduous forest class. RMSE = root mean square error, RMSE% = error relative to mean value of plots, n = number of plots used in validation.
| Modeling Area/Model | n | Basal Area, m2 ha−1 | Volume, m3 ha−1 | ||||||
|---|---|---|---|---|---|---|---|---|---|
| RMSE | RMSE% | Bias | Bias% | RMSE | RMSE% | Bias | Bias% | ||
| South | 130 | 3.5 | 20.4 | 0.0 | 0.0 | 32.2 | 23.0 | 0.9 | 0.6 |
| West | 5 | 5.5 | 24.9 | 1.2 | 5.3 | 59.3 | 35.0 | 8.9 | 5.2 |
| Southwest | 45 | 3.5 | 20.0 | 0.7 | 3.8 | 34.8 | 24.5 | 4.3 | 3.0 |
| East, model 11 | 38 | 5.5 | 20.2 | 1.3 | 4.9 | 46.5 | 22.1 | 12.6 | 6.0 |
| East, model 22 | 38 | 5.5 | 20.2 | 1.3 | 4.9 | 51.6 | 24.6 | 13.9 | 6.6 |
1 Model with leaf-on and leaf-off images.2 Model with leaf-on images.
Table 7. Summary of compartment-level validation results. RMSE = root mean square error, RMSE% = error relative to mean value of compartments, n = number of plots used in validation.
Table 7. Summary of compartment-level validation results. RMSE = root mean square error, RMSE% = error relative to mean value of compartments, n = number of plots used in validation.
| Modeling Area | n | Basal Area, m2 ha−1 | Volume, m3 ha−1 | ||||||
|---|---|---|---|---|---|---|---|---|---|
| RMSE | RMSE% | Bias | Bias% | RMSE | RMSE% | Bias | Bias% | ||
| South | 56 | 2.7 | 13.4 | 0.6 | 2.9 | 24.4 | 14.8 | 5.9 | 3.6 |
| West | 32 | 3.2 | 13.3 | 0.4 | 1.8 | 24.9 | 13.1 | 3.3 | 1.7 |
| Southwest | 80 | 3.9 | 17.7 | 0.8 | 3.7 | 36.3 | 20.3 | 6.6 | 3.7 |
| East | 43 | 2.4 | 9.3 | 0.6 | 2.3 | 22.4 | 10.1 | 7.1 | 3.2 |
Table 8. Land use/land cover error matrix for maximum likelihood classification. Reference values in horizontal and mapped values in vertical.
Table 8. Land use/land cover error matrix for maximum likelihood classification. Reference values in horizontal and mapped values in vertical.
| Land Use/Land Cover Class | Conifer | Mixed Conifer | Deciduous | Swamp Grass | Harvested | Regrowth | OPEN LAND | Water | Burnt Area |
|---|---|---|---|---|---|---|---|---|---|
| Conifer | 26 | 15 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| Mixed conifer | 13 | 51 | 5 | 1 | 0 | 0 | 1 | 0 | 0 |
| Deciduous | 0 | 13 | 39 | 0 | 0 | 0 | 1 | 0 | 0 |
| Swamp grass | 0 | 0 | 0 | 7 | 0 | 0 | 0 | 0 | 0 |
| Harvested | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 |
| Regrowth | 0 | 1 | 2 | 0 | 0 | 5 | 0 | 0 | 0 |
| Open land | 4 | 1 | 1 | 0 | 0 | 0 | 4 | 0 | 0 |
| Water | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
| Burnt area | 3 | 1 | 2 | 0 | 0 | 3 | 0 | 0 | 4 |
© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Share and Cite
MDPI and ACS Style
Rybakov, G.; Peuhkurinen, J.; Latva-Käyrä, P.; Villikka, M.; Sirparanta, S.; Kolesnikov, A.; Junttila, V.; Kauranne, T. Combining Camera Relascope-Measured Field Plots and Multi-Seasonal Landsat 8 Imagery for Enhancing the Forest Inventory of Boreal Forests in Central Russia.Remote Sens.2018,10, 1796. https://doi.org/10.3390/rs10111796
AMA Style
Rybakov G, Peuhkurinen J, Latva-Käyrä P, Villikka M, Sirparanta S, Kolesnikov A, Junttila V, Kauranne T. Combining Camera Relascope-Measured Field Plots and Multi-Seasonal Landsat 8 Imagery for Enhancing the Forest Inventory of Boreal Forests in Central Russia.Remote Sensing. 2018; 10(11):1796. https://doi.org/10.3390/rs10111796
Chicago/Turabian StyleRybakov, Georgy, Jussi Peuhkurinen, Petri Latva-Käyrä, Maria Villikka, Sanna Sirparanta, Alexander Kolesnikov, Virpi Junttila, and Tuomo Kauranne. 2018. "Combining Camera Relascope-Measured Field Plots and Multi-Seasonal Landsat 8 Imagery for Enhancing the Forest Inventory of Boreal Forests in Central Russia"Remote Sensing 10, no. 11: 1796. https://doi.org/10.3390/rs10111796
APA StyleRybakov, G., Peuhkurinen, J., Latva-Käyrä, P., Villikka, M., Sirparanta, S., Kolesnikov, A., Junttila, V., & Kauranne, T. (2018). Combining Camera Relascope-Measured Field Plots and Multi-Seasonal Landsat 8 Imagery for Enhancing the Forest Inventory of Boreal Forests in Central Russia.Remote Sensing,10(11), 1796. https://doi.org/10.3390/rs10111796
Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further detailshere.
Article Metrics
No
No
Article Access Statistics
For more information on the journal statistics, clickhere. Multiple requests from the same IP address are counted as one view.
