An Integrated GIS and Remote Sensing Approach for Monitoring Harvested Areas from Very High-Resolution, Low-Cost Satellite Images

Abstract

Advanced monitoring and mapping of forest areas using the latest technological advances in satellite imagery is an alternative solution for sustainable forest management compared to conventional ground measurements. Remote sensing products have been a key source of information and cost-effective options for monitoring changes in harvested areas. Despite recent advances in satellite technology with a broad variety of spectral and temporal resolutions, monitoring the areal extent of harvested forest areas in managed forests is still a challenge, primarily due to the highly dynamic spatiotemporal patterns of logging activities. Our goal was to introduce a plot-based method for monitoring harvested forest areas from very high-resolution (VHR), low-cost satellite images. Our method encompassed two data categories, which included vegetation indices (VIs) and texture analysis (TA). Each group of data was used to model the amount of harvested volume both independently and in combination. Our results indicated that the composition of all spectral bands can improve the accuracy of all models of average volume by 23.52 RMSE reduction and total volume by 33.57 RMSE reduction. This method demonstrated that monitoring and extrapolation of the calculated relation and results from smaller forested areas could be applied as an automatic remote-based supervised monitoring method over larger forest areas.

1. Introduction

Knowing the extent of forest cover disturbances related to anthropogenic interventions, particularly timber harvesting, over large managed areas is important to assess: (1) the construction or update of forest inventories [1,2]; (2) forest health and productivity [3,4]; (3) risks of climate change [5]; (4) carbon flux rates [6]; and (5) support required for forest management plans (FMPs) and decision making [7]. Despite the economic importance of the timber industry in many regions, logging activities on a global scale have resulted in significant ecological damage due to excessive timber volume extraction, which in some cases is considered to be irreversible [8].

It is evident that logging has multiple effects on forest composition and structure [9]. Logging not only decreases the tree density [10] and basal area [11], but also regulates canopy gaps dynamics with altered light conditions [12]. Even-aged forest stand management commonly leads to the final harvest, by removing the mature stand, either by applying several shelterwood cutting phases or in one-cut (clear-cut). On the contrary, selective logging, and uneven-aged stand management never decrease the growing stock to zero. Selective logging is a source of anthropogenic disturbance and therefore is regarded as forest degradation. Approximately, 50%–90% of forest cover is maintained after selective logging activities [13,14], crown openness can cause significant environmental consequences, like for instance semantic timber stock reduction [15].

More traditional methods of ground truth data collection to support forest resource management, such as field samples, are usually expensive and time-consuming. Also, there are challenges related to the size of study areas, accessibility, and the capacity to conduct repeated measurements over time (regular time-series datasets). However, the collection of ground truth data can help to calibrate other sources of data, such as remote sensing data using geographic information systems (GIS) and also support the interpretation and analysis of satellite imagery and aerial photography, among other methods of remote data collection.

To alleviate these challenges, satellite remote sensing (SRM) is an economical and practical solution to consistently monitor forest cover changes over large geographical areas [16]. The detection and monitoring of deforested areas using satellite images are essential for forest managers because they allow for the creation of regular time-series datasets with different spatial resolutions. These datasets can be used for multiple purposes, including forest cover and change detection [17,18] and estimation of structural (e.g., diameter at breast height, height, stem volume) and textural metrics, with acceptable accuracy [19,20,21,22,23]. In addition to volume estimation, satellite images can be used to estimate the area of forest classes, for which volume can be extrapolated using field-based measurements [24].

Some of the most significant challenges to mapping and monitoring forested areas are related to: (1) highly dynamic spatiotemporal patterns of logging activities, where SRM is able to detect changes in patterns for only short period of time, because of the quick canopy closure after silviculture activities [8]; (2) high levels of complexity inherent in forest ecosystems [25]; and (3) timber harvesting techniques of differing intensities and patterns [26]. Indeed, logging impacts can vary significantly depending on the intensity and the frequency of logging activities, since they can vary from very low, to more intensive and repeated logging activities.

Consistent or long-term remote-based monitoring is the key to dynamically discriminating and monitoring degraded forest areas either caused by biotic (e.g., pest infestation) or abiotic factors (e.g., detection of illegal logging activities, or as a control method for the evaluation of silvicultural treatment in larger regions. It is also suitable as support for terrestrial data collection for management plan elaboration, forest inventory and monitoring of carbon sequestration. Even though a single-date image [27] can identify the type of degradation, however, a time-series satellite is needed in order to better detect and assess forest cover dynamics [28].

To address the dynamics of monitoring vegetation by optical means, the best compromise is using sensors with a large field of view (FOV), such as MODIS or NOAA/AVHRR (these sensors combine moderate spatial resolution with high temporal resolution/revisit at a rate of approximately one per day); sensors with finer spatial resolution (e.g., Landsat ~30 m) or very high spatial resolution (e.g., WorldView-4 ~0.3 m, Pleiades’ (1A/1B) ~0.5 m, and Quick bird ~1m) tend to have revisit periods of several weeks. For example, the monitoring roles and mapping of timber harvested areas in Brazilian Amazon has largely relied on Landsat images [25]. The availability of large volumes of free images, sufficient computational capacity, and advances in image processing techniques have stimulated widespread use of Landsat imagery for a broad range of uses, including assessment of vegetation phenology [29] and the classification and assessment of forest changes [30]. In another study, Toutanova et al. [31], used Landsat-based archive to assess long-term regional forest dynamics in eastern Europe [31]. Their results clearly showed that long-term monitoring from remote sensing-based products that employ biophysical criteria of forest cover (defined according to tree canopy cover thresholds) is better suitable for monitoring of timber harvesting than forest inventory data. On the other hand, the uses of finer spatial resolution satellite images can significantly enhance the detection of harvested areas, but there may be computational constraints in many cases and they are typically costly [8]. In some cases, there are also geographical territories with limited or no access to satellite data for a particular time period.

Spectral readings are based on the interaction between objects and incoming radiation; this interaction allows for the detection and discrimination of target objects because of the characteristic signatures of reflected light, as all objects return different amounts of energy in different bands of the electromagnetic spectrum. These signatures are essentially functioning of vegetation light absorption properties at given wavelengths that can be used to estimate spectral indices [32].

Numerous vegetation indices (VIs) utilizing spectral information as a proxy for vegetation productivity, monitoring, and mapping have been developed; the normalized difference vegetation index (NDVI) is the most common index for high-throughput screening [33]. Using the visible and near-infrared (NIR) parts of the electromagnetic spectrum, NDVI estimates the fraction of the photosynthetically active radiation (FPAR) absorbed by the forest canopy [34]. The enhanced vegetation index (EVI) is an optimized index developed to reduce the saturation of the reflectance signal at increasing levels of green biomass, and diminish the noise produced by atmospheric aerosols and soil background [34]. Huene [35] intended to minimize the effects of soil background on the vegetation signal by incorporating a constant soil adjustment factor (L) into the denominator of the NDVI equation. The customized NDVI equation was named soil-adjusted vegetation index (SAVI). The adjustment factor (L) value varies based upon the soil reflectance characteristics, such as brightness and color, among others, and the appropriate value depends exclusively on vegetation density. For example, in the case of high densities, L = 0.25; for intermediate densities, L = 0.5; and for low vegetation densities, L = 1.0 [36]. However, Eastman [37] proposed that the best value for L is the point at which the difference between SAVI for light and dark soil is minimum. Therefore, for L = 1, SAVI approaches perpendicular vegetation index (PVI), which uses the perpendicular distance from each pixel coordinate to the soil line, while for L = 0, SAVI is equal to NDVI [37]. Many studies have identified the usefulness of long series of vegetation indices (mainly NDVI) derived from satellite data for monitoring purposes of forest degradation, but most have been based on low-to medium-resolution satellite images (e.g., [38,39,40]).

In this context, our main goal and novelty of this paper relies on integrating GIS and remote sensing, to utilize all potential information from very high-resolution (VHR) satellite images to help examine (i) if multi-spectral VHR, low-cost satellite images, can be used to successfully detect and estimate the amount of harvested trees vs living trees based on satellite image time series (SITS) from two different periods before and after the logging activities, (ii) the possibility to create a continuous map of harvested areas (iii) which composition of independent data can increase the model accuracy.

2. Materials and Methods

2.1. Characterization of the Study Area

Our study area was located in the city of Doksy on the shores of Lake Mácha in northern Bohemia in the Czech Republic; it is largely surrounded by dense forest area, covering approximately 300 km². The entire study area is five hectares in a flat terrain, extending geographically from 50°33’26.48”N; 14°43’33.72”E to 50°33’19.88”N; 14°43’13.95” E at 310 m a.s.l. We selected 48 circular permanent plots with a radius of 12.6 m representing different harvest intensities in mature forest stands (mean = 113 years; Figure 1). All plots were dominated by Scots pine (Pinus sylvestris L.) with admixtures of Norway spruce (Picea abies (L.) Karst.). The experimental plots were originally established for monitoring of natural regeneration success of Scots pine in different stand densities. For this purpose, stripes of shelterwood cutting of different intensity using harvester technology were established. Harvests were conducted in March 2017 and a total wood volume of 680.5 m³ was harvested from the plots.

Parent rock in the area was formed by late Cretaceous sandstone, with dominant soil types being Arenic Cambisol and Arenic Podzol [41]. The vegetative period is typically rather warm and dry. Mean annual air temperature is 7.3 °C and average maximum temperature is 31.5 °C. Mean annual precipitation is 635 mm, with only 354 mm during the growing season (1 April to 30 October). Mean annual number of cloudy days ranges from 120 to 150, of which 50 to 60 occur during the growing season [42].

2.2. Terrestrial Data Acquisition and Processing

To determine the center of each circular plot, we used a Leica TCRP 1201+ robotic total station using the polar coordinate method and the GNSS GPS Leica system 1200 (Leica Geosystems AG: Heerbrugg, Switzerland) with centimeter accuracy. Within each circular plot, we measured in two perpendicular directions the diameter at breast height (1.3 m; DBH_1.3) of standing trees and the stump diameter of harvested trees with a computer caliper DP II (accuracy 1 mm; Haglöf Inc., Sweden). Forty trees of each species were measured separately to create a regression equation for estimating the DBH_1.3 of harvested trees. Similarly, 25 trees of each species were also sampled to evaluate tree heights using the Vertex IV hypsometer (accuracy 0.1 m; Haglöf Inc., Sweden) (Table 1).

The stem volume was calculated using DBH_1.3 and tree height as predictors based on volumetric equations (Equations (1) and (2); Table 2) [43] by Petráš and Pajtík (Figure 2, Figure 3 and Figure 4). The field plot data were used as reference data to evaluate the extrapolation process of harvested volumes in the form of a continuous map using extracted spectral data for all plots.

All the processing and statistical analyses of the terrestrial data were conducted in Matlab R2017b professional edition (MathWorks©, Inc.; Massachusetts, U.S.A.).

(1)$\begin{array}{ll}{V}_{Scotspine}=& \left(A\ast DB{H}_{1.3}\left(\left(B-C\ast log10\left(DB{H}_{1.3}+1\right)\right)\ast H^E\right)\right)\\ & -\left(A_1\ast {\left(DB{H}_{1.3}+1\right)}^{B_1\ast H^{C_1}}\right)\end{array}$

(2)$V_{Norwayspruce}=\left(A\ast {\left(DB{H}_{1.3}+1\right)}^{B\ast H^C}\right)-\left(E\ast {\left(DB{H}_{1.3}+1\right)}^{F\ast H^G}\right)$

2.3. Acquisition and Pre-Processing of Satellite Images

We used two satellite images by Pléiades-HR 1A and Pléiades-HR 1B, which included an image from 2016 (pre-harvest) and one from 2019 (post-harvest); the images were taken on 27 March 2016 and 6 April 2019, respectively. The image consisted of one panchromatic (PAN) and four multi-spectral (red = R, green = G, blue = B, and NIR) bands (Table 3). To achieve accurate georeferencing of the satellite images, we used 20 ground control points (GCPs) randomly distributed around the study area; these points were measured using the Leica real-time kinematic (RTK) system, model RX1250XC (Leica Geosystems AG; Heerbrugg, Switzerland) with centimeter accuracy. The RTK correction was conducted using a base/rover set, which sends and receives fast-rate over-the-air RTK data corrections, using the PDLGFU15 radio module. The entire phase of satellite image pre-processing was conducted in ArcGIS desktop V.10.6.1 (ESRI Inc.; Redlands, CA, U.S.A.) and PCI geomatica 2016 (PCI Geomatics; Markham, ON, Canada).

2.4. Satellite Image Processing

To develop a deeper understanding of the spatial and temporal changes in patterns caused by timber harvesting, we applied several different VIs to evaluate the amount of vegetation cover change for the period between 2016 and 2019. We categorized the VIs into two major groups, denoted as the slope-based (Figure 5 and Figure 6) and distance-based (Figure 7 and Figure 8) groups (Table 4). These indices were calculated by simply dividing the reflectance values of the red band by those of the NIR. Since soil effect is found to have a considerable impact in sparsely-vegetated areas [44], we measured the soil line slope (Figure 9) to create the SAVI index using algebraic equations of the NIR and R bands (Table 4). The SAVI index was used to suppress the soil effect of the sandstone nature which was dominant in the study area.

We conducted one-way ANOVA analyses (p-value < 0.01) to test the differences in VI values in 2016 and 2019 using Matlab R2017b professional edition (Figure 6 and Figure 8; Table 5).

The PCI geomatica software (PCI Geomatics; Markham, ON, Canada) was used to create the texture analysis (TA) for each band separately. We used 13 TA variables, including mean, variance, entropy, contrast, heterogeneity, homogeneity, angular second moment, correlation, gray-level difference vector (GLDV), angular second moment, GLDV entropy, GLDV mean, GLDV contrast, and inverse difference with kernel window size 12 × 12 for RGB and NIR. For each year we included 52 TA data layers (Figure 10).

2.5. Modeling and Extrapolation of Harvested Parameters

We initially used a feature selection algorithm to identify the most important layers in the modeling process for each dependent variable separately (Supplementary Tables S1–S3). The most significant layers are presented in Table 6 which includes both VI and TA layers.

Using regression analysis as part of random forest (RF) algorithm in STATISTICA software V.12 (Dell Inc.; Round Rock, TX, U.S.A.), we modeled the variety of dependent variables (Table 7). For the application of the RF algorithm, we considered the following options, for random test data: 30%, subsample proportion: 50%, minimum n of cases: 5, minimum n in child node: 5, maximum n of levels: 20, maximum n of nodes: 100, and number of trees: 1000.

For evaluating the model’s accuracy, root mean square error (RMSE) was calculated as below

(3)$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(X_{obs,i}-X_{model,i}\right)}^2}{n}}$

