A Low-Cost and Robust Landsat-Based Approach to Study Forest Degradation and Carbon Emissions from Selective Logging in the Venezuelan Amazon

Selective logging in the tropics is a major driver of forest degradation by altering forest structure and function, including significant losses of aboveground carbon. In this study, we used a 30-year Landsat time series (1985–2015) to analyze forest degradation and carbon emissions due to selective logging in a Forest Reserve of the Venezuelan Amazon. Our work was conducted in two phases: the first, by means of a direct method we detected the infrastructure related to logging at the sub-pixel level, and for the second, we used an indirect approach using buffer areas applied to the results of the selective logging mapping. Pre- and post-logging forest inventory data, combined with the mapping analysis were used to quantify the effects of logging on aboveground carbon emissions for three different sources: hauling, skidding and tree felling. With an overall precision of 0.943, we demonstrate the potential of this method to efficiently map selective logging and forest degradation with commission and omission errors of +7.6 ± 4.5 (Mean ± SD %) and −7.5% ± 9.1 respectively. Forest degradation due to logging directly affected close to 24,480 ha, or about ~1% of the total area of the Imataca Forest Reserve. On average, with a relatively low harvest intensity of 2.8 ± 1.2 trees ha−1 or 10.5 ± 4.6 m3 ha−1, selective logging was responsible for the emission of 61 ± 21.9 Mg C ha−1. Lack of reduced impact logging guidelines contributed to pervasive effects reflected in a mean reduction of ~35% of the aboveground carbon compared to unlogged stands. This research contributes to further improve our understanding of the relationships between selective logging and forest degradation in tropical managed forests and serves as input for the potential implementation of projects for reducing emissions from deforestation and forest degradation (REDD+).


Introduction
More than 400 million hectares (ha) of natural tropical forests have been designated as production forests globally [1][2][3]. Moreover, about 40% of sawn wood traded annually in tropical regions has an origin in natural forests [4], often under a "selective logging" approach in which large trees of a relatively low number of tree species are harvested in rotation cycles of 30 years on average [2,5,6]. With some exceptions, one of the main features of selective logging across the tropics has been the insufficient adoption of reduced impact methods with negative environmental effects on forest structure and function [7,8].
Forest degradation is a change process caused by anthropogenic and/or environmental forces that result in alterations within any given forest, negatively affecting the structure and or function of the stand and site, and thereby lowering the capacity to sustain a continuous supply of products and/or services [9]. Selective logging, a major degradation driver in tropical forests, may account for at least half of the total carbon emissions coming from tropical forest degradation, representing~6% of total tropical greenhouse gas emissions [10,11]. Logging can cause significant losses in the aboveground carbon in tropical forests up to 10-50% [12,13], and in some cases can add up to 123% more forest-area damage above what has been reported for deforestation alone [14]. Moreover, in tropical regions of Latin America and Asia, selective logging has been responsible for more than 70% of the area of total forest degradation [15], and is often consider a preamble to deforestation processes [16][17][18]. Consequently, in recent years, forest degradation has become a major topic of discussion within the international scientific community as a factor of great importance in the global carbon cycle [11,19,20]. This process has been addressed in the United Nations Framework Convention on Climate Change (UNFCCC) since the Bali agreement in 2007 (CP. 13), when the concept of reducing emissions from deforestation (RED) was expanded to REDD+ to include forest management activities [21], and currently with the so called Nationally Determined Contributions to the UN Paris Climate Agreement [22].
Venezuela has a long history of natural forest management under long-term concession contracts that started in the early 1970s [23,24]. By 1992, almost 3.2 million ha were allocated to more than 30 forest management units (FMUs) and had management plans approved by the national government [25]. However, a remarkable decline in the use of forests as reliable sources of timber has been evident in the last 30 years. According to the last available official data from 2018, only about 2.5% of the wood legally consumed in the country came from FMUs in an estimated area of 246,313 ha of forests with formal management plans mostly in the Amazon region of the country [26].
Forests in the Venezuelan Amazon account for at least 83% of the national forest cover [27]. In this region, forests have been managed via various legal protection schemes, from strictly protected areas such as national parks to others with sustainable use objectives as the case of forest reserves, forest lots and forest areas under protection [28,29]. One of the most important reserves in the region is the Imataca Forest Reserve (IFR), which harbor important levels of biodiversity [30]. In addition, different studies have estimated that the region represents one of the most carbon-rich areas in Venezuela with an average of 205 ± 15 Mg C ha −1 in the aboveground biomass (AGB) [31,32]. At the same time, the IFR is considered a hot spot of deforestation [33] and was recently labelled as one of the deforestation fronts of the tropical belt [34]. In the last two decades, approximately 1144 ha of forests have been lost every year due to mining activities (55.5%), land use change for livestock (26.5%) and agriculture expansion (17.9%) [35]. Yet, from the close to 4 million ha that IFR covers, about 97% is still covered by different forest-types, mostly dominated by lowland "terra firme" forests [35].
Selective logging, via legal, yet mostly unplanned conventional timber harvesting operations was formally authorized in IFR around the 1980's decade and has become a major factor of anthropogenic disturbance since then [35]. In a standard logging operation in this region, large trees with diameters at breast height (dbh) > 40 cm and with an average height between 20−30 m from a few commercially valuable species are harvested, where for every tree harvested close to 11 additional trees can be severely affected [36]. Moreover, without adequate planning, the impacts of these operations can double the background rates of tree mortality when compared to unlogged stands while also causing a significant reduction in aboveground carbon [13].
To quantify carbon emissions caused by forest degradation as a result of selective logging, two analytical techniques can be used: the first combines logging rates, management plans, and high-resolution imagery for activity data (AD) and the gain/loss approach for emission factors (EFs) [10,37]. The second combines remote detection of medium-resolution images for AD and an assessment of the changes in carbon stocks for EFs [38,39]. For the second technique, the AD can be obtained by a direct method, identifying and mapping canopy damage [40][41][42][43][44][45], or by mapping canopy damage in combination with intact forests and regeneration patches [46,47]. Using an indirect method that implies identifying selective logging in the images, and the addition of buffer areas via geographic information system (GIS) tools, the area of degraded forests due to logging can be quantified [18,48,49]. With this method, the intact forest concept and evaluation criteria are applied to categorize intact and non-intact forests, which discriminate against forests with different carbon stocks [16,46]. The EFs can then be obtained using data collected from field inventories measured before and after logging, with the general change in the carbon reserves calculated using the difference between these two measurements [39].
In this study, we propose an analytical approach based on a Landsat time series [40,41,50,51], developed for the Monitoring System of Deforestation of the Amazon (TerraAmazon) [52]. One of the main goals of our approach is to produce a local-based computing analytical approach capable of functioning under conditions of limited connectivity. We believe this to be an important advantage to conduct an assessment of forest degradation produced by logging in Venezuela's Imataca Forest Reserve (IFR) for a 30-year period between 1985 and 2015 under different conditions. Overall, because of the general lack of reduced impact logging guidelines being applied [29], we expected that despite a relatively low logging intensity typically between 10 to 12 m 3 ha −1 [36,53,54], compared to an average of 20 m 3 ha −1 for the Amazon [55], carbon losses could be significant. We put our results in the context of other studies where the effects of logging on carbon have been addressed.

