Mapping When and Where Invasive Buffelgrass is Green at Saguaro National Park in Arizona.
Our research detects when and where invasive buffelgrass (Pennisetum ciliare) is photosynthetically active in the Tucson, Arizona region. The methods and results can alert land managers in a timely manner so they can implement condition-based management and treat the buffelgrass during optimal conditions. USGS scientists mapped the presence and phenological status of the invasive by integrating ground observations of buffelgrass phenology collected by citizen scientists and professionals, a map of buffelgrass distribution in Saguaro National Park (SNP), climate information, and Moderate-resolution Imaging Spectroradiometer (MODIS) imagery. Knowing where buffelgrass is located and when it is green enables the strategic application of herbicide, which is only effective when the plants are photosynthetically active and are at least 50% green. Climate Landscape Response (CLaRe) phenometrics, developed for this project, capture the strength of the landscape greenness response to climate and expose buffelgrass due to its more rapid response to precipitation events than native vegetation, i.e., areas with buffelgrass display strong correlations between MODIS greenness and lagged precipitation. These metrics reveal that precipitation data can predict green-up, whereas MODIS data can confirm green-up. The research was conducted to assist buffelgrass control activities of the National Park Service at SNP, but the results are expected to transfer easily to other invasive species because many display similar rapid response to environmental cues.
Buffelgrass is a perennial grass that is invasive to the Sonoran Desert. It was introduced to North America from Africa in the 1940s as cattle forage because it can thrive in arid environments. In contrast to the native Sonoran Desert vegetation, which typically consists of spatially dispersed individual plants, buffelgrass forms a continuous mat of highly flammable vegetation that carries fire quickly and broadly across the landscape. Buffelgrass has proven well adapted to the Sonoran Desert, threatening desert ecosystems by out-competing native plants and altering fire regimes. Repeated exposure of non-fire adapted native plants to fires fueled by buffelgrass ultimately transforms the Sonoran Desert ecosystem from a diverse assemblage of plants that supports a rich diversity of fauna to a desert grassland monoculture that burns regularly, threatening lives and property.

Buffelgrass was first observed in Saguaro National Park (SNP) in 1989; within 20 years the invasive grass covered an estimated 800 hectares, despite the mobilization of substantial park and volunteer resources to manually remove plants. The most recent estimates of the buffelgrass extent are 2000 acres. SNP managers consider the complete eradication of buffelgrass in the park to be effectively impossible, but aggressively work to constrain further encroachment given the plant’s destructive potential. Infestations in rugged, remote areas of SNP that cannot be readily reached by field crews are treated with aerial applications of herbicide. The timing of application is critical since the herbicide is most effective while the plant is actively growing and is at least 50% green. However, the appropriate timing can be difficult to determine given the heterogeneity of the landscape and variability of monsoon-driven greenness in the region. For example, in 2015 aerial spraying in SNP’s western Tucson Mountain District was cancelled because the buffelgrass in that area had already entered dormancy. To minimize the impact on surrounding native vegetation and ensure the optimal effectiveness of herbicide applications, park managers need to know not only where outlying buffelgrass patches are, but when they are most photosynthetically active.
The goal of our research was to map WHERE buffelgrass is located and detect WHEN buffelgrass is green. This will enable the strategic application of buffelgrass herbicide, which is most effective when vegetation is photosynthetically active, during optimal treatment windows.

Satellite Data
To compile landscape-wide indications of vegetation greenness, we downloaded all available 8-day composite Moderate-resolution Imaging Spectroradiometer (MODIS) reflectance images (46 per year) for 2011, 2012, and 2013. MODIS products are compiled at a range of spatial resolutions; we accessed the highest resolution dataset (250 m). The composite dataset is produced by selecting the pixel obtained during the eight-day period that contains the best possible observation as evaluated on the basis of favorable coverage and view angle, the absence of cloud contamination, and aerosol loading, resulting in data that are pre-filtered for quality. We retrieved MODIS bands 1 (red; 620–670 nm) and 2 (near infrared (NIR); 841–876 nm) and created Normalized Difference Vegetation Index (NDVI) images from the normalized ratio of the near-infrared (NIR) and red bands, as follows:
Vegetation exhibits a distinctive spectral reflectance pattern, with low reflectance of red wavelengths due to chlorophyll absorption and strong reflectance in the near-infrared (NIR) due to cellular structure. The NDVI exploits this distinctive signature to provide a measure of the photosynthetic activity of plants (Figure 3).

Citizen Science Observations-Pima Canyon
We accessed a record of high-frequency buffelgrass observations taken at a single location near Tucson (in lower Pima Canyon) for this study. Pima Canyon is located in the Santa Catalina foothills on the northern edge of the city (Figure 2) in geomorphic conditions similar to those of Saguaro National Park–East; i.e., mountainous terrain composed of felsic granitic and metamorphic rocks. The south-facing slope of the Santa Catalina Mountains, which includes the Pima Canyon area, has a particularly dense and widespread buffelgrass infestation. The observations were collected by a single observer on a near-weekly basis from July 2010 to December 2014 using standardized protocols; data were archived in the USA National Phenology Network (USA-NPN) database. The visual estimates of percent greenness of up to five individual buffelgrass plants located along a south-facing slope of the trail were assessed on each day of observation (Figure 2).
Citizen Science Observations-Tucson
Beginning in June 2013, we recruited 10 citizen scientist volunteers to participate in the USA-NPN project Nature’s Notebook Citizen and Professional Science project for buffelgrass observations in the Tucson region. The volunteers were trained through the USA-NPN Education Coordinator and the Southern Arizona Buffelgrass Coordination Center (SABCC) as part of an effort to engage the general public in long-term phenology monitoring. Those interested in collecting data were encouraged to select a buffelgrass patch or individual plant to follow. All proposed locations were field checked and on-site training was provided. As at the Pima Canyon site, observations followed a series of strict protocols designed to capture detailed information regarding the timing and proportion of plant green-up, flowering, fruiting, and senescence. More information about the project can be found on the USA-NPN Tucson Phenology Trail website. By the end of 2014, 6913 observation records and 896 site visits had been entered into the database by the volunteer participants. The broad geographical distribution of these observations provided a complementary dataset to the single-site, long-term records from Pima Canyon.
Map of Buffelgrass Extent and Density
The National Park Service supplied a map of buffelgrass locations and densities within both SNP districts produced from a September 2012 helicopter survey conducted to document the distribution and quantity of buffelgrass, estimate its growth and spread, and prioritize areas for treatment. Each area was tagged with a visual description of the assessed buffelgrass density: trace (1% to 25%), low (>25% to 50%), medium (>50% to 75%), or high (>75% to 100%). To transform the polygonal coverages to the raster-based format of the remote sensing data, we calculated the percent of buffelgrass within each MODIS pixel by assigning the midpoint value to each polygon for low to medium density values (i.e., 12.5%, 37.5% and 62.5% for trace, low and medium respectively) and 80% to the high density category. We then converted the polygon coverage to a 5-m raster attributed with the percent values, and averaged all pixel values within each MODIS footprint (Figure 5).

Precipitation Data
Two different sources of rainfall data were used for analyses conducted at a single site (Pima Canyon) and those conducted over a landscape extent (SNP East and West). For Pima Canyon we accessed the closest source of recorded precipitation information, which was MesoWest daily data ([40]) at station AT021 (32.33050°N, 110.9120°W). The station is within 4 km and at a similar elevation (890 m) as the observation site. For the landscape comparison, we downloaded PRISM (“Parameter-elevation Regressions on Independent Slopes Model”) daily rainfall data ([41]) over the SNP districts. Using nearest neighbor techniques, we resampled the 4 km PRISM data to 250 m to conform to the spatial resolution of the MODIS data. For both datasets, we summed daily precipitation values over eight-day intervals to coincide with the MODIS composite data
Methods
Relationships between field-observed greenness, MODIS satellite greenness and precipitation data were explored and characterized at the Pima Canyon site, given the temporal depth of these data. Simple linear regression analyses and Spearman’s correlation coefficients were used to quantify relationships among the variables. On the basis of the correlations determined at the Pima Canyon site, we created a novel suite of metrics (“Climate Landscape Response” (CLaRe) metrics) to predict buffelgrass locations in SNP from park-wide satellite NDVI measurements and gridded 4 km precipitation data. These metrics range from 0 to 1.0 and quantify the strength of the linear relationship between the satellite greenness value and various cumulative precipitation values. They are created by conducting a linear regression between the NDVI and a precipitation value at each pixel and extracting the correlation coefficient.
Since we wanted to use the 2012 helicopter mapping as validation, we restricted the calculation of metrics to the 2011–2013 time frame, thereby bracketing the year the “ground truth” mapping data were generated. Statistical tests, including Student’s t-tests, were conducted to evaluate if the metrics were effective at separating pixels of pure native vegetation from those with buffelgrass. The CLaRe metrics with a balance between spatial and temporal prediction were used to create models of buffelgrass presence/absence in the park. That is, we used metrics that spatially displayed strong overall correlations and separability, but temporally allowed for earlier prediction of buffelgrass presence and condition. Buffelgrass presence/absence models were validated by examining their spatial coincidence with confirmed buffelgrass locations.
The linear regression analyses with strong (>0.50) correlations between greenness observed on the ground and various cumulative and lagged precipitation values were solved for the desired 50 percent greenness condition, which is the minimum greenness for which herbicide application can be effective at killing the invasive grass.
Results
Pima Canyon
Strong relationships exist at Pima Canyon between buffelgrass greenness observed on the ground and both MODIS greenness and PRISM precipitation. These relationships are visually evident in Figure 6, which plots the variables compiled for eight-day intervals and uses MODIS NDVI greenness. In particular, the field observations of buffelgrass greenness and those extracted from MODIS data are temporally synchronous, while buffelgrass greenness distinctly lags precipitation events. Figure 7 explores these relationships at Pima Canyon further, and graphs the correlation coefficients for the linear relationships between a sequence of cumulative precipitation values and both eight-day MODIS NDVI values and field observed greenness values, revealing strong (>0.5) correlation for all cumulative precipitation totals except the current (ppt0) and the prior time period (ppt1).


Saguaro National Park (SNP)
Several CLaRe metrics were developed from correlations quantified at the Pima Canyon site and mapped across the SNP study area. These metrics capture the strength of the relationships between satellite-measured greenness and various lagged precipitation values. At each pixel for each year, the correlation coefficient describing the strength of the linear relationship between the 46 MODIS NDVI values and the cumulative precipitation for time periods of different lags was calculated. For example, the correlation between MODIS greenness and the cumulative precipitation for the prior three eight-day periods (ppt123) for 2011, 2012 and 2013 is shown in Figure 8. Overlaid on these figures are the polygons of buffelgrass mapped in 2012 via the helicopter survey. We observed an especially strong spatial coincidence of the buffelgrass patches and the high correlation values for all years using this CLaRe metric. Of particular note is that the pattern of high CLaRe values is the same in all three years, even though the climate varied, with 2011, 2012 and 2013 years of above, normal and below average monsoon precipitation, respectively.

To explore the response of native vegetation communities to precipitation, we isolated all pixels that had had no buffelgrass mapped in the 2012 aerial helicopter survey, and stratified them according to the major vegetation classes in the park; their precipitation/MODIS greenness correlations are presented in Figure 9. Several patterns are evident. In all cases but one, the native vegetation types show very low correlation to recent precipitation, with values less than 0.5 for current precipitation and up to three eight-day periods of cumulative lagged precipitation. In all cases, correlation strength increases through time, showing that native vegetation types are more increasingly correlated to longer-term cumulative precipitation. In addition, the CLaRe phenometrics effectively stratify the vegetation types by rooting depth, with deep-rooted woodland and forest showing the lowest correlation values and shallower-rooted grasslands and scrublands showing higher correlation.