where ${\mathit{X}}_{\mathit{o}\mathit{b}\mathit{s}}$ is observed values and ${\mathit{X}}_{\mathit{m}\mathit{o}\mathit{d}\mathit{e}\mathit{l}}$ is modeled values at time/place i.

Equation (4) was also used to convert the RMSE values to a non-dimensional or normalized form of RMSE

(4)$\mathit{R}\mathit{M}\mathit{S}\mathit{E}\left(\mathit{r}\mathit{e}\mathit{d}\mathit{u}\mathit{c}\mathit{t}\mathit{i}\mathit{o}\mathit{n}\right)=\frac{\mathit{RMSE}}{{\overline{\mathit{X}}}_{\mathit{obs}}}\times \mathbf{100}$

where ${\overline{\mathit{X}}}_{\mathit{obs}}$ is mean of observed values. After the evaluation of the modeled parameters, we created a continuous map of harvested parameters (Figure 11).

To show the differences between observed and estimated parameters we used bias and bias% as below

(5)$Bias=\frac{1}{n}{\displaystyle \sum }{x}_i-{\overline{x}}_j$

(6)$Bias\%=\frac{Bias}{\overline{x_j}}\times 100$

where $n$ is the number of samples, $x_i$ is the observed value and ${\overline{x}}_j$ is mean of estimated values.

