Exploiting Differential Vegetation Phenology for Satellite-Based Mapping of Semiarid Grass Vegetation in the Southwestern United States and Northern Mexico

: We developed and evaluated a methodology for subpixel discrimination and large-area mapping of the perennial warm-season (C 4 ) grass component of vegetation cover in mixed-composition landscapes of the southwestern United States and northern Mexico. We describe the methodology within a general, conceptual framework that we identify as the differential vegetation phenology (DVP) paradigm. We introduce a DVP index, the Normalized Difference Phenometric Index (NDPI) that provides vegetation type-speciﬁc information at the subpixel scale by exploiting differential patterns of vegetation phenology detectable in time-series spectral vegetation index (VI) data from multispectral land imagers. We used modiﬁed soil-adjusted vegetation index (MSAVI 2 ) data from Landsat to develop the NDPI, and MSAVI 2 data from MODIS to compare its performance relative to one alternate DVP metric (difference of spring average MSAVI 2 and summer maximum MSAVI 2 ), and two simple, conventional VI metrics (summer average MSAVI 2 , summer maximum MSAVI 2 ). The NDPI in a scaled form (NDPI s ) performed best in predicting variation in perennial C 4 grass cover as estimated from landscape photographs at 92 sites (R 2 = 0.76, p < 0.001), indicating improvement over the alternate DVP metric (R 2 = 0.73, p < 0.001) and substantial improvement over the two conventional VI metrics (R 2 = 0.62 and 0.56, p < 0.001). The results suggest DVP-based methods, and the NDPI in particular, can be effective for subpixel discrimination and mapping of exposed perennial C 4 grass cover within mixed-composition landscapes of the Southwest, and potentially for monitoring of its response to drought, climate change, grazing and other factors, including land management. With appropriate adjustments, the method could potentially be used for subpixel discrimination and mapping of grass or other vegetation types in other regions where the vegetation components of the landscape exhibit contrasting seasonal patterns of phenology.


Introduction
Across the southwestern and western United States, grazing of grass forage by livestock is the dominant land use, occurring in 73% and 83% of the land area in Arizona and New Mexico, respectively, and 62% of the combined land area of the 11 western-most conterminous U.S. states [1]. Several historical, ongoing or emerging stress factors, notably overgrazing, landscape fragmentation, climate change, persistent drought, woody plant encroachment, and invasive plant species, pose significant threats to the biological integrity and resilience of Southwest ecosystems and present major challenges to the design of effective strategies for land management, mitigation and conservation [2][3][4][5][6][7][8][9][10][11][12][13][14]. From the standpoint of agricultural and mixed-use land management, detrimental impacts on the current and future capacity of grass vegetation in rangelands of the Southwest and beyond to support grazing for livestock production is a paramount concern [15][16][17][18][19]. Knowledge and understanding of these issues from both scientific and land management perspectives, as well as the design of effective strategies for mitigation and adaptation, can benefit from improved, timely data and information about the distribution and surface cover of grass vegetation in the Southwest and their response to these stress factors over time. By enabling consistent observations over seasonal, annual, decadal and longer time periods, image-based, multispectral remote sensing can be an effective tool for mapping and monitoring of the distribution and relative productivity of grass and other vegetation types and their changes over time.

Multispectral Image-Based Approaches to Mapping of Semiarid Grass
A diverse set of analytical methods that exploit phenological information in time-series of multispectral images has been developed for distinguishing and mapping vegetation cover types [20][21][22][23][24][25], estimating and mapping biomass and vegetation fuel loads [26][27][28][29][30] and for investigating how vegetation and related ecosystem processes respond to management, disturbance, weather, drought and climate change [31][32][33][34][35][36][37]. These methods examine the seasonal growth dynamics of the vegetation through quantification and analysis of various metrics of the time signal observed within one or more years. Time-series observations of a spectral vegetation index (VI), such as the normalized difference vegetation index (NDVI), are commonly used for this purpose [38,39]. Typical metrics include, for example, the dates of the beginning and end of the growth period, the length of the growth period, and date and magnitude of the peak value of the NDVI or other VI [40].
Interpretation of such metrics is straightforward for land areas with a dominant vegetation type or a common phenology, but can be complex for areas with horizontal (surface cover) or vertical (understory and overstory) variation in vegetation type and phenology [41][42][43][44]. In the mixed-composition landscapes of the Southwest, the observed time signal may be composed of one or more discrete and potentially unique signals attributable to the individual grass, shrub, tree and succulent components of the landscape, with their respective contributions modulated by their fractional surface cover within the pixel. When seeking to discern the perennial warm-season grass component of the vegetation cover, such mixed time signals confound the extraction and interpretation of standard phenology metrics.
Many past studies have addressed the mixed signal problem in efforts to distinguish signal contributions from herbaceous and woody vegetation. For example, Moody and Johnson [45] showed that the first and second harmonics from a discrete Fourier transform of NDVI time-series data acquired by the Advanced Very High Resolution Radiometer (AVHRR) over a study area in southern California are sensitive to variations in grass cover, woody vegetation, and agriculture. Hirosawa et al. [46] applied principal components analysis to time-series NDVI data from the AVHRR. They showed that variations in the histograms for the first two principal components provided a measure of the homogeneity/heterogeneity of the cover types within the grasslands and other major vegetation communities of Arizona. Roderick et al. [47] and Lu et al. [48] used methods based on the seasonal-trend decomposition by LOESS (STL) procedure [49] to map fractional cover of evergreen (woody) vegetation by isolating its signal from the seasonally green (herbaceous) signal in time series NDVI data for the Australian continent acquired by AVHRR. Pisek et al. [50] used a semi-empirical approach [51] and a physically-based approach [52] to retrieve the understory (herbaceous) and overstory (woody) vegetation components of the seasonal NDVI signal observed at high latitude forest sites in Northern Europe. Scanlon et al. [53] used seasonally averaged NDVI and interannual variations in wet-season NDVI as state-space variables in a linear unmixing model to estimate fractional cover of trees, grass and bare soil along an aridity gradient in southern Africa. Johnson et al. [54] achieved improved accuracy in estimating fractional green vegetation cover with spectral mixture modeling by using spatial interpolation to account for spatial variation in the values of spectral end-members. Gessner et al. [55] applied decision tree regression to a set of satellite-derived phenology metrics for the wet and dry season to estimate fractional cover of woody, herbaceous and bare soil in savanna regions of southern Africa. Brandt et al. [25] estimated and mapped woody plant cover in the Sehelian drylands of Africa by using multiple regression modeling of the relationships between various satellite-derived phenological metrics and field data.
Helman et al. [56] and Chekroun et al. [57] used a distinctly different approach to distinguish the contributions of herbaceous and woody vegetation to the VI time signals in their Mediterranean study areas. In contrast to conventional signal decomposition, their approaches were designed a priori to take strategic advantage of two intra-annual periods, the dry season and the wet season, when the seasonal phenology of herbaceous and woody vegetation diverge. Because herbaceous vegetation is absent or senescent during the dry season, the VI signal during this period is attributed solely to woody vegetation and background soil. During the wet season when the herbaceous vegetation becomes photosynthetically active, the composite signal from the ecosystem includes contributions from all components (woody vegetation, herbaceous vegetation, and background soil). These studies demonstrated that the signal component associated with seasonal herbaceous vegetation could be isolated from the woody signal by exploiting this seasonal divergence in the VI response. Chekroun et al. [57] used this approach with several VI data types from the Moderate Resolution Imaging Spectroradiometer (MODIS) to achieve improved sensitivity to drought across landscapes of Northern Tunisia with large variations in the proportional surface cover of trees, shrubs and grass. Helman et al. [56] showed this approach is effective for estimating and mapping spatial variation in the subpixel fractional surface cover of woody vegetation and herbaceous vegetation in their study area in Israel.

A Paradigm for Exploiting Vegetation Phenology Information in Satellite-Based Vegetation Mapping
The method presented by Helman et al. [56] is a particular example of a general, conceptual approach that has strong potential for quantitative mapping and monitoring of vegetation in mixed-and variable-composition landscapes, yet to our knowledge remains relatively unexplored. We identify the conceptual approach by its emphasis on the strategic leveraging of a priori information and understanding of distinct, seasonal differences in the phenology of the vegetation components within a landscape. It is distinguished from other phenology-based (frequency-domain) approaches and those that emphasize spatial and spectral information by its use of metrics that exploit predetermined differences in phenology. The metrics are used strategically to extract from the composite time signal information on the fractional cover (and potentially, the condition in response to stress) of a specific vegetation component of the landscape. We identify this approach as a differential vegetation phenology (DVP) approach to mapping and monitoring of vegetation cover.

