Multi-Temporal Sentinel-2 Data Analysis for Smallholding Forest Cut Control

Abstract

Land fragmentation and small plots are the main features of the rural environment of Galicia (NW Spain). Smallholding limits land use management, representing a drawback in local forest planning. This study analyzes the potential use of multitemporal Sentinel-2 images to detect and control forest cuts in very small pine and eucalyptus plots located in southern Galicia. The proposed approach is based on the analysis of Sentinel-2 NDVI time series in 4231 plots smaller than 3 ha (average 0.46 ha). The methodology allowed us to detect cuts, allocate cut dates and quantify plot areas due to different cutting cycles in an uneven-aged stand. An accuracy of approximately 95% was achieved when the whole plot was cut, with an 81% accuracy for partial cuts. The main difficulty in detecting and dating cuts was related to cloud cover, which affected the multitemporal analysis. In conclusion, the proposed methodology provides an accurate estimation of cutting date and area, helping to improve the monitoring system in sustainable forest certifications to ensure compliance with forest management plans.

1. Introduction

Forest ecosystems provide multiple ecological, economic and social benefits. They harbor a portion of the world’s biodiversity, play a key role in regulating water flows, protect soils and contribute directly to national incomes and the local livelihoods of millions of people worldwide [1]. Furthermore, these ecosystems regulate key aspects of the global carbon cycle [2] and weather patterns via a number of different mechanisms, such as forest albedos [3], sensible heat and aerodynamic roughness [4,5]. Specifically, forest areas are among the land use types that capture and store more carbon than they use. For example, returning farmland to larch forests has been shown to increase the carbon concentration in mineral soils at a rate of 100× g m⁻² a⁻¹ [6]. Moreover, approximately 31% of the land surface of the planet is covered by forests [7], representing approximately 38% in the EU and 36.5% in Spain [8], the country where this study was carried out. Additionally, global forest cover has increased by 7% since the 1980s [9]; this increase is directly caused by the abandonment of agricultural land where it is difficult to mechanize processes to develop modern agriculture [10]. The abandonment of this agricultural land has mainly occurred in districts already dominated by forests [11]. Spain is not an exception, however: reforestation has decreased in the last 10 years, and in 2007, 54,000 ha of forested areas were repopulated, while in 2017, this increase was approximately 12,000 ha [12]. It is important to highlight the economic support that forestry generates in today’s society, either as a livelihood or as an economic complement. It is a broad sector that encompasses different activities and industries. European forest-based industries include woodworking, manufacturing pulp, paper and paper products, furniture, printing and bioenergy [13]. Together, these industries comprise approximately 420,000 companies with a total turnover of more than 520,000 million EUR [14]. For all these reasons stated, sustainable management and proper policies are essential for maintaining the ecological and socioeconomic functions of forest ecosystems [15]. To reach this goal, international forest certification processes and systems (FCSs) have been developed. The total area of certified forests worldwide steadily increased from less than 25 million hectares in 1998 to 405 million hectares in 2018, comprising 11% of the world’s forests [16]. In this sense, a requirement for the sustainable development of forest management is to ensure the ecological security of forests [17].

Therefore, multitemporal observations can be useful tools to detect changes in forest ecosystems. In the face of rapid and multiple forest changes, remote sensing has become the most practical and efficient means of extracting information with great temporal, spatial and thematic detail [18]. Consequently, multiple detection methods have been developed to monitor changes in forest areas. Most of these methods involve bitemporal images [19], and the newest method continuously records the dynamics of disturbances in a dense remotely sensed time series [20], making it easier to track rapid forest changes, such as fires. Regardless of the specific process, these remote sensing methods outperform traditional methods [21]. Hyperspectral sensors, which monitor the Earth’s surface in contiguous and narrow bands, allow us to capture the biochemical composition of vegetation [22], providing a significant level of detail. However, an optimal set must be selected from the broad set of wave bands, as most of the bands are highly correlated and require large amounts of computational power [23]. Generally, when very-high-spatial-resolution images are used, the spectral responses of individual trees are affected by differences in canopy illumination and background signals [24]. For this reason, in vegetation studies that cover large geographic areas, relatively dense and freely accessible multispectral imagery, such as Sentinel-2, appears to be the best solution [25,26]. Despite the wide availability of data, few studies have been conducted focusing on characteristic smallholding plots in which temporal analyses conducted using Sentinel-2 data to analyze forest vegetation and deal with the detection of cutting in a wide territory [27,28] were the objectives. Studies using satellite imagery have focused on the deforestation or degradation of ecosystems [18,29], especially in large areas such as forests where deforestation affects highly valuable ecosystems, and these images have also been used to analyze vegetation segmentation [30,31]. However, it is necessary to develop a methodology that, in the face of the rapid and multiple changes that occur in forests due to harvesting and forestry activities, analyzes information at the plot level. In addition, this method must contribute to and facilitate the planning and execution of the extraction of wood to control and guarantee good practices and guarantee the objectives of good forest use. In other words, the development of a methodology that is adapted to a given study area will enable forest managers to recognize when cutting has been carried out in that area. In this sense, the mapping of forest changes will provide information on the potential for obtaining wood, patterning forest growth and detecting risks associated with the multiple impacts that these ecosystems may suffer.