3. Results

3.1. Vegetation Indices

We created several VI maps to contribute in modeling the differences between the forest parameters pre- and post-harvest. To measure VIs, including the SAVI, MSAVI 1, and MSAVI 2, the slope of soil line, known as the L factor, was required (Table 4). As evident in Figure 8, the slope of the soil line differed for 2016 and 2019, thus we calculated the L factor for each year separately in the form of a raster layer.

The VI values calculated for pre- and post-harvest images were significantly different between years (Figure 6 and Figure 8; Table 5). Among all VIs, both the slope-based and distance-based, RVI and NRVI indices were the only indices that significantly increased post-harvest.

3.2. Texture Analysis

Figure 10 presents the results obtained from the processing of multi-spectral bands with 13 different descriptive statistic algorithms using a kernel size of 12 × 12 pixels. Since irrelevant or partially relevant features can negatively impact model performance, a feature selection algorithm was used to ensure that all independent variables were significantly related to the dependent variables (Supplementary Tables S1–S3). Some of the most significant TA and VI layers are presented in Table 6. Due to feature selection algorithm results, all the TA layers had a significant impact on modeling the input variables. As a result, non of the created layers were excluded from the modeling process.

3.3. Modeling and Extrapolation

After the preparation of all VIs and TA layers, the RF algorithm was used to model changes in several forest parameters (Table 7).

