Using Landsat and Sentinel Imagery to Monitor Alfalfa and Grass Pastures in Western Colorado

Crop phenology of alfalfa and grass pasture fields on the western slope of Colorado, USA, was tracked using an NDVI time series derived from Landsat 7 (ETM+), Landsat 8 (OLI) and Sentinel 2a MSI to obtain high temporal frequency information of crop growth stages. These perennial pastures, which occupy a majority of irrigated land in the region and typically have more than one growing/cutting cycle, are quite diverse in terms of agronomic conditions, and their growth stages do not match average reported, conventional values documented almost two decades ago. Divergence of crop growth from these conventional published values have been recognized for other crop types and regions in the past, and refinement of crop growth calendars to reflect actual regional patterns encouraged. Regional practitioners often turn to the conventional published values for guidance in crop and water management, however we document here that the conventional published values are not applicable for this study region. Besides the climatic differences due to latitudinal and elevation characteristics of the region, anthropogenic interventions like multiple cuttings or grazing in a growing season progressively compound the effect of divergence from the conventional published values in the growing season. Using the data available, we summarize regional crop growth stage lengths, which may be used as a reference in absence of real time remote sensing or other ground information. This study serves as an example of the application of a combination of freely available Landsat and Sentinel satellite data to monitor crop calendars in a region of interest, in order to enhance crop management. The results have an associated uncertainty due to limited temporal coverage of cloud-free satellite images, and will benefit from further improvements such as using increased sites, more years of data, and additional satellites (e.g., addition of Sentinel 2b) for increased temporal frequency.