A Differential Vegetation Phenology Approach to Mapping of Perennial C 4 Grass in the Southwest
The surface cover and productivity of the perennial grass component of the arid and semiarid Southwest landscapes can vary widely across the region and over time, reflecting the interacting influences of climatic and edaphic conditions (expressed primarily as water limitations on growth), inter-species competition for resources, disturbance (e.g., wildfire and pests), herbivory, and land use associated with human activities, particularly (but not exclusively) grazing [14,[58][59][60][61][62]. Under these particular, highly variable growth conditions, grasses (as well as other plant types) exhibit unique adaptive strategies that provide distinct competitive advantages, which are integral to their occurrence and sustained existence in the Southwest. Such adaptations are expressed in numerous ways, including not only unique physiognomy and physiology (such as shallow roots and C 4 photosynthesis), but also distinct patterns of seasonal growth, or phenology [58,[63][64][65].
The perennial C 4 grasses, which are the dominant grass type in our study area, exhibit a phenological pattern that is tightly coupled with the region's characteristic but highly variable seasonal cycle of temperature (winter minimum, summer maximum) and precipitation (winter and summer, with drought in spring and autumn) [62,66]. Their principal, active growth period occurs during the summer and early autumn when warm temperatures suitable for C 4 photosynthesis coincide with the arrival of precipitation from the North American Monsoon, and is typically followed by a period of senescence and dormancy in winter and spring [62]. This pattern contrasts starkly with the phenology of the other major vegetation types, primarily shrubs, trees and succulents, that coexist with perennial grasses in Southwest landscapes and remain green during a longer period (mesquite and riparian trees) or throughout the year (most shrubs and succulents) [67][68][69]. These contrasting seasonal patterns of phenology can be exploited to estimate and map the fractional cover of perennial grass through a DVP-based analysis of time-series, satellite-derived multispectral image data.
Here, we present the conceptual foundation and our development, application and evaluation of a DVP methodology that uses times-series observations from Landsat 5 and the Moderate-Resolution Imaging Spectroradiometer (MODIS) to estimate and map the fractional cover of perennial warm-season grass within a study area that encompasses parts of southern Arizona, USA and northern Sonora, Mexico. The methodology is based on a remote sensing-derived DVP metric, which we developed as a Normalized Difference Phenometric Index (NDPI). Our application in the Southwest suggests the NDPI can be an effective source of information on the distribution of grass and rangelands to support science and land management. With appropriate adjustments, the NDPI could potentially be used for subpixel discrimination and mapping of grass or other vegetation types in other regions where the vegetation components of the landscape exhibit contrasting seasonal patterns of phenology.

Study Area
The study area covers approximately 30,000 km 2 in southern Arizona, USA and northern Sonora, Mexico ( Figure 1). The area encompasses parts of two distinct ecoregions of the Southwest, the Apache Highlands and the Sonoran Desert, which dominate the southeastern and northwestern segments, respectively. The Apache Highlands ecoregion extends across approximately 12.1 million hectares in southeastern Arizona, southwestern New Mexico, northern Sonora, and northwestern Chihuahua [8,70,71]. Its physiography is characterized by isolated mountain ranges, or sky islands, and intervening valleys, which underlie the region's unique ecology and biodiversity. Topographic elevations within the study area range from roughly 450 m at the desert lowland in the northwest, to more than 2500 m at the highest peaks of the sky islands ( Figure 2).
The study area contains a diverse array of vegetation communities, with pine-oak woodlands and mixed conifer forests at higher elevations, grassland and shrubland at moderate elevations, and desert scrub at lower elevations [8]. The grass vegetation across the study area is dominated by perennial, warm-season (C 4 ) grass (both native and invasive species); annual grasses are a relatively minor constituent (Table 1). We therefore focus our effort on the mapping of perennial, C 4 grasses across the study area. The Apache Highlands ecoregion transitions to the Sonoran Desert ecoregion as elevation generally decreases and temperature and aridity increase toward the northwestern and western parts of the study area. This transition is characterized by a general decline in the presence and surface cover of perennial grasses until they become essentially absent in the Sonoran Desert's core regions, except for riparian zones and other locations where water limitations are less extreme.
The climate of the study area, as well as the broader region encompassing the full extent of the Apache Highlands and eastern portion of the Sonoran Desert ecoregions, is greatly influenced by the precipitation regime of the North American Monsoon [72,73]. Annual total precipitation is Remote Sens. 2016, 8, 889 5 of 33 characteristically low and highly variable in both space (non-uniform distribution at multiple spatial scales) and time (large fluctuations across seasons and among years). Despite inherent variability, the seasonal distribution of precipitation exhibits a distinct, bimodal pattern ( Figure 3): winter precipitation (primarily from frontal storms) followed by spring drought, then summer convective precipitation (primarily from the North American Monsoon), followed by autumn drought. The inner, white rectangle is the area of our Landsat-based analysis, which is located within the overlap of the two Landsat path-rows. The small Arizona map inset at the top-right shows the coverage of the satellite data sets: two adjacent, overlapping path-rows of Landsat, and a tiled MODIS product.        Figure 2. Topographic elevation across the study area with locations indicated for the Tucson metropolitan area, the communities of Picacho and Elgin, Arizona, and major sky island mountain ranges. The MODIS-based analysis encompassed the full extent of the map; the area of the Landsatbased analysis is indicated by the broken line.

Conceptual Basis
The major vegetation types within our study area, primarily perennial warm-season (C 4 ) grasses, shrubs, trees and succulents, represent successful evolutionary responses to climate, edaphic conditions, and interspecific competition for resources. For an individual vegetation type, this response has multiple expressions, including unique physiognomy, physiology and phenology [76,77] and variation in spatial distribution and surface cover. The phenology exhibited by the perennial C 4 grasses closely follows the seasonal precipitation pattern dictated by the North American Monsoon [60,62,68,72,73,78].
The perennial C 4 grass phenology in our study area is distinct from the phenology of the other major vegetation components of the landscape (trees, shrubs and succulents) [61]. The photosynthetic activity and green-leaf display of perennial C 4 grasses differ dramatically between the pre-Monsoon period (May and June) and the subsequent Monsoon-dominated period (July through September with influences extending through October), however remain relatively stable for trees, shrubs and succulents throughout both periods. This difference reflects distinct adaptive strategies to the water-limited growth environment [63,64,79]. The perennial C 4 grasses have relatively shallow but extensive root systems that facilitate efficient capture and rapid photosynthetic response to variations in soil water from ephemeral rain events of the North American Monsoon. In contrast, the deeper root systems of the trees and shrubs enable access to greater depths below ground where soil moisture is more consistently available and seasonal fluctuations are dampened [59][60][61]80].
Because seasonal changes in leaf display and photosynthetic pigments are the primary phenology-related signal detectable by multispectral optical remote sensing, such differential expressions of phenology among perennial C 4 grasses and the trees, shrubs and succulents (in contrast to their physiognomy and physiology) can be directly observed and analyzed with time-series observations from Earth-observing satellites such as Landsat and MODIS. We take advantage of this capability with a DVP approach to discriminate and map the perennial grass component within our study area.

Normalized Difference Phenometric Index
The seasonal pattern observed in time-series, satellite-based, multispectral (visible and near-infrared) images at an individual pixel location represents the integrated seasonal response of the various vegetation types of the landscape within the pixel, or overall landscape phenology. We introduce here a quantitative, phenology-based, temporal vegetation index for application to time-series of multispectral (visible and near-infrared) satellite images.
Our objective in developing this index is to provide information on the presence and surface cover of a specific vegetation type at the subpixel scale within a mixed-composition landscape and how that vegetation component responds to interannual and longer term changes in climate, and/or disturbance caused by grazing, fire, pests and invasive species. We pursued this objective by exploiting known differences in the phenology, or seasonal pattern of green leaf display and photosynthetic activity, among the major vegetation cover types within the landscape. We define and demonstrate the index in the context of our present research objective for satellite-based mapping of the perennial grass cover in our Southwest study area. We developed this phenology-based index as the Normalized Difference Phenometric Index (NDPI) and calculate it in raw form as where VI is the average spectral vegetation index for a specified, time period. The subscript parameters p1 and p2 are discrete, non-overlapping time periods within a year. Whereas the equation for NDPI has the same form as the equation for the well-established normalized difference vegetation index (NDVI), it is altogether distinct with respect to its input variables (vegetation index vs. spectral reflectance) and its purpose (exploiting temporal information vs. spectral information).
For mapping and analysis, we convert the raw NDPI to a scaled index, NDPI s : where NDPI obs is the raw NDPI value observed for an individual pixel, and NDPI ref is the raw NDPI value associated with the reference vegetation cover conditions specified for the target vegetation type (perennial C 4 grass in this study). We define the reference vegetation cover as nominally pure, continuous, and homogeneous. Equation (2) converts NDPI obs values to a consistent, convenient and intuitive measurement scale in which 0 indicates maximum (essentially continuous) surface cover and −NDPI ref (the additive inverse of NDPI ref ) indicates absence of the target vegetation type. We propose that the magnitude of NDPI s between 0 and −NDPI ref is a potential measure of the fractional surface cover of the target vegetation type in the observed pixel. We demonstrate this relationship in our present application to mapping of perennial C 4 grass in the Southwest. Although we consider the nominal range of output values from Equation (2) to be 0 to −NDPI ref , the equation itself does not preclude positive values of NDPI s . However, in practice positive values of NDPI s will occur only under certain, unique circumstances, and such occurrences are not necessarily problematic. The first circumstance is associated with our process for specifying the value of NDPI ref .
In the present application, we did not determine NDPI ref from the single, absolute maximum NDPI value observed within the study area for our target vegetation type (perennial C 4 grass), rather we used the mean value of a larger sample set of pixels. NDPI values that are greater than the mean of this sample would produce positive NDPI s values, however their deviation from 0 is typically small. Indeed, careful examination of positive NDPI s values in an iterative process to specify NDPI ref can provide valuable information and feedback for optimizing the selection of sample pixels. Another circumstance that could produce positive NDPI s values is when a component of the surface cover exhibits seasonal changes in the VI (between p1 and p2) that mimic the phenology of the target vegetation type, and the resulting NDPI value exceeds the specified value of NDPI ref . As we show subsequently in our results, such conditions may occur in areas of agricultural crops, where the seasonal cropping pattern is similar to the phenology of the target vegetation type (perennial C 4 grass in our study area).
Conceivably, other circumstances may occur that could also produce positive NDPI s values. For example, when the NDPI method is used for monitoring over multiple years, a positive deviation in NDPI obs from NDPI ref in any individual year would produce a positive value of NDPI s . In such a monitoring application, positive NDPI s values could be readily identified and examined with reference to climate data and management information, and provide meaningful insights about potential, causal factors, such as increased summer precipitation or reduced grazing pressure. Additional circumstances in either mapping or monitoring applications of the NDPI, besides those considered here, could potentially occur that would also produce positive NDPI s values. Regardless, occurrences of positive NDPI s values are not necessarily problematic and can often be meaningful and informative.

Specification of NDPI Parameters and Variables
The parameters p1 and p2 in the NDPI equation (Equation (1)) are selected to isolate seasonal intervals when the phenology of the target vegetation type diverges substantially from the phenology of other vegetation types within the landscape. The selections are made with reference to divergence in seasonal photosynthetic activity and green leaf display that occurs consistently, even if variably (in response to weather and climate) among years. We distinguish such divergence from episodic changes in seasonal photosynthesis resulting from disturbance events, for example, fire and land cover conversion [81]. For mapping of perennial C 4 grass in our Southwest study area, we specify p1 as day 121 to 176 (1 May through 25 June), and p2 as day 177 to 304 (26 June through 31 October). The break point between p1 and p2 (25/26 June) accounts for the end and start dates of the corresponding 8-day composite periods of the MODIS data product, described subsequently in Section 2.8.
The dates for p1 (1 May to 25 June) bracket a period when the perennial C 4 grasses in the study area exhibit minimal photosynthetic activity and greenness. This period immediately precedes the period when rainfall from the onset of the North American Monsoon typically becomes significant. The dates for p2 (26 June to 31 October) capture the dominant period of rainfall associated with the North American Monsoon (July through September) when the perennial C 4 grasses exhibit annual maximum photosynthesis and productivity, but also includes October when prior Monsoon (and October) rainfall can still support additional productivity before cooling temperatures become limiting for C 4 photosynthesis.
We purposely selected the p1 start-date (1 May) and the p2 end-date (31 October) to limit our analysis to the predominant leaf-on period of the deciduous vegetation in the study area (mesquite, cottonwood and other riparian trees). In comparison to perennial C 4 grass, the greenness of all other major vegetation components in the study area (both deciduous and evergreen) may therefore be reasonably assumed to be relatively stable throughout our specified p1 and p2 periods. This assumption is supported by observation-based spectral vegetation index data used in our subsequent numerical modeling of the NDPI s response to variation in fractional cover. Our specified time periods also help to minimize the signal contribution from annual grasses that are active in the cool season.
Our specification of p1 and p2 shares a common rationale with Helman et al. [56], who used MODIS-derived NDVI data to exploit similar differential patterns of vegetation phenology in wet and dry seasons to distinguish the herbaceous and woody components of vegetation cover in Mediterranean forests and woodlands in Israel. Helman et al. [56] used the average NDVI during the dry season to represent the signal contribution from woody vegetation and background soil (NDVI W ), and the difference of the maximum NDVI in the wet season and NDVI W to represent the herbaceous component (NDVI H ).
Our development and application of the NDPI for mapping of perennial warm-season grasses in the Southwest modifies the seasonal difference calculation used by Helman et al. [56] in four ways. Firstly, we purposely use time-series VI data that represent the long-term (15-year) average instead of one or more individual years to minimize weather-related differences in the VI response. While such interannual differences can potentially provide useful information for monitoring purposes, multiple-year average VI data (in our application, 15-year averages using 8-day MSAVI 2 data) are preferable for vegetation cover mapping because they emphasize spatial information by reducing confounding information associated with disturbance-or weather-induced interannual variations in the occurrence, greenness, or biomass of the perennial C 4 grass. Moreover, for monitoring purposes, an NDPI map produced from multiple-year average VI data is useful in defining reference conditions against which such fluctuations in yearly NDPI values can be quantitatively evaluated.
Secondly, whereas Helman et al. [56] used average and maximum VI (NDVI) data for the seasonal dry and wet periods, respectively, we use average VI (MSAVI 2 ) data for both periods (p1 and p2). By using the average VI value in p2 (wet season), we seek to reduce potential bias associated with the effects of interannual variability of the within-season temporal distribution of precipitation.
Thirdly, whereas Helman et al. [56] use the difference of the wet-and dry-season VI values, we use the NDPI formulation (Equation (1)) to transform the simple ratio of the wet-season (p2) and dry season (p1) VI values to a fixed scale with a finite range of −1.0 to +1.0. The response of VI(p1) and VI(p2) to changes in subpixel fractional cover of perennial C 4 grass is functionally (but not physically) analogous to the response of red and NIR reflectance to variation in subpixel fractional cover of healthy, green vegetation of any type. With increasing fractional cover of the surface component of interest (perennial C 4 grass for NDPI in this study; any green vegetation for NDVI), both VI(p1) and red reflectance decrease as both VI(p2) and NIR reflectance increase.
The purpose of the NDPI transformation is not improved quantification of the subpixel fractional coverage of target vegetation type (perennial C 4 grass in this study). Instead, our motivations parallel those for using the NDVI in the spectral domain [82]: to improve the dynamic range under low vegetation-signal conditions, improve the visualization and interpretation of the information contained in the data, and facilitate comparison and interpretation of results among different applications and studies. In the NDPI context, low vegetation-signal conditions manifest as small differences between the VI values of p1 and p2 (as expected when the fractional cover of the target vegetation type is small relative to other surface cover), whereas in the NDVI context they manifest as small differences in the NIR-and red-band reflectances (as expected when the land surface is sparsely vegetated). Because the NDPI and NDVI transformations have the same mathematical expression, their numerical results therefore have equivalent sensitivity to variation in the values of their corresponding input variables. We discuss this issue and related implications for both our present study and future, potential applications of the NDPI subsequently in Section 2.5. Fourthly, we convert the NDPI values (Equation (2)) to an intuitive scale in which 0 indicates maximum surface cover and −NDPI ref indicates absence of the target vegetation type.
By maximizing the seasonal contrast in perennial C 4 grass phenology, while avoiding or minimizing seasonal differences in the phenology of all other major vegetation types in the study area, our specification of p1 and p2 enables maximum sensitivity to the phenology exhibited by our target vegetation type, perennial C 4 grass. With this parameterization, the NDPI provides a means to isolate and extract information on specifically the perennial C 4 grass component of an integrated landscape phenology signal arising from variable combinations and cover fractions of vegetation types within an individual, satellite-observed pixel. Because satellite-based land imagers cannot directly detect grass cover that may be obscured by densely foliated over-story vegetation such as trees and large shrubs, the method is most suitable for discrimination and mapping of the exposed, non-understory grass component of vegetation surface cover. As the exposed areas of perennial C 4 grass cover across the study area decrease from a reference condition (characterized by essentially dense, continuous surface cover), it is necessarily replaced by other vegetation types, bare soil, or manmade materials. Decreasing perennial grass cover is manifested through systematic changes in the seasonal pattern of a satellite-derived VI, as displayed in an idealized form in Figure 4. Our specification of p1 and p2 in Equation (1) maximizes the sensitivity of the NDPI to these changes. In this research we investigated the NDPI s as a quantitative indicator of the divergence of the phenology of an observed landscape from the phenology of a reference landscape associated with continuous surface cover of perennial grass (Figure 4). other vegetation types, bare soil, or manmade materials. Decreasing perennial grass cover is manifested through systematic changes in the seasonal pattern of a satellite-derived VI, as displayed in an idealized form in Figure 4. Our specification of p1 and p2 in Equation (1) maximizes the sensitivity of the NDPI to these changes. In this research we investigated the NDPIs as a quantitative indicator of the divergence of the phenology of an observed landscape from the phenology of a reference landscape associated with continuous surface cover of perennial grass ( Figure 4). For VI in Equation (1), we used the modified soil-adjusted vegetation index (MSAVI2) [19]. MSAVI2 is calculated as For VI in Equation (1), we used the modified soil-adjusted vegetation index (MSAVI 2 ) [19]. MSAVI 2 is calculated as where ρ is reflectance in the near-infrared (NIR) or red band [83]. By normalizing the contribution of background soil signal to the integrated spectral reflectance, MSAVI 2 is a better indicator of the vegetation signal than the NDVI in land areas, such as much of our study area, where exposed soil can be a significant component of the observed surface [83,84]. For development, application and empirical analysis of the NDPI, we used time-series image data from two satellite sensors, Landsat 5 Thematic Mapper (TM, 30 m resolution, 16-day repeat cycle) and MODIS on the Terra satellite (250 m resolution, daily repeat cycle) to calculate the MSAVI 2 . The overlapping temporal records of Landsat 5 (1984Landsat 5 ( -2013 and MODIS Terra (2000-present) support comparable results. The red and NIR bands correspond to spectral wavelengths of 630-690 nm and 760-900 nm, respectively, for TM, and 620-670 nm and 841-876 nm, respectively, for MODIS.
The conceptual foundation and our initial development and confirmation of the DVP method and NDPI is based on our Landsat data analysis. Landsat's 30-m resolution facilitated confident identification of pixels to represent the reference vegetation cover conditions and NDPI ref . However, the geographic extent of our Landsat analysis was limited to the overlapping area of two adjacent Landsat path-rows ( Figure 1) where, for years 2009 and 2011, sufficient cloud-free observations were available to capture seasonal patterns of landscape phenology with reasonable fidelity. We used 15 years (2000-2014) of 250-m resolution, 8-day composite data from MODIS to calculate the NDPI over the full extent of our study area. We evaluated the performance of our DVP method by statistical analysis of the relationships between the NDPI s and categorical estimates of perennial C 4 grass cover from field observations (produced in this research) and field-based estimates of C 4 grass cover along an elevation gradient in a sky island mountain range [85]. We also used numerical modeling to examine the behavior of the NDPI s in responding to variations in perennial C 4 grass cover in hypothetical combinations with four major vegetation types that occur in the study area: shrub, mesquite, riparian cottonwood, and pinyon-juniper-evergreen oak.

Potential Variation in the NDPI-Vegetation Cover Relationship
In the case of the spectral vegetation index, NDVI, the relationship between the NDVI and the subpixel fraction of vegetation cover is sensitive to the brightness of the exposed soil (the sum of the red and NIR reflectance values), and the maximum sensitivity (greatest curvilinearity in the relationship) occurs when the soil component of the surface cover is very dark [86,87]. Likewise, in the temporal context of the NDPI, we expect the relationship between NDPI and the subpixel fractional cover of the target vegetation type varies with the sum of VI(p1) and VI(p2), and the maximum sensitivity occurs at the lowest cumulative values of VI(p1) and VI(p2).
We consider here two distinctly different, potential sources of variation in the satellite-derived VI that could affect the sensitivity (curvilinearity of the relationship) of the NDPI to subpixel fractional cover of the target vegetation (perennial C 4 grass in our study). The first source includes year-to-year differences in weather or disturbance that cause differences in vegetation surface cover, biomass, or greenness. It also includes potential noise in the satellite-derived reflectance data that are unrelated to surface and vegetation conditions (including atmospheric contamination and sensor calibration drift). The second source is relevant to applications of the NDPI that differ significantly from those of the present study with regard to the target vegetation cover type and/or the vegetation, exposed soil and phenology characteristics of the study area. In such applications, differences in the characteristic (for example, multiple-year average) magnitudes and dynamics of VI(p1) and VI(p2) of the vegetation components (both the target and non-target vegetation types and exposed soil) relative to those of our study could cause the NDPI relation to fractional cover of the target vegetation type to vary.
Following the method described by Hanan et al. [86] in their analysis of the NDVI, we used linear mixture modeling to perform a sensitivity analysis of how the NDPI responds to variation in the cumulative magnitude of MSAVI 2 (p1) and MSAVI 2 (p2). The response of NDPI to subpixel grass cover (not shown) becomes increasingly curvilinear as MSAVI 2 values of the non-grass surface components during p1 and p2 decline below the 15-year values observed in our study area and toward the lowest potential values. This systematic change in the NDPI relationship is similar to the change in the NDVI relationship in the spectral domain when the exposed soil reflectance becomes very dark [86]. However, for the actual, observed range of MSAVI 2 (p1) and MSAVI 2 (p2) average values in our study area, the sensitivity of the NDPI-grass cover relationship is modest: NDPI varies only within 0%-18% of the NDPI (0.40) associated with the average observed MSAVI 2 (p1) and MSAVI 2 (p2) in the study area. Because our sensitivity analysis does not represent all conditions of vegetation, exposed soil, and phenology that would be potentially encountered in other NDPI applications, we recommend incorporating a sensitivity analysis in applications that differ substantially from our present study in the Southwest.

Landsat Data and Analysis
Our initial application of the DVP-based method with the NDPI used Landsat 5 TM data for a sample area located in the southeastern part of the study area (Figures 1 and 5). We used surface reflectance data [88] and selected cloud-free or nearly cloud-free images over the Landsat sample area for analysis. We masked cloud and cloud shadow-affected pixels from each individual image by differencing the scene with clouds from the nearest cloud-free image date and manually editing the resultant image. All cloud or cloud shadow-affected pixel locations were masked from the entire temporal image stack. Simple linear interpolation from the temporally adjacent scenes was used when an image was removed because of excessive cloud contamination.
The sample area covers approximately 1772 km 2 (6%) of the study area, and is located entirely within the overlap of two horizontally adjacent Landsat World Reference System path-rows (paths 35 and 36, row 38). Landsat observes overlap areas twice as frequently as the non-overlap areas of a scene (8 days average vs. 16 days), which effectively diminishes the tradeoff between spatial resolution and observation frequency that is typically associated with Landsat-based applications [89,90]. Consequently, in the overlap areas landscape phenology can be captured at 30 m resolution with greater confidence than is possible elsewhere in a Landsat scene. By taking advantage of the overlap, we obtained cloud-free observations for 24 dates in 2009 and 25 dates in 2011. We selected years 2009 and 2011 for the Landsat analysis because they yielded the highest number of usable Landsat images from among the years since the MODIS data record began (2000). We deemed these years to be most effective for capturing seasonal patterns of landscape phenology and informing our MODIS-based analysis and mapping application, described henceforth. grassland to desert scrub. We evaluated the validity of this expectation, and its consistency with the idealized response in Figure 4, by examining the divergence of the phenology pattern of 10 NDPI classes from our reference phenology pattern.

Numerical Modeling of NDPIs Response to Variation in Perennial C4 Grass Cover
We used linear mixture modeling to simulate the behavior of NDPIs under a range of vegetation composition and surface cover conditions. The NDPIs was computed for a set of hypothetical Landsat pixels consisting of perennial C4 grass with fractional cover varying from 0.0 to 1.0 in increments of We used the 2009 and 2011 Landsat data to calculate and map the NDPI within the Landsat study area ( Figure 5). To specify the reference phenology and NDPI ref for this area, we selected pixel locations with the highest NDPI values as candidates to represent dense, essentially continuous, perennial C 4 grass cover. We verified that the candidate pixels exhibited the expected phenology of perennial C 4 grass by examining their time-series MSAVI 2 patterns and rejecting agricultural and other non-conforming pixels. The selected pixels were located within the Las Cienegas National Conservation Area (LCNCA) (Figure 5), which is recognized for its extensive perennial warm-season grasslands [91,92]. We further verified the suitability of these locations in our photographic field survey (described subsequently in Section 2.9.1). We used the time-series of average MSAVI 2 values for selected pixels to define the reference phenology patterns and NDPI ref values for 2009 and 2011, which indicate the seasonal growth response of a nominally pure, perennial C 4 grassland within our Landsat study area to the specific weather conditions in these years.
The next step in our Landsat analysis was designed to provide a graphical demonstration and confirmation of the systematic changes in seasonal phenology patterns that occur when the vegetation cover in an observed landscape deviates from the vegetation cover associated with our reference phenology, a nominally pure, dense, and continuous cover of perennial C 4 grass. To achieve this objective, we calculated the NDPI for all pixel locations within our Landsat analysis area, then applied an unsupervised classification to the resulting NDPI image to generate a set of 10 NDPI classes, ranging from the lowest NDPI values (class 1) to the highest NDPI values (class 10). Our DVP method is predicated on the expectation that the unsupervised classification of NDPI data would aggregate pixels with similar surface cover of perennial C 4 grasses, and thereby account for differences in perennial C 4 grass cover across a range of vegetation compositions, from pure grassland to desert scrub. We evaluated the validity of this expectation, and its consistency with the idealized response in Figure 4, by examining the divergence of the phenology pattern of 10 NDPI classes from our reference phenology pattern.

Numerical Modeling of NDPI s Response to Variation in Perennial C 4 Grass Cover
We used linear mixture modeling to simulate the behavior of NDPI s under a range of vegetation composition and surface cover conditions. The NDPI s was computed for a set of hypothetical Landsat pixels consisting of perennial C 4 grass with fractional cover varying from 0.0 to 1.0 in increments of 0.1. Four simulations were performed in which the remaining surface cover was represented by either shrubs, mesquite, riparian cottonwood, or pinyon-juniper-evergreen oak vegetation. Performing this simulation required MSAVI 2 values for p1 and p2 that represent pure samples (100% surface cover) of each vegetation type. We estimated these values by sampling from the 2009 Landsat 5-derived MSAVI 2 images at locations exhibiting maximum surface cover of each vegetation type. We were unable to identify confidently any locations with 100% cover of the individual vegetation types. To overcome this limitation, we identified within the Landsat study area pixel locations with the highest MSAVI 2 values associated with each of the five vegetation types based on our a priori knowledge of the study area. We verified the vegetation type in these locations by visual reference to high-resolution Google Earth imagery.
We estimated the cover fraction of the selected shrub site to be approximately 0.3, and approximately 0.8 for the selected sites of perennial C 4 grass, mesquite, cottonwood, and pinyon-juniper-evergreen oak. To convert the MSAVI 2 values in p1 and p2 for these locations to values that represent 100% cover, we divided the observed MSAVI 2 values by these estimated surface cover values to produce adjusted MSAVI 2 values (Table 2). We then used the adjusted MSAVI 2 values (representing 100% surface cover) in our simulation of how NDPI s responds to variation in perennial C 4 grass cover in combination with each of the four other vegetation types.

MODIS Data and Analysis
Whereas our initial development of the NDPI with Landsat was limited in both spatial extent (overlapping area of adjacent path-rows) and time (two years with abundant cloud-free observations), the wide swath (2330 km) and frequent (daily) revisit cycle of MODIS provided sufficient observations for us to apply the method over multiple years and across the full extent of our study area. We used a 15-year (2000-2014) time-series of 8-day composite MODIS images with 250 m resolution. Our purpose in using multiple years of MODIS data was to reduce the potential for bias in our NDPI-based mapping results that could result from interannual variations in weather (primarily rainfall) on the observed seasonal MSAVI 2 pattern, and hence on the NDPI [93]. Such year-to-year, weather-related variations are apparent in the Landsat results for 2009 and 2011, as discussed subsequently in Section 3.1. Sensitivity of the NDPI to weather effects on a specific vegetation component of mixed-composition landscapes provides an additional, unique source of information from the DVP method. This information may be exploited to detect and monitor at the subpixel scale the response of specifically the perennial C 4 grass component of mixed-vegetation landscapes in the Southwest to drought and climate change. Results from our investigation of the monitoring potential of the NDPI will be reported in a separate, forthcoming paper. For the present research in which our focus is mapping (rather than monitoring), we purposely removed weather-induced interannual variations in the MODIS-derived NDPI by computing a 15-year average NDPI image. We excluded from our analysis all image pixels located above 2134 m (7000 ft) to avoid effects of orographic cloud cover (particularly during the monsoon period) and snow cover during winter months. Although our initial processing and analysis of the MODIS time-series data included all months, we ultimately excluded snow-prone months in our specification of p1 and p2 for modeling of the NDPI.
The major steps in our processing of the MODIS 250-m, 8-day composite surface reflectance data (MOD09Q1 product) [94] are as follows. We excluded MODIS images with substantial, residual contamination of pixels by clouds and cloud shadows based on a visual comparison of false-color images for each weekly composite to the nearest cloud-free image in the time-series. This was performed by rapidly alternating the display of a reference cloud-free image and a target image. Although the temporal compositing procedure that generates the 8-day MODIS data product is intended to remove cloud-related artifacts that can occur in the individual, daily MODIS images, in our study area we found it was not consistently successful. Our visual-based approach to detecting such artifacts enabled us to identify and remove images with significant, residual cloud-related contamination.
We used data from MODIS band 1 (red) and band 2 (near-infrared) to generate a 15-year time-series of MSAVI 2 images representing remaining 8-day MODIS composite periods. To fill gaps in the MSAVI 2 image time-series associated with excluded MODIS images, we performed simple linear interpolation between corresponding pixels in the remaining, temporally adjacent MSAVI 2 images. We then calculated on a per-pixel basis the average MSAVI 2 values in each year for the seasonal periods specified for p1 and p2. With these yearly MSAVI 2 data, we calculated (Equation (1)) a set of 15 yearly NDPI images. We then computed a 15-year average NDPI image and used it to produce (Equation (2)) an NDPI s image that represents the 2000-2014 period. We produced a second pair of 15-year image data sets consisting of yearly average MSAVI 2 values for p1 and yearly maximum MSAVI 2 values for p2. We used these data to generate three additional 15-year average MSAVI 2 -based image products, described below, for comparison with the 15-year average NDPI s .
The remaining procedures for our MODIS-based application of the DVP method generally follow those of our Landsat-based application, however we extended the analysis to the spatial domain to evaluate the effectiveness of the NDPI s for mapping perennial C 4 grass cover across the full range of vegetation types, composition, and surface cover in our study area. We used the 15-year average NDPI image to identify and select an area to represent the reference seasonal phenology pattern associated with nominally pure, dense, continuous cover of perennial C 4 grasses. Based on our conceptual model of how NDPI changes in response to variation in perennial C 4 grass cover ( Figure 4) and results from our Landsat analysis (presented in Section 3.1) that confirm this model, we expected the highest observed NDPI values would be associated with the highest fractional cover. However, high average NDPI values also occurred in areas dominated by irrigated agriculture, primarily northwest of Tucson. Because the seasonal growth of irrigated crops is decoupled from precipitation-related constraints on moisture availability, their phenology can deviate substantially from the characteristic phenology of natural vegetation upon which our selections of p1 and p2 are based. In determining our reference phenology pattern, we therefore excluded these irrigated land areas from consideration, and limited the scope to areas with substantial rain-fed grasslands as determined from our a priori knowledge of the study area.
We selected the average phenology pattern and NDPI observed within the San Rafael Valley in southern Arizona to represent our reference phenology, and the associated NDPI ref value, for perennial C 4 grassland in our study area. The grasslands that cover the gentle, rolling hills of the San Rafael Valley are among the most expansive, pristine and productive grasslands of the Southwest [95]. As the setting for the 1955 musical film production, "Oklahoma!" [96], the San Rafael Valley grasslands became established as an iconic (if not a geographically accurate) representation, in widescreen Technicolor, of the North American shortgrass prairie. Following Equation (2), we generated an NDPI s image by subtracting our specified NDPI ref value from all pixels in the NDPI obs image.
To support a comparative analysis of the performance of our DVP method in estimating and mapping the perennial C 4 grass cover relative to other methods, we generated images for three additional VI-based metrics. These metrics represent potential, alternate approaches to grassland mapping that also exploit temporal information in time-series MODIS observations. The selected metrics are the 15 where, for our purposes, "spring" and "summer" refer respectively to May through June, and July through October. The July-October period corresponds to the principal growth period of perennial C 4 grasses (p2 in our parameterization of Equation (1)) and was selected for its potential to maximize sensitivity to specifically this vegetation type. The summer average MSAVI 2 and summer maximum MSAVI 2 are examples of conventional VI metrics for detecting differences in vegetation phenology, surface cover, green biomass, and production. The difference of the summer maximum MSAVI 2 and spring average MSAVI 2 was selected to represent the DVP-based method used by Helman et al. [56] to discriminate sub-pixel compositions of evergreen and herbaceous vegetation in Mediterranean forests. We evaluated the relative performance of these three satellite-based VI metrics, and our DVP-based method using the NDPI s , in estimating perennial C 4 grass cover by comparisons with field-based observations at sample sites within the study area.

Field Observations for Validation of Map Results
The preferred approach to validating the performance of our NDPI-based application of the DVP method would be through direct comparisons with an independent set of quantitative, spatial data on the distribution and fractional surface cover of specifically perennial C 4 grasses across the study area. However, we were unable to determine the existence and availability of such data from a search of the research literature. Producing a new data set with precise estimates of fractional cover with geographic coverage similar to our study area was beyond the scope of the present research. As an alternative, we used two independent, field-based data sets for validation. As part of this research, we produced a data set that identifies general categories of perennial C 4 grass cover at sample locations distributed throughout the study area. The other data set consists of estimates of the relative surface cover of C 4 grass species at field plots in the Huachuca Mountains reported by an earlier study [85].

Data from Photographic Field Survey
We used ground-based, natural-color, digital photography to record and qualitatively assess the presence and surface cover of perennial C 4 grass at selected field sites in the study area. During a field campaign in August of 2015, we traversed by vehicle all quadrants of the study area and collected natural-color digital photographs at 92 locations. We recorded the geographic coordinates of each site with a handheld GPS device (Oregon 650t, Garmin International, Inc., Olathe, KS, USA). Practical considerations, including site accessibility and time efficiency, precluded a random sample-based approach to site selection. Forest and woodland vegetation, which are limited to upper elevations of the sky island mountain ranges of our study area, were not represented in the survey. We assume the landscape conditions captured in our photographic survey of the 92 sites in 2015 provide a reasonable representation of the average, overall range of vegetation and grass cover conditions, from pure desert scrub to pure grassland, during our 2000-2014 study period.
We developed a portable camera system for efficient, mobile operation and used it to acquire landscape photographs at the field sites. The system included a digital single-lens reflex (DSLR) camera (Canon model EOS 50D with EF-S 18-55 mm lens) attached to the top of a 3.2 m vertical mast. The mast was manually deployed at each site and oriented to capture high-oblique, natural-color images of the adjacent landscape. The camera was connected to a portable computer located in the rear cargo area of the field vehicle and controlled remotely with a software application (DSLR Remote Pro for Windows, Breeze Systems Ltd., Camberley, UK). Image files were recorded in raw and JPEG formats.
Through visual examination of the natural-color images acquired by the DSLR camera, we assessed the vegetation cover characteristics of the landscape at each field site with regard to the presence and relative surface cover of perennial C 4 grass. We assigned each of the 92 sites to one of five qualitative categories: 1.
dense, continuous grass cover (n = 12) The locations of the 92 field sites and sample photographs representing each of the five categories are displayed in Figure 6. To enable direct, quantitative comparisons among the results from the four satellite-based metrics, we normalized the metric data values extracted from the image products at the 92 photographic survey locations.  Because our analysis of the landscape photographs could not yield reliable, precise, quantitative estimates of fractional cover, each set of sites assigned to a qualitative landscape category represents an undetermined amount of variation. Precise estimates of fractional cover could potentially be acquired through more extensive and detailed field studies than our photographic survey, or through analysis of extremely high-resolution (1-5 cm) multispectral or hyperspectral imagery from an airborne sensor, however these options could not be realized within the practical constraints of the present research project. Consequently, the uncertainty inherent in the five qualitative categories reduces the quantitative information available to support validation. We acknowledge this uncertainty, as well as potential uncertainties associated with using single-year observations (2015) to represent the 15-year average (2000-2014) conditions, and interpret our validation results with appropriate caution.
Our use of the categorical data for validation assumes that the individual sets of data values assigned to the landscape categories have distinct, although indeterminate, central tendencies that are separable in data space. On this basis, we use linear least-squares regression to quantify the strength of the relationships between the normalized satellite-based metrics and the landscape category numbers. We refer to the regression results to evaluate the relative performance of the satellite-based metrics in predicting variation in the central tendency among the landscape categories.
As an additional test of relative performance, we quantified and compared the overlapping ranges of the normalized data among the five landscape categories. Results from the overlap analysis enable us to compare the relative success of the satellite-based metrics in accounting for the unique, albeit quantitatively indeterminate information on perennial C4 grass cover represented by the set of sample locations associated with each landscape category. We expect the most reliable indicator of perennial C4 grass cover is the metric that exhibits the least amount of overlap among the landscape categories. Overlap calculations were performed for intervals of 1, 2 and 3 from each landscape category. In this analysis scheme, an interval of 1 accounts for overlap between categories 1 and 2, 2 and 3, 3 and 4, and 4 and 5. An interval of 2 accounts for overlap between categories 1 and 3, 2 and 4, Because our analysis of the landscape photographs could not yield reliable, precise, quantitative estimates of fractional cover, each set of sites assigned to a qualitative landscape category represents an undetermined amount of variation. Precise estimates of fractional cover could potentially be acquired through more extensive and detailed field studies than our photographic survey, or through analysis of extremely high-resolution (1-5 cm) multispectral or hyperspectral imagery from an airborne sensor, however these options could not be realized within the practical constraints of the present research project. Consequently, the uncertainty inherent in the five qualitative categories reduces the quantitative information available to support validation. We acknowledge this uncertainty, as well as potential uncertainties associated with using single-year observations (2015) to represent the 15-year average (2000-2014) conditions, and interpret our validation results with appropriate caution.
Our use of the categorical data for validation assumes that the individual sets of data values assigned to the landscape categories have distinct, although indeterminate, central tendencies that are separable in data space. On this basis, we use linear least-squares regression to quantify the strength of the relationships between the normalized satellite-based metrics and the landscape category numbers. We refer to the regression results to evaluate the relative performance of the satellite-based metrics in predicting variation in the central tendency among the landscape categories.
As an additional test of relative performance, we quantified and compared the overlapping ranges of the normalized data among the five landscape categories. Results from the overlap analysis enable us to compare the relative success of the satellite-based metrics in accounting for the unique, albeit quantitatively indeterminate information on perennial C 4 grass cover represented by the set of sample locations associated with each landscape category. We expect the most reliable indicator of perennial C 4 grass cover is the metric that exhibits the least amount of overlap among the landscape categories. Overlap calculations were performed for intervals of 1, 2 and 3 from each landscape category. In this analysis scheme, an interval of 1 accounts for overlap between categories 1 and 2, 2 and 3, 3 and 4, and 4 and 5. An interval of 2 accounts for overlap between categories 1 and 3, 2 and 4, and 3 and 5. An interval of 3 accounts for overlap between categories 1 and 4, and 2 and 5. None of the satellite-based metrics exhibited overlap between landscape categories 1 and 5, therefore we neglected the 4-category interval in our analysis of overlap.
The complete set of data and information from the field survey, including the location, category designation, and landscape image for each site are provided in the Supplementary Materials (S1).

Data from Field Plots in a Sky Island Mountain Range
The occurrence and relative surface cover of C 4 grass species at 68 sample plots across an elevation gradient on the northeastern flank of the Huachuca Mountains, a sky island mountain range in the southeastern part of our study area was reported for prior study by Wentworth [85]. The data are summarized as average values representing the 10 contiguous, 50-m intervals from 1600 to 2100 m elevation, and are reported in a line plot (Figure 2d in [85]). We used linear interpolation to estimate from visual inspection of the line plot the numerical values of relative cover of C 4 grasses. For comparison, we extracted NDPI s values at all grid cell locations on the northeastern side of the Huachuca Mountains and computed the average NDPI s values within the same elevation intervals as those used by Wentworth [85]. Because the Wentworth [85] data were collected nearly three decades before our MODIS observation period, we exercised caution by considering only elevation-related trends, rather than absolute values of relative cover.

Landsat-Based Results
The Landsat-derived reference phenology (seasonal MSAVI 2 curves) for 2009 and 2011 are consistent with expectations for landscape conditions characterized by nominally pure, continuous, perennial C 4 grass cover (Figure 7). The MSAVI 2 remains low and relatively stable throughout the winter and spring period when the perennial C 4 grasses are predominantly dormant. A significant feature of the general pattern evident in both years is the rapid rise in the MSAVI 2 after approximately day 176 (25 June). This abrupt increase reflects the increase in photosynthesis and productivity by perennial C 4 grasses with the onset of precipitation from the North American Monsoon and concurrent warm temperatures. Differences in the time patterns for 2009 and 2011 are attributable primarily to differences in the amount and timing of precipitation (data not shown).
Remote Sens. 2016, 8, 889 18 of 32 and 3 and 5. An interval of 3 accounts for overlap between categories 1 and 4, and 2 and 5. None of the satellite-based metrics exhibited overlap between landscape categories 1 and 5, therefore we neglected the 4-category interval in our analysis of overlap. The complete set of data and information from the field survey, including the location, category designation, and landscape image for each site are provided in the Supplementary Materials (S1).

Data from Field Plots in a Sky Island Mountain Range
The occurrence and relative surface cover of C4 grass species at 68 sample plots across an elevation gradient on the northeastern flank of the Huachuca Mountains, a sky island mountain range in the southeastern part of our study area was reported for prior study by Wentworth [85]. The data are summarized as average values representing the 10 contiguous, 50-m intervals from 1600 to 2100 m elevation, and are reported in a line plot (Figure 2d in [85]). We used linear interpolation to estimate from visual inspection of the line plot the numerical values of relative cover of C4 grasses. For comparison, we extracted NDPIs values at all grid cell locations on the northeastern side of the Huachuca Mountains and computed the average NDPIs values within the same elevation intervals as those used by Wentworth [85]. Because the Wentworth [85] data were collected nearly three decades before our MODIS observation period, we exercised caution by considering only elevationrelated trends, rather than absolute values of relative cover.

Landsat-Based Results
The Landsat-derived reference phenology (seasonal MSAVI2 curves) for 2009 and 2011 are consistent with expectations for landscape conditions characterized by nominally pure, continuous, perennial C4 grass cover (Figure 7). The MSAVI2 remains low and relatively stable throughout the winter and spring period when the perennial C4 grasses are predominantly dormant. A significant feature of the general pattern evident in both years is the rapid rise in the MSAVI2 after approximately day 176 (25 June). This abrupt increase reflects the increase in photosynthesis and productivity by perennial C4 grasses with the onset of precipitation from the North American Monsoon and concurrent warm temperatures. Differences in the time patterns for 2009 and 2011 are attributable primarily to differences in the amount and timing of precipitation (data not shown). Deviations in the Landsat-derived MSAVI2 values for NDPI classes labeled 3 through 10 from the corresponding MSAVI2 values that define our reference phenology patterns exhibit significant features that constitute the conceptual and practical underpinnings of the NDPI (Figure 8). These features, which we assert convey unique, quantitative information on the fractional surface cover of Deviations in the Landsat-derived MSAVI 2 values for NDPI classes labeled 3 through 10 from the corresponding MSAVI 2 values that define our reference phenology patterns exhibit significant features that constitute the conceptual and practical underpinnings of the NDPI (Figure 8). These features, which we assert convey unique, quantitative information on the fractional surface cover of perennial C 4 grasses in our study area, would not have been produced from an unsupervised classification of the image date layer-stack itself. The features and the information they contain emerge only from calculation and classification of an NDPI image or another way of initially stratifying the data by grass cover.
The deviations displayed in Figure 8 are positive in the winter/spring and negative in the summer. This pattern is consistent with our expectations from a priori knowledge of the phenology of the major vegetation components of the region as discussed previously and expressed in idealized form in Figure 4. We attribute the positive deviations in winter/spring to the MSAVI 2 of perennial C 4 grasses (which are dormant) being lower than the MSAVI 2 of other vegetation types in the landscape, which are green during part or all of this period (trees, shrubs, and/or succulents). Conversely, we attribute the negative deviations in summer to the MSAVI 2 values for perennial C 4 grasses being higher than the values of the other vegetation types during this period when the photosynthetic productivity and green leaf display of C 4 grasses are approaching or at their maximum.
Remote Sens. 2016, 8, 889 19 of 32 perennial C4 grasses in our study area, would not have been produced from an unsupervised classification of the image date layer-stack itself. The features and the information they contain emerge only from calculation and classification of an NDPI image or another way of initially stratifying the data by grass cover. The deviations displayed in Figure 8 are positive in the winter/spring and negative in the summer. This pattern is consistent with our expectations from a priori knowledge of the phenology of the major vegetation components of the region as discussed previously and expressed in idealized form in Figure 4. We attribute the positive deviations in winter/spring to the MSAVI2 of perennial C4 grasses (which are dormant) being lower than the MSAVI2 of other vegetation types in the landscape, which are green during part or all of this period (trees, shrubs, and/or succulents). Conversely, we attribute the negative deviations in summer to the MSAVI2 values for perennial C4 grasses being higher than the values of the other vegetation types during this period when the photosynthetic productivity and green leaf display of C4 grasses are approaching or at their maximum. Differences among the NDPI classes in the magnitude of their MSAVI2 deviations in both winter/spring and summer exhibit a systematic pattern. In this application, increasing positive deviations in winter/spring are associated with increasing negative deviations in summer. We assert that this pattern of change in seasonal deviations of the MSAVI2 from the reference phenology pattern, which corresponds to decreasing values of the NDPI, is indicative of decreasing surface cover of perennial C4 grass cover within the land area captured by an individual pixel. In our subsequent analysis, we evaluate the validity of this assertion on a theoretical basis by examining the results of our numerical modeling. We then evaluate it on an empirical basis by comparing 15-year averages of NDPIs from MODIS with results of categorical estimates of perennial C4 grass cover from field observations.

Numerical Modeling Results
The simulation results show mild, non-linear relationships between NDPIs and subpixel fractional cover for each vegetation-type combination (Figure 9). The regression lines diverge slightly from a common point corresponding to 100% perennial grass cover, with maximum divergence at values of fractional cover that vary with the curvilinearity of the regression lines. The differences among the relationships at their points of maximum divergence are 7%-12% of the average range of NDPIs, and have an average difference of 7% for all cases of simulated fractional cover. A more detailed and realistic sensitivity analysis could be performed in future research by using, for example, a geometric-optical canopy radiative transfer model instead of linear mixture modeling. Nevertheless, these results reinforce our expectation that the NDPIs can be effective in discriminating the presence and subpixel surface cover of perennial C4 grass throughout our study area. However, the absence of a single, universal relationship between NDPIs and perennial C4 grass cover suggests differences in the phenology of the non-grass vegetation components are a potential source of Differences among the NDPI classes in the magnitude of their MSAVI 2 deviations in both winter/spring and summer exhibit a systematic pattern. In this application, increasing positive deviations in winter/spring are associated with increasing negative deviations in summer. We assert that this pattern of change in seasonal deviations of the MSAVI 2 from the reference phenology pattern, which corresponds to decreasing values of the NDPI, is indicative of decreasing surface cover of perennial C 4 grass cover within the land area captured by an individual pixel. In our subsequent analysis, we evaluate the validity of this assertion on a theoretical basis by examining the results of our numerical modeling. We then evaluate it on an empirical basis by comparing 15-year averages of NDPI s from MODIS with results of categorical estimates of perennial C 4 grass cover from field observations.

Numerical Modeling Results
The simulation results show mild, non-linear relationships between NDPI s and subpixel fractional cover for each vegetation-type combination (Figure 9). The regression lines diverge slightly from a common point corresponding to 100% perennial grass cover, with maximum divergence at values of fractional cover that vary with the curvilinearity of the regression lines. The differences among the relationships at their points of maximum divergence are 7%-12% of the average range of NDPI s , and have an average difference of 7% for all cases of simulated fractional cover. A more detailed and realistic sensitivity analysis could be performed in future research by using, for example, a geometric-optical canopy radiative transfer model instead of linear mixture modeling. Nevertheless, these results reinforce our expectation that the NDPI s can be effective in discriminating the presence and subpixel surface cover of perennial C 4 grass throughout our study area. However, the absence of a single, universal relationship between NDPI s and perennial C 4 grass cover suggests differences in the phenology of the non-grass vegetation components are a potential source of uncertainty. A priori knowledge of the seasonal variation in the spectral reflectance properties of the major vegetation components of the landscape is therefore useful for reliable application and interpretation of the NDPI s . Remote Sens. 2016, 8, 889 20 of 32 uncertainty. A priori knowledge of the seasonal variation in the spectral reflectance properties of the major vegetation components of the landscape is therefore useful for reliable application and interpretation of the NDPIs.

MODIS-Based Map Results
The spatial distribution of the 15-year average summer maximum MSAVI2 ( Figure 10) generally corresponds to variation in topographic elevation (Figure 2). The lowest MSAVI2 values are observed in the most arid and sparsely vegetated lowlands. The highest MSAVI2 values are seen at the highest elevations within the sky island mountain ranges where woodland and forest vegetation are the dominant cover type. Significant exceptions are the anomalously high values apparent at low elevations where irrigated agriculture is the dominant land use (for example, northwest and immediately south of Tucson). Anomalous, high MSAVI2 values are also evident along riparian zones, where soil moisture is sufficient to support dense tree vegetation, and in some locations, agricultural crops. The riparian zone anomalies are evident throughout the study area, but are prominent along the Santa Cruz River (south of Tucson) and along the San Pedro River (extending north-south across the eastern part of the study area).
The overall spatial patterns displayed in Figure 10, including those of the anomalies, suggest the 15-year average maximum MSAVI2 is generally sensitive to variations in total aboveground biomass, green leaf area, and primary production. Such a relationship is consistent with expectations for a satellite-based, seasonal maximum, average, or integrated VI [97,98]. However, the patterns do not indicate the average summer maximum MSAVI2 is effective for distinguishing, in terms of either biomass, leaf area or surface cover, perennial C4 grass from the other vegetation types. This result highlights shortcomings of a basic, conventional VI-based approach to mapping specifically the perennial C4 grass component of mixed-composition landscapes in the Southwest.

MODIS-Based Map Results
The spatial distribution of the 15-year average summer maximum MSAVI 2 ( Figure 10) generally corresponds to variation in topographic elevation (Figure 2). The lowest MSAVI 2 values are observed in the most arid and sparsely vegetated lowlands. The highest MSAVI 2 values are seen at the highest elevations within the sky island mountain ranges where woodland and forest vegetation are the dominant cover type. Significant exceptions are the anomalously high values apparent at low elevations where irrigated agriculture is the dominant land use (for example, northwest and immediately south of Tucson). Anomalous, high MSAVI 2 values are also evident along riparian zones, where soil moisture is sufficient to support dense tree vegetation, and in some locations, agricultural crops. The riparian zone anomalies are evident throughout the study area, but are prominent along the Santa Cruz River (south of Tucson) and along the San Pedro River (extending north-south across the eastern part of the study area).
The overall spatial patterns displayed in Figure 10, including those of the anomalies, suggest the 15-year average maximum MSAVI 2 is generally sensitive to variations in total aboveground biomass, green leaf area, and primary production. Such a relationship is consistent with expectations for a satellite-based, seasonal maximum, average, or integrated VI [97,98]. However, the patterns do not indicate the average summer maximum MSAVI 2 is effective for distinguishing, in terms of either biomass, leaf area or surface cover, perennial C 4 grass from the other vegetation types. This result highlights shortcomings of a basic, conventional VI-based approach to mapping specifically the perennial C 4 grass component of mixed-composition landscapes in the Southwest. The spatial distribution of the 15-year average NDPIs ( Figure 11) differs significantly from the summer-maximum MSAVI2 pattern ( Figure 10) along the upper elevation gradients of the sky island mountains (c.f. Figure 2) and in riparian zones. Whereas the summer maximum MSAVI2 generally increases with increasing elevation within the sky island mountain ranges, the NDPIs decreases with elevation in these upland areas. These elevation-related patterns are illustrated in data from a transect across the Huachuca Mountains in the southeastern part of the study area ( Figure 12). The patterns along the transect are consistent with the influence of climatic and ecophysiological factors on vegetation distribution along elevation gradients [99,100]. With increasing elevation, cooler temperatures and reduced evaporative demands favor woodland and forest vegetation, while reducing the competitive advantage and viability of perennial C4 grasses. Thus, the decline in NDPIs along the uppermost elevation gradients suggests NDPIs, in contrast to the maximum summer MSAVI2, is sensitive to the expected elevation-related decline of perennial C4 grasses [85].
The spatial pattern of the 15-year average NDPIs also deviates significantly from the 15-year average maximum summer MSAVI2 in the riparian zones. Whereas the summer maximum MSAVI2 exhibits anomalously high values in those areas, the NDPIs does not. Although grasses occur in riparian ecosystems of the semiarid Southwest, woody vegetation, including gallery forests of cottonwood and willow, is typically the dominant cover type [101][102][103]. Whereas the high MSAVI2 values observed in the riparian zones suggest this index is a misleading indicator of the surface cover of perennial C4 grasses, the low NDPIs values observed in these areas conform to our expectations. The overall results from visual comparison of the MODIS-derived map products (Figures 10 and 11) The spatial distribution of the 15-year average NDPI s ( Figure 11) differs significantly from the summer-maximum MSAVI 2 pattern ( Figure 10) along the upper elevation gradients of the sky island mountains (c.f. Figure 2) and in riparian zones. Whereas the summer maximum MSAVI 2 generally increases with increasing elevation within the sky island mountain ranges, the NDPI s decreases with elevation in these upland areas. These elevation-related patterns are illustrated in data from a transect across the Huachuca Mountains in the southeastern part of the study area ( Figure 12). The patterns along the transect are consistent with the influence of climatic and ecophysiological factors on vegetation distribution along elevation gradients [99,100]. With increasing elevation, cooler temperatures and reduced evaporative demands favor woodland and forest vegetation, while reducing the competitive advantage and viability of perennial C 4 grasses. Thus, the decline in NDPI s along the uppermost elevation gradients suggests NDPI s , in contrast to the maximum summer MSAVI 2 , is sensitive to the expected elevation-related decline of perennial C 4 grasses [85].
The spatial pattern of the 15-year average NDPI s also deviates significantly from the 15-year average maximum summer MSAVI 2 in the riparian zones. Whereas the summer maximum MSAVI 2 exhibits anomalously high values in those areas, the NDPI s does not. Although grasses occur in riparian ecosystems of the semiarid Southwest, woody vegetation, including gallery forests of cottonwood and willow, is typically the dominant cover type [101][102][103] in the riparian zones suggest this index is a misleading indicator of the surface cover of perennial C 4 grasses, the low NDPI s values observed in these areas conform to our expectations. The overall results from visual comparison of the MODIS-derived map products (Figures 10 and 11) support our expectation that the DVP method using NDPI s can provide improved capabilities for discriminating the perennial C 4 grass component from other vegetation components of Southwest landscapes, relative to a simplified, conventional VI-based approach represented here by the summer maximum MSAVI 2 .
Images for both the NDPI s ( Figure 11) and summer-maximum MSAVI 2 ( Figure 10) exhibit high values in areas of irrigated cropland. Such areas appear as the distinct, dark green features (indicating high MSAVI 2 or NDPI s values) in the northwestern part of the study area, directly northwest of Tucson (c.f. Figure 1). Whereas these patterns are a reasonable result for MSAVI 2 with respect to green biomass and primary production, they represent an anomalous and misleading result for the NDPI s with respect to perennial C 4 grass cover. The NDPI s response in irrigated croplands may be attributed to a seasonal cropping pattern that effectively mimics the phenology of the perennial grasses in the region. While this artifact does not diminish the potential effectiveness of NDPI s and the overall DVP method for discriminating perennial C 4 grasses in non-irrigated land areas, it does indicate the need for a priori identification and masking of cropland areas. Conversely, the anomalous results for irrigated agriculture suggest the DVP method has potential for discrimination and mapping of agricultural crop vegetation within mixed-vegetation landscapes.
Remote Sens. 2016, 8,889 22 of 32 support our expectation that the DVP method using NDPIs can provide improved capabilities for discriminating the perennial C4 grass component from other vegetation components of Southwest landscapes, relative to a simplified, conventional VI-based approach represented here by the summer maximum MSAVI2. Images for both the NDPIs ( Figure 11) and summer-maximum MSAVI2 ( Figure 10) exhibit high values in areas of irrigated cropland. Such areas appear as the distinct, dark green features (indicating high MSAVI2 or NDPIs values) in the northwestern part of the study area, directly northwest of Tucson (c.f. Figure 1). Whereas these patterns are a reasonable result for MSAVI2 with respect to green biomass and primary production, they represent an anomalous and misleading result for the NDPIs with respect to perennial C4 grass cover. The NDPIs response in irrigated croplands may be attributed to a seasonal cropping pattern that effectively mimics the phenology of the perennial grasses in the region. While this artifact does not diminish the potential effectiveness of NDPIs and the overall DVP method for discriminating perennial C4 grasses in non-irrigated land areas, it does indicate the need for a priori identification and masking of cropland areas. Conversely, the anomalous results for irrigated agriculture suggest the DVP method has potential for discrimination and mapping of agricultural crop vegetation within mixed-vegetation landscapes. Gray color indicates the data mask applied to snow-prone areas (elevation > 2134 m). White color pixels that appear in sporadic, narrow bands indicate data masks applied to missing or unreliable MODIS data. The diagonal transect and the rectangular area, indicated by the broken lines across the Huachuca Mountains, are where the image data were sampled for comparative analysis.

Validation of MODIS-Based Map Results
We now present and discuss results from comparison of the four satellite-based metrics with data from our photographic survey. The scatter plots displayed in Figure 13 show results, including R 2 and associated p-values, from the regression of data sampled from the four image-based map products at the 92 photographic survey sites against the numerical value of the corresponding landscape category. As discussed earlier, the results require careful interpretation because the regressions were performed with numerical values that correspond to categories (indicating ranges of perennial C4 grass cover), rather than absolute values of cover. Accordingly, we avoid inferences to absolute values of cover, which would be fundamental for a rigorous quantitative validation. In our validation, we consider the regression results as general indicators of relative performance in accounting for the increasing amounts of perennial C4 grass cover represented by the sequence (1 to 5) of the qualitative landscape categories.

Validation of MODIS-Based Map Results
We now present and discuss results from comparison of the four satellite-based metrics with data from our photographic survey. The scatter plots displayed in Figure 13 show results, including R 2 and associated p-values, from the regression of data sampled from the four image-based map products at the 92 photographic survey sites against the numerical value of the corresponding landscape category. As discussed earlier, the results require careful interpretation because the regressions were performed with numerical values that correspond to categories (indicating ranges of perennial C 4 grass cover), rather than absolute values of cover. Accordingly, we avoid inferences to absolute values of cover, which would be fundamental for a rigorous quantitative validation. In our validation, we consider the regression results as general indicators of relative performance in accounting for the increasing amounts of perennial C 4 grass cover represented by the sequence (1 to 5) of the qualitative landscape categories. Results from least-squares linear regression indicate the 15-year averages of the two conventional VI-based metrics are moderately successful in predicting the variation in perennial C4 grass cover represented by the five landscape categories (Figure 13a,b). The coefficient of determination (R 2 ) is 0.56 (p < 0.001) for the summer average MSAVI2 (Figure 13a) and 0. 62 (p < 0.001) for the summer maximum MSAVI2 (Figure 13b). Relative to these results from simple, conventional VI metrics, regression results from the two DVP-based metrics suggest substantially greater sensitivity to the differences in perennial C4 grass cover among the landscape categories ( Figure  13c,d). The R 2 value for the difference of the summer maximum and spring average MSAVI2 is 0.73 (p < 0.001) (Figure 13c). With an R 2 value of 0.76 (p < 0.001), the results suggest the NDPIs is potentially the most successful predictor of perennial C4 grass cover among the four metrics that we compared (Figure 13d).
Our interpretation of the regression results is reinforced by quantitative analysis of the overlap in the ranges of the normalized data among the five landscape categories, which is apparent in the scatter plots in Figure 13. Differences in the amount of overlap among the four satellite-based metrics indicate their relative performance in distinguishing the variation in perennial C4 grass cover represented by the five landscape categories. Results from our overlap calculations (Table 3) are Results from least-squares linear regression indicate the 15-year averages of the two conventional VI-based metrics are moderately successful in predicting the variation in perennial C 4 grass cover represented by the five landscape categories (Figure 13a,b). The coefficient of determination (R 2 ) is 0.56 (p < 0.001) for the summer average MSAVI 2 (Figure 13a) and 0. 62 (p < 0.001) for the summer maximum MSAVI 2 (Figure 13b). Relative to these results from simple, conventional VI metrics, regression results from the two DVP-based metrics suggest substantially greater sensitivity to the differences in perennial C 4 grass cover among the landscape categories (Figure 13c,d). The R 2 value for the difference of the summer maximum and spring average MSAVI 2 is 0.73 (p < 0.001) (Figure 13c). With an R 2 value of 0.76 (p < 0.001), the results suggest the NDPI s is potentially the most successful predictor of perennial C 4 grass cover among the four metrics that we compared (Figure 13d).
Our interpretation of the regression results is reinforced by quantitative analysis of the overlap in the ranges of the normalized data among the five landscape categories, which is apparent in the scatter plots in Figure 13. Differences in the amount of overlap among the four satellite-based metrics indicate their relative performance in distinguishing the variation in perennial C 4 grass cover represented by the five landscape categories. Results from our overlap calculations (Table 3) are consistent with the results from regression analysis. At each interval, the largest overlaps are exhibited by the two conventional VI-based metrics ( Table 3). The average overlap for the combined set of intervals is 49% for the summer average MSAVI 2 , and 44% for the summer maximum MSAVI 2 . Consistent with the regression results, the two DVP-based metrics show substantially improved performance relative to the conventional metrics ( Table 3). The overlap for the combined set of intervals is 26% for the difference of the summer maximum and spring average MSAVI 2 , and 21% for the NDPI s (Table 3). Thus NDPI s exhibits 5% less overlap between categories than the other DVP-based metric and 23%-28% less overlap than the conventional VI metrics. These relative reductions in overlap potentially represent commensurate increases in unique information captured by the NDPI s . The results from the overlap analysis provide additional evidence that the NDPI s is potentially the most effective among the four satellite-based metrics in distinguishing differences in perennial C 4 grass cover among the field sites that represent the landscape categories. Results from comparison of the 15-year average NDPI s with estimates of the relative cover of C 4 grass species across an elevation gradient in the Huachuca Mountains [85] are displayed in Figure 14. The NDPI s decreases monotonically with increasing elevation across the sample range from 1600 to 2100 m. Although the decline in field-based estimates of the relative cover of C 4 grasses along the same elevation gradient is not monotonic, the general pattern corresponds well with the NDPI s pattern (Figure 14a). The additional variability apparent in the field plot-based C 4 grass data may be at least partially attributable to their representation of a limited set of conditions. The area from which the NDPI s data were sampled covers the entire northeastern side of the Huachuca Mountains ( Figure 11) and therefore captures all landscape conditions (including variations in surface slope and aspect). The 68 sites from which the C 4 grass cover data were collected represent a sample set of locations that may not represent the full range of landscape conditions observed within the same elevation intervals over the broader sample area represented by the NDPI s . Regardless, least-squares linear regression confirms a statistically significant relationship between the observed, elevation-related changes in NDPI s and C 4 grass cover (R 2 = 0.70, p < 0.01) (Figure 14b). These results support our interpretation, discussed earlier, of the divergence in the 15-year averages of the NDPI s and the summer maximum MSAVI 2 observed in the elevation transect across the Huachuca Mountains (Figures 10 and 11). In stark contrast to the conventional VI-based metric, the NDPI s captures the expected decline in C 4 grass cover with increasing elevation in the upland areas of the sky island mountain ranges.  (Figures 2 and 11), and displayed (a) as a function of elevation, and (b) as a scatter plot with results from least-squares linear regression. Both data sets represent average values for 50-m elevation intervals. The NDPIs data are from the rectangular area bounded by the dotted line in Figure 11, whereas the relative cover data are from field plots studied by Wentworth [85].

Conclusions
We developed a methodology for mapping of the exposed perennial, warm-season (C4) grass component of mixed-and variable-composition landscapes in the Southwest, and applied the method to a study area in southern Arizona, USA and northern Sonora, Mexico. Our application demonstrates and exemplifies a general, conceptual approach to mapping of an individual vegetation component of vegetation cover in mixed-and variable-composition landscapes, which heretofore has remained relatively unexplored [56]. Insomuch as the conceptual approach embodies a distinct category of remote sensing analytical methods that, to our knowledge, has not been explicitly identified in the remote sensing research literature, we have identified it as a differential vegetation phenology (DVP) approach to vegetation mapping and monitoring.
We distinguished the DVP approach from other vegetation mapping methods by its emphasis on strategic leveraging of a priori information on the seasonal differences in phenology of the vegetation cover components that comprise a landscape. We distinguished the DVP approach from other phenology-based (frequency-domain) methods, and from those that use only spatial and spectral information, by its use of metrics that strategically exploit such predetermined differences in phenology. In our DVP application to mapping of perennial C4 grass in the Southwest, we introduced  (Figures 2 and 11), and displayed (a) as a function of elevation, and (b) as a scatter plot with results from least-squares linear regression. Both data sets represent average values for 50-m elevation intervals. The NDPI s data are from the rectangular area bounded by the dotted line in Figure 11, whereas the relative cover data are from field plots studied by Wentworth [85].

Conclusions
We developed a methodology for mapping of the exposed perennial, warm-season (C 4 ) grass component of mixed-and variable-composition landscapes in the Southwest, and applied the method to a study area in southern Arizona, USA and northern Sonora, Mexico. Our application demonstrates and exemplifies a general, conceptual approach to mapping of an individual vegetation component of vegetation cover in mixed-and variable-composition landscapes, which heretofore has remained relatively unexplored [56]. Insomuch as the conceptual approach embodies a distinct category of remote sensing analytical methods that, to our knowledge, has not been explicitly identified in the remote sensing research literature, we have identified it as a differential vegetation phenology (DVP) approach to vegetation mapping and monitoring.
We distinguished the DVP approach from other vegetation mapping methods by its emphasis on strategic leveraging of a priori information on the seasonal differences in phenology of the vegetation cover components that comprise a landscape. We distinguished the DVP approach from other phenology-based (frequency-domain) methods, and from those that use only spatial and spectral information, by its use of metrics that strategically exploit such predetermined differences in phenology. In our DVP application to mapping of perennial C 4 grass in the Southwest, we introduced one such metric-a temporal vegetation index as opposed to a spectral vegetation index-that we refer to as a normalized difference phenometric index, or NDPI. The combined results from our analysis provide a consistent body of evidence that the DVP methods, and the NDPI s in particular, performs best in mapping perennial C 4 grass in our Southwest study area relative to two conventional VI-based metrics. Desirable goals for future research include additional validation and refinement of the NDPI (and other potential, DVP-based metrics), and quantitative comparisons against additional, existing methods for subpixel discrimination and mapping of the perennial C 4 grass component of mixed-composition landscapes in the Southwest. Because the relationship between NDPI and the subpixel fractional cover of the vegetation cover type of interest is not necessarily constant for all potential vegetation, exposed soil, and phenology conditions, we recommend performing a sensitivity analysis in future NDPI applications that differ substantially from our present study in the Southwest.
In this mapping application in the Southwest, we used 15-year (2000-2014) average data sets from MODIS to reduce the effects of interannual variability in vegetation phenology that are attributable to weather conditions, especially variability in precipitation. The resulting map of the 15-year average NDPI s ( Figure 11) therefore represents the average conditions of perennial C 4 grass cover during the observation period. Whereas the 15-year average NDPI s does not account for potential changes or variability in the spatial distribution or extent of perennial C 4 grass in the study area during 2000-2014, monitoring and analysis of the NDPI s on an annual basis could provide useful information about how perennial C 4 grass in mixed-composition landscapes responds to fire and other disturbances, weather, drought, climate change, and land management strategies (such as grazing regimes). Likewise, analysis of variations and trends in the NDPI s over multiple years may provide an integrated, ecosystem-based indicator of the strength and variability of the North American Monsoon. Our ongoing research is investigating this potential with an aim toward providing improved, timely, satellite-based data products on grassland and other vegetation components that support the information needs of the science and land management communities in the Southwest and beyond.
Acknowledgments: This research was supported by the USGS Land Remote Sensing Program, Climate and Land Use Change Mission Area.
Author Contributions: Dennis Dye and Barry Middleton conceived and designed the research, performed the field campaign, analyzed the data, and prepared the graphics; John Vogel and Dennis Dye configured the mobile imaging system; Miguel Velasco performed MODIS data preprocessing; John Vogel prepared the Supplementary Data File; Barry Middleton, Zhuoting Wu, and John Vogel provided editorial support; Dennis Dye wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: