Next Article in Journal
Detection of Pine Shoot Beetle (PSB) Stress on Pine Forests at Individual Tree Level using UAV-Based Hyperspectral Imagery and Lidar
Next Article in Special Issue
Sequential PCA-based Classification of Mediterranean Forest Plants using Airborne Hyperspectral Remote Sensing
Previous Article in Journal
Hyperspectral Anomaly Detection Based on Separability-Aware Sample Cascade
Previous Article in Special Issue
Changes in Forest Net Primary Productivity in the Yangtze River Basin and Its Relationship with Climate Change and Human Activities
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Azadeh Abdollahnejad
1,*,
Dimitrios Panagiotidis
1 and
Lukáš Bílek
2
1
Department of Forest Management, Faculty of Forestry and Wood Sciences, Czech University of Life Sciences (CULS), Kamýcká 129, 165 21 Prague, Czech Republic
2
Department of Silviculture, Faculty of Forestry and Wood Sciences, Czech University of Life Sciences (CULS), Kamýcká 129, 165 21 Prague, Czech Republic
*
Author to whom correspondence should be addressed.
Remote Sens. 2019, 11(21), 2539; https://0-doi-org.brum.beds.ac.uk/10.3390/rs11212539
Submission received: 20 September 2019 / Revised: 20 October 2019 / Accepted: 27 October 2019 / Published: 29 October 2019
(This article belongs to the Special Issue Monitoring Forest Change with Remote Sensing)

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.

Graphical Abstract

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 km2. 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 m3 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; DBH1.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 DBH1.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 DBH1.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.).
V S c o t s   p i n e = ( A D B H 1.3 ( ( B C l o g 10 ( D B H 1.3 + 1 ) ) H E ) ) ( A 1 ( D B H 1.3 + 1 ) B 1 H C 1 )
V N o r w a y   s p r u c e = ( A ( D B H 1.3 + 1 ) B H C ) ( E ( D B H 1.3 + 1 ) F H G )

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
R M S E = i = 1 n ( X o b s , i X m o d e l , i ) 2 n
where X o b s is observed values and X m o d e 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
R M S E ( r e d u c t i o n ) = RMSE X ¯ obs × 100
where X ¯ 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
B i a s = 1 n x i x ¯ j
B i a s % = B i a s x j ¯ × 100
where n is the number of samples, x i is the observed value and 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 (m3) 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://0-www-mdpi-com.brum.beds.ac.uk/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.

