Forest Degradation Assessment Based on Trend Analysis of MODIS-Leaf Area Index: A Case Study in Mexico

: Assessing forest degradation has been a challenging task due to the generally slow-changing nature of the process, which demands long periods of observation and high frequency of records. This research contributes to e ﬀ orts aimed at detecting forest degradation by analyzing the trend component of the time series of Leaf Area Index (LAI) collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) over Central Mexico from 2002 to 2017. The analysis of the trend component is proposed to overcome the challenge of identifying very subtle and gradual changes that can be undetected if only the raw time series is examined. Additionally, the use of LAI as an alternative to other widely used indexes (e.g., Normalize Di ﬀ erence Vegetation Index and Enhanced Vegetation Index) facilitates consideration of the structural changes evident from degradation though not necessarily observable with spectral indices. Overall, results indicate that 52% of the study area has experienced positive trends of vegetation change (i.e., increasing LAI), 37% has remained unchanged, and 11% exhibits some level of forest degradation. Particularly, the algorithm estimated that 0.6% (385 km 2 ) is highly degraded, 5.3% (3406 km 2 ) moderately degraded, and 5.1% (3245 km 2 ) slightly degraded. Most of the moderate and highly degraded areas are distributed over the east side of the study area and evergreen broadleaf appears to be the most a ﬀ ected forest type. Model validation resulted an accuracy of 63%. Some actions to improve this accuracy are suggested, but also a di ﬀ erent approach to validate this type of study is suggested as an area of opportunity for future research.


Introduction
Forest ecosystems play an essential role maintaining the carbon cycle balance because they are able to store carbon [1]. However, when these ecosystems are perturbed, stored carbon is released as carbon dioxide (CO2) into the atmosphere, where it contributes to global warming by absorbing heat and impeding transmission to space. The assessment of forest degradation depends on how "forest" and "forest degradation" are defined. Both concepts, which have been the subject of extensive debates [2,3], have been defined by many institutions based on particular perspectives and purposes. In this research, forests are "lands dominated by woody vegetation with a percent cover of more than 60 percent and height exceeding 2 meters" [4], which is in agreement with the definition used by the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover product (MCD12Q1). Regarding Mexico's National Forest and Soils Inventory (NFSI).