Focusing on the regional object of this study, Galicia is the most important forest region in Spain, with 8% of the total Spanish forest area [32]. This region contains 2,040,754 million ha, wherein almost 1,500,000 million ha belongs to wooded areas, with 28% conifers, 51% hardwoods and 21% mixed-use forests [33]. Eucalyptus globulus and Pinus pinaster are the most widespread forest species in this area due to commercial interest. Forest management and industry account for 12% of the final agricultural production in the rural economy of Galicia, compared to 3.5% in Spain overall. Forestry and wood transformation contribute 13% and 43%, respectively, to the gross domestic product (GDP) of Galicia, and together, these processes generate almost 3% of the total employment of this region [34]. On the other hand, this forested area is characterized by high fragmentation in terms of ownership, and more than two-thirds of the forestlands of Galicia belong to 670,000 owners. Individual forest holdings vary in size between 1.5 and 2.0 ha on average. Private forests comprise 97.3% of Galician forestland, and this land is often divided into up to 10 noncontiguous parcels, yielding a mean surface area of only 0.26 ha [35]. Two-thirds of the land under this private ownership belongs to individual properties, and the remaining third is collectively owned by residents through Communities of Communal Forests. These dispersed and very small properties are one of the main obstacles preventing the further development of the sector, although great efforts have been made in recent decades to overcome this structural problem through the consolidation of territories and the creation of forest societies [35]. In addition to this high dispersion problem, the low profitability caused by a lack of mechanization caused by the situation of plots on slopes, poor forestry or the volume of felling poses additional issues [36,37]. In addition, the low technical qualifications of owners, together with their advanced age, induce the low innovation and development of the sector [38]. Finally, it is necessary to highlight some of the impacts that these owners suffer the most, such as forest fires [39] and the presence of pests and insect diseases [40].

In this social, economic and environmental context, this work has the objective of verifying the potential of Sentinel-2 time-series images as key information in the process of the detection and temporary control of very-small-plot forest cuts within a management and control system. The priority of this study involves checking the level of effectiveness of these images and the influence that the size, shape, number and casuistry of the composition of the plots have on the detection of forest cuts. Thresholds are established to indicate, with high probability, the existence of felling, as well as the precision of this technology. At the same time, the influence of the irregularity of the applied time series resulting from the absence of information caused by the presence of clouds or other problems is analyzed. This innovative research focuses on the characteristics of the plots, and small parcels are designated to determine whether the proposed methodology is valid for smallholdings. The results achieved through this research could be integrated into the control process of a forest management group made up of thousands of owners and plots.

2. Materials and Methods

2.1. Study Area

This study was carried out in Galicia (Northwestern Spain) (Figure 1). This region comprises an area of 29,575 km² and has a population of almost three million people [41]. It has a Mediterranean oceanic climate with mild summers [42]. The predominant soils are developed on granitic rock and acid schist, have a loamy or sandy loam texture and are well-drained.