References

  1. Gillis, M.D.; Leckie, D.G. Forest Inventory Mapping Procedures Across Canada; Information Report PI-X-114; Forestry Canada, Petawawa National Forestry Institute: Chalk River, ON, Canada, 1993.
  2. Wulder, M.A.; White, J.C.; Fournier, R.A.; Luther, J.E.; Magnussen, S. Spatially explicit large area biomass estimation: Three approaches using forest inventory and remotely sensed imagery in a GIS. Sensors 2008, 8, 529–560. [Google Scholar] [CrossRef] [PubMed]
  3. Hanewinkel, M.; Cullmann, D.A.; Schelhaas, M.J.; Nabuurs, G.J.; Zimmermann, N.E. Climate change may cause severe loss in the economic value of European forest land. Nat. Clim. Chang. 2012, 3, 203–207. [Google Scholar] [CrossRef]
  4. Nabuurs, G.J.; Lindner, M.; Verkerk, P.J.; Gunia, K.; Deda, P.; Michalak, R.; Grassi, G. First signs of carbon sink saturation in European forest biomass. Nat. Clim. Chang. 2013, 3, 792–796. [Google Scholar] [CrossRef]
  5. He, H.S.; Mladenoff, D.J.; Crow, T.R. Linking an ecosystem model and alandscape model to study forest species response to climate warming. Ecol. Model. 1999, 114, 213–233. [Google Scholar] [CrossRef]
  6. Dixon, R.K.; Brown, S.; Houghton, R.A.; Solomon, A.M.; Trexler, M.C.; Wisniewski, J. Carbon pools and flux of global forest ecosystem. Science 1994, 263, 185–190. [Google Scholar] [CrossRef] [PubMed]
  7. Kašpar, J.; Hlavatý, R.; Kuželka, K.; Marušák, R. The Impact of Assumed Uncertainty on Long-Term Decisions in Forest Spatial Harvest Scheduling as a Part of Sustainable Development. Forests 2017, 8, 335. [Google Scholar] [CrossRef]
  8. Asner, G.P.; Rudel, T.K.; Aide, T.M.; Defries, R.; Emerson, R. A contemporary assessment of change in humid tropical forests. Conserv. Biol. 2009, 23, 1386–1395. [Google Scholar] [CrossRef]
  9. Oliver, C.D.; Larson, B.C. Forest Stand Dynamics; John Wiley and Sons: Hoboken, NJ, USA, 1996. [Google Scholar]
  10. Cannon, C.H.; Peart, D.R.; Leighton, M.; Keratinase, K. The structure of lowland rain-forest after selective logging in West Kalimantan, Indonesia. For. Ecol. Manag. 1994, 67, 49–68. [Google Scholar] [CrossRef]
  11. Perrotta, J.A.; Francis, J.K.; Knowles, O.H. Harvesting intensity affects forest structure and composition in an upland Amazonian forest. For. Ecol. Manag. 2002, 169, 243–255. [Google Scholar] [CrossRef]
  12. Jackson, S.M.; Frederiksen, T.S.; Malcolm, J.R. Area disturbed and residual stand damage following logging in a Bolivian tropical forest. For. Ecol. Manag. 2002, 166, 271–283. [Google Scholar] [CrossRef]
  13. Asner, G.P.; Keller, M.; Pereira, R.; Weeded, J.C. Remote sensing of selective logging in Amazonia. Remote Sens. Environ. 2002, 80, 483–496. [Google Scholar] [CrossRef]
  14. Laporte, N.T.; Staubach, J.A.; Grosch, R.; Lin, T.S.; Goetz, S.J. Expansion of industrial logging in Central Africa. Science 2007, 316, 1451. [Google Scholar] [CrossRef] [PubMed]
  15. Keller, M.; Palace, M.; Asner, G.P.; Pereira, R.; Silva, J.N.M. Coarse woody debris in undisturbed and logged forests in the eastern Brazilian Amazon. Glob. Chang. Biol. 2004, 10, 784–795. [Google Scholar] [CrossRef]
  16. Shimabukuro, Y.E.; Beachie, R.; Gracchi, R.C.; Achard, F. Assessment of forest degradation in Brazilian Amazon due to selective logging and fires using time series of fraction images derived from Landsat ETM+ images. Remote Sens. Lett. 2014, 5, 773–782. [Google Scholar] [CrossRef]
  17. Lunetta, R.S.; Ediriwickrema, J.; Johnson, D.M.; Lyon, J.G.; McKerrow, A. Impacts of vegetation dynamics on the identification of land-cover change in a biologically complex community in North Carolina, USA. Remote Sens. Environ. 2002, 82, 258–270. [Google Scholar] [CrossRef]
  18. Chiteculo, V.; Abdollahnejad, A.; Panagiotidis, D.; Surový, P.; Sharma, R.P. Defining Deforestation Patterns Using Satellite Images from 2000 and 2017: Assessment of Forest Management in Miombo Forests—A Case Study of Huambo Province in Angola. Sustainability 2019, 11, 98. [Google Scholar] [CrossRef]
  19. Kayitakire, F.; Hamel, C.; Defourny, P. Retrieving forest structure variables based on image texture analysis and IKONOS-2 imagery. Remote Sens. Environ. 2006, 102, 390–401. [Google Scholar] [CrossRef]
  20. St-Onge, B.; Hu, Y.; Vega, C. Mapping the height and above-ground biomass of a mixed forest using LiDAR and stereo Ikonos images. Int. J. Remote Sens. 2008, 29, 1277–1294. [Google Scholar] [CrossRef]
  21. Ozdemir, I.; Karnieli, A. Predicting forest structural parameters using the image texture derived from WorldView-2 multi-spectral imagery in a dryland forest, Israel. Int. J. Appl. Earth Obs. Geoinform. 2011, 13, 701–710. [Google Scholar] [CrossRef]
  22. Stibig, H.; Achard, F.; Carboni, S.; Raši, R.; Miettinen, J. Change in tropical forest cover of Southeast Asia from 1990 to 2010. Biogeosciences 2014, 11, 247–258. [Google Scholar] [CrossRef] [Green Version]
  23. Abdollahnejad, A.; Panagiotidis, D.; Shataee Joybari, S.; Surový, P. Prediction of Dominant Forest Tree Species Using QuickBird and Environmental Data. Forests 2017, 8, 42. [Google Scholar] [CrossRef]
  24. Abdollahnejad, A.; Panagiotidis, D.; Surový, P. Estimation and Extrapolation of Tree Parameters Using Spectral Correlation between UAV and Pléiades Data. Forests 2018, 9, 85. [Google Scholar] [CrossRef]
  25. Souza, C.M., Jr.; Siqueira, J.V.; Sales, M.H.; Fonseca, A.V.; Ribeiro, J.G.; Numata, I.; Cochrane, M.A.; Barber, C.P.; Roberts, D.A.; Barlow, J. Ten-year Landsat classification of deforestation and forest degradation in the Brazilian Amazon. Remote Sens. 2013, 5, 5493–5513. [Google Scholar] [CrossRef]
  26. Melaas, E.; Friedl, M.A.; Zhu, Z. Detecting interannual variation in deciduous broadleaf forest phenology using Landsat TM/ETM+ data. Remote Sens. Environ. 2013, 132, 176–185. [Google Scholar] [CrossRef]
  27. Souza, C., Jr.; Firestone, L.; Silva, L.M.; Roberts, D. Mapping forest degradation in the Eastern Amazon from SPOT4 through spectral mixture models. Remote Sens Environ. 2003, 87, 494–506. [Google Scholar] [CrossRef]
  28. Mitchell, A.L.; Rosenqvist, A.; Mora, B. Current remote sensing approaches to monitoring forest degradation in support of countries measurement, reporting and verification (MRV) systems for REDD+. Carbon Balance Manag. 2017, 12, 9. [Google Scholar] [CrossRef]
  29. Kennedy, C.; Cuddihy, J.; Engel-Yan, J. The changing metabolism of cities. J. Ind. Ecol. 2007, 11, 43–59. [Google Scholar] [CrossRef]
  30. Thompson, I.D.; Guariguata, M.R.; Okabe, k.; Bahamondez, C.; Nasi, R.; Heymell, V. An operational framework for defining and monitoring forest degradation. Ecol. Soc. 2013, 18, 20. [Google Scholar] [CrossRef]
  31. Turubanova, S.; Potapov, P.; Krylov, A.; Tyukavina, A.; McCarty, J.L.; Radeloff, V.C.; Hansen, M.C. Using the landsat data archive to assess long-term regional forest dynamics assessment in Eastern Europe, 1985–2012. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2015, 531–537. [Google Scholar] [CrossRef]
  32. Zarco-Tejada, P.J.; González-Dugo, V.; Williams, L.E.; Suárez, L.; Berni, J.A.J.; Goldhamer, D.; Fereres, E. A PRI-based water stress index combining structural and chlorophyll effects: Assessment using diurnal narrow-band airborne imagery and the CWSI thermal index. Remote Sens. Environ. 2013, 138, 38–50. [Google Scholar] [CrossRef]
  33. Glenn, E.; Nagler, P.; Huete, A. Vegetation index methods for estimating evapotranspiration by remote sensing. Surv. Geophys. 2010, 31, 531–555. [Google Scholar] [CrossRef]
  34. Tucker, C.J.; Townshend, J.R.G.; Goff, T.E. African land-cover classification using satellite data. Science 1985, 227, 369–375. [Google Scholar] [CrossRef] [PubMed]
  35. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213. [Google Scholar] [CrossRef]
  36. Huete, A.R. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 1988, 25, 295–309. [Google Scholar] [CrossRef]
  37. Eastman, J.R. IDRISI Kilimanjaro: Guide to GIS and Image Processing; Clark University: Worcester, MA, USA, 2003. [Google Scholar]
  38. Blackard, J.A.; Finco, M.V.; Helmer, E.H.; Holden, G.R.; Hoppus, M.L.; Jacobs, D.M.; Lister, A.J.; Moisen, G.G.; Nelson, M.D.; Riemann, R.; et al. Mapping U.S. forest biomass using nationwide forest inventory data and moderate resolution information. Remote Sens. Environ. 2008, 112, 1658–1677. [Google Scholar] [CrossRef]
  39. Hansen, J.; Makiko, S.; Pushker, K.; David, B.; Robert, B.; Valerie, M.D.; Mark, P.; Maureen, R.; Dana, L.R.; James, C.Z. Target atmosphere CO2: Where should humanity aim. Open Atmos. Sci. J. 2008, 2, 217–231. [Google Scholar] [CrossRef]
  40. Romero-Sanchez, M.E.; Ponce-Hernandez, R. Assessing and Monitoring Forest Degradation in a Deciduous Tropical Forest in Mexico via Remote Sensing Indicators. Forests 2017, 8, 302. [Google Scholar] [CrossRef]
  41. Soil Map. Available online: http://www.geology.cz/extranet-eng/services/web-applications/ (accessed on 19 August 2019).
  42. Tolasz, R. Atlas podnebí Česka (Climate atlas of Czechia), 1st ed.; Český hydrometeorologický ústav, Univerzita Palackého v Olomouci: Praha, Czech Republic, 2007; p. 255. [Google Scholar]
  43. Petráš, R.; Pajtík, J. Sústava česko-slovenských objemových tabuliek drevín. Lesnický Časopis 1991, 1, 49–56. [Google Scholar]
  44. Purevdorj, T.S.; Tateishi, R.; Ishiyama, T.; Honda, Y. Relationships between percent vegetation cover and vegetation indices. Int. J. Remote Sens. 1998, 19, 3519–3535. [Google Scholar] [CrossRef]
  45. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation system in the great plains with ERTS. In Proceedings of the Third Earth Resources Technology Satellite-1 Symposium, Greenbelt, MD, USA, 10–14 December 1974. pp. 309–3017. Available online: https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19740022614.pdf (accessed on 19 August 2019).
  46. Perry, C.R., Jr.; Lautenschlager, L.F. Functional equivalence of spectral vegetation indices. Remote Sens. Environ. 1984, 14, 169–182. [Google Scholar] [CrossRef]
  47. Baret, F.; Guyot, G. Potentials and limits of vegetation indices for LAI and APAR assessment. Remote Sens. Environ. 1991, 35, 161–173. [Google Scholar] [CrossRef]
  48. Birth, G.S.; McVey, G.R. Measuring the color of growing turf with a reflectance spectrophotometer. Agron. J. 1968, 60, 640–643. [Google Scholar] [CrossRef]
  49. Richardson, A.J.; Wiegand, C.L. Distinguishing vegetation from soil background information. Photogram. Eng. Remote Sens. 1977, 43, 15–41. [Google Scholar]
  50. Qi, J.; Chehbouni, A.; Huete, A.R.; Kerr, Y.H. Modified Soil Adjusted Vegetation Index (MSAVI). Remote Sens. Environ. 1994, 48, 119–126. [Google Scholar] [CrossRef]
  51. Sader, S.A.; Bertran, M.; Wilson, E.H. Satellite change detection of forest harvest patterns on an industrial forest landscape. For. Sci. 2003, 49, 341–353. [Google Scholar]
  52. Masiliūnas, D. Evaluating the Potential of Sentinel-2 and Landsat Image Time Series for Detecting Selective Logging in the Amazon; Wageningen University and Research Centre: Wageningen, The Netherlands, 2014; p. 54. [Google Scholar]
  53. Shchur, A.; Bragina, E.; Sieber, A.; Pidgeon, A.M.; Radelof, V.C. Monitoring selective logging with Landsat satellite imagery reveals that protected forests in Western Siberia experience greater harvest than non-protected forests. Environ. Conserv. 2017, 44, 191–199. [Google Scholar] [CrossRef] [Green Version]
  54. Read, J.M. Spatial analyses of logging impacts in Amazonia using remotely sensed data. Photogramm. Eng. Remote Sens. 2003, 69, 275–282. [Google Scholar] [CrossRef]