We evaluated the ability of various CLaRe phenometrics to locate buffelgrass by stratifying the landscape first by vegetation type and then by the density of buffelgrass within the pixels (Figure 10). The map shows the dominant vegetation types in each MODIS pixel footprint as mapped by the Southwest Regional Gap Analysis Project (SWreGAP). Note that even small amounts of invasive buffelgrass in the landscape have a strong effect on the correlation coefficients in the CLaRe phenometrics.

As seen in previous plots, the curves of all pure native vegetation types (the green curves in Figure 10) are monotonically increasing, whereas several of the curves for high density buffelgrass (the brown curves in Figure 10) peak and then decrease as the number of lagged time intervals increases. This is the same pattern seen at Pima Canyon (Figure 7). The buffelgrass curves are distinctly different than the curves for native Sonoran Desert vegetation communities. Not only do they display consistently higher correlations to all variations of cumulative precipitation, revealing that the invasive responds quickly to precipitation, but also these curves often peak, reflecting the invasive’s aggressive growth pattern.
Mapping Buffelgrass–Where?
Using the CLaRe metrics and the results for the individual vegetation types, we created several models of buffelgrass presence/absence. For modeling purposes, we focus on SNP-East, given its geologic and geomorphic similarity to the Pima Canyon site and its diversity of vegetation types and elevational gradients. We also focus on the shorter-term cumulative precipitation lags, given that we seek to provide operationally timely data and that these CLaRe metrics display strong differences between pure native and buffelgrass-invaded pixels. Our modeling was conducted by choosing options that seemed logical, useful and easily understood by practitioners; we did not apply statistical methods to optimize classification.
One of the best models we created is shown in Figure 11. This was developed to be easily created with no prior knowledge of buffelgrass distribution. Based on exploratory statistics and prior results, these models used the CLaRe ppt123 for each year. All CLaRe ppt123 values were extracted for each mapped vegetation type and the top 5th of the highest correlations were mapped as buffelgrass presence (i.e., the upper 20 percent of calculated correlation values). For modeling and validation purposes, we eliminated the higher elevation vegetation types that are unsuited for buffelgrass (e.g., Madrean pine-oak forest and woodland). The three sub-models show the results for each year (left panels, Figure 11). The composite model (right image, Figure 11) maps buffelgrass presence if 2 of 3 sub-models mapped buffelgrass presence.

Validation results suggest that the three sub-models do a reasonable job identifying known patches of buffelgrass, capturing 45 to 49 percent correctly. The composite models improve the accuracies: composite model 2 of 3 (shown in Figure 11) has an overall accuracy of 83 percent and correctly maps 46 percent of known buffelgrass. There are tradeoffs between overall accuracy and the total amount of buffelgrass correctly identified compared to the amount misclassified and missed. In all cases, the models show poor user’s accuracies (typically <25%) for mapping buffelgrass, with the model classifying areas as buffelgrass that were confirmed native vegetation in the field. Suggestions for improving the accuracies are presented in the discussion section below.
Mapping Buffelgrass—When?
Because optimal use of herbicide to control buffelgrass requires that the plant is photosynthetically active when treated, it is important to map not only where buffelgrass is located, but also when it is green. To predict the phenological status of buffelgrass, we analyzed the citizen science data and solved the linear equations relating field observed buffelgrass greenness to various cumulative precipitation totals for precipitation at specified amounts of buffelgrass greenness.
Results for the Santa Catalina Mountain sites are shown in Tables 1. To achieve the management threshold of greenness for effective herbicide application, results suggest buffelgrass will be 50% green 8 to 16 days after a 24 day period that totals over 46 mm (1.8 inches) of rain in the Santa Catalina Mountains (less is required in the Tucson Mountains, with different soils, etc.). In SNP-east, where we were unable to monitor buffelgrass on the ground, similar geology and geomorphology suggests that the thresholds observed in the Santa Catalina sites would be most appropriate to use. Of course, knowledge of these thresholds will evolve and become more accurate as practitioners interact with the landscape. In addition, the following section will show how MODIS data can be used to confirm green up in an area targeted for treatment.
Operational Protocol
The research results can be used to create a general protocol for operationally informing management activities related to herbicide treatment of buffelgrass at SNP and the Tucson region. Certain decisions, although produced from the data specific to this area, are likely to be generalizable. For example, the metric CLaRe ppt123 was used to model buffelgrass presence in the landscape for our region. Given the characteristic aggressive response of buffelgrass to precipitation, it is likely this metric will be useful, if not optimal, in other areas as well. Suggested steps for optimizing the timing of herbicide treatment for buffelgrass in the landscape are as follows:
- Map where buffelgrass is by using historic MODIS NDVI data and PRISM precipitation data to create CLaRe ppt123 metrics for several years of data. Include any model refinements in this step, as presented in the discussion section. Modeling results from this study suggest that an effective way to screen a region for buffelgrass presence is to extract CLaRe ppt123 metrics (i.e., cumulative precipitation for the three prior eight-day periods) for vulnerable vegetation types, and examine areas with the highest correlation metrics as candidate landscapes for hosting buffelgrass. In our study, good results were achieved by examining the top 5th (highest 20%) of the correlation values. The upper threshold selected can reflect expert knowledge of the extent of buffelgrass invasion in the management area, the size of the area that can be reasonably examined given resource constraints, and the relative vulnerability to invasion of landscapes of concern. Creating composite models from several years of data is shown to improve accuracy; composites are created using CLaRe ppt123 by mapping buffelgrass presence if half or more of the suite of individual annual models map buffelgrass.
- Select areas to target for treatment and collect the MODIS NDVI and PRISM precipitation data at that location. For example, eight-day MODIS (solid green line) and PRISM data (solid blue line) for 2011 at a pixel in SNP-East that contains a patch of high density buffelgrass is shown in Figure 12. Note that PRISM values are available daily and, as of this writing, eMODIS data, which is the composite NDVI value for the prior seven days, are also available daily. MODIS eight-day composite data used in this study are shown to compare very well with eMODIS data. This means that it is now possible to access for free PRISM 4-km precipitation data and 250-m eMODIS NDVI data every day that are current through yesterday.
- Predict when buffelgrass will be 50 percent green by calculating a running tally of the cumulative precipitation for the past 24 days (dashed magenta line). When the prior 24-day period totals over 46 mm (1.8 inches) of rain (vertical red bar), place crews and aircraft on notice for herbicide application 8 to 16 days from that time (shaded region of graph).
- At eight days prior to herbicide application, confirm that the target has greened up using current MODIS data (green line), adjust areas for treatment as needed.
On the day of treatment, confirm that the target has greened up using current MODIS data (open black circle), adjust areas for treatment as needed.