The selected plots come from the forest management database of the local company Asefor Ingeniería Forestal, S.L. From this database, the plots that are within the limit of Tile 29TNH and are smaller than 3 hectares with clear-cutting, cleaning or forest cutting in general occurring between 2018 and 2020 were selected as the study area. Of these plots, all the information regarding their characteristics and uses is known thanks to previous inventory work carried out by the company. Small parcels were designated to determine whether the proposed methodology is valid for smallholdings. In total, 4231 cadastral parcels (Figure 1), that is, individual and separate plots, were considered. This number of pilot plots was considered sufficient to verify the operability and adaptability of the developed control method. The plots have an average area of 0.47 ha, with a minimum area of 0.0034 ha and a maximum of 3 ha. The average perimeter of these plots is 324 m, with a minimum of 27 m and a maximum of 1545 m. An area of 1,157,000 ha, divided into pixels of 10 × 10 m, was analyzed and treated with the results presented in the subsequent sections and assigned to each of the 4231 cadastral parcels. All the plots contained forest species that were cut down at some point. Only species of the genera Pinus and Eucalyptus were cut. Some plots that were analyzed were not cut in their entire area; this may occur for multiple reasons, the most common of which are described below. Some plots contain several species (Figure 2b); therefore, two-aged stands (growing areas with trees of two distinct age classes [43]) or uneven-aged stands (stands of trees of three or more distinct age classes, either intimately mixed or in groups [43]) were considered. Other plots, despite having the same forest species, had different age classes (Figure 2a). Therefore, these plots had two-aged or uneven-aged silvicultural systems. Finally, some plots, or those that were also found in the two previous cases, contained “non-cuttable” forest species such as deciduous hardwoods, unique species or riparian vegetation. On the other hand, some of these plots were not entirely forested, sharing their areas with agricultural uses (meadows and crops), a condition that may affect the detection of cuts due to changes in vegetation that are typical of crop rotation.

2.2. Data Collection and Preprocessing

The dataset used in this study came from a selection of images obtained from the standard Sentinel-2 Level-2A. The Sentinel satellite carries an innovative wide-swath, high-resolution, multispectral instrument (MSI) with 13 spectral bands with wavelengths ranging from 440 to 2200 nm and from 10 to 60 m depending on the spectral band [44]; in this case, the resolution used was 10 m. The Sentinel-2 Level-2A product comprises 100 km × 100 km mosaics in the UTM/WGS84 projection [45] and provides orthorectified images with reflectance levels below the atmosphere (BOA) [46].

Sentinel-2 imagery collected between January 2018 and December 2020 was downloaded from the Copernicus Open Access Hub [47]. Initially, images with cloud percentages equal to or less than 5% in the selected time period were searched. On dates for which this condition was not met, due to the absence of images with less than 5% cloud cover, the condition was increased to a maximum of 15% cloudiness. The images were processed with the free software QGIS 3.8.2 using the SCP (semiautomatic classification) tool. An orthophoto was obtained from the Spanish National Geographic Institute’s PNOA (National Aerial Orthophotography Plan) [48]. Subsequently, each of the images was visually analyzed, and the images that exhibited cloudiness in the plots under analysis were discarded, as were those that were not within the coverage of the Tile 29TNH. The presence of elements in the images that hinder the correct analysis of data, such as aerosols, cloud shadows or the presence of fog in river valleys, must be analyzed and solved to guarantee the correct control of forest cutting. In total, 54 images were acquired for the selected period (

Table 1. Images acquired from the Copernicus Open Access Hub for the time period analyzed from 2018 to 2019.

2018	2019	2020
24/02/2018	14/02/2019	24/02/2020
No data	16/03/2019	18/03/2020 25/03/2020
18/04/2018	No data	14/04/2020 24/04/2020
05/05/2018 18/05/2018	13/05/2019	19/05/2020 24/05/2020 27/05/2020 29/05/2020
17/06/2018 19/06/2018 24/06/2018	No data	No data
09/07/2018	12/07/2019	No data
01/08/2018 21/08/2018 26/08/2018	21/08/2019	05/08/2020
07/9/2018	15/09/2019	29/09/2020
10/09/2018 10/10/2018 20/10/2018 22/10/2018	07/10/2019 10/10/2019 22/10/2019	16/10/2020
14/11/2018	No data	13/11/2020 20/11/2020 30/11/2020
26/12/2018	04/12/2019 26/12/2019	No data

). The data were visualized in QGIS and geospatially managed in the PostGIS database manager.