Using each group of independent values independently and in combination, we were able to model the relation between the living trees and spectral information for the years 2016 and 2019 separately. In most cases, the combination of VIs and TA data increased the accuracy of the models, and it increased by as much as 20% when modeling the changes due to harvesting. For example, the fusion of VIs and TA data increased the accuracy of the modeling of the total harvested volumes (m³) per plot by almost 20% (Table 7).

Based on the equations estimated using RF, we constructed maps of continuous changes of different forest parameters. Based on our plot-based study design, the output cell size was equal to the plot size (25 × 25 m) (Figure 11).

4. Discussion

We demonstrated the potential of low-cost satellite imagery to detect and model structural changes in forested areas [13,51]. Based on our analyses, all VIs, with the exceptions of RVI and NRVI, exhibited decreased values after the harvest operations [51]. Forest harvesting contributed to lower absorption capabilities of red light (less trees, less chlorophyll), thus resulting in higher values in the red band, and lower reflectance values in NIR (fewer trees, fewer leaves), thus resulting in lower values in the NIR band [52]. Consequently, VI values significantly decreased; our results clearly identified the contrasts between the red and NIR bands for vegetated pixels, with high VI values produced by the combinations of low red (due to the low absorption rates by chlorophyll) and high NIR (as a result of the decreased leaf density) reflectance. Because RVI is the inverse of the standard ratio VI, it was no surprise that RVI values increased post-harvest.

The results of the feature selection algorithm demonstrated that TA layers were significantly related to the dependent values (Figure 10, Supplementary Tables S1–S3), as evident in Table 6 Although the fusion of VIs and TA layers exhibited the highest model accuracy in most cases, using solely TA layers provided reliable accuracy in modeling the forest parameters. This proves that the harvesting operation caused significant changes in the crown canopy texture, which led to significant differences in TA values pre- and post-harvesting. The results from TA are in accordance with the result of our previous study [23], where TA had a significant performance on spatial detection of dominant tree species (i.e., Fagus orientalis L.), that had the highest influence on canopy structure. That verifies the capability of TA in modeling variables with significant influence on canopy structure.

The RF algorithm modeled the total volume and average volume (m3) per plot for the years 2016 and 2019 with high accuracy (Table 7). The RMSE reduction for the average volume total for 2016 and 2019 was 16.04% and 16.97%, respectively (Table 7). From all estimated parameters, the total harvested volume was the only one with the lowest accuracy (RMSE reduction 33.57%) comparing to the rest forest parameters (Table 7).

In addition, using a combination of VIs and TA values, RF could model changes in stem volume and crown density related to harvesting activities. Models were able to detect the harvested areas within plots (i.e., those that experienced selective logging) as well as outside the plots. Although it is beyond the scope of this study, it is worthwhile to note that the RF algorithm also successfully detected the clear-cut area [53] located south of the area of interest (Figure 11). The detection ability of our method allows for the successful detection of clear-cut areas, and, more importantly, it demonstrates an ability to detect areas that have been selectively logged. Our findings are consistent with those of two other studies [52,54], where the authors tried to compare IKONOS 1-m and 4-m resolution and Sentinel-2 MSI images with Landsat 7 ETM+ images to identify the significance of higher resolution sensors compared to lower or medium spatial resolution sensors to adequately identify areas of selective logging. Our study revealed the potential of our method to monitor larger forested areas with high reliability because it allows for the extrapolation of the results (local observations) to larger areas.