Study Area
Imataca Forest Reserve (IFR), located in southeastern Venezuela between Bolívar and Delta Amacuro states, was officially created in 1961 and has a total area of 3,821,900 ha, which represents 8.1% of the total area of the Venezuelan Amazon (~46.9 M ha) ( Figure 1a). The reserve has been divided into 23 different management units distributed across three major zones (north, central, and south) (Figure 1b). 90% of these units with extensions between 120,000 and 340,000 ha have been designated for permanent forest production, while the rest is officially allocated for mining activities and/or conservation of biodiversity. Between 1985 and 2012, about half of the timber-production areas at IFR was managed under a private concession model where national government granted management rights after an official management plan was approved with cutting cycles ranging from 25 to 40 years [29]. In recent years, with the enactment of the Forests and Forest Management Law of 2008, a policy shift began with regards to how forest management should be planned and applied in Venezuela. Along with government agencies and the newly created National Forest Company (ENAFOR), guidelines for developing new forest management plans were put in place to gradually shift from the model of private concessions to a more government-dominated approach. At present, the company supervises the management for all production forests in the country and is directly responsible for an active management operation in the Imataca Forest Reserve.
The IFR has a northeast-southwest pattern in the distribution of precipitation from 1000 mm to 3000 mm per year approximately. The average annual temperature is around 25 and 27 • C, evapotranspiration ranges from 1250 mm to 1400 mm per year. Overall, in this area we find lowland tropical humid forests, seasonal evergreen forests, deciduous forests and swamp forests. From the standpoint of species diversity around 2800 species of plants, 450 species of birds, 153 species of mammals, 90 species of reptiles, 62 species of amphibians and 242 species of fish have been identified [30].
logging inventories conducted at seven of these compartments namely: research plot I (RP I), research plot II (RP II), experimental development plot (EDP), compartment 1 (C1), compartment 2 (C2), compartment 3 (C3) and compartment 4 (C4). After operations were halted by the end of 1990's, the management was later transferred in 2012 to the Empresa Nacional Forestal (ENAFOR), where two additional compartments were allocated for timber harvest: Santa Maria I (STM I) and Santa Maria II (STM II) (Figure 1c). Depending on the characteristics of the overall planning process and the harvest intensity, selective logging at Unit V can be classified into three types [18]. First, unplanned conventional logging (CL) was characteristic of RP I, RP II and EDP between 1985 and 1988 without a formal management plan. Secondly, planned managed logging (ML) occurred in two ways: in the first (ML1), a pre-commercial inventory of trees was carried The management Unit V located in the central region of IFR was selected as a study case. The unit has a total area of 180,000 ha and was first granted management rights to the company Industria Técnica de Maderas C.A (INTECMACA) in 1982. The unit was originally divided into 25 logging compartments of an approximate equal area each according to a 25-year cutting cycle. One component of our work is based on pre-and post-logging inventories conducted at seven of these compartments namely: research plot I (RP I), research plot II (RP II), experimental development plot (EDP), compartment 1 (C1), compartment 2 (C2), compartment 3 (C3) and compartment 4 (C4). After operations were halted by the end of 1990's, the management was later transferred in 2012 to the Empresa Nacional Forestal (ENAFOR), where two additional compartments were allocated for timber harvest: Santa Maria I (STM I) and Santa Maria II (STM II) (Figure 1c).
Depending on the characteristics of the overall planning process and the harvest intensity, selective logging at Unit V can be classified into three types [18]. First, unplanned conventional logging (CL) was characteristic of RP I, RP II and EDP between 1985 and 1988 without a formal management plan. Secondly, planned managed logging (ML) occurred in two ways: in the first (ML1), a pre-commercial inventory of trees was carried out, followed Remote Sens. 2021, 13, 1435 5 of 25 by a general planning of logging roads and landing sites. These activities were part of a more detailed management plan in which each compartment should be-in theory-logged every year. This method was applied in C1, C2, C3 and C4 between 1990 and 1995 [56]. The second ML method (ML2) is in many ways similar to the ML1 case, especially with regards to how logging operations were applied. However, a major difference is that planning was organized at the landscape scale, thus watersheds and small-watersheds were used as management units. Two additional features that were also unique for this approach were that, on one hand, commercial trees were spatially mapped to facilitate planning of logging roads and landings. On the other, the minimum harvest diameter was modified from a common threshold of 40 cm for all species to dbh > 50 cm for high-wood density (WD) species, 60 cm for medium WD and 70 cm low WD species. This method was used in the case of STM I and STM II between 2012 and 2015 [57]. The mean area of each logging compartment was 2640 ± 848.5 ha (Mean ± SD).

Landsat Time Series
Fifty (50) Landsat 4, 5, 7 and 8-time series images were used, corresponding to route 233 and row 55 (Table 1). These were obtained from the collection of the US Geological Survey (http://glovis.usgs.gov/ (accessed on 10 March 2021)), with a processing level L1T. The time period for these datasets was selected approximately between one and two years after logging occurred, as rapid canopy closure after disturbance and lower understory revegetation may inhibit logging detection [41,58,59].

Field Data
Two independent datasets composed by information collected from temporary and permanent ground plots were used in support of the analysis:

•
The first group was obtained from INTECMACA inventories conducted between 1986 and 1995. Two groups of permanent plots were established in RP II: four plots of 0.5 ha (100 × 50 m) located in unlogged forests and four plots of 2 ha (500 × 40 m) in logged forests. The data include all living trees with diameters at breast height (dbh) > 10 cm [13,60].

•
The second set of data was obtained from ENAFOR inventories conducted between 2012 and 2015. A group of temporary and permanent plots was systematically established in the STM I and STM II logging subunits; a total of 65 plots of 1 ha (1000 × 10 m) with subplots of 0.01 ha (10 m × 10 m) were measured in the pre-and post-logging periods [61,62] ( Figure A1 in Appendix A). In all cases, a complete taxonomic identification was made to every individual to account for species composition and diversity.
In addition, official reports from the logging companies were used to collect data on the number and size of the harvested trees (diameter, height and species) [63]. Total volume per tree was estimated using the Smalian scale formula (cm 3 ), so total volume of wood harvested at the compartment level could be estimated.

Analytical Approach
Our analytical approach consisted of five different phases as follows: we first mapped selective logging, followed by a validation process of the resulting maps. The third phase consisted of the construction and validation of the forest degradation maps, followed by the estimation of aboveground biomass (AGB) and carbon to close with the estimation of committed carbon emissions (CCE) (Figure 2). in logged forests. The data include all living trees with diameters at breast height (dbh) > 10 cm [13,60].

•
The second set of data was obtained from ENAFOR inventories conducted between 2012 and 2015. A group of temporary and permanent plots was systematically established in the STM I and STM II logging subunits; a total of 65 plots of 1 ha (1000 × 10 m) with subplots of 0.01 ha (10 m × 10 m) were measured in the pre-and post-logging periods [61,62] ( Figure A1 in Appendix A). In all cases, a complete taxonomic identification was made to every individual to account for species composition and diversity.
In addition, official reports from the logging companies were used to collect data on the number and size of the harvested trees (diameter, height and species) [63]. Total volume per tree was estimated using the Smalian scale formula (cm 3 ), so total volume of wood harvested at the compartment level could be estimated.

Analytical Approach
Our analytical approach consisted of five different phases as follows: we first mapped selective logging, followed by a validation process of the resulting maps. The third phase consisted of the construction and validation of the forest degradation maps, followed by the estimation of aboveground biomass (AGB) and carbon to close with the estimation of committed carbon emissions (CCE) (Figure 2).

Mapping Selective Logging
The TerraAmazon system was used to map selective logging. This system involves a configuration that includes the creation of a PostgreSQL database, definition of the conceptual model, access control, phase control, project and control rules, definition of classes, definition of the control rules, and the definition of the control area [64]. The Landsat time series datasets were then exported to generate a linear spectral mixing model (LSMM) for Remote Sens. 2021, 13, 1435 7 of 25 each image to be used in the detection of selective logging [14,43,45,48,49]. The LSMM uses the red (0.63-0.69 µm), near infrared (0.76-0.90 µm) and mid-infrared (1.55-1.75 µm) bands from Landsat 4, 5, 7 and 8 [44], from which the samples of soil cover, vegetation and shade can be extracted to estimate the proportions in each pixel and in their respective images (Equation (1)). The soil fraction image was then used to estimate the area affected by selective logging [43,44,48,49,65].
where r i = is the response of the pixel in band i; a, b, and c = proportions of vegetation, soil, and shade, respectively; vege i , soil i and shadow i = spectral responses of the components of vegetation, soil and shade, for each band respectively; e i = is the error in band i.
A cloud mask and cloud shadow were applied to each image of the soil fraction, using thresholds of the minimum and maximum values of the blue band (0.45-0.51 µm) to detect shadows and the infrared thermal band (10.60-11.19 µm) to detect clouds [64,66]. Selective logging was detected by means of a binary classification of areas with and without selective logging based on the processed images of the soil fraction [44,49]. A value of zero was assigned to those pixels with soil fractions lower than 37% (i.e., areas without evidence of selective logging), and a value of one corresponded to the pixels with soil fractions between 37 and 100% (i.e., areas with signs of selective logging). A decision tree algorithm was then used, for which the 37% limit was statistically defined based on 150 points visually selected from a processed soil fraction image [49]. Once the binary images were obtained, these were added to generate the mapping of selective logging.

Validation of the Selective Logging Maps
To determine the quality and degree of agreement between the mapping of selective logging and field conditions, maps were validated via the comparison with an external source that is considered a realistic representation of the characteristics on the ground [67][68][69][70]. Thus, we applied a systematic sampling approach [50,71] to 36 blocks of 100 ha (1 km 2 ) [72], which represented approximately 11% of the study area. This sampling technique allowed us to precisely and quickly estimate the error of the analysis [67,73]. By using a visual on-screen interpretation of a minimum cartographic area of 1 ha in each sampling block [72] we can generate logged and unlogged forest datasets ( Figure A2 in Appendix A). These were considered the ground-truth data and were then used to analyze the thematic quality of the selective logging map. A confusion matrix was generated and the errors of omission and commission with the level of global accuracy were also calculated.
To confirm the logged classification, a spatially non-localized analysis was used by comparing the proportions of the area in each sample block in the selective logging map and in the ground-truth data. These proportions were compared via simple linear regression, with the rationale that if the mapping of selective logging and ground-truth data were similar, adjustment values would be high and the coefficient (R 2 ), would be close to 1 [67,74]. This method has been widely used in the mapping of forest fires [75,76], analyses of land use [74] and deforestation [77].

Construction and Validation of the Forest Degradation Maps
Using an indirect method, mapping of forest degradation was performed for each logging compartment. In doing so, we estimated the average radius between log landings in each of the soil fraction images as proposed by Monteiro et al. [48]. In our case, this value was 600 meters (m), so 300 m was used as a threshold to estimate the approximate area of forest degradation caused by logging, and a square buffer was applied to the mapping of Remote Sens. 2021, 13, 1435 8 of 25 selective logging using GIS tools. Maps of forest degradation were validated by comparing these with the area logged reported in the management plans, which allowed for the calculation of commission and omission errors [68].

Estimation of Aboveground Biomass (AGB) and Carbon
Forest inventory data for pre-and post-logging conditions were used to estimate AGB and carbon, following the approach by the Global Observation of Forest and Land Cover Dynamics panel (GOFC-GOLD) [46] for the establishment of REDD+ projects. The AGB per tree in each ground plot was estimated using the pantropical allometric regression from Chave et al. [78] (Equation (2)). All the estimates were generated for each plot and the estimations were scaled to 1 ha when necessary. Values of aboveground carbon (Mg C ha −1 ) were assumed to be 50% of AGB [79].
where: AGB = is the aboveground biomass of the individual trees expressed in kilograms (kg); E = is a water stress factor that shows an important covariance with the diameterheight ratio in tropical trees and includes information on seasonal temperature (ST) and the climatic water deficit (CWD). Based on the geographic location of each plot, E and CWD were derived from a 2.5 arc-minute resolution raster file available at http://chave.ups-tlse. fr/pantropical_allometry.htm. (accessed on 10 March 2021); ρ = is the density of the wood in g cm −3 , with data assigned for each taxonomic group from the pantropical database of Zanne et al. [80] and Chave et al. [81]; D = is the diameter of each tree in cm.

Estimation of Committed Carbon Emissions (CCE)
A stock-difference method was used to estimate emissions related to selective logging following the 2019 refinement of the 2006 IPCC guidelines [39]. In doing this, the following assumptions were considered:

•
To simplify the carbon accounting process, the committed emissions approach was used, in which all carbon removed is assumed to be emitted at the time of its removal via logging [37,39]. • Emissions were estimated in each compartment by multiplying the area affected by degradation (activity data) with the difference of the carbon content in the pre-and post-logging period (emission factor) [39].

•
The different harvesting activities were classified in the selective logging map into: log landings, caused when the forest is cleared for the purpose of temporary log storage before final transportation; logging roads, built to transport timber from log landings to sawmills; and logging gaps, created by tree felling and skid trails, resulting in damage or death to other standing trees [18,46]. These categories were associated with the emissions in each compartment to determine the overall emission contribution for each activity.

•
Using the reported values of timber extracted from each compartment, carbon losses from logging were estimated by calculating the equivalent carbon of the volume of extracted roundwood, which considered the wood specific density to obtain AGB. A factor of 0.5 was used to estimate the amount of carbon [79,82]. impact factor (CIF) (Mg Mg −1 ) also called "mean carbon export ratio" [82], represents the emissions of each compartment relative to the emissions of the total volume of extracted roundwood.

Accuracy of Maps of Selective Logging and Forest Degradation
The results from the error matrix are presented in terms of area proportions. The row totals of the error matrix (Sum = 0.113) provide what was represented in the ground-truth samples in each class and constitutes the total proportion of selective logging. Conversely, the column totals provide the estimated proportions according to the ground-truth data, and this was estimated in 0.013 for the logged class. Multiplying this value by the total number of pixels on the map (i.e., 330,386), the result is 4158 pixels, or approximately 374 ha. The area of agreement between the mapping of selective logging and the ground-truth data was 2452 pixels or close to 221 ha. There was an underestimation of 1706 pixels or around 154 ha that was confused with the unlogged class. On the other hand, the logged class had the highest user error (commission), with a proportion of 0.146, and a producer error (omission) of 0.4103 because greater proportions of areas were included and excluded in this class, respectively ( Table 2). The proportion of these errors was observed in nine of the 36 blocks that were used in the validation, three for the case of ENAFOR and six for INTECMACA ( Figure A3 in Appendix A). The global precision of the selective logging map was 0.943. The values of the ground-truth proportions and those of the mapping of selective logging showed an overall good fit. The proportions presented high similarity, as indicated by a high value in the coefficient of determination after a simple linear regression (R 2 = 0.82) ( Figure A4 in Appendix A). In relation to the area of forest degradation obtained from the map, this represented approximately 13.6% (24,484 ha) of the entire area of Unit V. Analyzing this area for each compartment, six had commission errors, ranging from +1.6% in compartment C2 (71 ha) to +14.7% in compartment RP II (325 ha), while three compartments had errors of omission, ranging from −0.3% in compartment EDP (5 ha) to −17.7% in compartment C1 (292 ha) ( Figure A5 in Appendix A). Figure 3 shows the full process to detect and outline selective logging activities for the case of the EDP compartment (i.e., experimental development plot), logged during the first three months of 1988. m ha −1 in length, and a disturbed area of roads close to 583 ± 109 m 2 ha −1 . Finally, the average area affected in logging gaps was estimated at 174 ± 111 m 2 ha −1 (Table 3). Spatial patterns of the different areas affected by logging were highly heterogeneous across all compartments. More regular arrangements, especially in terms of road building and distribution, seem clearer in the case of the oldest logging operations (e.g., RP1, RP2, EDP) in contrast with the most recent logged areas (i.e., STM I, STM II) (Figure 4). Stepwise process for the mapping of selective logging for the case of compartment EDP. Composition 453 (a) which was used to generate the soil fraction of the linear spectral mixing model (LSMM), a cloud mask (b) was then applied, and later a decision tree classification was used to obtain a binary image of the logged and unlogged areas (c). This image was then categorized to analyze the three main logging activities (d). Stepwise process for the mapping of selective logging for the case of compartment EDP. Composition 453 (a) which was used to generate the soil fraction of the linear spectral mixing model (LSMM), a cloud mask (b) was then applied, and later a decision tree classification was used to obtain a binary image of the logged and unlogged areas (c). This image was then categorized to analyze the three main logging activities (d).

Mapping Selective Logging at Imataca Forest Reserve
The total area directly affected by selective logging-related activities was around 2535 ha, from which compartment C2 had the largest area with 539 ha (21.2% of total area) and compartment C1 the lowest with 89 ha (3.5%), with an overall mean area of forest degradation of 282 ± 124 ha (±SD). The number of log landings built was, on average, 2.1 ± 0.6 (±SD) per 100 ha, with a total area of 2926 ± 496 m 2 and an area affected by logging of 59 ± 16 m 2 ha −1 . In terms of the extension of logging roads, we estimated an average of 16 ± 2.3 m ha −1 in length, and a disturbed area of roads close to 583 ± 109 m 2 ha −1 . Finally, the average area affected in logging gaps was estimated at 174 ± 111 m 2 ha −1 (Table 3). Spatial patterns of the different areas affected by logging were highly heterogeneous across all compartments. More regular arrangements, especially in terms of road building and distribution, seem clearer in the case of the oldest logging operations (e.g., RP1, RP2, EDP) in contrast with the most recent logged areas (i.e., STM I, STM II) (Figure 4).

Analysis of Forest Degradation
The total area of forest degradation caused by selective logging activities was estimated in 24,484 ha. The average area of forest degradation by compartment was 2720 ± 911 ha (Mean ± SD), with a maximum value of 3457 ha for compartment C4 and a minimum of 1360 ha in compartment C1 (Figure 4; Table 4). In total, 65,036 trees were harvested in all compartments with a total volume of 247,034 m 3 , for a mean logging intensity of 2.8 ± 1.2 trees ha −1 or 10.5 ± 4.6 m 3 ha −1 . Average ground area damaged per tree logged was estimated in 334 ± 121.9 m 2 , equivalent to 91 ± 41.9 m 2 for each cubic meter of timber harvested. The mean of the total area affected by logging represented 8.2 ± 1.8 % of the entire degraded area, 0.6 ± 0.2 % of this area corresponded to the construction of log landings, 5.8 ± 1.1 % to roads, and 1.7 ± 1.1 % to logging gaps (Table 4).

Aboveground Carbon Density
The values of forest carbon density for pre-and post-logging varied in the different compartments, with the pre-logging ranging from 268 Mg C ha −1 (C4 compartment) to 312 Mg C ha −1 (C2 compartment), with an average of 308 ± 76 (± SD) Mg C ha −1 . After logging, aboveground carbon density varied from 222 Mg C ha −1 (RP I compartment) to 263 Mg C ha −1 (STM I compartment), with an average of 229 ± 30 (± SD) Mg C ha −1 . The average difference between the two conditions was estimated in 35% (Table 5).

Committed Carbon Emissions (CCEs) from Selective Logging
Mean CCEs per volume harvested was 2.9 ± 1.1 (mean ± SD) Mg C m −3 , those from construction of logging roads were 1.9 ± 0.   Figure 5B).  CCEs for all nine compartments were highly variable. The CCEs expressed in terms of the Carbon Impact Factor (CIF), ranged between a minimum of 5.0 and a maximum of 13.2 Mg C Mg C −1 for compartments C1 and EDP respectively, and their baseline was 8.6 ± 3.7 Mg C Mg C −1 . On a volume basis, CCEs ranged between 1.8 and 4.4 Mg m −3 for the case of C1, C4, STM I and EDP compartments, and their baseline was 2.9 ± 1.1 Mg m −3 . On an area basis, CCEs ranged between 21.9 and 86.1 Mg ha −1 for STM I and RP II compartments respectively, and their baseline was 61 ± 22 Mg ha −1 . The effect of the intensity of logging on CCEs was very similar in all three ways of reporting emissions ( Figure 6). In terms of the logging approach, logging intensity values for the EDP and RP I compartments, where a conventional unplanned (CL) approach was common, were relatively low (between 6 and 8 m 3 ha −1 ), but the effects on carbon were above the baselines of the three ways of reporting emissions. In the case of the RP II compartment, logging intensity was high (17.5 m 3 ha −1 ) but below the CCEs baselines expressed either by volume or CIF. Planned selective logging 1 (ML1) used in compartments C1, C3 and C 4 had values below the baseline for all three ways of reporting emissions. However, the intensities increased from 8.3 m 3 ha −1 for compartment C3, 13.6 m 3 ha −1 for compartment C4 and 18.1 m 3 ha −1 for compartment C1, while in compartment C2 this intensity was low (7.5 m 3 ha −1 ) but its CCE values were still the highest. For planned selective logging 2 (ML2) applied in the compartment STM I, all values when reporting emissions were below the baseline, with an intensity of 5.2 m 3 ha −1 , while compartment STM II the CCEs values by volume and CIF were above the baseline with an intensity of 8.4 m 3 ha −1 (Figure 6). CCEs for all nine compartments were highly variable. The CCEs expressed in terms of the Carbon Impact Factor (CIF), ranged between a minimum of 5.0 and a maximum of 13.2 Mg C Mg C −1 for compartments C1 and EDP respectively, and their baseline was 8.6 ± 3.7 Mg C Mg C −1 . On a volume basis, CCEs ranged between 1.8 and 4.4 Mg m −3 for the case of C1, C4, STM I and EDP compartments, and their baseline was 2.9 ± 1.1 Mg m −3 . On an area basis, CCEs ranged between 21.9 and 86.1 Mg ha −1 for STM I and RP II compartments respectively, and their baseline was 61 ± 22 Mg ha −1 . The effect of the intensity of logging on CCEs was very similar in all three ways of reporting emissions ( Figure 6). In terms of the logging approach, logging intensity values for the EDP and RP I compartments, where a conventional unplanned (CL) approach was common, were relatively low (between 6 and 8 m 3 ha −1 ), but the effects on carbon were above the baselines of the three  (Figure 6).

