Remote Sensing of Pasture Degradation in the Highlands of the Kyrgyz Republic: Finer-Scale Analysis Reveals Complicating Factors

Abstract

Degradation in the highland pastures of the Kyrgyz Republic, a small country in Central Asia, has been reported in several studies relying on coarse spatial resolution imagery, primarily MODIS. We used the results of land surface phenology modeling at higher spatial resolution to characterize spatial and temporal patterns of phenometrics indicative of the seasonal peak in herbaceous vegetation. In particular, we explored whether proximity to villages was associated with substantial decreases in the seasonal peak values. We found that terrain features—elevation and aspect—modulated the strength of the influence of village proximity on the phenometrics. Moreover, using contrasting hotter/drier and cooler/wetter years, we discovered that the growing season weather can interact with aspect to attenuate the negative influences of dry conditions on seasonal peak values. As these multiple contingent and interactive factors that shape the land surface phenology of the highland pastures may be blurred and obscured in coarser spatial resolution imagery, we discuss some limitations with prior and recent studies of pasture degradation.

1. Introduction

The Kyrgyz Republic (Kyrgyzstan) is a small, highly mountainous nation of ~6.5 million in 2019, where more than 60% live in rural areas and where agropastoralism is the predominant land use [1]. Degradation of pasture resources in highly mountainous Kyrgyz Republic has been both widely reported [2,3,4,5,6,7,8] and disputed [9,10]. Many studies have attempted to detect land degradation in Central Asia by using remote sensing products based on finer temporal but coarser spatial resolution imagery [4,7,8,11,12,13,14,15,16]. While the finer temporal resolution products are better able to obtain clear views of the land surface, coarser spatial resolution products blend together heterogeneous surfaces. This blending or mixing is of particular concern in mountainous terrain, where differential insolation regimes arising from the interaction of aspect and slope generates microclimates that can accommodate distinct vegetation communities exhibiting different phenologies [3,17,18,19].

Untangling the differential influences of climate and human activity is complicated by coarser spatial resolution [7,14,20]. Pastoralism in Kyrgyzstan is characterized by transhumance, the seasonal movement of livestock to fresh pastures, and vertical transhumance in particular, where the herds are moved from low-elevation winter pastures near villages through transitional pastures in spring (and again in fall) to higher elevation summer pastures distant from settlements. This combination of spatial heterogeneity, terrain effects, and seasonality of pasture use presents challenges for detecting pasture degradation from coarser resolution remote sensing imagery. A further complication is presented by the high interannual variation in weather arising from the land-locked location of the country, its relatively high elevation, and regional climate change [12,21,22,23,24].

Here, we explore the pasture degradation question across rural Kyrgyz Republic using phenological metrics (phenometrics) derived from finer spatial resolution, but less frequent imagery (Landsat 30 m) linked to more frequent, but coarser spatial resolution products (MODIS 1 km) across 17 years (2001–2017). The phenometrics that we use here indicate the first seasonal peak in NDVI and the amount of thermal time (measured as accumulated growing degree-days calculated from land surface temperature time series) required to reach that peak NDVI. These phenometrics can be calculated for each pixel in each year for which there are sufficient, high-quality data and for which the joint time series exhibit sufficient seasonality to enable the modeling of the land surface phenology (LSP). Pasture degradation can be evaluated with phenometrics in a variety of ways, including trends, variation, extremes, and abrupt spatial and/or temporal shifts [12,20,25,26]. However, we are interested here in three questions about patterns that can complicate interpretation of patterns and trends: (1) What are the spatial patterns of the phenometrics as a function of distance from village center? (2) How do these patterns change as a function of elevation? and (3) How do these patterns change as a function of aspect?

Although we focus on the patterns of the temporal mean phenometrics during the study period, we also examine the patterns during two years with contrasting weather: hotter, drier 2007 versus cooler, wetter 2009. This study leverages the findings from prior studies [25,26] using the same approach to modeling land surface phenology and advances our understanding of the spatio-temporal dynamics of vegetation in the socio-ecological landscapes of montane Central Asia, specifically in the Kyrgyz Republic.