Introduction
Monitoring the phenological dynamics and growth of vegetation is important for understanding climate-biosphere interactions, estimating net primary productivity, managing natural and agricultural ecosystems, assessing accurate yields, and quantifying water use (Gao et al., 2017;Kimball et al., 2004;Zhang et al., 2003;Duchemin et al., 2006). Remote sensing offers opportunities for spatial monitoring of vegetation growth and development over large areas in a consistent, standardized manner. Evidence of strong similarity between satellite and ground phenology has been observed in the past (Reed at al., 1994;Gao et al., 2017). In recent years, satellite vegetation index data, such as the Normalized Difference Vegetation Index (NDVI) developed from the Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) at spatial resolutions of 1 km (both AVHRR and MODIS), 500 m (MODIS), and 250 m (MODIS), have been utilized to distinguish vegetation phenology and its changes (e.g., Xin et al., 2002;Zhang et al., 2003;). However, these resolutions are too coarse for field-scale detection of crop phenology in many parts of the world (Gao et al., 2017).
The Landsat suite of satellites provides multispectral data at 30 m resolution, but their low temporal resolution (16 days for each platform), which is further deteriorated by clouds and shadows, is usually not enough to monitor rapidly developing agricultural crops, whose phenology can be complex, especially when the crop exhibits multiple growth cycles in a single climatic growing season (e.g., perennial alfalfa and grass hay/pastures). The recent launch of the Sentinel-2 series of satellites, which match, and have been calibrated to, Landsat 8 characteristics, improve the temporal resolution of these sensors for monitoring high-frequency changes in vegetation. The minimal time gap between their overpasses reduces radiometric differences due to atmospheric and illumination conditions (Mandanici and Bitelli, 2016). Although some co-registration errors (~6 m (Yan et al., 2016); between the Landsat 8 OLI (30 m multispectral resolution) and Sentinel-2A MSI (10 m resolution for Red and Near Infra Red (NIR) bands) were observed by Yan et al. (2016), the effect can be reduced by downsampling and spatial averaging. Korhonen et al. (2017) observed highly accurate registration of images from both satellites. Linear correlations tested preliminarily between different bands of Landsat 8 and Sentinel-2a MSI were generally good (Pearson coefficient >0.99 for Red and >0.985 for NIR), as long as spatial heterogeneity is avoided (Mandanici and Bitelli, 2016). Flood, 2017 in their evaluation over Australia concluded that although further improvements to exactly match the two satellites may be required, the errors computed for biophysical quantities (NDVI RMSE of 0.036 and fractional cover RMSE of 0.040-0.054) are relatively low (the errors may vary with region). The shorter revisit period when Landsat 7 and 8 and Sentinel-2a satellites are used can significantly enhance cropland monitoring, which was previously not successful using only Landsat (Navarro et al., 2016).
For agricultural crops, it is a common practice to employ crop growth stages defined by the Food and Agriculture Organization (FAO) paper 56 (FAO56, Allen et al., 1998), where the growing period is divided into four distinct stages in order to manage the crop. The four stages are: 1) Initial (in perennials: from 'greenup date' when initiation of new leaves occurs to approximately 10% green cover); 2) Development (10% green cover to effective full cover (fractional cover=0.7-0.8)); 3) Mid-season (maximum growth, relatively constant growth conditions); and 4) Late season (maturity to harvest/end of season). Using the information from these crop growth stages and assuming linearity for the four phenological stages, a crop coefficient curve is constructed for "standard crop conditions" which refers to excellent growing conditions where no limitations due to water or nutrient shortage, disease, insect, weed, or salinity incidence are prevalent on crop growth or evapotranspiration (ET). The curve is constructed using time (length of growth season) on the x-axis, and the corresponding crop coefficients (Kc) on the y-axis (see figure 25, pg 100 of FAO56). The first and third stages have a slope of zero, while the second stage has a constant increasing and the last stage a constant decreasing slope. The linearity is an approximate average trend over time and has proven to be within acceptable uncertainty (Pereira et al., 2015). The information is used by crop and water managers to determine appropriate times for intervening in the crop growth with specific management techniques, such as irrigation; and diagnosing seasonal water use and yield.
The FAO56 crop coefficient curve is a simple crop-specific model but is accurate for both research and practice (Pereira et al., 2015). It has been adopted into various computer operational management models including Cropwat (Martin, 1992) and Aquacrop ) to be automated for applications like irrigation scheduling and yield assessment (Vashisht and Satpute, 2015;Todorovic et al., 2009). The authors of the original paper (Allen et al., 1998) as well as others have pointed out that the general tabulated crop growth calendar represents average standard conditions and encourage the use of local and regional observations of crop development including the use of tools like remote sensing . Researchers have made modifications to crop coefficient curves for both crop growth calendars as well as crop coefficient values for different crops, crop varieties, climates and management practices (e.g. Tasumi et al., 2005;Er-Raki et al., 2008;Howell et al., 2015). Tasumi et al. (2005) advised that satellite data of regional crop populations should be used for improving crop coefficient curves. Since this simple crop model will continue to be used extensively by agricultural producers, water managers, modelers and policy makers (Pereira et al., 2015), it is necessary to investigate the validity of its important aspects using actual observations in local and regional settings and possibly make refinements.
The objective of this paper is to detect alfalfa and grass hay/pasture crops' phenological development at a field-scale using the enhanced temporal frequency obtained by drawing on data from Landsat 7 ETM+, Landsat 8 OLI and Sentinel-2A MSI satellites/sensors, and the adherence or deviation of this crop development from generalized values commonly used in the agroclimatic setting of Western Colorado, USA. Demonstrating the utility of integrating these three data sources is important because information derived from this data can be immediately available and integrated into yearly cropland management practices. This paper also serves as an example of how data from the Landsat 7, Landsat 8 and Sentinel-2A, as well as the new Sentinel-2B, satellites can be integrated for multi-cycle crop phenology detection as well as future validation or further refinement of crop progress calendars.

Study Area
The study area comprises the Uncompahgre region of the Western Slope of Colorado, USA (figure 1). Situated west of the continental divide in a semi-arid climate, it spans approximately 1,000 square kilometers (100,000 ha). The elevation ranges from about 1500 m to about 2000 m. Climatology is summarized here: the average annual precipitation is 72 cm (22 cm without snowfall), and the average growing season precipitation is 15 cm; the average maximum temperature occurring in July-August months is 33˚C, and the average minimum temperature occurring in December-January months is -6˚C (source: www.usclimatedata.com). The growing season usually lasts from early April to late October. Alfalfa and grass occupy 90% of the irrigated land and consume 92% of the total irrigation water in the western Colorado (MWH, 2012). Being hay/pasture crops that are cut or grazed periodically, both crops have more than one growth cycle in a growing season. Lastly, surface irrigation is the most prevalent form of irrigation, limiting the size of agricultural parcels.