After the selection of images that were as cloud-free as possible, preprocessing was carried out to eliminate the presence of clouds or cloud shadows that were in the range of the studied plots. Pixels with blue band reflectances below 0.01 were considered to be predominantly related to undetected cloud shadows [49]. For this reason, the blue band (B2) of each image was used to eliminate pixels above the 1150 spectral profile. The objective of this process was that for each date analyzed, a cloud mask could be created that was omitted from the processing; in this way, atmospheric pollution that affected the calculation of the NDVI index was eliminated. Although this process is not perfect, as it does not totally eliminate all types of cloudiness or fog, it was valid for eliminating a large portion of the subsequent anomalies. In this study, only forested plots were worked on, and the elimination of other land use types (urban, water, etc.) did not affect the object of this methodology. The threshold of 1150 was obtained by performing systematic sampling in places where clouds were detected, and this value was chosen after several visual tests of the obtained results were performed.

2.3. Normalized Difference Vegetation Index (NDVI)

The spectral index used in this study was the normalized difference vegetation index (NDVI) [50]. If NDVI is compared with other vegetation indices, it has the best correlation with canopy cover, especially in the dry season [51]. NDVI is the most reliable vegetation index used to estimate forest density [52,53]. For these reasons, it was selected for use in this study. The description of NDVI is as follows:

$$\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}=\frac{B8-B4}{B8+B4}$$

(1)

where NIR refers to the spectral reflectance at near-infrared wavelengths (0.85–0.88 µm), and RED refers to the spectral reflectance at red wavelengths (0.64–0.67 µm) .

The aim of applying this index was to study how NDVI variations occur in different forest plots and to analyze the capacity of the proposed method to detect abrupt changes in forest cover, such as clear-cuts, cleaning and forest cuts in general. The mean NDVI value was obtained for each plot in each image through RStudio [54]; the same software was used for the generation and presentation of raster images. This process allows the spatial and temporal characteristics of these changes to be detected. Through this analysis with RStudio, each plot was temporally studied on each date (Table 1), and the phenology present in each plot was statistically studied, as were the changes caused by cuts. In addition, the processes were automated through the use of scripts with RStudio. The 4231 felling plots were reforested by the Pinus and Eucalyptus genera, although in many of the plots, other species or intercropped covers (mainly deciduous leafy plants, scrub, pastures, etc.) existed that were not cut, as was mentioned already; these plots experienced differing evolution over the course of the annuity and seasonality. Pinus and Eucalyptus have no marked seasonality, so the differences will be taken as periods. Therefore, abrupt changes in NDVI in one pixel are expected to be due to cutting in most cases. NDVI can help identify and quantify changes in the analyzed forest plots. To accomplish this, the differences in NDVI between different periods were calculated using the following equation:

$$NDVI\left(t\right)-NDVI\left(t-k\right)$$

(2)

where “NDVI” is the value of the raster at each point, “t” is the initial period of analysis and “k” is the number of previous periods with which the most recent image is compared. Generally, in this case study, k is equal to 1, as the most recent image is always compared with the immediately previous period. An exception occurred when, in cases where no image was available for the previous period, it was necessary to compare the most recent image with images from previous dates.

All the locations of the selected pixels and their values were stored in tables managed in PostGIS for later treatment in RStudio and visualization in QGIS.

2.4. Detection of the Forest Cutting Threshold

To detect the values that allow the probability of felling to be discerned, a statistical analysis was carried out. The packages used to carry out the classification tree were: library rpart, rpart plot and caret. In this process, NDVI differences were established for each satellite information period. All trend data were extracted in 53 time periods with 197,398 pixels analyzed. A statistical boxplot analysis was performed for all pixel trends for each plot over the period analyzed in this study (Figure 3). The lower extreme was then used to identify the cut plots. This extreme represents the global minimum value used for the detection. The value used was −0.1254; that is, a pixel was “cut” when the trend drop between periods was ≤−0.1254. The detected threshold values for the trend drops for which cuts could be detected were obtained for the set of all pixels of all the plots regardless of the coverage of the individual plots. This analysis also provided a median value (0.0009), lower extreme value (−0.1254) and highest (0.1279) and lowest extreme values (−0.1254).

