Monitoring Approach for Tropical Coniferous Forest Degradation Using Remote Sensing and Field Data

: Current estimates of CO 2 emissions from forest degradation are generally based on insu ﬃ cient information and are characterized by high uncertainty, while a global deﬁnition of ‘forest degradation’ is currently being discussed in the scientiﬁc arena. This study proposes an automated approach to monitor degradation using a Landsat time series. The methodology was developed using the Google Earth Engine (GEE) and applied in a pine forest area of the Dominican Republic. Land cover change mapping was conducted using the random forest (RF) algorithm and resulted in a cumulative overall accuracy of 92.8%. Forest degradation was mapped with a 70.7% user accuracy and a 91.3% producer accuracy. Estimates of the degraded area had a margin of error of 10.8%. A number of 344 Landsat collections, corresponding to the period from 1990 to 2018, were used in the analysis. Additionally, 51 sample plots from a forest inventory were used. The carbon stocks and emissions from forest degradation were estimated using the RF algorithm with an R 2 of 0.78. GEE proved to be an appropriate tool to monitor the degradation of tropical forests, and the methodology developed herein is a robust, reliable, and replicable tool that could be used to estimate forest degradation and improve monitoring, reporting, and veriﬁcation (MRV) systems under the reducing emissions from deforestation and forest degradation (REDD + ) mechanism.


Introduction
Forest monitoring has been an important scientific objective mainly due to the large number of ecosystem services that serve humanity. One of the most efficient methods to monitor them is through geographically explicit and consistent mapping over time [1,2].
Currently, forest monitoring has been focused on quantifying deforestation; spatial representation and the monitoring of forest degradation are poorly studied, mainly because there is no clear, standardized, and recognized definition of 'forest degradation' globally [3]. Additionally, international distribution of C stocks in the forest, and the amount of CO2 emissions from forest degradation per unit area for a given period of analysis.

Study Area
Our research area is located in the Dominican Republic, mainly in two thirds of the so-called Hispaniola Island in the east, which is the second largest island in the Greater Antilles. The territory of the country covers 48,198 km 2 (18°28′35″ N and 69°53′36″ W) ( Figure 1). The Dominican Republic has a diverse bioclimatic and topographic zones, ranging from dry regions where precipitation reaches 450 mm yr −1 to humid regions where precipitation reaches 2500 mm yr −1 , at altitudes over 3,000 m.a.s.l. The northwest-southeast trending mountain range includes the highest peaks in the Caribbean, Pico Duarte (3098 m.a.s.l) [27]. This wide variety of geographic conditions has given rise to diverse ecosystems and habitats, including arid, semi-arid, humid, and tropical sub-humid zones [28].
The study was carried out in the pine forests of the Dominican Republic, which cover 3287 km 2 . Most of this area lies in the Cordillera Central (the highest elevation mountain range on the island) and comprises four large protected areas that were declared national parks (NPs): (i) NP Armando Bermúdez, (ii) NP José del Carmen Ramírez, (iii) NP Valle Nuevo, and (iv) NP Sierra de Bahoruco. The water service for human consumption and agricultural use in most of the country is one of the main environmental services offered by the Central Cordillera to the Dominican Republic. The vegetation patterns vary mainly due to large-scale climatic factors, such as the direction of the northeast or southeast winds. This mountain range features the principal pine forests of the country. The higher sites of the mountain range include moist broadleaf forests, while the windward side includes forests of West Indies pine (Pinus occidentalis Swartz) [29]. This pine is endemic to the island of Hispaniola (19° N, 71° W), although it has been known to scientists for more than 200 years and still covers extensive areas of the Dominican Republic and Haiti [30].
Recent studies conducted under the REDD+ program by the Ministry of the Environment and Natural Resources (MARN) have indicated that the illegal extraction of wood for firewood and charcoal to be used as fuel, timber, and weak management are the principal drivers of pine forest degradation in the Dominican Republic [31]. The water service for human consumption and agricultural use in most of the country is one of the main environmental services offered by the Central Cordillera to the Dominican Republic. The vegetation patterns vary mainly due to large-scale climatic factors, such as the direction of the northeast or southeast winds. This mountain range features the principal pine forests of the country. The higher sites of the mountain range include moist broadleaf forests, while the windward side includes forests of West Indies pine (Pinus occidentalis Swartz) [29]. This pine is endemic to the island of Hispaniola (19 • N, 71 • W), although it has been known to scientists for more than 200 years and still covers extensive areas of the Dominican Republic and Haiti [30].
Recent studies conducted under the REDD+ program by the Ministry of the Environment and Natural Resources (MARN) have indicated that the illegal extraction of wood for firewood and charcoal Remote Sens. 2020, 12, 2531 4 of 25 to be used as fuel, timber, and weak management are the principal drivers of pine forest degradation in the Dominican Republic [31].

Forest, Deforestation, and Degradation Definitions
In recent times, the definition of 'forest' has taken on particular relevance due to the challenges of countries to monitor the CO 2 emissions from the forest sector as part of the objectives to establish robust MRV systems for REDD+. In general, the definitions of 'forest' include references to threshold parameters that include the minimum area of land, minimum tree height, and minimum canopy cover. Many countries are aligned with the minimum thresholds described by the United Nations Food and Agriculture Organization's (FAO) Global Forest Resource Assessment (FRA). Through the FRA, since 2000, all countries have aligned themselves to adopt a definition of 'forest' with common parameters, such as (i) a canopy coverage of more than 10%, (ii) trees of 5 m, and (iii) land of at least 0.5 ha [32]. The present study subscribes to the definition of 'forest' adopted by the Dominican Republic, with a focus on pine forests in accordance with the Reference Emission Levels/Forest Reference Levels (FREL/FRL): "land of at least 0.5 ha covered by pine trees higher than 5 m and with a canopy cover of more than 30%, or by trees able to reach these thresholds, and predominantly under forest land use, this excludes land that is mainly under agricultural or urban land uses" [33].
Based on the definition of forest described above, in the current study, forest degradation is defined as "the loss of carbon content in forest lands that remain as forest lands with a decrease in canopy cover that does not qualify as deforestation and that can be caused by anthropogenic activities". Forest degradation has a human-induced negative impact on carbon stock changes; our operational definition for measuring forest degradation is based on indicators, such as forest structure (changes in canopy cover) [34] that affect the ability of the forest to store carbon under natural conditions [35]. These definitions demonstrate that deforestation and forest degradation involve different conditions, processes, and concepts. Deforestation suggests a change in land use from forest to non-forest land use, altering the original structure and environment of the forest, while degradation occurs in forest lands that are maintained as forest lands but suffer losses in their forest ecosystem functions ( Figure 2

Forest, Deforestation, and Degradation Definitions
In recent times, the definition of 'forest' has taken on particular relevance due to the challenges of countries to monitor the CO2 emissions from the forest sector as part of the objectives to establish robust MRV systems for REDD+. In general, the definitions of 'forest' include references to threshold parameters that include the minimum area of land, minimum tree height, and minimum canopy cover. Many countries are aligned with the minimum thresholds described by the United Nations Food and Agriculture Organization's (FAO) Global Forest Resource Assessment (FRA). Through the FRA, since 2000, all countries have aligned themselves to adopt a definition of 'forest' with common parameters, such as (i) a canopy coverage of more than 10%, (ii) trees of 5 m, and (iii) land of at least 0.5 ha [32]. The present study subscribes to the definition of 'forest' adopted by the Dominican Republic, with a focus on pine forests in accordance with the Reference Emission Levels/Forest Reference Levels (FREL/FRL): "land of at least 0.5 ha covered by pine trees higher than 5 m and with a canopy cover of more than 30%, or by trees able to reach these thresholds, and predominantly under forest land use, this excludes land that is mainly under agricultural or urban land uses" [33].
Based on the definition of forest described above, in the current study, forest degradation is defined as "the loss of carbon content in forest lands that remain as forest lands with a decrease in canopy cover that does not qualify as deforestation and that can be caused by anthropogenic activities". Forest degradation has a human-induced negative impact on carbon stock changes; our operational definition for measuring forest degradation is based on indicators, such as forest structure (changes in canopy cover) [34] that affect the ability of the forest to store carbon under natural conditions [35]. These definitions demonstrate that deforestation and forest degradation involve different conditions, processes, and concepts. Deforestation suggests a change in land use from forest to non-forest land use, altering the original structure and environment of the forest, while degradation occurs in forest lands that are maintained as forest lands but suffer losses in their forest ecosystem functions ( Figure 2).

Method
Once the baseline data were collected and the key concepts (forest, deforestation, and forest degradation) were defined, a method to quantify the degradation of pine forests in the Dominican Republic was developed. The overall structure of the method (Figure 3) consists of four stages: (i) the preprocessing and selection of Landsat images, (ii) the computation of the spectral indices to map land cover for the 1990-2018 period, (iii) changing the magnitude of mapping, and (iv) mapping the carbon stocks in pine forests. The entire process was accompanied by an accuracy analysis for each

Method
Once the baseline data were collected and the key concepts (forest, deforestation, and forest degradation) were defined, a method to quantify the degradation of pine forests in the Dominican Republic was developed. The overall structure of the method (Figure 3) consists of four stages: (i) the Remote Sens. 2020, 12, 2531 5 of 25 preprocessing and selection of Landsat images, (ii) the computation of the spectral indices to map land cover for the 1990-2018 period, (iii) changing the magnitude of mapping, and (iv) mapping the carbon stocks in pine forests. The entire process was accompanied by an accuracy analysis for each step in which a classification, probability, and regression model was applied using the Smile random forest (RF) algorithm. The GEE cloud-based platform was also used in our research.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 29 step in which a classification, probability, and regression model was applied using the Smile random forest (RF) algorithm. The GEE cloud-based platform was also used in our research. Our models were developed using the Landsat Thematic Mapper (TM) and an operational land imager (OLI) by applying the RF algorithm in GEE to generate a dynamic land cover change map, degraded forest map, and carbon forest map. The model was trained and validated using sample plots from the forest inventory and satellite images.

Landsat TM/OLI Data Processing
We used Landsat-5 TM and Landsat 8 OLI surface reflectance data with 16 day and 30 m resolutions (available in the GEE computing platform) [21]. All Landsat-5TM surface reflectance data from year 1990 ± 0.5 (a total of 22 images) and Landsat-8 OLI surface reflectance data from year 2018 ± 0.5 (a total of 322 images) available in GEE were used in this study.
The Landsat surface reflectance data in GEE were atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (TM) and Landsat Surface Reflectance Corrected (LaSRC) (OLI) algorithms [36,37]. The CFmask algorithm was used to mask the clouds and cloud shadows [38,39]. Landsat 5 TM Top of Atmosphere (TOA) and OLI TOA collections were also used [40].
Digital elevation data were obtained from the Shuttle Radar Topography Mission (SRTM) [41] in GEE. These data have a 30 m spatial resolution. SRTM data were used to calculate the topographic slope and elevation. In addition, the empirical Earth rotation model (ERM) was used as a basis to apply a terrain illumination correction algorithm [42], which allowed us to topographically correct each image. For reflectance images, we used the medoid method [43] (Figure 4).
Once the images were preprocessed, a composite mosaic was developed. This mosaic was formed by combining spatially overlapping images into a single image based on a function of multiple spectral and temporal aggregation ranges [44]. This mosaic (multiband and multidate) was built with the images of the years 1990 and 2018. Our models were developed using the Landsat Thematic Mapper (TM) and an operational land imager (OLI) by applying the RF algorithm in GEE to generate a dynamic land cover change map, degraded forest map, and carbon forest map. The model was trained and validated using sample plots from the forest inventory and satellite images.

Landsat TM/OLI Data Processing
We used Landsat-5 TM and Landsat 8 OLI surface reflectance data with 16 day and 30 m resolutions (available in the GEE computing platform) [21]. All Landsat-5TM surface reflectance data from year 1990 ± 0.5 (a total of 22 images) and Landsat-8 OLI surface reflectance data from year 2018 ± 0.5 (a total of 322 images) available in GEE were used in this study.
The Landsat surface reflectance data in GEE were atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (TM) and Landsat Surface Reflectance Corrected (LaSRC) (OLI) algorithms [36,37]. The CFmask algorithm was used to mask the clouds and cloud shadows [38,39]. Landsat 5 TM Top of Atmosphere (TOA) and OLI TOA collections were also used [40].
Digital elevation data were obtained from the Shuttle Radar Topography Mission (SRTM) [41] in GEE. These data have a 30 m spatial resolution. SRTM data were used to calculate the topographic slope and elevation. In addition, the empirical Earth rotation model (ERM) was used as a basis to apply a terrain illumination correction algorithm [42], which allowed us to topographically correct each image. For reflectance images, we used the medoid method [43] (Figure 4).
Once the images were preprocessed, a composite mosaic was developed. This mosaic was formed by combining spatially overlapping images into a single image based on a function of multiple spectral and temporal aggregation ranges [44]. This mosaic (multiband and multidate) was built with the images of the years 1990 and 2018.

Field Inventory Data
The reference field data included in this study are based on the National Forest Inventory (NFI) collected by the Ministry of the Environment and Natural Resources (MARN) of the Dominican Republic with the support of the REDD/CCAD-GIZ program and the World Bank's Forest Carbon Partnership Facility (FCPF) (https://www.forestcarbonpartnership.org/country/dominican-republic (ERPD document, September, 2019)) [45]. In 2012, the Dominican Republic designed its MRV strategy. This strategy proposed two major lines of monitoring forest resources: (i) satellite monitoring and (ii) terrestrial monitoring. For terrestrial monitoring, the country executed an NFI between 2017 and 2018 with a plan to develop permanent sampling plots to be measured every 5 years according to the action plan of the country's MRV System. The NFI of the Dominican Republic contains 404 sampling units located in the different forest classes, such as moist broadleaf forests (204 plots), subdivided into semi-humid broadleaf forests (117 plots), humid broadleaf forests (76 plots), and broadleaf cloud forest (11 plots); pine forests (59 plots), subdivided into high canopy cover density (19 plots) and low canopy cover density (40 plots); dry forests (71 plots); and mangrove forests (70 plots).
The plot is rectangular with a size of 0.125 hectares (ha) (25 m × 50 m). Different forest characteristics and topographical factors were measured at all plots (tree species, height, diameter at breast height (DBH), soil organic matter, number of trees, geographical coordinates, elevation, and slope). To design the NFI, the methodology proposed by the REDD/CCAD-GIZ program was used [46,47]. For our study, we used the data from 51 plots located in the pine forest areas. More information about the methodology and the results of the NFI is available in the FREL/FRL submission of the Dominican Republic [33,48] to REDD+ UNFCCC ( Figure 5) (Appendix D Table A2 and Appendix E Figure A1).

Field Inventory Data
The reference field data included in this study are based on the National Forest Inventory (NFI) collected by the Ministry of the Environment and Natural Resources (MARN) of the Dominican Republic with the support of the REDD/CCAD-GIZ program and the World Bank's Forest Carbon Partnership Facility (FCPF) (https://www.forestcarbonpartnership.org/country/dominican-republic (ERPD document, September, 2019)) [45]. In 2012, the Dominican Republic designed its MRV strategy. This strategy proposed two major lines of monitoring forest resources: (i) satellite monitoring and (ii) terrestrial monitoring. For terrestrial monitoring, the country executed an NFI between 2017 and 2018 with a plan to develop permanent sampling plots to be measured every 5 years according to the action plan of the country's MRV System. The NFI of the Dominican Republic contains 404 sampling units located in the different forest classes, such as moist broadleaf forests (204 plots), subdivided into semi-humid broadleaf forests (117 plots), humid broadleaf forests (76 plots), and broadleaf cloud forest (11 plots); pine forests (59 plots), subdivided into high canopy cover density (19 plots) and low canopy cover density (40 plots); dry forests (71 plots); and mangrove forests (70 plots).
The plot is rectangular with a size of 0.125 hectares (ha) (25 m × 50 m). Different forest characteristics and topographical factors were measured at all plots (tree species, height, diameter at breast height (DBH), soil organic matter, number of trees, geographical coordinates, elevation, and slope). To design the NFI, the methodology proposed by the REDD/CCAD-GIZ program was used [46,47]. For our study, we used the data from 51 plots located in the pine forest areas. More information about the methodology and the results of the NFI is available in the FREL/FRL submission of the Dominican Republic [33,48] to REDD+ UNFCCC ( Figure 5) (Appendix D Table A2 and Appendix E Figure A1).

Classification of Dynamic Land Cover Change
One of the most relevant tasks for RS is land cover mapping. Different spectral indices are used to improve such mapping techniques [49]. For land cover change mapping, we generated a composite mosaic for the 1990-2018 period, and spectral indices were selected based on the known characteristics of land cover classes. Landsat time series and spectral analyses were used to detect deforestation and degradation forests.

Classification of Dynamic Land Cover Change
One of the most relevant tasks for RS is land cover mapping. Different spectral indices are used to improve such mapping techniques [49]. For land cover change mapping, we generated a composite mosaic for the 1990-2018 period, and spectral indices were selected based on the known characteristics of land cover classes. Landsat time series and spectral analyses were used to detect deforestation and degradation forests.
Classes that were determined in the dynamic land cover map correspond to a stable forest, stable non-forest, degradation, deforestation, and restored forest. Seven different vegetation indices were used to monitor the dynamics of forest change during the 1990-2018 period. Among the most relevant indices used are the enhanced vegetation index (EVI), which was designed to enhance the vegetation signal with improved sensitivity in high biomass areas [50]; the soil adjust vegetation index (SAVI) [51]; and the normalized difference fraction index (NDFI), which was constructed to highlight degraded or cleared forest areas. The NDFI values in intact forests are expected to be high (i.e., approximately 1) due to the combination of high "green vegetation" (GV) and low nonphotosynthetic vegetation (NPV) and soil values [52]. The spectral indices used are closely related to the land cover defined in our research; Appendix C, Table A1 details each index used.

Land Cover Change Samples
Training samples for land cover change classification were derived from a visual analysis using Landsat from GEE and high spatial resolution images from Google Earth (GE) ( Figure 6).
First, we established four cover change classes: stable forest, stable non-forest, deforestation, and forest restoration. We allocated 97 samples to areas that appeared to be stable forest or stable nonforest, and the remaining 15 samples were allocated to land cover change. The second step was to detect forest degradation based on the stable forest class, commonly called the 'forest mask', following the criteria established in the definition of a forest in our study. We focused on two classes: non-degraded forest and degraded forest. A total of 90 samples were allocated to areas with disturbances observed in the 1990-2018 period. Reference samples were randomly distributed over each land cover class in the study area, while a single pixel was used as the sample unit ( Figure 7).
The training dataset was used to improve the supervised classification, the per-band pixel values of the stacked composite images were extracted from the training samples, and the resulting data were used to train the RF classifiers [53]. We used the RF algorithm, because it is a built-in classifier Classes that were determined in the dynamic land cover map correspond to a stable forest, stable non-forest, degradation, deforestation, and restored forest. Seven different vegetation indices were used to monitor the dynamics of forest change during the 1990-2018 period. Among the most relevant indices used are the enhanced vegetation index (EVI), which was designed to enhance the vegetation signal with improved sensitivity in high biomass areas [50]; the soil adjust vegetation index (SAVI) [51]; and the normalized difference fraction index (NDFI), which was constructed to highlight degraded or cleared forest areas. The NDFI values in intact forests are expected to be high (i.e., approximately 1) due to the combination of high "green vegetation" (GV) and low non-photosynthetic vegetation (NPV) and soil values [52]. The spectral indices used are closely related to the land cover defined in our research; Appendix C, Table A1 details each index used.

Land Cover Change Samples
Training samples for land cover change classification were derived from a visual analysis using Landsat from GEE and high spatial resolution images from Google Earth (GE) ( Figure 6).
First, we established four cover change classes: stable forest, stable non-forest, deforestation, and forest restoration. We allocated 97 samples to areas that appeared to be stable forest or stable non-forest, and the remaining 15 samples were allocated to land cover change. The second step was to detect forest degradation based on the stable forest class, commonly called the 'forest mask', following the criteria established in the definition of a forest in our study. We focused on two classes: non-degraded forest and degraded forest. A total of 90 samples were allocated to areas with disturbances observed in the 1990-2018 period. Reference samples were randomly distributed over each land cover class in the study area, while a single pixel was used as the sample unit ( Figure 7).
The training dataset was used to improve the supervised classification, the per-band pixel values of the stacked composite images were extracted from the training samples, and the resulting data were used to train the RF classifiers [53]. We used the RF algorithm, because it is a built-in classifier in GEE and has been widely demonstrated to improve the accuracy of maps by combining random subsets of trees to classify the training samples. In GEE, the RF algorithm is applied through the following function: (ee.Classifier.smileRamdomForest). Moreover, the algorithm can be configured in three ways (ee.setOutputMode) based on the classification mode (class/type maps), regression mode (maps with continuous values predicted), and the probability mode (map with rescaled values between 0 and 1). In the current study, the RF algorithm was applied in the classification mode using GEE to obtain the land cover, in the regression mode to estimate the carbon maps, and in the probability mode to estimate the change magnitude maps.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 29 in GEE and has been widely demonstrated to improve the accuracy of maps by combining random subsets of trees to classify the training samples. In GEE, the RF algorithm is applied through the following function: (ee.Classifier.smileRamdomForest). Moreover, the algorithm can be configured in three ways (ee.setOutputMode) based on the classification mode (class/type maps), regression mode (maps with continuous values predicted), and the probability mode (map with rescaled values between 0 and 1). In the current study, the RF algorithm was applied in the classification mode using GEE to obtain the land cover, in the regression mode to estimate the carbon maps, and in the probability mode to estimate the change magnitude maps.

Carbon Stock and Change Magnitude
Estimating the spatial-temporal distributions of forest carbon stocks subject to land cover changes is critical for estimating and reporting GHG emissions [54]. To spatially represent explicit in GEE and has been widely demonstrated to improve the accuracy of maps by combining random subsets of trees to classify the training samples. In GEE, the RF algorithm is applied through the following function: (ee.Classifier.smileRamdomForest). Moreover, the algorithm can be configured in three ways (ee.setOutputMode) based on the classification mode (class/type maps), regression mode (maps with continuous values predicted), and the probability mode (map with rescaled values between 0 and 1). In the current study, the RF algorithm was applied in the classification mode using GEE to obtain the land cover, in the regression mode to estimate the carbon maps, and in the probability mode to estimate the change magnitude maps.

Carbon Stock and Change Magnitude
Estimating the spatial-temporal distributions of forest carbon stocks subject to land cover changes is critical for estimating and reporting GHG emissions [54]. To spatially represent explicit

Carbon Stock and Change Magnitude
Estimating the spatial-temporal distributions of forest carbon stocks subject to land cover changes is critical for estimating and reporting GHG emissions [54]. To spatially represent explicit forest degradation along with degraded carbon, we generated a carbon map of pine forests using 51 sampled plots from the NFI data and Landsat 2018. The size of each parcel is 0.125 hectares (ha) (25 m × 50 m). As the size of each plot and the Landsat pixels do not match, we used the geographic coordinates of the centroid of the plot. Thus, we applied the carbon values in units per hectare (t ha −1 ) to each pixel.
For each plot, four of the five C pools were measured [55]: aboveground biomass (AGB), belowground biomass (BGB), deadwood (DW), and leaf litter (Table 1) [47,56,57] (See Section 2.4.2 Field Inventory Data for more details). To convert the biomass to carbon, the IPCC default carbon factor value (0.47) was used. Each pool was modeled independently using spectral responses (vegetation indices) from Landsat, applying a regression model with the RF algorithm in GEE (ee.Classifier.randomForest (30).setOutputMode('REGRESSION'). Through this process, each spatial pixel acquired an AGB (t C ha −1 ) value and its related standard deviation as a measure of uncertainty.
The code generated to model the C stock using the GEE is available in Appendix B.
The change magnitude is assumed to be an approximate indicator of the amount of tree removal or canopy damage that occurred due to disturbances [58]. The change magnitude was estimated from the spectral indices for the stable forest class using the composite mosaic for the 1990-2018 period. We fitted an RF probability model in GEE to represent the structural forest changes in each of the spectral indices described previously (ee.Classifier.randomForest (50).setOutputMode('PROBABILITY').
The disturbance monitoring algorithm used to identify the forest changes in each pixel location is closely related to the Continuous Change Detection and Classification algorithm (CCDC) [20,58,59] but was adapted using RF models to predict the change magnitude probabilities for bitemporal observations. Several studies have successfully applied this algorithm to different sensors and different spectral indices to detect changes [60,61]. This change detection algorithm operates on the time series of each pixel in the study area.
The decrease in a spectral index caused by a disturbance is recognized by a certain change magnitude. For example, values close to 1 in the NDFI spectral indices indicate high proportions of green vegetation (GV) or stable forests, while values close to 0 in the NDFI imply higher proportions of soil (So). The code generated to estimate the change magnitude using the GEE is available in Appendix B.
Pixel-based mapping facilitates comparisons and evaluations of changes with direct algebraic calculations [62]. Once the carbon stored in the four pools of the forest was estimated for the year 2018 along with the magnitude of change between 1990-2018, we used (Equation (1)) to determine the C stocks in each pixel for the year 1990. Finally, the degraded carbon was estimated as the difference of the C stored in the pine forests between 1990 and 2018 combined with the change magnitude observed in the same period (Equation (2)). To estimate forest degradation, only disturbances occurring in the stable forest class were considered for the period analyzed: where C t1 is the C stock in 1990 (Mg ha −1 ), C t2 is the C stock in 2018, and CM is the change magnitude (value > 0 and < 1).
where CD is the carbon degraded for the 1990-2018 period (Mg ha −1 ), C t1 is the C stock in 1990, and C t2 is the C stock in 2018. Finally, to calculate the annual rate of degraded carbon, the stock-difference method (SDM) was used, where changes in carbon stock (∆Carbon) represent the difference between carbon stocks for a given forest area estimated at two time points (Equation (3)): where ∆Carbon is the annual change in C stocks (Mg·C·ha −1 ·yr −1 ), Carbon t1 is the C stock in 1990 (Mg·C·ha −1 ), and Carbon t2 is the C stock in 2018 (Mg C·ha −1 ). Non-woody biomass is recorded, which includes dead leaves (dead biomass) and herbaceous vegetation (living non-woody biomass on the ground). The maximum diameter for woody material to be considered is 2 cm.

Model Evaluation: Carbon Stock
Cross-validation (CV) is one of the most commonly used techniques to evaluate the efficiency of a machine learning (ML) technique; this is due to its wide application in the scientific arena and its efficiency in detecting a model's overfitting problems [63]. To evaluate the performance of the machine learning model applied to map the change magnitude and forest carbon, the following functions were used: coefficient of determination (R 2 ), mean square error (MSE), root mean square error (RMSE), mean absolute deviation (MAD), cumulated forecast error (CFE), and mean absolute percentage error (MAPE): where n (i = 1, 2, . . . , n) is the number of samples used for the machine learning model, y i is the value observed (C stock), y i is the corresponding mean value, f i is the predicted value (C stock), and f i is the mean value.

Accuracy Assessment and Analysis
We used the confusion matrix statistical accuracy assessment method to evaluate the dynamic land cover change classification. The overall accuracy (OA), user's accuracy (UA), and producer's accuracy (PA) were applied to each class, and the Kappa coefficient was used to assess the class map and determine the level of agreement between two raters. The standard deviations and the confidence intervals (at a 95% significance level) were also estimated.
Since land change classes (degradation, deforestation, and forest restoration) tend to cover only a small portion of the study objectives compared to stable areas (stable forest and stable non-forest), it is recommended to stratify the study based on a map that represents the classes of principal interest to ensure an effective statistical sample representation in land change classes, such as degradation, deforestation, and forest restoration [64].
An accuracy assessment of the dynamic land cover change map (for the 1990-2018 period), generated through a sampling-based approach to estimate the area of forest degradation in the Dominican Republic (Figure 6), was performed on the following land cover change classes: Stable forest: pine forests that remain pine forests without disturbance; this forest contains over 30% canopy cover.
Stable non-forest: other non-forest lands, such as agriculture, wetlands, grasslands. Deforestation: elimination of the forest canopy cover that exceeds 30%; results in a land-use change.
Degradation: this entails any disturbance that changes the canopy cover density between 100% and 30% and does not result in a land-use change.
Forest restoration: conversion of non-forested land to forest; this includes forest restoration with a canopy cover greater than 30% (through natural and artificial means) on deforested land.
For the accuracy assessment, a total of 1124 spatial sampling points ( Table 2) were established for the study area using the stratified random sampling approach following best practices [64] (Equation (10)). It was necessary to modify the minimum sample size to determine the objective standard error of the degradation area, rather than the OA of the map, and thereby ensure that sample size would be large enough to produce sufficiently accurate estimates [65]. The stratified area estimator-design tools hosted in the System for Earth Observation Data Access, Processing and Analysis for Land (SEPAL) were used to generate random spatial points ( Table 2). SEPAL is a cloud-based computing platform developed by FAO, which uses the GEE and OpenForis Geospatial Toolkits [66].
where n = number of points in the study area, S Ô is the standard error of the estimated OA, W i is the mapped proportion of the class area i, and S i is the standard deviation of land cover classes i. We performed an analysis following best practices to assess the accuracy of the map classification, and the area of change was estimated using a classification error matrix. For details on the matrix nomenclature, refer to Olofsson et al. (2014) [64]. Collecting reference observations of forest degradation is a complex task, primarily because degradation is a continuous process that must be observed over a long period. In this sense, satellite images with high spectral resolutions have become a key tool, but they are not sufficient to reconstruct the landscape's historical dynamics. Therefore, this study required the use of a Landsat observation time series supported by very high spatial resolution (VHRS) imagery. Independent stratified validation samples were visually interpreted from the VHRS time-series images of Collect Earth (CE). We built a survey in CE that helped us access multiple satellite images, including archives including VHSR imagery (Google Earth, Bing Maps) and a set of satellite images from the GEE catalog, along with their derived spectral indices [67] (Figure 8).
To facilitate the historical collection of reference data, other GEE assessment tools were adapted and used, such as the Accuracy and Area Estimation Toolbox (AREA2) developed by Bullock and Olofsson (2018) (see github.com/bullocke/AREA2). The sample interpretation tool allowed us to determine reference labels for the 1124 samples collected. This algorithm helped estimate the map accuracy and disturbance area and visualize the time series trends of each sample using a dataset created using the Continuous Degradation Detection (CODED) methodology in GEE [20] (see Appendix B for the code developed in this study using the GEE).
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 29 stratified validation samples were visually interpreted from the VHRS time-series images of Collect Earth (CE). We built a survey in CE that helped us access multiple satellite images, including archives including VHSR imagery (Google Earth, Bing Maps) and a set of satellite images from the GEE catalog, along with their derived spectral indices [67] (Figure 8).
To facilitate the historical collection of reference data, other GEE assessment tools were adapted and used, such as the Accuracy and Area Estimation Toolbox (AREA2) developed by Bullock and Olofsson (2018) (see github.com/bullocke/AREA2). The sample interpretation tool allowed us to determine reference labels for the 1,124 samples collected. This algorithm helped estimate the map accuracy and disturbance area and visualize the time series trends of each sample using a dataset created using the Continuous Degradation Detection (CODED) methodology in GEE [20] (see Appendix B for the code developed in this study using the GEE). To review and assign reference labels to each of the 1124 selected special sample units, three trained interpreters were delegated. These interpreters did not know the classes of the assigned samples. An additional interpreter reviewed all samples with low or medium confidence. At least three interpreters reviewed all units labeled as degradation. The decision matrix and labels assigned to the samples evaluated in the time series are shown in Table 3.  To review and assign reference labels to each of the 1124 selected special sample units, three trained interpreters were delegated. These interpreters did not know the classes of the assigned samples. An additional interpreter reviewed all samples with low or medium confidence. At least three interpreters reviewed all units labeled as degradation. The decision matrix and labels assigned to the samples evaluated in the time series are shown in Table 3.

Dynamic Land Cover Changes from 1990 to 2018
The total study area corresponded to 328,777 ha of pine forests in the Dominican Republic. The results showed that degraded forests accounted for 11% ± 1.21% (95% confidence interval) of the total study area between 1990 and 2018, while 79% ± 1.28% remained stable and did not suffer any disturbances; further, 2% ± 0.61% were deforested. In total, we estimated that 36,808 ± 446 ha of pine forests was degraded.
The margin of error of the area estimate for forest degradation was 10.8% (95% CI). The user's accuracy was 70.7%, and the producer's accuracy for forest degradation was 91.3%. The overall accuracy of the dynamic land cover change map was 92.8%. The main results corresponding to the accuracy assessment are shown in Table 4. We analyzed the dynamics of land cover change, as they relate to the country's protected areas, and identified that 71% of the degraded pine forests exist within the protected area. Among these, the main forest belongs to Sierra de Bahoruco National Park (NP), with 14,166 ha (30% of the total degraded area), followed by Valle Nuevo NP and José del Carmen Ramírez NP, with 8,736 ha (18% of the total degraded area) and 6,462 ha (16% of the total degraded area), respectively. Table 5 shows the locations of the protected areas with the degraded pine forests from 1990 to 2018. A geographic representation of the dynamic land cover change map obtained in our study is also provided (Figure 9). Detailed results on land cover change are available via a dashboard called "Accuracy assessment and analysis tools" available in Appendix A.

Carbon Stock
The total carbon stock in the pine forest area analyzed was composed of AGB C, BGB C, DW C, and litter C pools. The results for the total carbon analysis are presented in

Carbon Stock
The total carbon stock in the pine forest area analyzed was composed of AGB C, BGB C, DW C, and litter C pools. The results for the total carbon analysis are presented in Table 6,  (Table 7). Detailed results on the carbon stored from the different pools and protected area categories are available via a dashboard provided in Appendix A. C) was stored in natural reserves (Table 7). Detailed results on the carbon stored from the different pools and protected area categories are available via a dashboard provided in Appendix A.

Carbon Degraded in the 1990-2018 Period
The total carbon degraded in the pine forest area analyzed was 3,479,159 Mg C. Converting this degraded carbon into emission-and climate-related units of CO2-equivalent emissions (metric tons CO2-equivalent units), the emissions caused by the degradation of the pine forests in the period 1990-2018 were 12,756,916 tCO2eq, with an annual average of 2.6 Mg C ha −1 yr −1 (9.5 tCO2eq ha −1 yr −1 ). Of the total degraded C stock, 73.9% (2,570,081 Mg C) was found in national parks, while 2.9% (102,401

Carbon Degraded in the 1990-2018 Period
The total carbon degraded in the pine forest area analyzed was 3,479,159 Mg C. Converting this degraded carbon into emission-and climate-related units of CO 2 -equivalent emissions (metric tons CO 2 -equivalent units), the emissions caused by the degradation of the pine forests in the period 1990-2018 were 12,756,916 tCO 2 eq, with an annual average of 2.6 Mg C ha −1 yr −1 (9.5 tCO 2 eq ha −1 yr −1 ). Of the total degraded C stock, 73.9% (2,570,081 Mg C) was found in national parks, while 2.9% (102,401 Mg C) and 22.8% (792,048 Mg C) of C were degraded in natural reserves and non-protected areas, respectively (Table 8). Detailed results on the degraded carbon in pine forests in the 1990-2018 period for the different pools and protected area categories are provided in Appendix A.  (Table 8). Detailed results on the degraded carbon in pine forests in the 1990-2018 period for the different pools and protected area categories are provided in Appendix A.

Validation of Dynamic Land Cover Change Map from the 1990-2018 Period
Most of the countries that are part of the REDD+ mechanism do not quantify or report their emissions caused by forest degradation [16]. Efforts to find such a method have been great, and the challenge of obtaining accurate estimates remains under investigation and debate in the scientific arena. One of the main agreements in the measurement and monitoring of forest degradation is that a long series of temporal and spatiotemporal observations are required to detect disturbances in forest cover, which is why satellite images, such as those from Landsat, are key inputs to establishing more robust MRV systems for REDD+.
Using Landsat data in GEE, we developed a methodology to monitor pine forest degradation in the tropics. This approach was determined to be precise, with an overall 92.9% and 91% producer's

Validation of Dynamic Land Cover Change Map from the 1990-2018 Period
Most of the countries that are part of the REDD+ mechanism do not quantify or report their emissions caused by forest degradation [16]. Efforts to find such a method have been great, and the challenge of obtaining accurate estimates remains under investigation and debate in the scientific arena. One of the main agreements in the measurement and monitoring of forest degradation is that a long series of temporal and spatiotemporal observations are required to detect disturbances in forest cover, which is why satellite images, such as those from Landsat, are key inputs to establishing more robust MRV systems for REDD+.
Using Landsat data in GEE, we developed a methodology to monitor pine forest degradation in the tropics. This approach was determined to be precise, with an overall 92.9% and 91% producer's accuracy in the degraded forest class of the dynamic land use change map of pine forests in the Dominican Republic. Our estimates of forest degradation are compatible with those of other studies on a sub-national scale. For example, the OA obtained in degradation and deforestation mapping in Rondônia, Brazil, was 91%, while the producer's accuracy reached 68% in the forest degradation class [58]. In the forests of the Brazilian Amazon, the OA in degradation and deforestation mapping was 92%, while the producer's accuracy was 80% in the forest degradation class [68]. Another study using SPOT images with spectral mixing models in the eastern Amazon showed results that also indicated good agreement (86% OA) [69].
The dynamic mapping analysis determined an efficient stratification in the study area and allowed for an impartial estimation. Margins of error of 10.8% were obtained when mapping forest degradation at a 95% confidence level. Although the mapping of forest degradation in tropical forests is scarce in the literature, we observed some consistency between our results and those of other studies in the temporal scale, spatial and spectral resolution of the images used, accuracy, and the use of vegetation index analysis as a method to evaluate and map tropical forest disturbances.
Historical data collection on forest changes is a challenging task because such data are not readily available everywhere, and temporal change data are not detailed enough for the validation of time series maps. The current study used the AREA2 algorithm developed by Bullock and Olofsson (2018) by applying the Time Series Viewer. This is a sample interpretation tool used to determine the reference labels derived from a mapped dataset. However, a new challenge involved assessing the changes detected by the model. For this process, independent datasets were selected and assessed using CE [70]. The combination of these tools proved to be efficient in our study.
Using spectral index measures to validate the dynamic land cover change map in our study allowed us to extend tele-interpretation techniques and facilitated the visual detection of historical change processes, especially for degraded forest detection. In other research, the NDFI was used to map degradation; ultimately, the NDFI was found to be more sensitive to disturbances from tropical forests than other spectral indices [71]. We use different spectral indices to improve our estimates of forest degradation and performed a regression analysis with random forest to determine the importance of the indices in the constructed model. We found that the NDFI and EVI were the main variables able to explain the model and thereby classify areas with forest degradation for the 1990-2018 period (Figure 13a,b).
Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 29 accuracy in the degraded forest class of the dynamic land use change map of pine forests in the Dominican Republic. Our estimates of forest degradation are compatible with those of other studies on a sub-national scale. For example, the OA obtained in degradation and deforestation mapping in Rondônia, Brazil, was 91%, while the producer's accuracy reached 68% in the forest degradation class [58]. In the forests of the Brazilian Amazon, the OA in degradation and deforestation mapping was 92%, while the producer's accuracy was 80% in the forest degradation class [68]. Another study using SPOT images with spectral mixing models in the eastern Amazon showed results that also indicated good agreement (86% OA) [69]. The dynamic mapping analysis determined an efficient stratification in the study area and allowed for an impartial estimation. Margins of error of 10.8% were obtained when mapping forest degradation at a 95% confidence level. Although the mapping of forest degradation in tropical forests is scarce in the literature, we observed some consistency between our results and those of other studies in the temporal scale, spatial and spectral resolution of the images used, accuracy, and the use of vegetation index analysis as a method to evaluate and map tropical forest disturbances.
Historical data collection on forest changes is a challenging task because such data are not readily available everywhere, and temporal change data are not detailed enough for the validation of time series maps. The current study used the AREA2 algorithm developed by Bullock and Olofsson (2018) by applying the Time Series Viewer. This is a sample interpretation tool used to determine the reference labels derived from a mapped dataset. However, a new challenge involved assessing the changes detected by the model. For this process, independent datasets were selected and assessed using CE [70]. The combination of these tools proved to be efficient in our study.
Using spectral index measures to validate the dynamic land cover change map in our study allowed us to extend tele-interpretation techniques and facilitated the visual detection of historical change processes, especially for degraded forest detection. In other research, the NDFI was used to map degradation; ultimately, the NDFI was found to be more sensitive to disturbances from tropical forests than other spectral indices [71]. We use different spectral indices to improve our estimates of forest degradation and performed a regression analysis with random forest to determine the importance of the indices in the constructed model. We found that the NDFI and EVI were the main variables able to explain the model and thereby classify areas with forest degradation for the 1990-2018 period (Figure 13a,b). While the Dominican Republic has experienced a decrease in deforestation in recent times, forest degradation has been on the rise. In our study, we found that an area equivalent to 14% of pine forests of the Dominican Republic was degraded between 1990 and 2018. This is a critical element that must be considered for the development of a REDD+ program in the country. The Dominican Republic established a FREL/FRL to obtain results-based payments under the REDD+ mechanism supported While the Dominican Republic has experienced a decrease in deforestation in recent times, forest degradation has been on the rise. In our study, we found that an area equivalent to 14% of pine forests of the Dominican Republic was degraded between 1990 and 2018. This is a critical element that must be considered for the development of a REDD+ program in the country. The Dominican Republic established a FREL/FRL to obtain results-based payments under the REDD+ mechanism supported by the FCPF. This includes emissions and removal in the remaining forest land (emissions from forest degradation) for the 2006-2015 period [45]. In this sense, the methodology proposed and the results obtained in the current study contribute directly to monitoring and quantifying CO 2, emissions and removal. Furthermore, this methodology is a valuable tool that can be used in other tropical countries to monitor forest degradation.

Carbon Model Assessment
Forest inventory is an important source of information for a variety of strategic purposes in forest management. Based on 51 carbon samples obtained in pine forests, we generated a predictive regression model of the carbon stored in AGB, BGB, DW, and litter. Satisfactory results were obtained when applying the RF algorithm to estimate the total carbon stock in the study area, obtaining an R 2 of 78.1% and an RMSE of 13.38 Mg C ha −1 . It was determined that 19,002,000 Mg C is stored in the pine forests of which 3,479,159 Mg C was degraded (18% of C stock), which is equivalent to 124,255 Mg C yr −1 for the 1990-2018 study period.
We found some differences between our results and previous estimates; for example, the MARN of the Dominican Republic estimated the local emissions from forest degradation at a level of 182,937 Mg C yr −1 , including all forest classes, and at a level of 46,591 Mg C yr −1 specifically for pine forests for the 2006-2015 reference period [33]. The differences between both estimates are mainly due to the scale of monitoring (sub-national vs. national), the methodology, the reference period, the algorithms used for image processing and classification, and especially the differences in the accuracy obtained between both estimates.
These differences also suggest that between 1990 and 2006, the degradation rate in the pine forest was much higher compared to that during the period 2006-2015. This provides a new research opportunity to understand the drivers of degradation that have decreased in the country. However, we believe that the degradation estimates found in both studies differ from each other due to the previously mentioned technical and methodological factors.
A robust and transparent national forest monitoring system to monitor and reporting five REDD+ activities is required [9]. Often, capabilities and national circumstances prevent the monitoring and reporting of CO2e emissions and their removal under the five REDD+ activities. In our study, the forest carbon stock increases in forest lands that remain forest lands were not estimated for the reference period due to the absence of information on the average annual increase in biomass in the studied forest. This offers a new research opportunity for a "step-wise" approach that will allow estimate net CO2e emissions as part of one of the five REDD+ activities called "enhancement of forest carbon stocks".

Google Earth Engine Platform
The processing and analysis of the data from our study were accomplished using the JavaScript API via the code editor of the GEE platform. Using GEE, the processing time efficiency increased, and satisfactory results were obtained (see the GEE code developed for this study in Appendix B). GEE is considered a multidisciplinary tool; since 2017, its application has increased notably, especially in land-use mapping and water resources. Landsat images (82%), the RF algorithm (52%), and the NDVI spectral indices are the most frequently used methods in recent studies related to vegetation [72].
In our study, we faced a series of complexities in the detection of forest degradation, especially since a very detailed approach is required. This approach involves detecting the reduction or modifications in the forest structure that are not considered a total loss over an extended period of time and at a large scale. Using the characteristics of the GEE platform by extracting spectral characteristics from satellite images, it was possible to satisfactorily solve the main challenges we encountered. The methodology proposed here was able to detect the dynamics of change in land use and forest degradation. While there is room to improve the methodology, the use of GEE with its computing capacity and the availability of free satellite images provide powerful support for mapping forest degradation. Additionally, GEE allowed us to quantify the carbon stored until 2018 and the degraded carbon in four pine forest pools. Thus, this method could become a key monitoring tool (at the sub-national and national levels). Moreover, this tool will allow authorities to monitor forest degradation according to indicator 15.3.1 of the Sustainable Development Goals (SDGs), which defines the area of degraded land and will serve as an essential element for developing national MRV systems for REDD+ strategy implementation and supporting the report on GHG emissions of the United Nations Framework Convention on Climate Change (UNFCCC) in an efficient, robust, and transparent manner.

Conclusions
The current study developed and applied a methodology for forest degradation mapping based on available data from the Open Access Landsat and the GEE platforms. Additionally, GEE was used in combination with tools, such as SEPAL and CE by the FAO. The main objective was to estimate the degraded pine forest area and quantify the degraded carbon for a period of 28 years (1990-2018) by applying machine learning models, such as random forest.
The methodology applied in this study shows new possibilities for forest degradation monitoring and estimating CO 2 emissions from forest degradation using spectral information derived from Landsat archives and data from the forest inventory; combining both sources of information can also help improve the MRV systems for the REDD+ mechanism.
The model assessment revealed a dynamic land change map with a cumulative overall accuracy of 92%, in relevant classes (such as forest degradation) with a UA of 70.7%, and a PA of 91.2%. A carbon stock model was also developed with an R 2 of 79% to estimate degradation in terms of the Mg C ha −1 . The applied models were built, trained, and validated to demonstrate the efficiency of the methodology. The results obtained indicate that this methodology can be an especially useful tool for time series processing to map forest degradation by applying technologies, such as GEE.
GEE has excellent potential for the "wall to wall" forest degradation mapping of tropical pine forest ecosystems. More research is still required to assess the ability of GEE to map degradation in broadleaf forests and dry forest ecosystems by applying machine learning techniques combined with spatial data and field measurements. The approaches presented herein could become a key tool for measuring and monitoring emissions from forest degradation in the tropics.