5. Conclusions

Our main objective was to determine if multi-spectral, VHR, time-series satellite images could be used to successfully detect and estimate the amount of harvested and living trees, and also create a continuous map of the harvested areas. In addition, we aimed to define which composition of independent data could be used to increase mapping accuracy. Our results demonstrated how low-cost satellite images can be successfully used for the detection of forest structural changes (i.e., volume, canopy density, stand density, etc.). Consequently, this method can be used to monitor volume dynamics caused by anthropogenic forest disturbances. Our findings also showed that estimated variables can be applied as an input for remote sensing-based supervised extrapolation over large areas, without the need for terrestrial interventions. For example, non-accessible areas such as national parks or private forest landowners could benefit by using this method for monitoring purposes. Because the proposed method is also based on the spectral and structural changes that mainly occur in forest canopies, it can also be used for the assessment of abiotic stress factors, such as windthrow and post-fire monitoring. Thus, when it comes to cost-effective monitoring of forest ecosystems over larger areas, the proposed method is considered to be a practical, data sufficient, and low-cost solution for the structural detection and estimation of the amount of harvested areas.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/11/21/2539/s1. Table S1: Correlation between spectral and terrestrial data pre-harvest (2016). Table S2: Correlation between spectral and terrestrial data post-harvest (2019). Table S3: Most correlated layers to changes in vegetation structure between 2016 and 2019 based on a feature extraction algorithm.

Author Contributions

A.A. conceived the idea for this manuscript; D.P. and A.A. designed and wrote the manuscript; A.A. processed the remote sensing and terrestrial data. D.P. and A.A. designed and performed the statistical analyses; L.B. provided all the terrestrial data and helped in finalizing the manuscript.

Funding

We acknowledge that this work was financially supported by the following projects: (1) “Advanced research supporting the forestry and wood-processing sector’s adaptation to global change and the 4th industrial evolution”—shortened project name: EVA4.0 [Grant No. CZ.02.1.01/0.0/0.0/16_019/0000803] of the Faculty of Forestry and Wood Sciences from the Czech University of Life Sciences in Prague and (2) the Ministry of Agriculture of the Czech Republic [No. QJ1520037].

Acknowledgments

Authors would like to thank Peter Surový from the department of forest management for the satellite data acquisition support and the Military Forests and Farms of the Czech Republic, state enterprise (VLS ČR, s. p.) for permission and assistance in establishing the study area. Also, the authors would like to thank Jakub Roztočil and Martin Pros for their contribution to the collection of terrestrial data.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Geographic locations of the study area and plots.

Figure 2. Shows the average amount of volume for harvested and living trees.

Figure 3. Visual comparison of the average and total volume before (2016) and after harvest (2019) using box-plots.

Figure 4. Shows the different harvesting intensities applied in each plot in the form of percentages. The red plots inside the dotted line frame are not part of the study and are being used for illustration purposes only.

Figure 5. Application of several slope-based vegetation indices (VIs) obtained from remote sensing-based canopies.

Figure 6. Visual comparisons of the group location parameters using boxplots for slope-based vegetation indices (VIs) for 2016 and 2019. The letters with red font indicate that there were significant differences between 2016 and 2019 (p-value < 0.01).

Figure 7. Application of distance-based vegetation indices obtained from remote sensing-based canopies.

Figure 8. Visual comparisons of the group location parameters using boxplots for distance-based VIs for 2016 and 2019. The letters with red font indicate that there were significant differences between 2016 and 2019 (p-value < 0.01).

Figure 9. Bi-spectral plot of NIR reflectance versus red reflectance of all pixels between 2016 and 2019, as described in Section 2.4. The straight red and black lines corresponds to the bare soil for this particular dataset. Both datasets were plotted together using the same x-axis and y-axis values so that the differences are clearly evident.

Figure 10. Texture analysis applied in NIR band for the years 2016 and 2019.

Figure 11. Continuous maps of forest parameters changes related to harvesting operation.

Table 1

Shows the used equations and their accuracies for estimation of height and DBH_1.3 for each species separately. The below equations were estimated using the relation between (i) stump diameter and DBH_1.3 (ii) DBH_1.3 and height.

Estimated Parameter	Tree Species	Number of Trees	Trendline	Equation	R2
Height	Scots pine	25	linear	y = 0.1581x + 21.141	0.48
	Norway spruce	25	linear	y = 0.2354x + 17.145	0.58
$DB{H}_{1.3}$	Scots pine	40	linear	y = 1.0601x + 5.7583	0.92
	Norway spruce	40	linear	y = 1.29x - 0.165	0.96

Table 2

Shows the parameters used for the stem volume equations for both tree species.

Scots Pine	A	B	C	E	${\mathit{A}}_1$	${\mathit{B}}_1$	${\mathit{C}}_1$
Coefficients	3.03E-05	2.075238	0.012492	0.961028	0.071975	−2.12449	1.372591
Norway spruce	A	B	C	E	F	G
Coefficients	4.01E-05	1.821816	1.132062	9.29E-03	−1.02037	0.896101

Table 3

Characteristics of Pléiades-HR 1A and Pléiades-HR 1B satellite sensors.

Pléiades HR 1A-1B
Imagery products	2 m multi-spectral resolution
	Bundle: 50 cm panchromatic—2 m multi-spectral
	Panchromatic resolution 50 cm
Spectral bands	Red: 600–720 nm
	Green: 490–610 nm
	Blue: 430–550 nm
	Near Infrared: 750–950 nm
	Panchromatic: 480–830 nm
Image location accuracy	With * GCPs is 1 m
	Without GCPs is 3 m (CE90)

* GCPs: ground control points.

Table 4

Equations of applied vegetation indices (VIs).

	VIs	References	Explanation of Symbols
Slope-based	$\ast NDVI=\frac{NIR-R}{NIR+R}$	Rouse et al. [45]	R = red bandB = blue bandL = soilNIR = near infrared band
	$\ast CTVI=\frac{NDVI+0.5}{ABS\left(NDVI+0.5\right)}$ × $\sqrt{ABS\left(NDVI+0.5\right)}$	Perry and Lautenschlager [46]
	$\ast NRVI=\frac{RVI-1}{RVI+1}$	Baret and Guyot [47]
	$\ast RATIO=\frac{NIR}{R}$	Birth and McVey [48]
	$\ast RVI=\frac{R}{NIR}$	Richardson and Wiegand [49]
	$\ast SAVI=\frac{NIR-R}{\left(NIR+R\right)}\times \left(1+L\right)$	Huete [35]
Distance-based	$\ast MSAVI1=\frac{NIR-R}{NIR+R+L}\times \left(1+L\right)$ Qi et al. [39]	Qi et al. [50]	R = reflectance in the visible redL = 1 − (2 × NDVI × WDVI)WDVI = NIR − (slope × R)
	$\ast MSAVI2=\frac{2pNIR+1-\sqrt{{\left(2pNIR+1\right)}^2-8\left(pNIR-pR\right)}}{2}$	Qi et al. [50]	pR = reflectance of the redpNIR = reflectance of the near infrared

* NDVI: normalized difference vegetation index; CTVI: corrected transformed vegetation index; NRVI: normalized ratio vegetation index; RATIO: ratio vegetation; RVI: ratio vegetation index; SAVI: soil-adjusted vegetation index; MSAVI-1 and MSAVI2: modified soil-adjusted vegetation indices; WDVI: weighted difference vegetation index.

Table 5

Analysis of variance (ANOVA) for slope- and distance-based categories of vegetation indices.

	* VIs	S ource	*SS	* df	* MS	F-Statistic	* Prob>F
Slope-based	CTVI	* Groups	0.02	1	0.02	120.57	0.00
		* Error	0.02	94	0.00
		* Total	0.04	95
	NDVI	Groups	0.06	1	0.06	132.78	0.00
		Error	0.04	94	0.00
		Total	0.10	95
	NRVI	Groups	0.06	1	0.06	132.78	0.00
		Error	0.04	94	0.00
		Total	0.10	95
	RATIO	Groups	0.61	1	0.61	169.74	0.00
		Error	0.34	94	0.00
		Total	0.95	95
	RVI	Groups	0.11	1	0.11	105.84	0.00
		Error	0.10	94	0.00
		Total	0.21	95
	SAVI	Groups	0.01	1	0.01	16.60	0.00
		Error	0.08	94	0.00
		Total	0.09	95