The locations of the pixels that met the threshold conditions defined in each period were related to the surfaces. A database was obtained, including the period, identification of the plot and fall surface and its percentage in relation to the surface of each plot. This information was compared with RGB images. A total of 135 plots were selected and statistically analyzed). When felling is detected, statistical analyses are carried out for pixels below −0.12: mean of the values, twenty-fifth percentile, fiftieth percentile, seventy-fifth percentile, the area with variation and the relationship of this area with that of the total area of the plot. The same values were also studied for the same plot in the periods in which no short is detected. These analyses (Appendix A) are used to create the database that will form the classification tree. A decision tree or classification tree is a supervised matching-learning algorithm that will show us a series of sequential decisions [55]. This classification tree was constructed (Figure 4) with all the data collected as a prediction model by obtaining an algorithm whose objective (dependent variable) was the detection of forest cutting; specifically, 1 = “Cutting” or 2 = “No cutting”. The aim of the decision tree is to divide a complex decision into several simple decisions. Using this approach, the statistics defined from the data are predictor variables, while the short one would be the target variable. Since the class proportions of a mixed pixel are measured on a continuous scale from 0 to 1, it is necessary to apply decision tree regression to estimate them. The data set was split in two, with 70% being training data and the remaining 30% being test data, to check the effectiveness of the model. From this 70% training data, an algorithm was obtained, which allowed us to predict the value of the target variable using the independent variables. The rpart function is used to train the model, and the target variable, in this case, “clearcut”, is formulated, while the rest of the variables are used as predictors. The algorithm finds the independent variable that best separated the data into groups; this separation corresponds to a rule, and each rule has a node. The data were then separated (partitioned) into groups based on the obtained rules. Then, for each of the resulting groups, the same process was repeated. It first checks the percentile value, then the ratio between the occupied pixels and finally, the mean pixel values. This process was carried out repeatedly until it was no longer possible to obtain a better separation. When this occurred, the algorithm ended. The model is then tested with the remaining 30% of data, from which the confusion matrix is obtained (

Table 2. Statistical analysis including equation and description of each indicator.

Accuracy	$\frac{TP+TN}{TP+FN+TN+FP}$	TP = true positive, TN = true negative, FP = false positive, FN = false negative
95% CI	$\overline{x}\pm 1.96\frac{Sx}{\sqrt{n}}$	$\overline{x}$ = the sample mean $Sx$= margin of error N = standard error
p-value (ACC < NIR)	$z=\frac{\widehat{p}-p_0}{\sqrt{\frac{p_0\left(1-p_0\right)}{n}}}$	$\widehat{p}$= the sample proportion p0 = the hypothesized proportion n = the sample size ACC = accuracy NIR = No information rate p-value for ACC > NIR
Kappa	$\frac{p_0-p_e}{1-p_e}$	P0 = proportion of trials in which judges agree Pe = proportion of trials in which agreement would be expected due to chance
Mcnemar’s Test	$x^2=\frac{{\left(b-c\right)}^2}{b+c}$	-
Sensitivity	$\frac{TN}{TP+FN}$	TP = true positive, FN = false negative
Specificity	$\frac{TN}{TN+FP}$	TN = true negative, FP = false positive
Pos. Pred Value	$\frac{TP}{TP+FP}$×100	TP = true positive, FP = false positive
Neg. Pred Value	$\frac{TN}{TN+FN}$×100	TN = true negative, FN = false negative
Prevalence	$I\times D$	I = incidence, D = duration
Detection Rete	$\frac{TP}{TP+FN}$	TP = true positive, FP = false positive
Detection Prevalence	$\frac{T_{detection}}{Total}\times 100$	-
Balanced Accuracy	$\frac{Sensitivity+Specificity }{2}$	-

). Then, the data contained in the RGB images were visually reviewed to verify and validate the detection of cuts based on the applied methodology (Figure 5).

3. Results

3.1. Variation in NDVI

NDVI was calculated for each forest parcel (the two examples in Figure 6 and Figure 7 reflect the variability in the characteristics of the plots). The aim was to analyze the sudden changes observed through these values and to study whether they allow the automatic detection of spatially and temporally limited cuts. Therefore, of the 4231 total plots, there were 197,398 pixels in which NDVI was obtained in 53 different time periods.

The temporal information obtained for each plot allowed the correct overall analysis of the biomass of each plot. The different felled plots were mainly reforested by Pinus and Eucalyptus, although in many plots, other species (mainly deciduous hardwoods, scrub, grasses, etc.) were mixed that were not felled and that had a different evolution over the course of the year and a different seasonality. In addition to displaying the felling characteristics of plots, obtaining the temporal analysis made it possible to reveal the phenological changes of the species in each plot. This phenomenon can be clearly observed in Figure 5a,b, where along the temporal analysis, slight fluctuations were caused by phenological changes. Finally, the last images of the time series show the felling of the plots. Plot N°2 shows a partial cutting in October 2019 and the recovery of the vegetation in the following months.