2. Materials and Methods

2.1. Study Area

The study area focuses on pasture lands within the territory of the Kyrgyz Republic that neighbors Uzbekistan (west), Kazakhstan (north), China (east and southeast), and Tajikistan (southwest) (Figure 1). The total area of the country is shy of 200,000 km² (96% in land), and the 2019 population was about 6.5 million, according to the World Bank (https://data.worldbank.org/country/kyrgyz-republic; accessed on 20 June 2021). It is a highly mountainous country, where more than 56% of the territory lies above 2500 m and where the mountain ranges of the Tien Shan, Pamir, and Alatau cover more than 90% of the total land area [27]. Pastoral rangelands constitute 87% [1] of the agricultural lands in the Kyrgyz Republic. Less than 10% of the land is used for crops, while forests cover only about 5%. Our study period extends from 2001 through 2017. The Kyrgyz Republic is divided into seven provinces (oblasts)—Talas, Chuy (including the capital city, Bishkek), and Issyk-Kul in the northern part as well as Jalal-Abad, Naryn, Osh, and Batken in the southern part—and 40 districts (rayons).

2.2. Geospatial Data

The land surface phenology metrics (or phenometrics) and terrain information used for this study were supplied from our previous work [25], where a detailed description of data, its processing, and phenometrics calculation methodology can be found. We provide an overview here.

To calculate LSP metrics, we used two products: MODIS land surface temperature (LST) and Landsat surface reflectance NDVI.

We downloaded two tiles (h23v04 and h23v05) of 8-day MODIS Terra and MODIS Aqua Land Surface Temperature (MOD11A2/MYD11A2 V006) products at 1 km spatial resolution [31] from 2001 (MODIS/Terra) and from 2002 (MODIS/Aqua) up to the end of 2017. We merged tiles, removed poor quality pixels, converted units from Kelvin to °C, and reprojected data to Albers Conic Equal Area at 30 m spatial resolution using bilinear resampling.

The surface reflectance NDVI dataset from Landsat Collection 1 Tier 1 Level-1 Precision and Terrain (L1TP) of Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) was acquired for years 2001 to the end of 2017 from the USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On-Demand Interface (https://espa.cr.usgs.gov/). We downloaded 13,285 images across 33 unique tiles (WRS-2 Paths 147–155 and Rows 30–33) which were already projected into Albers Conic Equal Area. We masked poor-quality pixels and applied an inter-calibration equation to adjust Landsat 5 TM surface NDVI and Landsat 7 ETM+ surface NDVI to the surface Landsat 8 OLI NDVI, which on average was shown to have higher values (cf., Table 3 in [32], Surface NDVI from OLI = 0.0235 + 0.9723 ETM+). Because of the small differences between the Landsat 5 TM and Landsat 7 ETM+ data [33,34], we used the same equation for both datasets.

For the analyses conducted in this study, we additionally used three other geospatial datasets: (1) digital elevation model, (2) pasture land-use mask, and (3) point coverage of settlements. We downloaded 133 tiles of SRTMGL1, the NASA Shuttle Radar Topography Mission Global 1 arc second (~30 m) V003 elevation product [30] from USGS Earth Explorer (https://earthexplorer.usgs.gov/). Tiles were merged and reprojected into the Albers Conic Equal Area at 30 m spatial resolution using bilinear resampling. We then generated aspect and slope layers. For this study, we masked out pixels from all layers where slope was greater than 30 degrees. Additionally, we created sets of three terrain layers where we left only pixels on the contrasting aspects (northern: NW-N-NE-E, southern: SE-S-SW-W) at four elevation ranges: 1800–2400 m, 2400–2900 m, 2900–3400 m, and 3400–4000 m, and elevation−aspect interactions.

The pasture land-use class (122,405 km²) was obtained from a Soviet-era land use map that was updated in 2008 using Landsat 7 ETM+ and MODIS datasets for the CACILM project [28,29]. We used those data to mask out non-pasture pixels in the terrain and LSP layers to focus only on the pasture land-use pixels.

The settlement point layer was also obtained from the CACILM project where the dataset was collected and consolidated by ECONET WWF Project [35] national teams. Over Kyrgyzstan, the layer has 703 points representing detailed administration levels, including small villages. We conducted a quality check on this dataset and eliminated points that were clearly mislocated, yielding 617 point locations (

Table 1. Number of settlement points after dataset revision and filtering per province (oblast).

Province (Oblast)	Total Points	Points after 10 km Filter
Batken	63	36
Chuy	67	33
Issyk-Kul	100	45
Jalal-Abad	113	58
Naryn	96	44
Osh	137	56
Talas	41	21
TOTAL	617	293

) for analysis.

2.3. Methods

2.3.1. Land Surface Phenology

We used a downward-arching convex quadratic (CxQ) function to characterize LSP [36,37,38]. The model uses a vegetation index—here the NDVI from a Landsat surface reflectance time-series—as proxy for active green vegetation, and thermal time—here accumulated growing degree-days (AGDD) from MODIS LST—as proxy for insolation.

To obtain AGDD, we first transformed two diurnal and nocturnal observations from the MODIS on Terra and on Aqua into a mean LST using the following Equation (1):

mean LSTt = [max(LSTt_TERRA1030, LSTt_AQUA1330) + min(LSTt_TERRA2230, LSTt_AQUA0130)]/2

(1)

where LSTt_TERRA1030 is the LST for period t at the Terra daytime overpass, LSTt_AQUA1330 is the LSTt at the Aqua daytime overpass, LSTt_TERRA2230 is the LSTt at the Terra nighttime overpass, and LSTt_AQUA0130 is the LSTt at the Aqua nighttime overpass.

We filled gaps that resulted from missing or filtering excluded pixels using the Seasonally Decomposed Missing Value Imputation method [39], replaced all negative values with 0 °C, and calculated growing degree-days GDD (2) at compositing period t as the maximum of mean LST and T_base, where T_base was set to 0 °C [37,40].

GDD_t = max(mean LSTt − T_base, 0)

(2)

For each year, we produced 46 GDD composites that were multiplied by 8 to account for the 8-day MODIS product composite period and accumulated each year into time series of AGDD (3), with an annual reset in January to 0 °C:

AGDD_t = AGDD_t-1 + (GDD_t × 8)

(3)

Prior the CxQ LSP model fitting, it was necessary to further clean and filter NDVI and AGDD datasets to reduce noise and spurious data. Therefore, we filtered out observations with NDVI < 0.1 and AGDD < 100 to avoid non-vegetated or snow-covered pixels. To account for cloud contamination that may have slipped through the masking process, we looked for unusual abrupt dips in the NDVI time-series. We first calculated the simple average of NDVI observations on either side of the focal observation. We then calculated the percentage difference between the average NDVI and the focal observation and excluded observations that were ≥15% than the average of the two neighboring observations [25,26]. Having NDVI and AGDD datasets prepared for each pixel and each year (from 2001 to 2017), we applied the CxQ LSP model shown in (4):

NDVI = α + β × AGDD + γ × AGDD²

(4)

Using the fitted coefficients (4) for the intercept (α), slope (β), and quadratic (γ) parameters, we calculated two LSP phenometrics: Peak Height [PH=α − (β²/4 × γ)], which is the maximum modeled NDVI; and Thermal Time to Peak [TTP = −β/2 × γ], which is the quantity of AGDD required to reach PH and corresponds to the duration of green-up phase. To control CxQ LSP model fitting performance, we used a suite of six quality criteria: (i) the fitted quadratic parameter coefficient was less than zero (γ < 0); (ii) TTP greater than the AGDD at the first observation; (iii) adjusted R² greater than 0.7; (iv) Root Mean Square Difference (RMSD) less than 0.08; (v) at least seven observations in the time series to be fit, where at least three observations were distributed before and at least three after the PH; and (vi) the PH less than or equal to 1.0.

If any criterion was not fulfilled during the fitting process, then the last observation in the time series was removed and the model fitting procedure was rerun over the reduced time series. We repeated this fitting procedure until either the fitted model passed all criteria or the length of the time series was fewer than seven observations. In the latter case, the model fit for that pixel was labeled as failed and no phenometrics were calculated for that pixel.

Next, for each pixel where model fit was successful and phenometrics were obtained, we calculated the mean values of PH and TTP across 17 years. We also highlighted PH and TTP from the years 2007 and 2009, which were drier and wetter weather conditions, respectively. As mentioned in Section 2.2, we used the pasture land-use layer to mask data.

2.3.2. Settlement Ring Buffer Analyses

From the settlement point coverage, we randomly selected (ArcMap Software 10.6.0.8321, Random Generator Type: default ACM599, seed:0) points with the spatial constraint of 10 km from each other, which resulted in 293 focal points for analysis (Table 1).

Then, we created 10 ring buffers around each settlement focal point from 500 m to 5000 m distance in 500 m intervals (0–500 m, 500–1000 m, 1000–1500 m, 1500–2000 m, 2000–2500 m, 2500–3000 m, 3000–3500 m, 3500–4000 m, 4000–4500, and 4500–5000 m). The spatial constraint of 10 km ensured that the ring buffers would not intersect. Finally, we prepared nine datasets (raster stacks) at four elevation classes—(1) 1800–2400 m, (2) 2400–2900 m, (3) 2900–3400 m, and (4) 3400–4000 m—yielding 36 in total, as follows:

• Elevation classes {1–4}, all aspects, slopes < 30°, PH_mean, TTP_mean;
• Elevation classes {1–4}, all aspects, slopes < 30°, PH₂₀₀₇, TTP₂₀₀₇;
• Elevation classes {1–4}, all aspects, slopes < 30°, PH₂₀₀₉, TTP₂₀₀₉;
• Elevation classes {1–4}, northern aspects, slopes < 30°, PH_mean, TTP_mean;
• Elevation classes {1–4}, northern aspects, slopes < 30°, PH₂₀₀₇, TTP₂₀₀₇;
• Elevation classes {1–4}, northern aspects, slopes < 30°, PH₂₀₀₉, TTP₂₀₀₉;
• Elevation classes {1–4}, southern aspects, slopes < 30°, PH_mean, TTP_mean;
• Elevation classes {1–4}, southern aspects, slopes < 30°, PH₂₀₀₇, TTP₂₀₀₇;
• Elevation classes {1–4}, southern aspects, slopes < 30°, PH₂₀₀₉, TTP₂₀₀₉.

Once these data were prepared, we extracted pixels from each set of settlement point ring buffers across all 36 datasets. During the extraction process, we ensured that only those pixels whose center point was within the ring buffer polygon were included so the same pixel would not spill into the neighboring ring buffers. Because of the prior pixel filtering by pasture land-use map, elevation classes, slope threshold value, and the fact that LSP CxQ modeling did not always succeed (i.e., no phenometrics calculated), the number of extracted pixels varied within each ring buffer for each dataset, and in some instances, there were no pixels to extract due to the multiple filtering steps. When there were extracted pixels, we calculated the mean values for each ring buffer, which were used for the complementary analyses of the influence of elevation, aspect, growing season weather, and distance from village on the two phenometrics.

3. Results

3.1. Buffer Mean Values of the Peak NDVI as a Function of Distance and Elevation

Figure 2 displays the spatial mean ±2SE values of the fitted NDVI Peak Height within each ring buffer calculated from the mean PH values across all years. Not surprisingly, there is considerable variation, but some patterns are evident. First, the PH decreased with increasing elevation. Second, the variation in the PH increased with elevation. Third, PH increased with distance from villages at 2400–2900 m and 2900–3400 m.

3.2. Contrasting Mean Values of Phenometrics Nearby and Far from Villages

Figure 3 subsets the distributions to focus on either end of the spatial series: the 0–500 m ring buffer captured those pasture areas closest to villages and the 4500–5000 m ring buffer captured pastures that were far from the focal village and not intruding on the ring buffer of another village. For PH (Figure 3, left panel), there was no difference in the lowest elevation class (1800–2400 m), but the differences between PH values nearby and far from villages were strong at both 2400–2900 m and 2900–3400 m, where the PH values in the distant ring buffer were higher than in the nearby ring buffer. However, the distributions of the mean PH values were significantly different between the 0–500 m and 4500–5000 m buffers only at the 2900–3400 m elevation range, according to the Kolmogorov–Smirnov two-sample test with the Dunn-Šidák correction for post-hoc multiple comparisons [41]. Note that in the highest elevation class (3400–4000 m), there were no pasture areas in the ring buffer nearest (0–500 m) to villages. It is also evident that the variation of PH in the nearby ring buffer appeared larger at higher elevations relative to that in the lowest elevation class, but this difference may result from fewer samples at high elevations.

The Thermal Time to Peak phenometric (Figure 3, right panel) decreased with elevation, as expected, and the distant TTP ring buffer means were consistently—but not significantly—lower than the nearby values.

3.3. Influence of Weather on Phenometrics

Figure 4 shows the PH (left) and TTP (right) means for the ring buffers nearby versus far from the villages. For PH, the contrast between years and distance from villages was significant only for the 2900–3400 m elevation class (Figure 4, left). At the lowest elevation class, the nearby ring buffer means were not significantly different, but the distant means were. TTP was clearly lower during 2009, the wetter year, at all elevations, and lower at higher elevations, as expected. Evaluation by the Kolmogorov-Smirnov two-sample test (with the Dunn-Šidák correction) confirms what Figure 4 (right panel) shows: the TTP distributions were significantly different between 2007 and 2009 at every elevation class, but the difference between 0–500 m and 4500–5000 m ring buffer distributions was significant only at the 2900–3400 m elevation class and only in the drier year of 2007.

3.4. Influence of Aspect on Phenometrics

We tested for differences in the distributions of buffer means as a function of aspect (NW-N-NE-E vs. SE-S-SW-W) by distance from settlement point and elevation class using the Kolmogorov-Smirnov two-sample test with the Dunn-Šidák correction for post-hoc multiple comparisons. Northerly aspects consistently showed higher Peak Height means than in southerly aspects, but these differences were not statistically different at every distance or elevation class (Figure 5). In the lowest elevation class (1800–2400 m), the effect of aspect on PH was significant starting at the 2500–3000 m ring buffer; in the 2400–2900 m elevation class, the difference in PH between aspects appeared significant 500 m closer to the settlement point (i.e., at the 2000–2500 m ring buffer). In the two higher elevation classes, the aspect differences were not significant, but the variation about the means appears lower in the more distant ring buffers.

3.5. Interaction of Weather, Aspect, and Distance on Phenometrics

Figure 6 shows the distributions of phenometrics for northerly and southerly aspect slopes from pasture areas nearby (0–500 m) and far from (4500–5000 m) villages in the drier year of 2007 and the wetter year of 2009. The Kolmogorov-Smirnov two-sample test with the Dunn-Šidák correction reveals four patterns of interest in the Peak Height distributions: (1) the aspect effect on PH is significant at the two lower elevation classes only far from villages; (2) at the far ring buffer at each elevation class between 1800 m and 3400 m, the distributions of PH means from southerly aspects in the drier year of 2007 are significantly different from the northerly aspect PH means in the wetter year of 2009; (3) at each distance and elevation, distributions of PH means are not significantly different by year for the same aspect, i.e., the distributions of southerly (or northerly) aspect PH means are not different between years; and (4) the distributions of northerly aspect PH means in the drier year of 2007 are not significantly different from the southerly aspect PH means in the wetter year of 2009. This last pattern is clearer in Figure 7 (top panels) and illustrates how the northerly aspect slopes moderate the impact of dry years on pasture LSP.

The patterns in the TTP distributions were not surprising (Figure 6, right): (1) decreasing TTP with increasing elevation; (2) lower TTP values in the cooler, wetter year of 2009; and (3) no apparent aspect effect in the TTP buffer mean distributions. The bottom panels in Figure 7 show more clearly how strongly the TTP mean distributions diverge between years.

4. Discussion

The results of these linked analyses address the questions posed earlier. The Peak Height phenometric tends to increase in pasture areas as distance from the village increases. One interpretation of this pattern is the pasture areas near villages tend to be degraded. Winter pastures are closer to villages than either summer or transitional pastures and have been shown to exhibit significant differences in biogeochemistry and vegetation composition [5]. However, this general pattern can be modulated by the terrain, both elevation and aspect. Of the four elevation classes, significantly larger PH differences were seen between nearby and distant buffer rings in the 2900–3400 m elevation class (Figure 3, left). TTP decreased with elevation class but showed no clear difference between nearby and distant values (Figure 3, right), due in part to the much coarser spatial resolution of the LST data.

For the same reason, there was no clear influence of aspect on TTP (Figure 6, right panels; Figure 7, bottom panels). In contrast, the influence of aspect on PH was strongest in the lower elevation classes (Figure 5, bottom panels). When distant from a settlement, the PH on a northern aspect slope during a dry year appeared very similar to the PH on a southern aspect slope during a wet year (Figure 7, top right). However, in pasture areas closest to settlements, this aspect influence vanishes (Figure 7, top left).

The advantages of our study are several. First, we integrated long time series of two independent remote sensing products through a simple biometeorological model of land surface phenology that responds to the progress of growing season temperature rather than to the mere passing of days [37,38]. Second, the CxQ model has been shown to capture the initial seasonal peak well in grassland LSP [25,37,38]. Third, the finer spatial resolution of our analysis captured the influence of terrain features on vegetation growth and development as captured by the phenometrics [25].

Four limitations of our study are key. First, even the 30 m Landsat data can miss important landscape features influencing LSP [42]. Second, the 1 km spatial resolution of the MODIS LST product is coarser than would be optimal for use in rugged terrain. However, there are no effective alternatives at this point. Third, the CxQ model captures the initial seasonal peak but not necessarily the post-peak decay, particularly if there is heavy grazing post-peak [25]. This limitation is not as serious as it may first appear: we expect clear evidence of pasture degradation to appear only after several years of overgrazing, particularly if coupled with a severe drought. Fourth, while the pasture land-use mask was critical for omitting non-pastoral land-uses, there is substantial uncertainty associated with its development and accuracy. We have no accuracy assessment associated with it, and we expect that commission error is more likely than omission error, which would have the effect of increasing variance and decreasing the significance of the patterns. Finally, a key caveat for interpretation of these results should be noted: the 5 km extent of analysis was not meant to capture the summer pastures associated with villages. Summer pastures can be 10–50 km or more from a given village, and the spatial allocation and arrangement of summer pastures can be quite complicated.

Given our findings, to what extent do they raise questions about previous remote sensing studies on widespread degradation of pasture resources? The relationship between the scale of observation and the scale of the phenomenon of interest is crucial to the understanding and interpretability of the remote sensing data [43]. It is clear from these results that terrain effects can be obscured by coarser spatial resolution data, yet there is another scale to consider as well. Prior studies have demonstrated significant effects of climate oscillation modes on LSP in Central Asia [12,44,45]. However, more recent work at higher spatial resolution was not able to discern a significant influence from climate modes because the influence on LSP was overridden by local landscape structure [26].

Another facet of the degradation monitoring problem relates to the impact of spatial heterogeneity on the characterization of LSP. An important empirical study [46] explored how spatial heterogeneity in the land surface phenology observed at finer spatial resolution influences the timing of phenophase detection when observed at coarser spatial resolution. They found, when the land cover was homogeneous, that greater than 60% of Landsat 8 OLI 30 m pixels exhibiting start of season (SOS) detections were within 1 day of the SOS detection within a VIIRS pixel of 500 m spatial resolution. In contrast, less than 20% of the 30 m SOS detections were within 1 day if the land cover within the 500 m pixel was heterogeneous. They further reported that the SOS detection timing at the coarser spatial resolution was controlled by the timing when roughly 30% of the 30 m pixels within the 500 m pixel transitioned to SOS [46]. Thus, phenophase detections at coarser spatial resolutions are biased toward the earlier components in the vegetated land surface if the target landscape exhibits spatial heterogeneity. The study shows that the scaling effect on LSP is not resolvable with simple averaging. The mountain pastures of the Kyrgyz Republic and elsewhere in Central Asia exhibit a higher spatial heterogeneity than the croplands of central Iowa examined in [46]. Using a nonstandard MODIS product at 250 m spatial resolution, a study of pasture degradation in western Kyrgyz Republic found that pastures with a higher abundance of non-palatable vegetation exhibited later timing of peak NDVI during dry years in sub-alpine and mountain steppe ecozones but negligible impact in normal years [8]. The study is notable in that it included an extensive field component and demonstrated that higher NDVI values can accompany pasture degradation.

Several prior studies that have used MODIS data to detect trends in browning or greening in Kyrgyzstan and elsewhere in montane Central Asia may have been biased towards detecting browning trends and interpreting them as evidence of degradation due to using Collection 5 (C5) of the MODIS products. The differences between C5 and Collection 6 (C6 aka V006) are pronounced due to the loss of sensitivity of the red channel in the Terra MODIS [47,48]. Trend comparisons between the two collections have shown that significant negative trends in C5 were no longer significant in C6, and many nonsignificant trends in C5 appeared now as significant positive trends in C6 [49,50]. Two earlier studies [4,14] finding degradation in Kyrgyz pastures clearly used C5 products. Another study finding degradation [8] used data that was unlikely to have included the C6 corrections. A recent study [7] does not state which collection was used, but as the analysis period extended only until 2014, it is likely that it used C5 instead of C6 as well. The key point here is that as remote sensing products undergo upgrades that significantly change observed patterns, it is important to revisit those studies relying on earlier product versions to re-evaluate their findings in light of the new information.

A further concern is the use of simple linear slopes to detect significant trends. Applying ordinary least squares regression to a time series to fit a slope and declare the slope to be the trend has serious limitations from a statistical standpoint [51,52]. Positive autocorrelation reduces residual variance and inflates significance, increasing the risk of a Type I inferential error. Alternatives exist, such as the nonparametric Seasonal Kendall trend test, but it is not resistant to strong interannual autocorrelation. Fitting time series with autoregressive models has not been common in remote sensing studies.

Our findings suggest potential degradation in pasture areas as indicated by lower PHs closer to villages, but our analysis cannot rule out functional degradation arising from an increase in the coverage of nonpalatable species that can enhance NDVI without providing accessible forage [8,13,19]. Furthermore, what constitutes degradation in a particular socio-ecological system is frequently more subjective and culture-bound than is typically acknowledged [13,53,54,55].

5. Conclusions

We analyzed the landscape patterns of phenometrics based on fitted parameter coefficients of the land surface phenology model applied to 17 years of Landsat and MODIS seasonal time series across the montane pastures of the Kyrgyz Republic. We found that the Peak Height of NDVI generally decreased closer to villages, but the patterns were modulated—sometimes strongly—by elevation, aspect, and growing season weather. These findings raise questions about reports of pasture degradation based on coarser spatial resolution image time series. Our goal here was not to differentiate the relative contributions of these factors, because they are not susceptible to linear “unmixing”. Rather, we sought to recognize their potential influence on the mixed signal observed at coarser spatial resolutions and to urge caution in interpreting, for instance, declines in NDVI trends with pasture degradation and increases with pasture remediation. The situation is more complicated. Due to the spatial heterogeneous distribution of pastures and pasture usage in mountainous landscapes, contextual information should be used to interpret remotely sensed patterns and trends in an appropriately nuanced manner.