Distance-based	MSAVI1	Groups	0.23	1	0.23	132.42	0.00
		Error	0.17	94	0.00
		Total	0.40	95
	MSAVI2	Groups	0.11	1	0.11	105.7	0.00
		Error	0.10	94	0.00
		Total	0.21	95

* VIs: vegetation indices; Groups: the variability between groups; Error: the variability within groups; Total: the total variability; SS: sum of squares; df: the degrees of freedom; MS: the mean squares for each source; F-statistic: the ratio of the mean squares; Prob>F: the returned p-value.

Table 6

Shows the summary of the most correlated layers to changes in vegetation structure between 2016–2019 based on the feature selection algorithm.

Dependent Variables	Variable Category	Variable	Loading
-Total Harvested Volume-Average Volume-Changes in Openness Based on Volume-Changes in Openness Based on Basal Area-Changes in Openness Based on Tree Density	* TA	GLDVContrast-NIR band	0.927
	TA	Contrast-NIR band	0.927
	TA	GLDVContrast-NIR band	0.927
	TA	Dissimilarity-blue band	0.924
	TA	GLDVMean-blue band	0.924
	TA	Entropy-NIR band	0.92
	TA	STANDDEV-blue band	0.918
	* VI	NDVI	0.975
	VI	NRVI	−0.975
	VI	MSAVI1	0.975
	VI	CTVI	0.975
	VI	MSAVI2	0.972
	VI	RVI	−0.972
	VI	SAVI	0.971
	VI	RATIO	0.965

* NDVI: normalized difference vegetation index; CTVI: corrected transformed vegetation index; NRVI: normalized ratio.vegetation index; RATIO: ratio vegetation; RVI: ratio vegetation index; SAVI: soil-adjusted vegetation index; MSAVI-1 and MSAVI2: modified soil-adjusted vegetation indices; GLDV: gray-level difference vector; TA: texture analysis; VI: vegetation index.

Table 7

Modeling several forest parameters using random forest (RF) algorithm.

	Parameters (Per Plot)	Independent Variables	RMSE	RMSE Reduction	Bias	Bias%
2016	Total volume (m3)	* VIs	4.49	16.98	0.66	2.43
		* TA	5.4	19.37	−1.63	−6.24
		Total	4.51	16.53	−0.27	−1.01
	Average volume (m3)	VIs	4.52	17.11	0.73	2.71
		TA	0.24	17.06	−0.002	−0.13
		Total	0.23	16.04	−0.01	−0.68
2019	Total volume (m3)	VIs	3.53	27.24	−1.62	−14.3
		TA	3.43	28.87	−0.05	−0.47
		Total	2.58	25.59	1.88	15.74
	Average volume (m3)	VIs	0.22	17.47	−0.036	2.89
		TA	0.2	15.74	−0.02	−2.25
		Total	0.21	16.97	0.03	2.67
Harvested or changes	Total volume (m3)	VIs	6.88	53.42	3.38	20.78
		TA	8.16	52.63	−0.63	−4.23
		Total	5.18	33.57	−0.53	−3.57
	Average volume (m3)	VIs	0.45	26.74	−0.15	−9.7
		TA	0.4	24.84	−0.023	−1.46
		Total	0.35	23.52	0.15	9.5
	* COV%	VIs	20.65	43.13	8.85	15.6
		TA	18.49	34.66	0.24	0.45
		Total	11.48	20.97	0.48	0.88
	* COB%	VIs	20.03	42.29	8.73	15.56
		TA	18.12	34.29	−0.04	−0.08
		Total	11.15	20.54	−0.09	−0.18
	* COD%	VIs	21.55	51.87	11.29	21.37
		TA	16.47	32.9	−1.68	−3.47
		Total	9.07	17.79	−1.91	−3.9

* COV: changes in crown openness based on volume; COB: changes in crown openness based on basal area (BA); COD: changes in tree density; VIs: vegetation indices; and TA: texture analysis.