Data
Both Landsat 7 and 8 imagery (atmospherically corrected reflectance products and spectral indices, downloaded from https://espa.cr.usgs.gov/) was utilized to decrease long time gaps due to cloud cover and to increase temporal frequency. Even though Landsat 7 ETM+ scenes have gaps (due to the failure of the Scan Line Corrector (SLC)) that cover a part of the scene (path 35/row 33), the study area was mostly free of any gaps. Sentinel-2a Multispectral Instrument (MSI) imagery available in 2016, was downloaded from Earthexplorer (www.earthexplorer.usgs.gov). All cloud-free images utilized are listed in Table 1. Crop ground truth information was collected from field visits. Ground truthing showed that the crop cover classification layer from United States Department of Agriculture's National Agricultural Statistics Service (USDA NASS) Cropscape (http://nassgeodata.gmu.edu/CropScape/ ) misclassifies grass as alfalfa or other hay/pastures, therefore it was not used. Only the crop ground truth data with a total of 26 study sites, including 16 grass and 10 alfalfa parcels spread out spatially in the study area, was utilized. June 26 June 20 (L 7) July 9 June 28 (L 8) July 16 July 6 (L 7) July 29 July 14 (L 8)

Vegetation Index (VI)
NDVI was chosen because of its simplicity and low sensitivity to non-uniform, topographically varying terrain. While indices like SAVI and EVI may be less sensitive to soil background effects, they have been found to be more sensitive to topographic effects as opposed to indices with a band ratio format (e.g., the NDVI) (Matsushita, 2007). NDVI has been utilized in the past to depict crop phenology. It has been shown to be directly linked to crop coefficients (Senay, 2008) and was able to explain the majority of deviance in crop coefficients for both alfalfa and grass pastures (Vashisht et al., submitted), thus befitting its use for tracing crop growth stages for a crop coefficient curve. Field boundaries were buffered in to avoid any edge effects, and average field as well as within-field standard deviation NDVI values were calculated for each field for each satellite overpass. Given the high number of growing cycles (also called cutting cycles in some literature) and limited temporal frequency of imagery available for each cycle, the unequally spaced NDVI time series was used with its original values, without adjustment to a pre-defined mathematical curve approximation as has been done for some previous studies (e.g. Zhang et al., 2003). This helped in retaining natural variations of multiple growing cycles caused by various management practices, without losing any crucial information from 'strict fitting'.

Selection and profiling
Fields with high within-field standard deviation (>0.05) for most of the growing season were initially removed. Homogeneity requirement is especially essential when two different spatial scales (here 30 m Landsat and 10 m Sentinel-2a MSI) are used (Zhang et al., 2017). Spatial averaging for homogenous fields minimizes effects of co-registration errors between Landsat and Sentinel satellites (Korhonen et al., 2017). Heterogeneity within fields could be caused due to grazing by animals on a part of the field, improper water or nutrient application etc. Next, through a scrutiny of growing season NDVI variations, fields that did not show significant variation over time or had non-optimum (or stressed) growing conditions were also removed. These fields possibly had extensive grazing (animals grazing throughout the growing season) or were under rather poor (or stressed) management conditions. This left a total of 17 nearly "standard (without any shortage of water or nutrients) growing condition" study sites for temporal investigation of crop growth stages, including 10 grass and 7 alfalfa fields. It is noted that initial filtering already highlights the diversity of growing patterns and management conditions; but we chose nearly optimal conditions as in FAO56 to track "standard" crop growth stages in the region. This also highlights the importance of utilizing fine spatial resolution imagery because of high occurrence of sub-field variations, perhaps due to uneven application of water/nutrients/fertilizers and management conditions.

Crop growth stages
Crop growth patterns were tracked with continuous NDVI time-series curves. These curves are similar to crop coefficient curves, except that the absolute values of the y-axis are not crop coefficients, but a function of the crop coefficients. The reader is referred to figure 25 (pg 100) of FAO56 for a basic visualization of the crop coefficient curve, showing various growth stages: initial, development, mid and late as approximated by linear segments. Since our goal was to obtain crop growth calendars, we focused on the y-axis of the time-series curves. The temporal crop development patterns were then qualitatively compared with FAO56 values (from table 11 for Idaho, USA; FAO56). For grass, FAO guidelines present general information about a cutting cycle without any specifications for different (individual) cutting cycles. The tabulated values are 10 and 20 days for initial and development growth stages, respectively (no specification for mid and late season stages, because they may vary with management conditions). Unlike grass pastures, FAO56 separates rotated grazing alfalfa pastures (extensive grazing, as opposed to rotated grazing sites were few and removed in filtering) into first cutting cycle and other cutting cycles. The number of days in initial, development, mid, and late stages are 10, 30, 25, and 10 respectively for the first cutting cycle; and 5, 20, 10, and 10 respectively for the other cutting cycles.
The FAO general guidelines are extremely valuable and adequate for a detailed document with several contributions, but need further expansion and incorporation of regional information from the study area. To obtain crop growth stage lengths (CGSLs), pinpointing the key transitions from one stage to the next is a non-trivial task. This is true especially because of high spatial diversity in the study region, even with a limited number of ground-truthed sites. Similar to Tasumi et al. (2005), who adjusted emergence, cover and termination dates for crop coefficients for multiple single cycle crops in Idaho, we identified the beginning of the four stages in order to determine average growth stage calendars for each cutting cycle under different management conditions of grass and alfalfa pastures. In perennials, since there is no planting date, the start of the season (initial stage for first cycle or greenup date) for the very first cycle was derived from temperature thresholds mentioned in FAO56 (which, for grass, were originally identified from a report based on high altitude grass pastures in Western Colorado itself; Kruse and Haise, 1974). For the rest of the growing (or cutting) cycles, the initial stage was only considered if there was a time period with negligible slope right after a late growth stage minima. The start of development stage (or the end of initial stage) was approximated as the date right before significant increasing growth was perceptible via slope change. Effective cover date, a.k.a. start of midseason (or the end of development stage) was taken to be the time when the slope of the curve starts to flatten out roughly. The start of late season (or the end of development stage) was approximated similarly to the start of development stage, but with a decreasing slope after a local maxima. The CGSLs were finally determined from approximate days between the transition dates of successive growth stages, separately for varied agronomic practices of grass and alfalfa.

Grass
The 2016 growing season for grass pastures was estimated to begin on April 1, i.e., 7 days before the last -4 ˚C day in spring (according to FAO56-Allen et al., 1998;Kruse and Haise, 1974).
Considerable variation among the fields was observed, but a few distinct patterns can be detected. To support easy visualization, the fields were separated into three categories, as defined by the total number of growing (or cutting) cycles observed in a growing season. Crop temporal growth from fields with 2, 3 and 4 cutting cycles are shown in figures 4, 5 and 6, respectively. "Idealized" FAO56 values from the tabulations are graphed in figure 7, without any mid and late season stages (because of absence of information on them). It is noted that unlike figures 7-9, the y-axis of figure 10 is the FAO56 crop coefficient, but for the purpose of assessing crop growth stages we are concerned with x-axis relative to the pattern (e.g. slope) of y-axis, not its absolute value. It is seen that not only is there a noticeable deviation of actual growth cycles and CGSLs from FAO-tabulations, but also a significant difference amongst the individual fields in the study area. In general, the length of a given crop growth stage for a given cutting cycle depends on the total number of cycles in the growing season and where that cycle lies in the growing season. The first development stage is similar for the three categories (but twice as long as tabulated values), while the mid and late growth stages vary widely because they are determined by management decisions, such as when to (and how much to) cut/graze the field, and that in turn determines the length of the next cycle(s).

Alfalfa
A graphical visualization of variation in population of fields is given in figures 8, 9, and 10. Among the 7 alfalfa plots tracked, a variety of management conditions were observed. For easy visualization, the fields were divided into three categories, i.e., 3, 4 and unconventional cutting cycles (unconventional being the one that seemed not to follow any number of cycles). Again, it is seen that not only is there a noticeable deviation of actual growth cycles and individual growth stage lengths from FAO tabulations, but also a significant difference between fields. Unlike FAO56, the initial stage was not observed in any growth cycle, but a gradual green-up period (NDVI increase) was observed in March. Like grass, the first development stage is similar for the first two categories (unconventional category is not discussed hereafter), while the mid and late growth stages vary widely because they are determined by management decisions like when to (and how much to) cut/graze the field, and they in turn determine the length of the next cycle(s). From this it is clear that an implementation of FAO56 would not be successful. As is known from direct field experience, four cutting cycles dominate the region. Tasumi et al. (2005) for Kimberly, Idaho, also observed the majority of alfalfa fields with four cycles, and some with three cycles. They also observed that the NDVI distribution showed times/dates being relatively uniform/constant in space (among or across fields) for the first cutting, but uniformity decreased in subsequent cycles due to increasingly compounding effect of different cutting dates and management.

CGSLs
A direct implication from the results above is that utilization of FAO56 tabulated CGSLs is not valid for Western Slope of Colorado. An indirect implication is that any tabulation or standardization of CGSLs is difficult because of the diversity in management conditions. Nonetheless, an attempt is made to summarize the CGSLs for each category (i.e., number of cycles) with the limited data available (table  10 for grass; table 11 for alfalfa. This is because for practical purposes like operational evapotranspiration and yield assessment, even an approximate knowledge of the regional growing season calendar might be better than using conventional values that are not valid. The growing season is taken to be from April to November, although a green-up period (NDVI increase) is seen in the month of March and some vegetation decrease is discernable in early November. Because a clear distinction could not be drawn between the initial and development stages for cutting cycles beyond the first cutting cycle, these two growth stages were combined to determine CGSLs. Similarly, the mid and late stages were combined. It is noted that these CGSLs should only be used as a guide if remote sensing data (preferably for near-time phenology tracking) or local crop cutting and management data is not available for purposes like estimation of crop evapotranspiration and yield.   Mid+Late=40-50 -Crop growth rate is strongly impacted by weather variables including mean air temperature, solar radiation, and relative humidity and their spatial variations manifested from latitudinal and topographic differences. Crop type and variety also has a role to play in diversity of management in the region, given that there are various types of cool-season grasses and sometimes even mixtures of grass and alfalfa grown together. Bounds for each phenological stage are somewhat broad, which is not uncommon. Although we do not have ground estimates of vegetation for validation, previous studies have found uncertainties that range similar to the broad approximations here. Sakamoto et al. (2005) found that the error of remotely sensed (from MODIS) growing period and each phenological stage for rice was about 12 days, which they considered adequate for discrimination of phenology in absence of other better information (albeit not sufficient for use in crop yield models). Average 3-year root mean square error for You et al. (2013) was approximately 17 days for both start and end of season estimation. Tasumi and Allen (2007) used 12 images in the growing season and determined midseason lengths, which had broad width (e.g. 10-15 days for bean, 45-50 for corn, 45-60 for sugar beet), and stated that true averages for a region are better than following inaccurate mean guidelines.

Limitations and future work
The diversity of field conditions encountered with the limited ground truth data calls for future validation. We assumed the start of season according to pre-defined thresholds, but it is possible that this may have changed due to increasing temperatures as climate change is recognized to cause early vegetation greening. As Tasumi and Allen (2007) also reported, accuracy of stage transition dates depends on temporal sampling, thus the results will have an associated error determined by temporal coverage of cloud-free satellite images. It is acknowledged that Landsat and Sentinel are different satellites with different overpass times and sensor characteristics, and there is scope to improve their combined application vis-a-vis co-registration and radiometric normalization (Korhonen et al., 2017). Further investigation should also include the upcoming Sentinel 2b satellite in order to decrease the uncertainty and improve the crop growth and management calendar for better farm and resource monitoring, decision-making and management.

Conclusions
Crop development of ground truthed grass and alfalfa pasture fields was tracked using NDVI time series from Landsat 7 ETM+, Landsat 8 OLI and Sentinel 2a MSI in order to assess the crop calendar. Through this, we investigated an important aspect of the FAO56 crop coefficient model that is used by managers to guide application in water use (through reference evapotranspiration) and yield estimation. The Crop coefficient model presented in FAO56 consists of four simplified categories of crop growth stages and their respective crop coefficient values. Regional remotely sensed data reveals marked differences of crop growth and length of various growth stages, with respect to the conventional standard FAO56 guidelines. This is because it has contributions from both generic, average crop agronomy (which the standard guidelines could be more indicative of) as well as regional and local management conditions. Even though the farm practitioners turn to these conventional guidelines in the absence of any other information source, they are not suggested for future use.
Our findings show that for both grass and alfalfa, the length of a given crop growth stage for a given cutting cycle depends on the total number of cycles in the growing season and where that cycle lies in the growing season. For grass, the development stage of the first cycle is found to be twice the tabulated values, while the mid and late growth stages vary widely because they are more likely to be determined by management decisions. For alfalfa, the first cycle development and initial stages are somewhat similar to FAO56, but the stages for latter cycles are longer. Like grass, mid and late stage lengths of neither the first nor consequent cycles match the FAO56 stages.
These results are not surprising as previous regional studies in other regions have found similar results and led to modifications to crop coefficient curves in those regions. We summarize crop growth stage lengths for various categorizable conditions for both grass and alfalfa pastures in the study region and suggest that these can be used as a guide in the absence of any real time remote sensing or ground information. The results have an associated uncertainty due to limited temporal coverage of cloud-free satellite images. Notably, the diversity of field conditions encountered with the limited ground truth data calls for future validation. Future work should facilitate an improvement of the combined application of Landsat and Sentinel satellites, including the upcoming Sentinel 2B to decrease uncertainty.

Acknowledgements
We appreciate Dr. Perry Cabot and Carter Stoudt for providing ground truthed data. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.