3.2. Detection of the Forest Cutting Threshold

From all the images analyzed, 53 time periods were analyzed for each plot, with a total of 197,398 analyzed pixels. With the aim of explaining the analysis carried out, one pilot plot was selected to show the detection capacity of the applied methodology (Figure 8a). A histogram representing the frequency of each of the pixel values along the analyzed time series (Figure 8b) confirmed the previously selected cut detection value (−0.1254). In addition, with the boxplot analysis (Figure 8c) of each pixel found within the plot, the variations in their values along the time series allowed the detection of the time at which the felling took place.

From the 135 selected plots, the decreases in NDVI in each period due to a cut or a “not cut” (caused by phenological changes) were detected. Moreover, if sufficient contiguous images were available, it was possible to estimate the duration of forest action. Estimates of the cut area could even be made daily, as could estimates of the total final area of the felling action. In this example (Figure 9), four different plots were analyzed, and a total felling area of 4.42 ha was estimated to have occurred between the end of December 2019 and February 2020. In this example, the cut began on 16 July 2019 (Figure 9a), with 58 points showing felling due to differences in NDVI values between two periods, denoting the detection of a cut with an estimated surface of 0.57 ha, comprising the entire plot. The forest action continued on 4 December 2019 (Figure 9b), affecting 0.5 ha (92.6% with 50 points with decreased differences in NDVI). Another cut was detected on 26 December 2020 (Figure 9c), with a cut of almost 52% over the total surface area (154 points and 1.54 ha), ending on 24 February 2020 (Figure 9d), with 61% of the area cut (181 points and 1.81 ha).

The statistical values of the falls (trends) of the plots were considered and correlated for each period in the selected pilot plot (

Table 3. Confusion matrix used for the classification of forest cutting on the analyzed plots in the study area.

Confusion Matrix and Statistics
Accuracy	0.8862
95% CI	(0.8164, 0.9364)
No Information Rate	0.7236
p-value (ACC < NIR)	1.044 × 10−5
Kappa	0.7047
Mcnemar’s Test	0.4227
Sensitivity	0.7353
Specificity	0.9438
Pos. Pred Value	0.8333
Neg. Pred Value	0.9032
Prevalence	0.2764
Detection Rete	0.2033
Detection Prevalence	0.2439
Balanced Accuracy	0.8396

). In total, cutting was detected on 111 plots, and 150 falls from periods identified as “not cut” were also detected. In the obtained classification tree (Figure 4), each rectangle represents a node with a classification rule. The rectangle of each node shows the proportion of cases belonging to each category and the proportion of the total data that are grouped together. For example, the rectangle at the bottom left of the graph shows 96% of cases in Type 1 (Cutting) and 4% in Type 2 (Not cutting), representing 24% of all data. These proportions provide the accuracy of the model when making predictions. Thus, the rules leading to the rectangle described above provide 88% correct classifications (Table 3). The classification tree indicates that the value of the percentile decreases and that the values of the means are important as detection variables. This corroborates the hypothesis that the trend percentage (relation between the surface of the felling points and the surface of the plot) can be used for the detection of cuts. The confusion matrix indicates that the model correctly predicts 25 of the 34 total cuts, while 9 are classified as “not cut”. On the other hand, the model correctly predicts 84 cuts, while 5 are predicted as cut but are not.

4. Discussion

One of the main challenges faced when mapping land cover using time series images is the lack of continuity caused by cloud cover. Sentinel-2 is vulnerable to cloudy weather, making it difficult to obtain sufficient clear images for monitoring in temperate climate zones. However, as pointed out by Vuolo et al. [56], the widespread use of NDVI in remote sensing as a tool for monitoring forest areas has stimulated the development of models that reduce the noise caused by clouds. Puletti and Bascietto [57] conclude that due to the high resolution of the data provided by Sentinel-2, spectral variability over short time periods (5 days) could be considered negligible. Moreover, this information would be more accurate than that provided by Landsat, whose resolution is lower (30 m) and temporal resolution is higher (16 days) [58]. Taking into account that this study was carried out on a type of property characteristic of smallholding, the lack of more frequent images due to cloudiness results in a small decrease in the value of the technology. On such small plots, forest felling can be completed in a single day; the felling of contiguous plots can be completed in several days (Figure 6). Due to the cloudiness problem, the first cut could be detected days later, or it could be assumed that the cut was made in the same time range when more than one plot is cut. In this case, to reduce the noise caused by cloudiness, additional masking was performed in this study using the threshold values obtained for Band 2, as tested in other studies for the additional removal of dense clouds [59]. In this way, interference and cloud contamination were reduced, although it should be highlighted that in a real control system, it is not always possible to eliminate cloud contamination in images. Considering this possibility, this study analyses a long period of time in the absence of images, allowing us to evaluate the behavior of the method in a real control period. In any case, it was verified that a drop in NDVI was detected despite the existence of very small plots or the image being taken days after cutting, implying an acceptable temporal detection margin for sustainable forest management.

This result would also imply the need for high-capacity data storage and analysis. Therefore, whether it is necessary to create large databases and use processing systems that cover these bases must be resolved, taking into account the objective that is being pursued. The use of multitemporal data may have inadvertently increased the noise level in the classification process, as increasing the dimensionality of the data leads to greater redundancy [60]. It would not be implausible that only the use of Sentinel-2 scenes with optimal timing and atmospheric correction provides good results [61]. However, the use of these images would depend on the objectives of the study. In this case, the multitemporal factor would be the most effective way to address the monitoring of such small plots.

On the one hand, the use of multitemporal data allows phenological information to be obtained with greater precision and allows more effective monitoring in which invalid images are discarded [62]. Furthermore, the use of these images on managed land would allow anthropogenic changes to be detected more effectively [63]. On the other hand, the applicability and use of vegetation indices to detect changes in land use must be considered. These indices are very effective in large areas where changes in land use are analyzed, such as deforestation in the Amazon [64]. However, in this study, the objective was to analyze the use of these technologies in sustainable forest management in a large number of small plots (with different stands or different felling periods), specifically for the detection of the logging or cutting of forest stands. In the latter case, NDVI products have shown good results in multiple studies for different types of forest stands [65,66]. Although it has not been the objective of this work, the described method could also be used for the detection of changes in land use as well as for the detection of forest fires, issues that can be raised as possible future lines of research. In this work, since these are parcels that are part of a certification system, changes in land use are not possible, and forest fires have not interfered in the study since no parcel registered this type of event.

In the present study, this index provided data that were consistent with the forest stands of each plot; in all systems, the index showed very similar patterns within the forested areas (Figure 4). The multitemporal NDVI evolution describes the intra-annual patterns in temperate forest climates [67], showing slight intra-annual and interannual variations, consistent with other studies. At the start of the spring season around Apri–May, NDVI increases until it reaches a maximum in July–August, even extending until October. The use of images from the four seasons facilitates the distinction of vegetation classes with different phenological cycles (for example, conifers and deciduous trees) [68]. These interannual variations were characterized by a rapid increase during the initial growth of each stand, from a minimum postharvest value (a few weeks after harvest) to a first maximum value that occurred during the first or second year after planting [69]. Figure 4b shows the results of a partial cutting that occurred around October 2019, decreasing the NDVI value, which recovered in the following months at different speeds depending on the vegetation present in the plots. This rapidly increasing NDVI index in plots where cutting occurred could be due to the presence of shrubs, which would explain why the values remained lower at the end than at the beginning of the study, reaching the highest point in spring 2020. On the other hand, Figure 4—N°1 shows the NDVI index for a Eucalyptus monoculture plot. This tree is evergreen, so it is normal that NDVI remained constant in this stand, with small phenological variations, until the time of cutting. Hua Lu et al. [70] found that the average NDVI values varied from 0.6 to 0.8 for Eucalyptus forests, similar to the values obtained in this study (1.0 to 0.8). Bare soil represented NDVI values close to 0, and water bodies represented negative NDVI values [71]. In this case, both sample plots had minimum values of 0.3–0.4, indicating the presence of undergrowth on the plots. The high variabilities in shrubland species induce higher variability in NDVI and facilitate the acquisition of greater variations among seasons [72].