Discussion
The composite model for mapping buffelgrass (Figures 11) is ~80 percent accurate overall and effective at capturing ~50 percent of known buffelgrass patches, but has poor user’s accuracy (<26 percent). Inspection of the geography of misclassified pixels suggests environmental factors could be affecting the accuracy. In the model and sub-models, the producer’s error tends to occur in an east-west band near the southern boundary of the park whereas the user’s error tends to occur in the northwest (Figure 11). This observation suggests that model results may be improved if areas with distinct aspect are isolated before extracting the highest correlations for composite models.
Other factors may contribute to over-mapping buffelgrass patches and poor user’s accuracies: the coarseness of the PRISM precipitation data used (4 km), the accuracy and relevance of the vegetation map, and the presence of other invasive plants, such as Schismus and fountain grass (Cenchrus macrouru) that display similarly aggressive green-up behavior. The model may be improved using finer resolution precipitation data and more accurate maps of vegetation and invasive species. The SWreGAP map was produced using data collected in 2000 through 2002; an update may be beneficial. Although SNP is protected, vegetation changes may have resulted from fire and drought over the 10+ years before the 2012 helicopter mapping survey was conducted.
This study pooled annual MODIS NDVI composites and PRISM precipitation data to produce CLaRe metrics. Additional insights will likely be gleaned by not only splitting the data into a winter (cool) season and a summer (warm) season, but also by creating CLaRe metrics using temperature. During exploratory analyses, when observed buffelgrass greenness at the Pima Canyon site was plotted against daily temperature variables, a bimodal distribution was: Buffelgrass greened up in response to rainfall in both cool and warm seasons, unlike native C4 grasses that respond primarily to warm-season precipitation. The cool-season greenup of buffelgrass is likely predicated on reaching a certain threshold for temperature; such relationships can be readily explored and mapped spatially using CLaRe phenometrics.
This study demonstrates the usefulness of CLaRe phenometrics in a specific management application related to control of invasive buffelgrass, but they hold promise for many issues related to dryland vegetation. The CLaRe phenometrics used here explicitly incorporate precipitation and capture the strength of the correlation between the satellite greenness data and precipitation data (both lagged and cumulative). These phenometrics are effective in locating populations of invasive buffelgrass since the grass responds more rapidly and strongly to precipitation events than native vegetation. Even small amounts of invasive buffelgrass in the landscape have a strong effect on the values of these phenometrics. Given that many invasive species are successful because they can out-compete native species for limited resources, it is likely the CLaRe phenometrics will prove valuable for mapping other invasive species of regional extent in dryland ecosystems.
Conclusions
Buffelgrass is a noxious invasive plant that threatens to degrade the iconic Sonoran Desert ecosystems by out-competing native species and altering fire regimes. Effective treatment of buffelgrass using herbicide requires knowledge of where buffelgrass is located in the landscape and when it is at least 50 percent green. Optimizing herbicide treatment is crucial for many reasons, including cost, efficacy, potential collateral damage to the ecosystem, public opinion, and social concerns. Our research enables the targeted and strategic application of buffelgrass herbicide by using remote sensing data to detect when and where buffelgrass is photosynthetically active and predicting when it will be at least 50 percent green.
CLaRe metrics capture the aggressive response of buffelgrass to recent precipitation: buffelgrass-invaded pixels typically display statistically higher values (i.e., stronger correlations, at ρ < 0.05) than pixels of pure native vegetation for all vegetation types in SNP. Composite models created from CLaRe metrics produced validation results with overall accuracies of ~80 percent and correctly mapped ~50 percent of known buffelgrass patches. Results from the study suggest accuracy can be improved by further stratifying the landscape into large-scale topographic facets of similar aspect. Refinements may also be realized by focusing analyses on specific seasons (winter vs. monsoon) and by creating CLaRe phenometrics with temperature data in addition to precipitation, especially in the winter rainy season.
Citizen science observations at Pima Canyon were key in linking the greenness observed on the ground to the satellite-based MODIS-NDVI data. In addition, the citizen science observations from across the region allowed for the prediction of when buffelgrass will be at least 50 percent green and thus suitable for effective herbicide treatment. In the general Tucson Area, a conservative estimate is that buffelgrass will be at least 50 percent green 8 to 16 days after a 24-day period with cumulative rainfall totals of 46 mm (1.8 inches).
The phenometrics introduced here are a new and innovative suite of metrics that are proven effective in mapping when and where invasive buffelgrass is green in Saguaro National Park near Tucson, Arizona. The ability to monitor and interpret trends in land surface phenology informs science-based land management. In contrast to temperate regions, which respond to the annual solar cycle, producing predictable greenness curves readily captured by traditional phenometrics, drylands respond to precipitation, which can be erratic. As such, the “start of season” in drylands is more correctly seen to be the onset of a significant precipitation event, such as the monsoon. The CLaRe phenometrics capture the landscape response to climate and will complement the current temperate-zone phenometrics; they will prove broadly useful for many applications, including monitoring dryland ecosystems, mapping invasive species, habitat modeling, and vegetation classification.
Acknowledgments: We wish to thank the National Park Service and Dana Backer, in particular, for supplying us with data and critical feedback on the methods used here. Data were provided by the USA National Phenology Network and the many participants who contribute to its Nature’s Notebook program. The USGS (Land Change Science, Land Remote Sensing and National Park Monitoring Project within the Status & Trends Program) provided funding for this research. Jessica Walker was supported by a USGS Mendenhall Research Fellowship during this work. We appreciate the thoughtful and constructive comments from our three anonymous reviewers and USGS internal reviewer Joel Sankey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Mapping When and Where Invasive Buffelgrass is Green at Saguaro National Park in Arizona.
Our research detects when and where invasive buffelgrass (Pennisetum ciliare) is photosynthetically active in the Tucson, Arizona region. The methods and results can alert land managers in a timely manner so they can implement condition-based management and treat the buffelgrass during optimal conditions. USGS scientists mapped the presence and phenological status of the invasive by integrating ground observations of buffelgrass phenology collected by citizen scientists and professionals, a map of buffelgrass distribution in Saguaro National Park (SNP), climate information, and Moderate-resolution Imaging Spectroradiometer (MODIS) imagery. Knowing where buffelgrass is located and when it is green enables the strategic application of herbicide, which is only effective when the plants are photosynthetically active and are at least 50% green. Climate Landscape Response (CLaRe) phenometrics, developed for this project, capture the strength of the landscape greenness response to climate and expose buffelgrass due to its more rapid response to precipitation events than native vegetation, i.e., areas with buffelgrass display strong correlations between MODIS greenness and lagged precipitation. These metrics reveal that precipitation data can predict green-up, whereas MODIS data can confirm green-up. The research was conducted to assist buffelgrass control activities of the National Park Service at SNP, but the results are expected to transfer easily to other invasive species because many display similar rapid response to environmental cues.
Buffelgrass is a perennial grass that is invasive to the Sonoran Desert. It was introduced to North America from Africa in the 1940s as cattle forage because it can thrive in arid environments. In contrast to the native Sonoran Desert vegetation, which typically consists of spatially dispersed individual plants, buffelgrass forms a continuous mat of highly flammable vegetation that carries fire quickly and broadly across the landscape. Buffelgrass has proven well adapted to the Sonoran Desert, threatening desert ecosystems by out-competing native plants and altering fire regimes. Repeated exposure of non-fire adapted native plants to fires fueled by buffelgrass ultimately transforms the Sonoran Desert ecosystem from a diverse assemblage of plants that supports a rich diversity of fauna to a desert grassland monoculture that burns regularly, threatening lives and property.

Buffelgrass was first observed in Saguaro National Park (SNP) in 1989; within 20 years the invasive grass covered an estimated 800 hectares, despite the mobilization of substantial park and volunteer resources to manually remove plants. The most recent estimates of the buffelgrass extent are 2000 acres. SNP managers consider the complete eradication of buffelgrass in the park to be effectively impossible, but aggressively work to constrain further encroachment given the plant’s destructive potential. Infestations in rugged, remote areas of SNP that cannot be readily reached by field crews are treated with aerial applications of herbicide. The timing of application is critical since the herbicide is most effective while the plant is actively growing and is at least 50% green. However, the appropriate timing can be difficult to determine given the heterogeneity of the landscape and variability of monsoon-driven greenness in the region. For example, in 2015 aerial spraying in SNP’s western Tucson Mountain District was cancelled because the buffelgrass in that area had already entered dormancy. To minimize the impact on surrounding native vegetation and ensure the optimal effectiveness of herbicide applications, park managers need to know not only where outlying buffelgrass patches are, but when they are most photosynthetically active.
The goal of our research was to map WHERE buffelgrass is located and detect WHEN buffelgrass is green. This will enable the strategic application of buffelgrass herbicide, which is most effective when vegetation is photosynthetically active, during optimal treatment windows.

Satellite Data
To compile landscape-wide indications of vegetation greenness, we downloaded all available 8-day composite Moderate-resolution Imaging Spectroradiometer (MODIS) reflectance images (46 per year) for 2011, 2012, and 2013. MODIS products are compiled at a range of spatial resolutions; we accessed the highest resolution dataset (250 m). The composite dataset is produced by selecting the pixel obtained during the eight-day period that contains the best possible observation as evaluated on the basis of favorable coverage and view angle, the absence of cloud contamination, and aerosol loading, resulting in data that are pre-filtered for quality. We retrieved MODIS bands 1 (red; 620–670 nm) and 2 (near infrared (NIR); 841–876 nm) and created Normalized Difference Vegetation Index (NDVI) images from the normalized ratio of the near-infrared (NIR) and red bands, as follows:
Vegetation exhibits a distinctive spectral reflectance pattern, with low reflectance of red wavelengths due to chlorophyll absorption and strong reflectance in the near-infrared (NIR) due to cellular structure. The NDVI exploits this distinctive signature to provide a measure of the photosynthetic activity of plants (Figure 3).

Citizen Science Observations-Pima Canyon
We accessed a record of high-frequency buffelgrass observations taken at a single location near Tucson (in lower Pima Canyon) for this study. Pima Canyon is located in the Santa Catalina foothills on the northern edge of the city (Figure 2) in geomorphic conditions similar to those of Saguaro National Park–East; i.e., mountainous terrain composed of felsic granitic and metamorphic rocks. The south-facing slope of the Santa Catalina Mountains, which includes the Pima Canyon area, has a particularly dense and widespread buffelgrass infestation. The observations were collected by a single observer on a near-weekly basis from July 2010 to December 2014 using standardized protocols; data were archived in the USA National Phenology Network (USA-NPN) database. The visual estimates of percent greenness of up to five individual buffelgrass plants located along a south-facing slope of the trail were assessed on each day of observation (Figure 2).
Citizen Science Observations-Tucson
Beginning in June 2013, we recruited 10 citizen scientist volunteers to participate in the USA-NPN project Nature’s Notebook Citizen and Professional Science project for buffelgrass observations in the Tucson region. The volunteers were trained through the USA-NPN Education Coordinator and the Southern Arizona Buffelgrass Coordination Center (SABCC) as part of an effort to engage the general public in long-term phenology monitoring. Those interested in collecting data were encouraged to select a buffelgrass patch or individual plant to follow. All proposed locations were field checked and on-site training was provided. As at the Pima Canyon site, observations followed a series of strict protocols designed to capture detailed information regarding the timing and proportion of plant green-up, flowering, fruiting, and senescence. More information about the project can be found on the USA-NPN Tucson Phenology Trail website. By the end of 2014, 6913 observation records and 896 site visits had been entered into the database by the volunteer participants. The broad geographical distribution of these observations provided a complementary dataset to the single-site, long-term records from Pima Canyon.
Map of Buffelgrass Extent and Density
The National Park Service supplied a map of buffelgrass locations and densities within both SNP districts produced from a September 2012 helicopter survey conducted to document the distribution and quantity of buffelgrass, estimate its growth and spread, and prioritize areas for treatment. Each area was tagged with a visual description of the assessed buffelgrass density: trace (1% to 25%), low (>25% to 50%), medium (>50% to 75%), or high (>75% to 100%). To transform the polygonal coverages to the raster-based format of the remote sensing data, we calculated the percent of buffelgrass within each MODIS pixel by assigning the midpoint value to each polygon for low to medium density values (i.e., 12.5%, 37.5% and 62.5% for trace, low and medium respectively) and 80% to the high density category. We then converted the polygon coverage to a 5-m raster attributed with the percent values, and averaged all pixel values within each MODIS footprint (Figure 5).

Precipitation Data
Two different sources of rainfall data were used for analyses conducted at a single site (Pima Canyon) and those conducted over a landscape extent (SNP East and West). For Pima Canyon we accessed the closest source of recorded precipitation information, which was MesoWest daily data ([40]) at station AT021 (32.33050°N, 110.9120°W). The station is within 4 km and at a similar elevation (890 m) as the observation site. For the landscape comparison, we downloaded PRISM (“Parameter-elevation Regressions on Independent Slopes Model”) daily rainfall data ([41]) over the SNP districts. Using nearest neighbor techniques, we resampled the 4 km PRISM data to 250 m to conform to the spatial resolution of the MODIS data. For both datasets, we summed daily precipitation values over eight-day intervals to coincide with the MODIS composite data
Methods
Relationships between field-observed greenness, MODIS satellite greenness and precipitation data were explored and characterized at the Pima Canyon site, given the temporal depth of these data. Simple linear regression analyses and Spearman’s correlation coefficients were used to quantify relationships among the variables. On the basis of the correlations determined at the Pima Canyon site, we created a novel suite of metrics (“Climate Landscape Response” (CLaRe) metrics) to predict buffelgrass locations in SNP from park-wide satellite NDVI measurements and gridded 4 km precipitation data. These metrics range from 0 to 1.0 and quantify the strength of the linear relationship between the satellite greenness value and various cumulative precipitation values. They are created by conducting a linear regression between the NDVI and a precipitation value at each pixel and extracting the correlation coefficient.
Since we wanted to use the 2012 helicopter mapping as validation, we restricted the calculation of metrics to the 2011–2013 time frame, thereby bracketing the year the “ground truth” mapping data were generated. Statistical tests, including Student’s t-tests, were conducted to evaluate if the metrics were effective at separating pixels of pure native vegetation from those with buffelgrass. The CLaRe metrics with a balance between spatial and temporal prediction were used to create models of buffelgrass presence/absence in the park. That is, we used metrics that spatially displayed strong overall correlations and separability, but temporally allowed for earlier prediction of buffelgrass presence and condition. Buffelgrass presence/absence models were validated by examining their spatial coincidence with confirmed buffelgrass locations.
The linear regression analyses with strong (>0.50) correlations between greenness observed on the ground and various cumulative and lagged precipitation values were solved for the desired 50 percent greenness condition, which is the minimum greenness for which herbicide application can be effective at killing the invasive grass.
Results
Pima Canyon
Strong relationships exist at Pima Canyon between buffelgrass greenness observed on the ground and both MODIS greenness and PRISM precipitation. These relationships are visually evident in Figure 6, which plots the variables compiled for eight-day intervals and uses MODIS NDVI greenness. In particular, the field observations of buffelgrass greenness and those extracted from MODIS data are temporally synchronous, while buffelgrass greenness distinctly lags precipitation events. Figure 7 explores these relationships at Pima Canyon further, and graphs the correlation coefficients for the linear relationships between a sequence of cumulative precipitation values and both eight-day MODIS NDVI values and field observed greenness values, revealing strong (>0.5) correlation for all cumulative precipitation totals except the current (ppt0) and the prior time period (ppt1).


Saguaro National Park (SNP)
Several CLaRe metrics were developed from correlations quantified at the Pima Canyon site and mapped across the SNP study area. These metrics capture the strength of the relationships between satellite-measured greenness and various lagged precipitation values. At each pixel for each year, the correlation coefficient describing the strength of the linear relationship between the 46 MODIS NDVI values and the cumulative precipitation for time periods of different lags was calculated. For example, the correlation between MODIS greenness and the cumulative precipitation for the prior three eight-day periods (ppt123) for 2011, 2012 and 2013 is shown in Figure 8. Overlaid on these figures are the polygons of buffelgrass mapped in 2012 via the helicopter survey. We observed an especially strong spatial coincidence of the buffelgrass patches and the high correlation values for all years using this CLaRe metric. Of particular note is that the pattern of high CLaRe values is the same in all three years, even though the climate varied, with 2011, 2012 and 2013 years of above, normal and below average monsoon precipitation, respectively.

To explore the response of native vegetation communities to precipitation, we isolated all pixels that had had no buffelgrass mapped in the 2012 aerial helicopter survey, and stratified them according to the major vegetation classes in the park; their precipitation/MODIS greenness correlations are presented in Figure 9. Several patterns are evident. In all cases but one, the native vegetation types show very low correlation to recent precipitation, with values less than 0.5 for current precipitation and up to three eight-day periods of cumulative lagged precipitation. In all cases, correlation strength increases through time, showing that native vegetation types are more increasingly correlated to longer-term cumulative precipitation. In addition, the CLaRe phenometrics effectively stratify the vegetation types by rooting depth, with deep-rooted woodland and forest showing the lowest correlation values and shallower-rooted grasslands and scrublands showing higher correlation.

We evaluated the ability of various CLaRe phenometrics to locate buffelgrass by stratifying the landscape first by vegetation type and then by the density of buffelgrass within the pixels (Figure 10). The map shows the dominant vegetation types in each MODIS pixel footprint as mapped by the Southwest Regional Gap Analysis Project (SWreGAP). Note that even small amounts of invasive buffelgrass in the landscape have a strong effect on the correlation coefficients in the CLaRe phenometrics.

As seen in previous plots, the curves of all pure native vegetation types (the green curves in Figure 10) are monotonically increasing, whereas several of the curves for high density buffelgrass (the brown curves in Figure 10) peak and then decrease as the number of lagged time intervals increases. This is the same pattern seen at Pima Canyon (Figure 7). The buffelgrass curves are distinctly different than the curves for native Sonoran Desert vegetation communities. Not only do they display consistently higher correlations to all variations of cumulative precipitation, revealing that the invasive responds quickly to precipitation, but also these curves often peak, reflecting the invasive’s aggressive growth pattern.
Mapping Buffelgrass–Where?
Using the CLaRe metrics and the results for the individual vegetation types, we created several models of buffelgrass presence/absence. For modeling purposes, we focus on SNP-East, given its geologic and geomorphic similarity to the Pima Canyon site and its diversity of vegetation types and elevational gradients. We also focus on the shorter-term cumulative precipitation lags, given that we seek to provide operationally timely data and that these CLaRe metrics display strong differences between pure native and buffelgrass-invaded pixels. Our modeling was conducted by choosing options that seemed logical, useful and easily understood by practitioners; we did not apply statistical methods to optimize classification.
One of the best models we created is shown in Figure 11. This was developed to be easily created with no prior knowledge of buffelgrass distribution. Based on exploratory statistics and prior results, these models used the CLaRe ppt123 for each year. All CLaRe ppt123 values were extracted for each mapped vegetation type and the top 5th of the highest correlations were mapped as buffelgrass presence (i.e., the upper 20 percent of calculated correlation values). For modeling and validation purposes, we eliminated the higher elevation vegetation types that are unsuited for buffelgrass (e.g., Madrean pine-oak forest and woodland). The three sub-models show the results for each year (left panels, Figure 11). The composite model (right image, Figure 11) maps buffelgrass presence if 2 of 3 sub-models mapped buffelgrass presence.

Validation results suggest that the three sub-models do a reasonable job identifying known patches of buffelgrass, capturing 45 to 49 percent correctly. The composite models improve the accuracies: composite model 2 of 3 (shown in Figure 11) has an overall accuracy of 83 percent and correctly maps 46 percent of known buffelgrass. There are tradeoffs between overall accuracy and the total amount of buffelgrass correctly identified compared to the amount misclassified and missed. In all cases, the models show poor user’s accuracies (typically <25%) for mapping buffelgrass, with the model classifying areas as buffelgrass that were confirmed native vegetation in the field. Suggestions for improving the accuracies are presented in the discussion section below.
Mapping Buffelgrass—When?
Because optimal use of herbicide to control buffelgrass requires that the plant is photosynthetically active when treated, it is important to map not only where buffelgrass is located, but also when it is green. To predict the phenological status of buffelgrass, we analyzed the citizen science data and solved the linear equations relating field observed buffelgrass greenness to various cumulative precipitation totals for precipitation at specified amounts of buffelgrass greenness.
Results for the Santa Catalina Mountain sites are shown in Tables 1. To achieve the management threshold of greenness for effective herbicide application, results suggest buffelgrass will be 50% green 8 to 16 days after a 24 day period that totals over 46 mm (1.8 inches) of rain in the Santa Catalina Mountains (less is required in the Tucson Mountains, with different soils, etc.). In SNP-east, where we were unable to monitor buffelgrass on the ground, similar geology and geomorphology suggests that the thresholds observed in the Santa Catalina sites would be most appropriate to use. Of course, knowledge of these thresholds will evolve and become more accurate as practitioners interact with the landscape. In addition, the following section will show how MODIS data can be used to confirm green up in an area targeted for treatment.
Operational Protocol
The research results can be used to create a general protocol for operationally informing management activities related to herbicide treatment of buffelgrass at SNP and the Tucson region. Certain decisions, although produced from the data specific to this area, are likely to be generalizable. For example, the metric CLaRe ppt123 was used to model buffelgrass presence in the landscape for our region. Given the characteristic aggressive response of buffelgrass to precipitation, it is likely this metric will be useful, if not optimal, in other areas as well. Suggested steps for optimizing the timing of herbicide treatment for buffelgrass in the landscape are as follows:
- Map where buffelgrass is by using historic MODIS NDVI data and PRISM precipitation data to create CLaRe ppt123 metrics for several years of data. Include any model refinements in this step, as presented in the discussion section. Modeling results from this study suggest that an effective way to screen a region for buffelgrass presence is to extract CLaRe ppt123 metrics (i.e., cumulative precipitation for the three prior eight-day periods) for vulnerable vegetation types, and examine areas with the highest correlation metrics as candidate landscapes for hosting buffelgrass. In our study, good results were achieved by examining the top 5th (highest 20%) of the correlation values. The upper threshold selected can reflect expert knowledge of the extent of buffelgrass invasion in the management area, the size of the area that can be reasonably examined given resource constraints, and the relative vulnerability to invasion of landscapes of concern. Creating composite models from several years of data is shown to improve accuracy; composites are created using CLaRe ppt123 by mapping buffelgrass presence if half or more of the suite of individual annual models map buffelgrass.
- Select areas to target for treatment and collect the MODIS NDVI and PRISM precipitation data at that location. For example, eight-day MODIS (solid green line) and PRISM data (solid blue line) for 2011 at a pixel in SNP-East that contains a patch of high density buffelgrass is shown in Figure 12. Note that PRISM values are available daily and, as of this writing, eMODIS data, which is the composite NDVI value for the prior seven days, are also available daily. MODIS eight-day composite data used in this study are shown to compare very well with eMODIS data. This means that it is now possible to access for free PRISM 4-km precipitation data and 250-m eMODIS NDVI data every day that are current through yesterday.
- Predict when buffelgrass will be 50 percent green by calculating a running tally of the cumulative precipitation for the past 24 days (dashed magenta line). When the prior 24-day period totals over 46 mm (1.8 inches) of rain (vertical red bar), place crews and aircraft on notice for herbicide application 8 to 16 days from that time (shaded region of graph).
- At eight days prior to herbicide application, confirm that the target has greened up using current MODIS data (green line), adjust areas for treatment as needed.
On the day of treatment, confirm that the target has greened up using current MODIS data (open black circle), adjust areas for treatment as needed.

Discussion
The composite model for mapping buffelgrass (Figures 11) is ~80 percent accurate overall and effective at capturing ~50 percent of known buffelgrass patches, but has poor user’s accuracy (<26 percent). Inspection of the geography of misclassified pixels suggests environmental factors could be affecting the accuracy. In the model and sub-models, the producer’s error tends to occur in an east-west band near the southern boundary of the park whereas the user’s error tends to occur in the northwest (Figure 11). This observation suggests that model results may be improved if areas with distinct aspect are isolated before extracting the highest correlations for composite models.
Other factors may contribute to over-mapping buffelgrass patches and poor user’s accuracies: the coarseness of the PRISM precipitation data used (4 km), the accuracy and relevance of the vegetation map, and the presence of other invasive plants, such as Schismus and fountain grass (Cenchrus macrouru) that display similarly aggressive green-up behavior. The model may be improved using finer resolution precipitation data and more accurate maps of vegetation and invasive species. The SWreGAP map was produced using data collected in 2000 through 2002; an update may be beneficial. Although SNP is protected, vegetation changes may have resulted from fire and drought over the 10+ years before the 2012 helicopter mapping survey was conducted.
This study pooled annual MODIS NDVI composites and PRISM precipitation data to produce CLaRe metrics. Additional insights will likely be gleaned by not only splitting the data into a winter (cool) season and a summer (warm) season, but also by creating CLaRe metrics using temperature. During exploratory analyses, when observed buffelgrass greenness at the Pima Canyon site was plotted against daily temperature variables, a bimodal distribution was: Buffelgrass greened up in response to rainfall in both cool and warm seasons, unlike native C4 grasses that respond primarily to warm-season precipitation. The cool-season greenup of buffelgrass is likely predicated on reaching a certain threshold for temperature; such relationships can be readily explored and mapped spatially using CLaRe phenometrics.
This study demonstrates the usefulness of CLaRe phenometrics in a specific management application related to control of invasive buffelgrass, but they hold promise for many issues related to dryland vegetation. The CLaRe phenometrics used here explicitly incorporate precipitation and capture the strength of the correlation between the satellite greenness data and precipitation data (both lagged and cumulative). These phenometrics are effective in locating populations of invasive buffelgrass since the grass responds more rapidly and strongly to precipitation events than native vegetation. Even small amounts of invasive buffelgrass in the landscape have a strong effect on the values of these phenometrics. Given that many invasive species are successful because they can out-compete native species for limited resources, it is likely the CLaRe phenometrics will prove valuable for mapping other invasive species of regional extent in dryland ecosystems.
Conclusions
Buffelgrass is a noxious invasive plant that threatens to degrade the iconic Sonoran Desert ecosystems by out-competing native species and altering fire regimes. Effective treatment of buffelgrass using herbicide requires knowledge of where buffelgrass is located in the landscape and when it is at least 50 percent green. Optimizing herbicide treatment is crucial for many reasons, including cost, efficacy, potential collateral damage to the ecosystem, public opinion, and social concerns. Our research enables the targeted and strategic application of buffelgrass herbicide by using remote sensing data to detect when and where buffelgrass is photosynthetically active and predicting when it will be at least 50 percent green.
CLaRe metrics capture the aggressive response of buffelgrass to recent precipitation: buffelgrass-invaded pixels typically display statistically higher values (i.e., stronger correlations, at ρ < 0.05) than pixels of pure native vegetation for all vegetation types in SNP. Composite models created from CLaRe metrics produced validation results with overall accuracies of ~80 percent and correctly mapped ~50 percent of known buffelgrass patches. Results from the study suggest accuracy can be improved by further stratifying the landscape into large-scale topographic facets of similar aspect. Refinements may also be realized by focusing analyses on specific seasons (winter vs. monsoon) and by creating CLaRe phenometrics with temperature data in addition to precipitation, especially in the winter rainy season.
Citizen science observations at Pima Canyon were key in linking the greenness observed on the ground to the satellite-based MODIS-NDVI data. In addition, the citizen science observations from across the region allowed for the prediction of when buffelgrass will be at least 50 percent green and thus suitable for effective herbicide treatment. In the general Tucson Area, a conservative estimate is that buffelgrass will be at least 50 percent green 8 to 16 days after a 24-day period with cumulative rainfall totals of 46 mm (1.8 inches).
The phenometrics introduced here are a new and innovative suite of metrics that are proven effective in mapping when and where invasive buffelgrass is green in Saguaro National Park near Tucson, Arizona. The ability to monitor and interpret trends in land surface phenology informs science-based land management. In contrast to temperate regions, which respond to the annual solar cycle, producing predictable greenness curves readily captured by traditional phenometrics, drylands respond to precipitation, which can be erratic. As such, the “start of season” in drylands is more correctly seen to be the onset of a significant precipitation event, such as the monsoon. The CLaRe phenometrics capture the landscape response to climate and will complement the current temperate-zone phenometrics; they will prove broadly useful for many applications, including monitoring dryland ecosystems, mapping invasive species, habitat modeling, and vegetation classification.
Acknowledgments: We wish to thank the National Park Service and Dana Backer, in particular, for supplying us with data and critical feedback on the methods used here. Data were provided by the USA National Phenology Network and the many participants who contribute to its Nature’s Notebook program. The USGS (Land Change Science, Land Remote Sensing and National Park Monitoring Project within the Status & Trends Program) provided funding for this research. Jessica Walker was supported by a USGS Mendenhall Research Fellowship during this work. We appreciate the thoughtful and constructive comments from our three anonymous reviewers and USGS internal reviewer Joel Sankey. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.