Potential of the Analytical Approach
In this study we were able to develop a reliable and useful method to analyze selective logging via a local network using the Amazon Deforestation Monitoring System (TerraAmazon), with an architecture developed in a client-server environment, and with a TerraLib database which uses PostgreSQL as a Database Management System [52]. This analytical approach allowed us to estimate the area of forest degradation caused by selective logging with a high level of accuracy that subsequently help estimating the amount of carbon emissions produced.
Our results show the potential for this method to efficiently map selective logging automatically as shown in other studies [44,48,49,84,85]. With a relatively high global precision (GP = 0.943) that is within the range recommended by the GOFC-GOLD [46] for the development of forest cover monitoring maps, and slightly higher than that reported in the Brazilian Amazon (0.92) [49], and a high determination coefficient (0.82) that shows a close fit between the model and real proportions of the logging classes [67,74,77], we believe this is a promising and powerful tool to study forest degradation in tropical countries with severe connectivity limitations as the case of Venezuela.
Nonetheless, our study also confirms that detecting and mapping forest degradation with optical remote sensing data is a complex task, because the pixels that indicate forest degradation are an intricated mix of different land cover with diverse signals (i.e., vegetation, dead trees, bark, tree branches, soil, shadows, etc.) [45]. In addition, evidence of logging can rapidly disappear in less than two years after logging due to canopy closure and understory revegetation potentially limiting the overall accuracy of the method [18,41,58].
Mapping of forest degradation caused by selective logging yielded acceptable results, with average commission and omission errors of +7.6 ± 4.5 % (mean ± SD) and −7.5% ± 9.1 respectively, well below the uncertainty threshold of ±20% of the estimate for area [46], and of the results reported for the Brazilian Amazon (18% and 20%) [18] and the entire Amazon region (12% and 32%) [40].
The use of a relatively simple GIS model where the average of the log yarding radius was used, as in other studies [48,49], obeys to the fact that the Landsat time series for this area of the Amazon was not robust enough. Large proportions of clouds and cloud shadows covering the images [86], and the fairly low number of images available during the 1980s and 1990s [87] limited the use of spectral indices specialized in forest degradation such as the Normalized Degradation Fraction Index (NDFI) [18,41] and Continuous Degradation Detection (CODED) [40], or other more specialized statistical indices such as the Forest Degradation Index (FDI), based on a Multi Criteria Decision Analysis (MCDA) approach using the Analytic Hierarchy Process (AHP) technique [88,89]. Indeed, a similar result was found in a recently study on forest degradation in the Amazon using the NDFI [40], where forest degradation was marginally detected in our study area, sometimes overlapping with deforested areas.
The radius of the storage yards used by this indirect method (300 meters) to determine the buffer areas is slightly lower than in the case of other forest-types with lower densities of commercial tree species in the central and southern Amazon where 350 m was used as a threshold [48]. Our number of 300 m is, however, higher than the 180 m buffer zone used in a study carried out in dense tropical forests in the south-central Amazon [49]. Although stands of Unit V at IFR are often dense to moderately dense with an average volume of commercial trees of 33 m 3 ha −1 [57], compared to 20 m 3 ha −1 for transitional forests [48] and the 38 m 3 ha −1 for dense forests [49], we argue that an intermediate threshold value can reflect the fact that logging intensity was overall low in this unit of the Imataca Forest Reserve (IFR).

Selective Logging Detection
The use of Landsat images with a spatial resolution of 30-m has been a common approach in other studies analyzing the effects of selective logging [41,45,51,[90][91][92], so these results are useful to compare the degree of agreement with our estimates. For instance, the mean mapped size of the log landings in our study area was 2926 ± 496 m 2 (mean ± SD), 17% higher than the originally planned area of 2500 m 2 (50 × 50 m) [56,61]. Despite this, the mean size mapped is within the range (1-4 pixels) of log landings detection in the soil fraction image [48], since approximately 3.3 pixels were detected.
In our study area, log landings were approximately nine times larger than those reported in southern Brazilian Amazonia (339 ± 31 m 2 ) for reduced impact logging (RIL) operations. Regarding the number of log landings, we found that for every 100 ha, a mean of 2.1 ± 0.6 landings was created, equivalent to a disturbed area of 59 ± 16 m 2 ha −1 . These estimates largely differ from those reported for RIL in the Brazilian Amazon by Feldpausch et al. [82], where for every 100 ha, a higher density of log landings was created (6.2 ± 0.4), but with a much lower area disturbed of 20.8 ± 1.2 m 2 ha −1 . We interpret these results as a reflection of the overall poor planning and practical operationalization of logging activities in the case of IFR in comparison with RIL operations in other parts of the Amazon.
With regards to the effects of construction of logging roads, the management plans indicate a maximum width of 10 meters for main roads and 5 meters in the case of secondary roads plus an extra 10 m portion at each side for shoulder and ditch purposes in both type of roads [56,61,93,94]. Our mapping analysis shows that the mean length of all roads mapped was 16 ± 2.3 m ha −1 and the disturbed area was 583 ± 109 m 2 ha −1 . If we consider that all the roads were established with an average of 30 m, the equivalent disturbed area that should have been expected would be around 483 ± 69 m 2 ha −1 , indicating a potential overestimation of~21%. This overestimation can be interpreted as a result of the inaccuracies in the measurements of these areas due to the effect of the pixel size [18].
The average disturbance of logging gaps indicates a large variability in the total area disturbed (i.e., 174 ± 111 m 2 ) likely a direct response of the different densities of commercial trees that can be found in this type of forests and the subsequent effects that felling and hauling one or more logged trees can have on the overall structure of the unlogged portion of the forest stands [36,94].

Relationship between Logging Intensity and Degradation
We found no significant differences between the CL and ML1 harvesting modalities (3.1 ± 1.3 trees ha −1 and 3.0 ± 1.2 trees ha −1 respectively). However, the ML2 modality was significantly different (1.8 ± 1.2 trees ha −1 ), likely a consequence of the increase made to the minimum harvest diameters (MHD) for this compartment [57]. Overall, logging intensity was relatively low and similar to other management units in the IFR [53,54]. Compared to other areas in the Amazon, logging intensity at IFR is lower (e.g., 4.4 trees ha −1 -Jackson et al. [95]; 4.5 trees ha −1 -Johns et al. [96]; 6.4 trees ha −1 -Verissimo et al. [97]), which can be explained by the differences in species composition and the abundance of commercial species among these areas.
Considering the different approaches of selective logging applied in our study case, the proportion of aboveground biomass (AGB) affected by logging was 35 ± 1.7% for the CL modality, below the 60% reported for other conventional logging cases in Amazonian forests [18,59]. In the case of planned logging, the average AGB damaged was 24 ± 8.5% for ML1 and 14.5 ± 9.7% for the ML2 case, above and below respectively compared with an overall 20% reported for this type of logging [18,59]. In general, all three modalities of logging were within the ranges of effects to living biomass (10-46%) found in other studies [97][98][99].
Of the total area of forest degradation (24,484 ha), the mean area affected by selective logging that was detected by the direct method was 8.2 ± 1.8%, slightly lower than the 10.2 ± 1.2% [100] and 13 ± 4.5 (SD) % [82], reported for the eastern and southern Brazilian Amazonia respectively, where reduced impact techniques were applied.
Logging disturbances were not homogeneous in each compartment and within the different types of logging. In some cases, however, by increasing the minimum harvest diameter (MHD), the logging intensity declined but not the proportion of area disturbed, such as the case of ML2 where MHDs were stratified by specific wood density. A potential explanation for these discrepancies can be found in the fact that while a lower number of trees were harvested, their average size also increased. Thus, without an adequate planning for liana removal prior to logging or the sufficient application of directional feeling during harvesting there is potential for a much larger disturbance effect [101].

Pre and Post Logging Aboveground Biomass (AGB) and Committed Emissions (CCE)
If the AGB is averaged prior to logging, the value obtained is 308 ± 76.3 Mg ha −1 , which coincides with several studies in humid tropical lowland forests of the Venezuelan Amazon [12,13,31,32,102]. Logging across all compartments reduced pre-disturbance AGB, on average, by 35% with higher values found for the unplanned logging approaches as expected (Table 3; Figure 6). While the characteristics of the logging at each compartment are unique, it is encouraging to see a lower reduction in carbon losses when selective logging includes a better preparation. However, it must be noted that a lower damage is not always a response of better planning. Instead, it can also be a direct response to the spatial aggregation of commercially valuable timber trees along with topographical conditions and other biophysical/economic factors that are particular at each site and at the time when logging occurred (e.g., 1985 vs. 2015). For instance, in a recent review conducted across the tropics, Putz et al. [101] found that an average of 57% (range 22-97%) of the area in logging blocks was not directly affected by timber harvests, with more forests being left intact in areas farther from roads, on slopes >40%, and within 25 m of perennial streams. In addition, our study, as many others focusing on the impacts of selective logging in the tropics (e.g., Verissimo et al. [97]; Gerwing, [98]; Veríssimo et al. [99]) is mostly based on aggregated means of logging intensity, and that often can be a relatively weak reflection of the conditions on the ground [101].
In our analysis, CCEs both by sources ( Figure 5) and those expressed by area, volume and impact factor ( Figure 6) were higher compared to other parts of the tropics (e.g., [98][99][100]). This corroborates what has been discussed regarding the intensity and the logging methods used across these areas of the IFR, highlighting the widespread low efficiency compared to other cases where reduced impact logging is formally applied. Although beyond the scope of our work, the absence of formal criteria and indicators for monitoring forest management practices in Venezuelan managed forests that has been demonstrated previously [29], our results add new evidence about the inadequate planning of logging activities, which may further limit the potential of forest management to serve as a climate change mitigation tool.

Conclusions
Compared to other regions of the Amazon basin and the tropics in general, our research reveals that in the northeastern Venezuelan Amazon, while the overall harvesting intensity has remained low over long periods of time, the disturbances associated to logging were considerably high. Selective logging activities showed rather poor planning and low efficiency, reflected in the fact that, regardless of the metric used (area-based, volume-based or CIF), carbon emissions were higher than most studies focusing on similar questions.
Digital image processing and GIS techniques used within the TerraAmazon system, and in conjunction with the 2019 Refinement IPCC Guidelines, enabled us to develop, for the first time in Venezuela, a low-cost and robust analytical approach to study the relationships between selective logging, forest degradation and carbon emissions. Our work also reveals that the use of spectral contrast along with Landsat time series, and ground-based data are excellent tools for the analysis of forest degradation, the evaluation of the impact at the canopy level due to the different activities and modalities of selective logging, and for the estimation of carbon emissions. This study is a step forward to improve and plan other more detailed analytical techniques (e.g., Lidar-Light Detection and Ranging) and for establishing a baseline of carbon emissions in the context of sustainable forest management. Acknowledgments: The authors would like to thank to the coordinators of the project "Sust Forest Lands Management and Conservation Under an Eco-Social Approach (GCP/VEN/011/ funded by the Global Environment Facility (GEF), for logistical support in field trips; and EN for sharing information from management plans and data from forest inventory plots. W thank Jose Ignacio Azuaje from Ministerio del Poder Popular para el Ecosocialismo for his s providing data and information about the management process for Imataca Forest Reserve.

Conflicts of Interest:
The authors declare no conflict of interest. Appendix A Figure A1. Location of field data.