From a holistic perspective in which the type of felling is not taken into account, but the small-scale management of a large number of small plots in a large forest area is facilitated, an accuracy of 88% was achieved. These results were similar to those obtained by [73], wherein the authors obtained a detection accuracy of approximately 90% for clear cuts in the size category above 15 ha, similar to studies that employ the Random Forest algorithm, 90% [74] and 98% [75]. In the case of this study, the methodology was specifically designed for small plots for which the average area was 0.47 ha and the maximum surface area was 3 ha; therefore, this difficulty found by other investigations was corrected. On the other hand, the Kappa index, which is used to assess concordance or reproducibility, provided a value of 0.70. Persson et al. [76], using Sentinel-2, obtained a maximum value of 0.84 and a minimum of 0.63, depending on the months used, from which it can be deduced that each season affects the analysis due to phenological variables. In our case, wherein the methodology covered 3 years and each season, it can be deduced that an optimum value was obtained. Furthermore, this study obtained similar values to Forstmaier et al. [77], where a model was used to predict Eucalyptus cover, with a sensitivity of 75.7% and a specificity of 95.8%.

The methodology developed herein offers a new opportunity that would improve efficient forest management as well as land-use planning. The studies focus on large plot sizes, demonstrating the great potential of using medium resolution temporal imagery [78]. However, there is no literature analyzing small plot sizes and their significance. Some studies analyzing medium-sized plots (28.4 ha) report overestimation [79]. Moreover, the methodologies are tedious and complicated, even combining three types of remote sensing data [80], and are hardly applicable by the average forest user. This study has selected the most effective and straightforward approach for small plots (<3 ha), where its control and management are more difficult. Therefore, this method can be assumed and applied for forest management by companies, owners or certifiers. Currently, in the study region, the data on the felling date came from loggers. This information is usually incomplete or wrong. Therefore, offering an automated system, which can save time and decrease error, can improve the wood handling and control system. One of the advantages of the methodology employed in this study is that it is based on pixels and the use of spectral signatures and allows the application of indices and vegetation. In this case, the index used was NDVI due to its highly demonstrated good performance. This study did not have the general objective of detecting phenological changes in the studied plots, but in this sense, NDVI demonstrated good effectiveness in correctly monitoring forested areas. This method provides the opportunity to detect problems caused by pests or other anthropogenic pressures by analyzing declines in the index value. According to the articles reviewed in which Sentinel-2 images were used, high efficiencies were obtained, but the adaptation of this method to the small plots that make up the study area was improved in this study. Therefore, this methodology would allow the creation of automation that detects felling assigned to a specific management system. In this way, resources would be saved when field data are collected through “in situ checks”. In addition to speeding up the detection capacity, this process would facilitate the management system’s ability to react. This application would allow monitoring to improve the achievement of forest policy objectives as well as the implementation of these detection systems in decision making for forest planning.

Regarding the required processing time, it is important to note that this time depends on the performance of the computer, the software solutions and the specific knowledge and skills of the operator. In addition, the operator’s knowledge and skills in (visual and digital) image interpretation and analysis, as well as knowledge of the characteristics and conditions of the observed environment, are necessary. On the other hand, the applied procedure is considered to be fast, facilitating the rapid monitoring of a large study area consisting of multiple plots with different compositions and characteristics.

5. Conclusions

Based on the results obtained herein, the use of the Sentinel-2 satellite can be considered a valuable tool for forest monitoring, especially considering that these data will be routinely available for many years and that the images are frequent and free of charge.

On the other hand, many of the abrupt changes obtained in the NDVI series of a given pixel are due to cutting in most cases. NDVI could identify and quantify changes in the analyzed forest parcels. The approach presented in this study could be used in other forest stands to monitor anthropogenic changes in small forest plots. An accuracy of approximately 95% was achieved when the whole plot was cut, with an accuracy of 81% obtained for partial cuts. At the same time, the methodology used herein also facilitated correct monitoring by the managers of these plots and effectively detected phenological changes, making it possible to assess the effects of changes in forest stands.

Finally, the presented workflow demonstrated the efficiency of the performed approach. This methodology would allow future works to improve and streamline, from an operational perspective, the processing of these images, allowing the reduction of processing time and storage capacity. This methodology is particularly useful in large research areas, in areas of multiple plots with different forest species and in areas where the property system is characterized by smallholdings.