Study Area
Due to the multitemporal nature of this work, the study area encompassed all 500 m pixels contained in the MODIS tile h08v06 that were covered by forest during at least one year between 2002 and 2017 ( Figure 1). According to the MODIS Land Cover product (MCD12Q1, type 3), these pixels comprised 63,964 km 2 , 60% of which corresponded to deciduous broadleaf forest, 36% to evergreen broadleaf, and 4% to evergreen needleleaf. These forests are home to a wide variety of species (in the order of hundreds), some of the most commonly found include Quercus magnoliifolia, Quercus rugosa, Quercus laeta, Quercus sideroxyla, Pinus durangensis, Abies sp., Lysiloma divaricatum, Laguncularia racemosa, Byrsonima crassifolia, Beilschmiedia mexicana, and Guazuma ulmifolia [34].
This vegetation is primarily distributed across the Sierra Madre Occidental (west mountain range) and Sierra Madre Oriental (east mountain range). In the western region, elevation ranges from sea level to 3104 m with a mean of 972 m and a standard deviation of 742 m, while slopes range from 0 to 30° with a mean of 6.9 o and a standard deviation of 4.5°. In the eastern region, the altitudinal range is wider, it goes from sea level to 3300 m with a mean of 1067 m and a standard deviation of 684 m, while slopes range from 0 to 35° with a mean of 6.9° and a standard deviation of 5°.
The climatic conditions range from tropical to semi-cold. Tropical and semi-tropical climates, whose mean annual temperature exceeds 18 °C, extend over most of the areas facing the coastal plains and cover 41% and 35% of the study area, respectively. Temperate climate, whose mean annual temperature oscillates between 12 and 18 °C, extends over higher elevations and covers 13% of the study area. Arid climates are generally found towards continental areas and represent 8% of the study area. The highest elevations are occupied by semi-cold climates, which cover the remaining 3% of the study area. This vegetation is primarily distributed across the Sierra Madre Occidental (west mountain range) and Sierra Madre Oriental (east mountain range). In the western region, elevation ranges from sea level to 3104 m with a mean of 972 m and a standard deviation of 742 m, while slopes range from 0 to 30 • with a mean of 6.9 • and a standard deviation of 4.5 • . In the eastern region, the altitudinal range is wider, it goes from sea level to 3300 m with a mean of 1067 m and a standard deviation of 684 m, while slopes range from 0 to 35 • with a mean of 6.9 • and a standard deviation of 5 • . The climatic conditions range from tropical to semi-cold. Tropical and semi-tropical climates, whose mean annual temperature exceeds 18 • C, extend over most of the areas facing the coastal plains and cover 41% and 35% of the study area, respectively. Temperate climate, whose mean annual temperature oscillates between 12 and 18 • C, extends over higher elevations and covers 13% of the study area. Arid climates are generally found towards continental areas and represent 8% of the study area. The highest elevations are occupied by semi-cold climates, which cover the remaining 3% of the study area.
In Mexico, nearly 60% of the forest lands are owned by ejidos (i.e., a tenurial system in which land is managed collectively) or agrarian communities [35], where forests are often managed for timber production. The remaining forested areas are managed by private property owners and, in a smaller proportion, by the government in the form of natural protected areas [36,37]. In the study area, about 50% of the land is managed by ejidos. Although this proportion is lower than at national level, ejidos are the dominant land tenure units.
Central Mexico was chosen as the study area for three reasons. First, studies focused on detecting forest disturbances or degradation in Mexico have been conducted primarily at local scales [16,[38][39][40][41] and rarely at a regional scale [22] like this study. Second, according to [37], the levels of human disturbance in Mexican forest are significantly high. Third, according to the Land Use and Vegetation Map produced by the Instituto Nacional de Geografía y Estadística (INEGI), 14,900 km 2 (23%) of the study area was covered by perturbed secondary vegetation in 2014 [42].

Data and Preprocessing
LAI, the ratio of half of the total leaf surface area per unit ground area, is a dimensionless measure that ranges from 0 (e.g., bare ground) to over 10 (e.g., dense forest) [43]. In this research, LAI estimates were acquired from the MODIS sensor onboard the Terra satellite. The morning overpass time of Terra (at approximately 10:30 at the equator) provides more reliable LAI estimates than the afternoon overpass time of Aqua (at approximately 13:30) due to lower cloud presence [44]. A visual comparison of the Terra and Aqua cloud fractions was performed by [45], who found a generally cloudier pattern over land in the Aqua products.
MODIS LAI is released in conjunction with Fraction of Photosynthetically Active Radiation (FPAR) in a single product. Both variables are estimated from an algorithm based on a tridimensional (i.e., it accounts for vegetation structure) radiative transfer model that selects the best LAI and FPAR estimates by comparing measured spectral values (from the red and near-infrared regions) with a look-up table that contains observed spectral values of six biome types and their corresponding LAI/FPAR values. This main algorithm may fail to obtain a solution due to large uncertainties in the input spectral values, cloud effects, or too low sun/view zenith angles. In this case, the retrievals, with relatively poor quality, are produced from a back-up algorithm based on the biome-specific empirical relationship between NDVI and LAI/FPAR [25,44,[46][47][48]. According to [44], Collection 6 (C6) of the MODIS LAI/FPAR product is considerably better than Collection 5 (C5). C6 properly captures the interannual variation of LAI, as well as the general seasonality of all biomes, except for evergreen broadleaf forests, where poor quality retrievals are produced. However, this shortcoming is also observed in other satellite products, such as the Carbon cYcle and Change in Land Observational Products from an Ensemble of Satellites (CYCLOPES) and those from the Global Land Surface Satellite (GLASS), and the SPOT-VEGETATION GEOV1 system. MODIS data were acquired through the Land Processes Distributed Active Archive Center (LP DAAC) managed by the NASA Earth Science Data and Information System (ESDIS) project. Based on the characteristics mentioned, MODIS Terra 8-day LAI/FPAR composites (MOD15A2H, C6) of the tile h08v06 were downloaded from 2002 to 2017 at 500 meters spatial resolution. To extract forest land covers, MODIS Terra + Aqua yearly Land Cover composites (MCD12Q1, C6) were also obtained for the same area and time period, and at the same spatial resolution. The preprocessing steps that are described below were conducted by means of ArcPy, a site package that allows for the use of ArcGIS functionalities through the Python programming language.
All datasets were projected to Lambert Conformal Conic projection in ITRF2008 datum and clipped to forested areas, which were delimited by identifying all pixels covered by forest during at least one year between 2002 and 2017. In the few cases for which the type of forest changed between years, the most frequent forest type was assigned. Only good quality LAI retrievals (i.e., obtained from the main algorithm) were selected to perform the analysis. The Land Data Operational Product Evaluation (LDOPE) tool was used to unpack the LAI Quality Control layers, which specify the type of algorithm used to derive the retrievals. Once the good quality values were extracted, monthly LAI means were calculated on a pixel-by-pixel basis for each month from January 2002 to December 2017. Finally, in order to manipulate the data in an R environment, all LAI means were extracted from the raster files and rearranged into a comma-separated values (CVS) file, where each row contained the monthly LAI time series of a single pixel. Due to the exclusion of poor-quality retrievals, 0.5% of the pixels (1505) were discarded from the analysis.

Algorithm to Analyse LAI Trends
A seasonal-trend decomposition procedure based on regression (STR; see [33] for a complete description of the procedure) was used to break down the time series (y t , where t is time) into trend (µ t ), seasonal (γ t ), and remainder ( t ) components. In a given time series (Figure 2a), the simplest STR model considers these components to be additive: The trend component ( Figure 2b) describes low frequency non-stationary long-term changes. The seasonal component (Figure 2c) shows a stable repeating pattern of change; thus, γ t represents only one element of the seasonal pattern, while γ sn(t),t represents a k x n matrix of seasonal shapes, where k is the number of seasons and n is the length of the time series. The remainder (Figure 2d) depicts variations considered as a normally distributed white noise process with mean zero and variance σ 2 [33,[49][50][51][52].
The trend (Figure 2b), which is the component of interest in this research, can be modeled as a local approximation to a linear trend: where η t and ζ t are normally distributed independent white noise processes with mean zero and variance σ 2 η and σ 2 ζ , respectively. The level and slope change slowly over time based on a random walk mechanism [49]. The reason behind analyzing the trend lies on the nature of the component. Its low frequency, non-stationary behavior, and small rate of change in a long-term basis allows for the identification of very subtle and gradual changes that may otherwise be obscured.
Using the "stR" package [53] in an R environment [54], trend components were extracted in the form of a new time series with as many values as the original time series. These new sets of values may or may not exhibit a general increase or decrease. This means that, in some cases, the trend component may actually show that the time series does not exhibit a trend. Therefore, to evaluate the presence of a monotonic trend in the trend component, a modified Mann-Kendall test for serially-correlated data [55][56][57] was applied. This non-parametric test evaluates the null hypothesis of independent and randomly ordered observations (i.e., no trend). The trend component (Figure 2b) describes low frequency non-stationary long-term changes.
The seasonal component ( Figure 2c) shows a stable repeating pattern of change; thus, represents only one element of the seasonal pattern, while ( ), represents a k x n matrix of seasonal shapes, where k is the number of seasons and n is the length of the time series. The remainder (Figure 2d) depicts variations considered as a normally distributed white noise process with mean zero and variance [33,[49][50][51][52].  The original Mann-Kendall [55,56] test determines if the Y variable (in this case LAI) tends to increase or decrease as the X variable (time) increases. To achieve this, data were ordered chronologically and a Kendall's τ correlation coefficient was computed as follows: where n is the number of observations and S is the test statistic, which is calculated by subtracting D, the number of "discordant pairs of observations", from C, the number of "concordant pairs of observations". In D, the later-in-time observation has a smaller value than the previous observation (i.e., Y decreases as X increases); while in C, the later-in-time observation has a larger value that the previous observation (i.e., Y increases as X increases). A total of n(n − 1)/2 pairs of observations are possible. If half of the pairs are concordant and half discordant, so that S = 0, τ will equal 0 (i.e., no trend). If all pairs are concordant, so that S = n(n − 1)/2, τ will equal 1 (i.e., strong increasing trend).
In contrast, if all pairs are discordant, so that S = −(n(n − 1)/2), τ will equal −1 (i.e., strong decreasing trend). Therefore, τ produces values, from −1 to 1, that show the strength of the monotonic association between the variable of interest (Y) and time (X) [58]. The significance of the trends is evaluated by comparing the standardized S statistic, Z = S/[var(S)] 0.5 , with the standard normal variate at the desired level or significance. Under the null hypothesis of trend absence, the S statistic tends to show zero mean and variance n(n − 1)(2n + 5)/18. However, the variance of S and, as a consequence, the p values are affected when serial correlation is detected. A positive serial correlation increases the probability of detecting a significant trend when it actually does not exist, while a negative serial correlation decreases the probability of detecting a significant trend when it actually exists [57]. A rank von Neumann ratio test [59], which assesses the null hypothesis of randomness, corroborated the presence of serial correlation in the extracted trend components of this study. Therefore, a S variance correction approach, proposed by [57], was implemented to address the issue of serial correlation and, consequently, generate the correct p values.
As explained, the modified Mann-Kendall test detects significant monotonic trends and provides a measure of their strength (i.e., how consistently decreasing or consistently increasing they are) through the Kendall's τ correlation coefficient. However, it does not specify their magnitude (i.e., how much they change). To account for the magnitude, a Theil-Sen's slope [60,61], which estimates a monthly rate of change that can be used to calculate the absolute change over the study period, was obtained by computing a set of N = n(n − 1)/2 slope values [62]: and then finding the median of those slopes: The "modifiedmk" R package [63] was used to calculate the Mann-Kendall p value, the Kendall's τ correlation coefficient, and the Theil-Sen slope for each of the trend components. After exploring the overall results and computing the occurrence of positive versus negative significant trends, a subset of negative trends was created for a deeper analysis. Negative trends are potential indicators of forest degradation. However, consistently decreasing trends alone (i.e., negative Kendall's τ correlation coefficients) cannot be considered processes of forest degradation because their magnitude of change may be insignificant. Similarly, negative magnitudes alone (i.e., Theil-Sen slopes) cannot be considered processes of forest degradation because their rate of change may be weakly monotonic and, therefore, not representative of a typical process of forest degradation. Therefore, the product of negative Kendall's τ correlation coefficients and negative Theil-Sen slopes is here proposed as a measure of forest degradation. This measure facilitates the identification of consistently decreasing trends with high magnitudes of change, which can be interpreted as severe processes of forest degradation. As a rule, the higher the product, the greater the degradation severity. The outcomes of the multiplication were reclassified into three categories of severity (i.e., high, moderate, and low) based on a geometrical interval classification method, which performs well with non-normally distributed data.

Validation of the Algorithm
To validate the results, the health conditions of 52,807 trees collected in 177 conglomerates were used as reference data [34]. This information was gathered between 2004 and 2017 by the Comisión Nacional Forestal (CONAFOR) as part of Mexico's National Forest and Soils Inventory (NFSI). A simple trend analysis of the tree conditions was performed and compared with the proposed measure of forest degradation. The Mexico's NFSI is completely updated every five years by revisiting 20% of the sampling conglomerates every year. Within the study area, there are approximately 1540 conglomerates separated from each other by roughly 5 km. These conglomerates have been measured two or, in a few cases, three times during the study period. Based on these spatio-temporal characteristics, only conglomerates that have been visited three times and whose centroids fall within a buffer of 150 m around the center of the MODIS pixels were considered for the analysis. The first criterion responds to the necessity of at least three observations to perform a trend analysis. The second criterion ensures that the conglomerates are representative of what the Instantaneous Field of View (IFOV) of the MODIS pixels captures. In other words, it warrants spatial coincidence. The 177 conglomerates that meet the two aforementioned criteria mimic the distribution of the original systematically collected NFSI data. Therefore, they are representative of the overall characteristics of the area (i.e., forest type, altitude, slope, climate). Each conglomerate contains four sampling sites that are distributed in an inverted "Y" shape and whose centroids are separated by 45.14 m. Each site has an area of 400 m 2 and a radius of 11.28 m. The health conditions of all trees falling within this area and having diameter greater than 7.5 cm were recorded. As a result, 13,905 trees were measured during the NFSI 2004-2009 in the 177 conglomerates, 15,730 during the NFSI 2009-2014, and 23,172 during the first 3 years of the NFSI 2014-2019. Two classification systems were used to record the tree conditions. The first system, employed for the NFSI 2004-2009 and the NFSI 2009-2014, comprises the six following classes: stump, dead-standing, very poor vigor, poor vigor, good vigor, and maximum vigor. The second system, used for the NFSI 2014-2019, comprises the five following classes: stump, dead-standing, low vigor, moderate vigor, and high vigor. This incompatibility between classification schemes implies a limitation that is addressed in Section 4.3. For the purposes of this validation, the conditions from the two classification schemes were unified as follows: stump, dead-standing, very poor or low, poor, moderate, good, and maximum or high.
To obtain a single measure representative of all trees within each conglomerate, a health index was developed. First, the number and proportion of trees within each class was calculated. Then, the proportion was multiplied by the weight of the class (i.e., stump = 1, dead-standing = 2, very poor or low vigor = 3, poor vigor = 4, moderate vigor = 4.5, good vigor = 5, and maximum or high vigor= 6).
The resulting values were summed to obtain an overall index that ranges from 1 to 6, where 6 is the optimum vigor condition (Table 1). Once three index values were computed per conglomerate (i.e., one per observation date), the presence of trends was evaluated by looking at the sign of the differences between two consecutive index values (later date minus earlier date). As there are three observations, two differences were calculated. When these differences share the same sign, a given location exhibits a trend, which is negative if both signs are negative and positive otherwise. The magnitude of change of all negative trends was calculated by subtracting the last observation from the first one and dividing the outcome by two. The resulting values were reclassified into three categories of severity (i.e., low, moderate, and high) based on a geometrical interval classification method. According to this methodology, each conglomerate can take one of the following classes: no trend, positive trend, low degradation, moderate degradation, or high degradation ( Table 2). By means of a confusion matrix, these results were compared with the results of the trend analysis proposed in this research.

LAI Trends
To achieve the objective of detecting forest degradation in Mexican forested areas, the STR procedure was used to decompose 295,298 monthly LAI pixel-based time series comprised of 192 values into trend, seasonal, and remainder components. Only the trend component was extracted from each of the time series to evaluate the presence and strength of significant monotonic trends through a modified Mann-Kendall test proposed by [57]. Additionally, a Theil-Sen slope was calculated to estimate the magnitude of change of the significant trends.
The trend analysis resulted in 63% of the study area (187,424 pixels) exhibiting trends at a significance level of 0.05 or higher. Of these significant trends, 52% are positive (i.e., they show a positive Kendall's τ correlation coefficient), while only 11% are negative ( Figure 3). Particularly, evergreen broadleaf forest presents the highest proportion of negative trends (18%), followed by evergreen needleleaf (12%), and deciduous broadleaf (7%). Significant trends do not cluster in a specific area, their distribution is regular across all the study area (Figure 4a). The majority of negative trends are located over the eastern portion of the study area (Figure 4b). The overall magnitude of change, measured by the Theil-Sen slope, ranges from −0.026 to 0.021 LAI/month and presents a median of 0.001 LAI/month. The steepest negative slopes are predominantly located in the state of Tamaulipas (eastern region), which, together with the state of Nayarit (western region), is also scenario of the steepest positive slopes. Lower magnitudes are located along the highest elevations of the Sierra Madre Oriental (east mountain range), as well as over the western states of Sinaloa and Durango (Figure 4c).  [57]. Additionally, a Theil-Sen slope was calculated to estimate the magnitude of change of the significant trends. The trend analysis resulted in 63% of the study area (187,424 pixels) exhibiting trends at a significance level of 0.05 or higher. Of these significant trends, 52% are positive (i.e., they show a positive Kendall's correlation coefficient), while only 11% are negative ( Figure 3). Particularly, evergreen broadleaf forest presents the highest proportion of negative trends (18%), followed by evergreen needleleaf (12%), and deciduous broadleaf (7%). Significant trends do not cluster in a specific area, their distribution is regular across all the study area (Figure 4a). The majority of negative trends are located over the eastern portion of the study area (Figure 4b). The overall magnitude of change, measured by the Theil-Sen slope, ranges from −0.026 to 0.021 LAI/month and presents a median of 0.001 LAI/month. The steepest negative slopes are predominantly located in the state of Tamaulipas (eastern region), which, together with the state of Nayarit (western region), is also scenario of the steepest positive slopes. Lower magnitudes are located along the highest elevations of the Sierra Madre Oriental (east mountain range), as well as over the western states of Sinaloa and Durango (Figure 4c).     Forest degradation is presented as the dimensionless product of the strength and magnitude of negative trends ( Figure 6). In this study area, the product ranges from 0 to 0.02 and was reclassified into three categories of severity (i.e., high, moderate, and low) based on a geometrical interval classification method. This classification method indicates that 385 km 2 (0.6% of the study area) are highly degraded, while 3406 km 2 (5.4%) and 3245 km 2 (5.1%) are moderately and slightly degraded, respectively. Forest degradation is mostly low in the west and generally higher in the eastern region of Mexico, which is characterized by consistently decreasing trends with high magnitudes of change. Although the western portion also displays consistently decreasing trends (Figure 5a), the magnitudes of change are generally not high (Figure 5b). The most degraded areas are irregularly distributed along the Sierra Madre Oriental and, to a lesser extent, over southwest Chihuahua, the western fringe of Durango and Jalisco, and some isolated areas of Nayarit.
Regarding forest types, evergreen broadleaf has the highest median severity score of forest degradation (0.0011), followed by evergreen needleleaf (0.0008), and deciduous broadleaf (0.0006). In general, evergreen broadleaf forest is characterized by a combination of consistently decreasing trends with high magnitudes of change, evergreen needleleaf by a combination of irregularly decreasing trends with high magnitudes of change, and deciduous broadleaf by a combination of consistently decreasing trends with low magnitudes of change. Forest degradation is presented as the dimensionless product of the strength and magnitude of negative trends (Figure 6). In this study area, the product ranges from 0 to 0.02 and was reclassified into three categories of severity (i.e., high, moderate, and low) based on a geometrical interval classification method. This classification method indicates that 385 km 2 (0.6% of the study area) are highly degraded, while 3406 km 2 (5.4%) and 3245 km 2 (5.1%) are moderately and slightly degraded, respectively. Forest degradation is mostly low in the west and generally higher in the eastern region of Mexico, which is characterized by consistently decreasing trends with high magnitudes of change. Although the western portion also displays consistently decreasing trends (Figure 5a), the magnitudes of change are generally not high (Figure 5b). The most degraded areas are irregularly distributed along the Sierra Madre Oriental and, to a lesser extent, over southwest Chihuahua, the western fringe of Durango and Jalisco, and some isolated areas of Nayarit. Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 21 Evergreen needleleaf shows the highest proportion of highly degraded areas (8%), followed by evergreen broadleaf (7%), and deciduous broadleaf (3%, Figure 6). Deciduous broadleaf is almost free of highly degraded areas in the west part of the study area, but it presents a concentration of highly degraded pixels in the west portion, specifically over the state of Tamaulipas (Figure 7a). Evergreen broadleaf exhibits highly degraded pixels distributed throughout its range as well as a high proportion of moderately degraded pixels (Figure 7b). Unlike the other forest types, evergreen needleleaf displays more degraded areas in the east portion of the study area (Figure 7c). Regarding forest types, evergreen broadleaf has the highest median severity score of forest degradation (0.0011), followed by evergreen needleleaf (0.0008), and deciduous broadleaf (0.0006). In general, evergreen broadleaf forest is characterized by a combination of consistently decreasing trends with high magnitudes of change, evergreen needleleaf by a combination of irregularly decreasing trends with high magnitudes of change, and deciduous broadleaf by a combination of consistently decreasing trends with low magnitudes of change.
Evergreen needleleaf shows the highest proportion of highly degraded areas (8%), followed by evergreen broadleaf (7%), and deciduous broadleaf (3%, Figure 6). Deciduous broadleaf is almost free of highly degraded areas in the west part of the study area, but it presents a concentration of highly degraded pixels in the west portion, specifically over the state of Tamaulipas (Figure 7a). Evergreen broadleaf exhibits highly degraded pixels distributed throughout its range as well as a high proportion of moderately degraded pixels (Figure 7b). Unlike the other forest types, evergreen needleleaf displays more degraded areas in the east portion of the study area (Figure 7c

Validation
The comparison between the trend analysis of the MODIS LAI data and the reference data shows an overall agreement of 63%. This concordance between the remotely sensed products and the field observations is slightly better in the west side (69%) than in the east side (56%) of the study area ( Figure 8). According to the confusion matrix (Table 3), the "high degradation" and "no trend" classes have the highest accuracy. All of the MODIS pixels classified as "high degradation" are actually highly degraded, while 81% of the MODIS pixels classified as "no trend" do not actually exhibit a significant trend. On the other hand, half of the reference data identified as highly degraded were correctly classified as "high degradation", whereas 72% of the reference data without any trend were correctly classified as "no trend". The comparison between the trend analysis of the MODIS LAI data and the reference data shows an overall agreement of 63%. This concordance between the remotely sensed products and the field observations is slightly better in the west side (69%) than in the east side (56%) of the study area ( Figure 8). According to the confusion matrix (Table 3), the "high degradation" and "no trend" classes have the highest accuracy. All of the MODIS pixels classified as "high degradation" are actually highly degraded, while 81% of the MODIS pixels classified as "no trend" do not actually exhibit a significant trend. On the other hand, half of the reference data identified as highly degraded were correctly classified as "high degradation", whereas 72% of the reference data without any trend were correctly classified as "no trend". In contrast, the "moderate degradation" class exhibits the poorest accuracy. Forty five percent of the MODIS pixels classified as "moderate degradation" are actually moderately degraded and only 19% of the reference data identified as moderately degraded were correctly classified as "moderate degradation". These pixels were erroneously classified as "no trend" and "positive trend". Finally, the class with the highest contrast between accuracy types is "positive trend". Although 78% of the reference positive trends were correctly classified as "positive trend", only 27% of the MODIS pixels classified as "positive trend" exhibited positive trends. In contrast, the "moderate degradation" class exhibits the poorest accuracy. Forty five percent of the MODIS pixels classified as "moderate degradation" are actually moderately degraded and only 19% of the reference data identified as moderately degraded were correctly classified as "moderate degradation". These pixels were erroneously classified as "no trend" and "positive trend". Finally, the class with the highest contrast between accuracy types is "positive trend". Although 78% of the reference positive trends were correctly classified as "positive trend", only 27% of the MODIS pixels classified as "positive trend" exhibited positive trends.

Potential Drivers of Forest Degradation
The objective of this research was to detect forest degradation through a trend analysis of monthly MODIS LAI collected over Central Mexico from 2002 to 2017. The overall results indicate that about half of the study area (52%, Figure 3) has exhibited positive trends of forest change, while the other half has either remained unchanged (37%) or experienced negative trends (11%). Of the 7036 km 2 exhibiting negative trends, the proposed measure of forest degradation estimates that 385 km 2 are highly degraded, 3406 km 2 moderately degraded, and 3245 km 2 slightly degraded. Most of the moderate and highly degraded areas are distributed over the east side of the study area ( Figure 6) and evergreen broadleaf seems to be the most affected forest type (Figure 7).
The authors are not aware of other studies that examine vegetation trends at the regional scale in Mexico, which precludes a direct comparison of the results of this research. However, the set of INEGI's Land Use and Vegetation Maps provides points of comparison. For example, this research reveals that 7036 km 2 (11% of the study area) of the study area have experienced forest degradation between 2002 and 2017. According to the most recent INEGI's cartography, 14,900 km 2 (23%) were covered by perturbed secondary vegetation in 2014. Of that area, 12,787 km 2 were already perturbed in 2002. As the level of disturbance of these areas is unknown for both dates, they could have experienced no changes, or slight negative or positive trends. In contrast, 1686 km 2 were covered by undisturbed primary vegetation in 2002, which clearly translates into negative trends since in 2014 were classified as perturbed. The remaining 427 km 2 were agriculture or water [42,64]. Thus, although this research and the INEGI's cartography are not directly comparable, both show that important processes of forest degradation are of the order of a few thousands of km 2 .
In terms of spatial distribution, this research shows more severe levels of forest degradation in the eastern portion. The most common disturbance recorded by Mexico's NFSI on this side of the study area is the presence of epiphytes, followed by drought and defoliating insects. On the western portion, the most frequent disturbance is fire, followed by logging and drought [34]. Although not all forest degradation can be attributed to these disturbances, they provide an indication of the factors affecting forest health in Central Mexico. For instance, changes in microclimate or vegetation structure may affect the abundance and distribution of epiphytes [65]. This seems to be a major problem because, although epiphytes do not obtain nutrients from their host as parasitic plants do, they use the host as support, which suffocates the branches and eventually kills the trees [66]. Drought has also been recorded as an important driver of forest disturbance. Indeed, in 2009 and 2011, Central Mexico suffered the most severe droughts in seven decades [67]. The effects of these droughts continue to affect the region and, according to [34], the impacts on forest vegetation have been most pronounced on the east side of the study area, where higher levels of forest degradation are encountered based on this research.

Methodological Approach
Assessing forest degradation has been a challenging task due to the generally slow-changing nature of the process, which demands long periods of observation with high frequency of records. In this research, the analysis of the trend component of monthly LAI time series helped to overcome the challenge of identifying subtle and gradual changes. This component extracts low frequency changes that occur throughout the time series. Figure 9 demonstrates how the proposed algorithm successfully classifies the trend components based on their strength and magnitude of change. Figure 9a is a good example of a very subtle consistently decreasing change that may pass unnoticed if only the original time series is analyzed. This was confirmed by a Seasonal Kendall test [58], which did not detect a significant trend in the raw time series. Figure 9b illustrates a very high magnitude of change, but a non-consistently decreasing trend. Although it is classified as highly degraded, it does not exhibit one of the highest scores of forest degradation because its trajectory is not typical of a forest degradation process. It is worth mentioning that many pixels show this abrupt change around 2009, which may correspond to the severe droughts that Central Mexico experienced in 2009 and 2011. Figure 9c  Assessing forest degradation has been a challenging task due to the generally slow-changing nature of the process, which demands long periods of observation with high frequency of records. In this research, the analysis of the trend component of monthly LAI time series helped to overcome the challenge of identifying subtle and gradual changes. This component extracts low frequency changes that occur throughout the time series. Figure 9 demonstrates how the proposed algorithm successfully classifies the trend components based on their strength and magnitude of change. Figure  9a is a good example of a very subtle consistently decreasing change that may pass unnoticed if only the original time series is analyzed. This was confirmed by a Seasonal Kendall test [58], which did not detect a significant trend in the raw time series. Figure 9b illustrates a very high magnitude of change, but a non-consistently decreasing trend. Although it is classified as highly degraded, it does not exhibit one of the highest scores of forest degradation because its trajectory is not typical of a forest degradation process. It is worth mentioning that many pixels show this abrupt change around 2009, which may correspond to the severe droughts that Central Mexico experienced in 2009 and 2011. Figure 9c

Potential Sources of Validation Uncertainty
The difficulty of validating degradation models has been recognized and attributed to the considerable amount of economic and human resources needed to collect reference data during long periods of study [68]. This research made use of reference data that were not collected for the specific purpose of validating this study. Therefore, although these data are considered the best possible

Potential Sources of Validation Uncertainty
The difficulty of validating degradation models has been recognized and attributed to the considerable amount of economic and human resources needed to collect reference data during long periods of study [68]. This research made use of reference data that were not collected for the specific purpose of validating this study. Therefore, although these data are considered the best possible information available, they present some limitations. In this context, there are at least four reasons that may be associated with the moderate agreement (i.e., 63%) between the remotely sensed LAI and the field observations of tree conditions. First, the spatial extent of the reference data is not ideally compatible with the spatial resolution of the MODIS pixels. Each sampling conglomerate extends over 1600 m 2 , while each MODIS pixel covers 216,516 m 2 . This means that the reference data are representative of less than 1% of the MODIS pixel with which they are compared. Second, the classification schemes of tree conditions are slightly different between forest inventories. The NFSI 2004-2009 and the NFIS 2009-2014 recognize four vigor conditions (i.e., very poor, poor, good, and maximum), while the NFSI 2014-2019 recognizes only three (i.e., low, moderate, and high). This incompatibility hindered an appropriate comparison of the vigor conditions through time. Third, the number of field observation dates within the study period is limited for the purposes of a trend analysis. Fourth, the collection of tree conditions may be subjective since it depends on the criterion and expertise of the data collector.
The accuracy of future validation efforts may be improved by including all field observations regardless of the position within a pixel, by establishing a unique classification scheme of tree conditions, by considering a longer study period that allows for the incorporation of more field observation dates, and by testing other reference data such as crown density or proportion of live crown. The latter two variables were not considered in this work due the high proportion of missing records in the study area and period of interest. On the other hand, forest degradation studies may need an entirely different validation approach. Although an index that aimed to identify the general condition of the forest in each conglomerate was created (see Section 2.4 for details), forest degradation studies may need field observations that evaluate forest as a whole in a given area rather than as single individuals (i.e., trees) or plots. Collecting data to properly validate forest degradation studies at a regional scale represents an area of opportunity for future research.

Conclusions
This research detected forest degradation over Central Mexico through a trend analysis of monthly MODIS LAI collected from 2002 to 2017. The results indicate that 52% of the study area has experienced positive trends of vegetation change, 37% has remained unchanged, and 11% has suffered some level of forest degradation. The proposed algorithm estimated that 385 km 2 are highly degraded, 3406 km 2 moderately degraded, and 3245 km 2 slightly degraded. Most of the moderate and highly degraded areas are distributed over the east side of the study area and evergreen broadleaf appears to be the most affected forest type.
This research overcame the problem of detecting very slow and gradual changes associated with forest degradation by measuring the strength and magnitude of the trend component of pixel-based LAI time series. None of the algorithms cited in this paper have made use of LAI as vegetation greenness indicator, nor have performed this type of trend analysis to detect forest degradation. The main limitation of this study is the moderate accuracy of the validation effort (i.e., 63%). Some suggestions to improve this accuracy are expressed in Section 4.3, but also a different approach to validate this type of studies is suggested as an area of opportunity for future research.
This study is expected to be a contribution to the regional efforts to measuring and monitoring forest degradation in Mexico. It could serve as one of the first steps towards the mitigation of the serious impacts that forest degradation has on the environmental services that these ecosystems provide. Thus, this work suggests that this methodology can be applied at the national level and be used as a first look into hot spots of forest degradation. Funding: This research was funded by the Department of Geography, the Graduate College, and the International Office of Texas State University, as well as by the Consejo Nacional de Ciencia y Tecnología (CONACYT).