Figure 1. Geographic locations of the study area and plots.
Figure 1. Geographic locations of the study area and plots.
Remotesensing 11 02539 g001
Figure 2. Shows the average amount of volume for harvested and living trees.
Figure 2. Shows the average amount of volume for harvested and living trees.
Remotesensing 11 02539 g002
Figure 3. Visual comparison of the average and total volume before (2016) and after harvest (2019) using box-plots.
Figure 3. Visual comparison of the average and total volume before (2016) and after harvest (2019) using box-plots.
Remotesensing 11 02539 g003
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 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.
Remotesensing 11 02539 g004
Figure 5. Application of several slope-based vegetation indices (VIs) obtained from remote sensing-based canopies.
Figure 5. Application of several slope-based vegetation indices (VIs) obtained from remote sensing-based canopies.
Remotesensing 11 02539 g005
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 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).
Remotesensing 11 02539 g006
Figure 7. Application of distance-based vegetation indices obtained from remote sensing-based canopies.
Figure 7. Application of distance-based vegetation indices obtained from remote sensing-based canopies.
Remotesensing 11 02539 g007
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 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).
Remotesensing 11 02539 g008
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 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.
Remotesensing 11 02539 g009
Figure 10. Texture analysis applied in NIR band for the years 2016 and 2019.
Figure 10. Texture analysis applied in NIR band for the years 2016 and 2019.
Remotesensing 11 02539 g010
Figure 11. Continuous maps of forest parameters changes related to harvesting operation.
Figure 11. Continuous maps of forest parameters changes related to harvesting operation.
Remotesensing 11 02539 g011
Table 1. Shows the used equations and their accuracies for estimation of height and DBH1.3 for each species separately. The below equations were estimated using the relation between (i) stump diameter and DBH1.3 (ii) DBH1.3 and height.
Table 1. Shows the used equations and their accuracies for estimation of height and DBH1.3 for each species separately. The below equations were estimated using the relation between (i) stump diameter and DBH1.3 (ii) DBH1.3 and height.
Estimated ParameterTree SpeciesNumber of TreesTrendlineEquationR2
HeightScots pine25lineary = 0.1581x + 21.1410.48
Norway spruce25lineary = 0.2354x + 17.1450.58
D B H 1.3 Scots pine40lineary = 1.0601x + 5.75830.92
Norway spruce40lineary = 1.29x - 0.1650.96
Table 2. Shows the parameters used for the stem volume equations for both tree species.
Table 2. Shows the parameters used for the stem volume equations for both tree species.
Scots PineABCE A 1 B 1 C 1
Coefficients3.03E-052.0752380.0124920.9610280.071975−2.124491.372591
Norway spruceABCEFG
Coefficients4.01E-051.8218161.1320629.29E-03−1.020370.896101
Table 3. Characteristics of Pléiades-HR 1A and Pléiades-HR 1B satellite sensors.
Table 3. Characteristics of Pléiades-HR 1A and Pléiades-HR 1B satellite sensors.
Pléiades HR 1A-1B
Imagery products2 m multi-spectral resolution
Bundle: 50 cm panchromatic—2 m multi-spectral
Panchromatic resolution 50 cm
Spectral bandsRed: 600–720 nm
Green: 490–610 nm
Blue: 430–550 nm
Near Infrared: 750–950 nm
Panchromatic: 480–830 nm
Image location accuracyWith * GCPs is 1 m
Without GCPs is 3 m (CE90)
* GCPs: ground control points.
Table 4. Equations of applied vegetation indices (VIs).
Table 4. Equations of applied vegetation indices (VIs).
VIsReferencesExplanation of Symbols
Slope-based N D V I = N I R R N I R + R Rouse et al. [45]R = red band
B = blue band
L = soil
NIR = near infrared band
C T V I = N D V I + 0.5 A B S ( N D V I + 0.5 ) × A B S ( N D V I + 0.5 ) Perry and Lautenschlager [46]
N R V I = R V I 1 R V I + 1 Baret and Guyot [47]
R A T I O = N I R R Birth and McVey [48]
R V I = R N I R Richardson and Wiegand [49]
S A V I = N I R R ( N I R + R ) × ( 1 + L ) Huete [35]
Distance-based M S A V I 1 = N I R R N I R + R + L × ( 1 + L ) Qi et al. [39]Qi et al. [50]R = reflectance in the visible red
L = 1 − (2 × NDVI × WDVI)
WDVI = NIR − (slope × R)
M S A V I 2 = 2 p N I R + 1 ( 2 p N I R + 1 ) 2 8 ( p N I R p R ) 2 Qi et al. [50]pR = reflectance of the red
pNIR = 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.
Table 5. Analysis of variance (ANOVA) for slope- and distance-based categories of vegetation indices.
* VIsS ource*SS* df* MSF-Statistic* Prob>F
Slope-basedCTVI* Groups0.0210.02120.570.00
* Error0.02940.00
* Total0.0495
NDVIGroups0.0610.06132.780.00
Error0.04940.00
Total0.1095
NRVIGroups0.0610.06132.780.00
Error0.04940.00
Total0.1095
RATIOGroups0.6110.61169.740.00
Error0.34940.00
Total0.9595
RVIGroups0.1110.11105.840.00
Error0.10940.00
Total0.2195
SAVIGroups0.0110.0116.600.00
Error0.08940.00
Total0.0995
Distance-basedMSAVI1Groups0.2310.23132.420.00
Error0.17940.00
Total0.4095
MSAVI2Groups0.1110.11105.70.00
Error0.10940.00
Total0.2195
* 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.
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 VariablesVariable CategoryVariableLoading
-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
* TAGLDVContrast-NIR band0.927
TAContrast-NIR band0.927
TAGLDVContrast-NIR band0.927
TADissimilarity-blue band0.924
TAGLDVMean-blue band0.924
TAEntropy-NIR band0.92
TASTANDDEV-blue band0.918
* VINDVI0.975
VINRVI−0.975
VIMSAVI10.975
VICTVI0.975
VIMSAVI20.972
VIRVI−0.972
VISAVI0.971
VIRATIO0.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.
Table 7. Modeling several forest parameters using random forest (RF) algorithm.
Parameters (Per Plot)Independent VariablesRMSERMSE ReductionBiasBias%
2016Total volume (m3)* VIs4.4916.980.662.43
* TA5.419.37−1.63−6.24
Total 4.5116.53−0.27−1.01
Average volume (m3)VIs4.5217.110.732.71
TA0.2417.06−0.002−0.13
Total 0.2316.04−0.01−0.68
2019Total volume (m3)VIs3.5327.24−1.62−14.3
TA3.4328.87−0.05−0.47
Total 2.5825.591.8815.74
Average volume (m3)VIs0.2217.47−0.0362.89
TA0.215.74−0.02−2.25
Total 0.2116.970.032.67
Harvested or changesTotal volume (m3)VIs6.8853.423.3820.78
TA8.1652.63−0.63−4.23
Total 5.1833.57−0.53−3.57
Average volume (m3)VIs0.4526.74−0.15−9.7
TA0.424.84−0.023−1.46
Total 0.3523.520.159.5
* COV%VIs20.6543.138.8515.6
TA18.4934.660.240.45
Total 11.4820.970.480.88
* COB%VIs20.0342.298.7315.56
TA18.1234.29−0.04−0.08
Total 11.1520.54−0.09−0.18
* COD%VIs21.5551.8711.2921.37
TA16.4732.9−1.68−3.47
Total 9.0717.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.

Share and Cite

MDPI and ACS Style

Abdollahnejad, A.; Panagiotidis, D.; Bílek, L. An Integrated GIS and Remote Sensing Approach for Monitoring Harvested Areas from Very High-Resolution, Low-Cost Satellite Images. Remote Sens. 2019, 11, 2539. https://0-doi-org.brum.beds.ac.uk/10.3390/rs11212539

AMA Style

Abdollahnejad A, Panagiotidis D, Bílek L. An Integrated GIS and Remote Sensing Approach for Monitoring Harvested Areas from Very High-Resolution, Low-Cost Satellite Images. Remote Sensing. 2019; 11(21):2539. https://0-doi-org.brum.beds.ac.uk/10.3390/rs11212539

Chicago/Turabian Style

Abdollahnejad, Azadeh, Dimitrios Panagiotidis, and Lukáš Bílek. 2019. "An Integrated GIS and Remote Sensing Approach for Monitoring Harvested Areas from Very High-Resolution, Low-Cost Satellite Images" Remote Sensing 11, no. 21: 2539. https://0-doi-org.brum.beds.ac.uk/10.3390/rs11212539

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop