Monitoring Broadscale Vegetational Diversity and Change across North American Landscapes Using Land Surface Phenology

Abstract

We describe a polar coordinate transformation of vegetation index profiles which permits a broad-scale comparison of location-specific phenological variability influenced by climate, topography, land use, and other factors. We apply statistical data reduction techniques to identify fundamental dimensions of phenological variability and to classify phenological types with intuitive ecological interpretation. Remote sensing-based land surface phenology can reveal ecologically meaningful vegetational diversity and dynamics across broad landscapes. Land surface phenology is inherently complex at regional to continental scales, varying with latitude, elevation, and multiple biophysical factors. Quantifying phenological change across ecological gradients at these scales is a potentially powerful way to monitor ecological development, disturbance, and diversity. Polar coordinate transformation was applied to Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) time series spanning 2000-2018 across North America. In a first step, 46 NDVI values per year were reduced to 11 intuitive annual metrics, such as the midpoint of the growing season and degree of seasonality, measured relative to location-specific annual phenological cycles. Second, factor analysis further reduced these metrics to fundamental phenology dimensions corresponding to annual timing, productivity, and seasonality. The factor analysis explained over 95% of the variability in the metrics and represented a more than ten-fold reduction in data volume from the original time series. In a final step, phenological classes (‘phenoclasses’) based on the statistical clustering of the factor data, were computed to describe the phenological state of each pixel during each year, which facilitated the tracking of year-to-year dynamics. Collectively the phenology metrics, factors, and phenoclasses provide a system for characterizing land surface phenology and for monitoring phenological change that is indicative of ecological gradients, development, disturbance, and other aspects of landscape-scale diversity and dynamics.

1. Introduction

The use of remote sensing to characterize vegetation, map spatial gradients, and monitor temporal change has significantly advanced the field of landscape ecology [1,2,3]. This importantly includes applications using measures of surface reflectance (“greenness”) through time to infer properties of the vegetative land cover. Remotely sensed seasonal variation in vegetation (land surface phenology (LSP)) based on such measures is strongly correlated with annual pulses of leaf-out through senescence [4,5,6]. Given appropriate choices of spatial and temporal resolutions, LSP approaches can be used to systematically track vegetation dynamics including disturbance, recovery, and development [7,8]. Beyond providing phenological information, LSP linkages to a wide variety of ecosystem and landscape properties and processes provide insight into large-scale ecological diversity, distributions, and change [4,9,10,11]. Ground-based approaches exist for observing landscape-level vegetation phenology over local areas, but satellite-based LSP is much broader in scale, often spanning biomes with strongly differing phenologies. These differences are useful for indicating ecological gradients and change, but they also present challenges for systematic and consistent continental-scale assessment and monitoring.

Satellite-based LSP from sources such as the moderate resolution imaging spectroradiometer (MODIS) allows large landscape assessment and change diagnosis [12,13]. These measurements cannot distinguish individual tree-crowns. At medium spatial resolutions (250 m in the case of MODIS), a phenological measurement typically reflects aggregate behaviors of multiple plant assemblages, which can be affected by different environmental factors (e.g., harvested and non-harvested patches), or differential responses within a plant assemblage to a common factor (e.g., species-specific drought responses). If environmental factors are strong enough to affect the aggregate surface reflectance signal, temporal and spatial variation in land surface phenology reflects the underlying dynamics of the ecosystem. Examples include biotic disturbances, such as defoliation from insect outbreaks and mortality [14,15,16], and abiotic disturbances tied to regional climate changes [17].

Due in part to complications from scale representativeness and asynchrony, much of the work in LSP has addressed thematically coarse land cover mapping and disturbance monitoring rather than variability along nuanced ecological gradients such as those associated with within-ecosystem plant species compositional and structural change. By contrast, classification schemes which group image pixels into phenologically self-similar clusters allow for quantifying and mapping phenological variation at arbitrarily fine levels of distinction [7,18,19,20]. Successive satellite measurements of a reliable vegetation indicator such as the normalized difference vegetation index (NDVI) are used to build a phenological profile across a year, and annual profiles are clustered. Clustering may be performed on the successive index values themselves or on derived metrics such as growing season start and length [19]. The process of generating metrics and phenological classes (phenoclasses) can be automated for each year to assess change. This approach has been used as a basis for monitoring vegetation dynamics [7].

The capacity to undertake broadscale monitoring of vegetation as it develops and responds to environmental conditions has been expanded by high-performance computing [21,22] and LSP algorithm development [23], but more can be done to boost abilities to interpret phenology across broad scales in a standard way. Further, while clustering approaches elegantly summarize phenological similarity and difference, their interpretation in intuitive phenological terms (e.g., growing season or seasonality) can be challenging when the number of starting variables or metrics is large. An optimal system should be complex enough to provide rich phenological characterization, yet simple enough to allow rapid processing and ready interpretation.

We introduce an LSP method that reduces data volume through a series of annual metric computations to isolate the fundamentally important dimensions of phenological variability. This approach provides intuitive metrics and phenoclass memberships with clear quantitative phenological descriptions for each pixel and year. This quantifies variability over space and time along multiple phenological dimensions amenable to comparison with various ecological, biophysical, and climatological data sets. Phenology data products resulting from our analysis can be explored through our online landscape phenology monitoring system (https://landat.org).

Our process is three-fold and begins by generating phenology metrics (e.g., length of season, average greenness) for each pixel using polar coordinate transformation (PCT). The phenology metrics can be tailored to summarize aspects of biophysical change that are of primary interest. Circular plots of time series data resulting from PCT visualize annual phenology as an explicitly cyclical phenomenon [24,25]. This allows for the efficient calculation of the normal phenological year timing (as distinct from the calendar year), extraction of phenology metrics relative to the phenological year, and enhanced the appearance of anomalous phenology years because departures appear as deviations from the normal annual ellipse.

Next, the metrics are transformed through factor analysis to reduce the data to a set of latent variables (factor scores), reducing the variables needed to represent the annual phenology profile while preserving nearly all of the variability. The factor analysis reveals fundamental dimensions of phenology which are intuitive because they correspond to common interpretations of annual phenology, and because their variability has straightforward ecological relevance. Last, the factor scores of each pixel in each year are clustered to yield discrete phenoclasses. Each phenoclass is defined by its mean factor score values, thereby quantifying similarities and differences among phenoclasses. The phenology metrics, factor scores, and phenoclasses can all be mapped to reveal the geography of the various dimensions of vegetation phenological behavior. Tracking year-to-year changes either in factor scores or phenoclass membership provides a quantitative means of assessing landscape dynamics.

An important advantage of polar coordinate transformation is that it facilitates the comparison of locations across large domains. Through PCT, each pixel is standardized to its own phenology calendar, regardless of when its phenological minimum occurs during the year. As a result PCT standardizes each pixel for its own latitudinal, elevational, or climatic characteristics. PCT provides a new capacity to measure phenology relative to the calendar year or the phenology year of each pixel beginning at its own minimum [25]. This has substantial advantages for phenological analysis at continental scales, especially with respect to climate change. This work focuses on developing a system for monitoring landscapes as they respond to ecological and climatic change.

2. Materials and Methods

2.1. Study Area and Data

The study area data are derived from normalized difference vegetation index time-series data set from the MODIS platform. This is a gap-filled and smoothed product with a 250-m nominal resolution, between 20° and 50° latitude [26]. The full data set consists of 188.6 million pixels spanning the Conterminous US and significant parts of Canada and Mexico, with 46 regularly spaced values per year (an 8-day interval) from 2000 through 2018. Processing procedures for this data set differ from the MYD13Q1 and MOD13Q1 16-day MODIS NDVI products. First, NDVI measurements from Aqua and Terra satellites taken on the same day were merged to enhance representation according to the methodology of [26]. Second, artifacts—primarily clouds—were filtered using maximum value compositing, as described by [27,28]. Last, additional temporal processing was performed to correct for other artifacts including snow cover [29].

While MODIS MCD12Q2 provides metrics that identify phenophase transitions, we use PCT to calculate phenological metrics directly from NDVI time series. An important goal of this approach was to rely on as few assumptions as possible, particularly about trends and trajectories of the phenology cycle. This precluded standard LSP products that identify phenophase transitions by first fitting models to the data and then looking for departures. In addition, our approach allows for the period of least phenological activity to be excluded from the analysis based on user-defined thresholds for the growing season. This aids a comparison across latitudes, as high latitudes tend to have snow cover artifacts during the winter that artificially reduce NDVI. PCT also provides an elegant solution for representing temporal patterns that span across the beginning or end of calendar years, as explained below.

2.2. Polar Coordinate Transformation and Phenology Metrics

Polar coordinates are used in atmospheric sciences, engineering, astronomy, and elsewhere to describe two-dimensional measurements relative to a central reference point (e.g., wind speed and direction). Less commonly, polar transformations have been applied to time series data, treating regular cycles as passes around the polar circle. Surprisingly, PCT has rarely been applied to temporal environmental data (but see [24,25,30,31]).

We applied PCT to land surface phenology as the first step in data analysis. We then calculated a set of milestones and descriptive statistics using transformed data to collectively characterize the phenological year and seasonal differences in terms that were relevant to our interest in landscape analysis (

Table 1. Description of phenological variables based on polar coordinate transformation (PCT).

Type	Variable Name	Descriptive Name	Units	Polar Description
Timing Variables	GSbegin	Beginning of growing season	Days	Number of days (or radial angle) corresponding to 15% of cumulative annual NDVI
	GSmid_early	Middle of early growing season	Days	Number of days (or radial angle) corresponding to 32.5% of cumulative annual NDVI
	GSmid	Middle of entire growing season	Days	Number of days (or radial angle) corresponding to 50% of cumulative annual NDVI
	GSmid_late	Middle of late growing season	Days	Number of days (or radial angle) corresponding to 65% of cumulative annual NDVI
	GSend	End of growing season	Days	Number of days (or radial angle) corresponding to 80% of cumulative annual NDVI
Greenness & Seasonality Variables	LOS	Length of growing season	Days	Number of days between early and late growing season thresholds
	mean_NDVI_grw	Average growing season greenness	NDVI	Average NDVI during the growing season (GSbegin to GSend)
	std_NDVI_grw	Variability in growing season greenness	NDVI	Standard deviation of NDVI during the growing season
	AVearly	Magnitude of early growing season seasonality	NDVI	Length of the average vector during early growing season (GSbegin to GSmid)
	AVgrw	Magnitude of entire growing season seasonality	NDVI	Length of the average vector during entire growing season (GSbegin to GSend)
	AVlate	Magnitude of late growing season seasonality	NDVI	Length of the average vector during late growing season (GSmid to GSend)
	Theta (Offset) 1	Offset between calendar year and start of phenological year	Days	Number of days between the beginning of the calendar year (1 January) and the start of the phenological year (defined by when the average minimum in NDVI occurs)

¹ Offset was not a variable used in analysis but was a timing point used to define the start of the phenological year within which all PCT variables were measured.

). The procedure used for PCT is fundamentally the same as commonly found in atmospheric science, e.g., to determine vector averaged wind speed and direction, except here, time is converted into radial coordinates. In the first step the day of the year, $d=1\dots 365$, was converted to radians $r$ using:

$$r=\left(\frac{d}{365}\right)\times 2\pi$$

(1)

Using $r$ and measured NDVI, $v$, the time series can be graphed using polar coordinate pairs $\left[r,v\right]$. Irregularities in the overlapping orbits around the polar plot reflect variations in the eight-day measurements from year to year (Figure 1b,c).

Summarizing polar coordinate data requires projection onto two-dimensional Cartesian coordinates. Each $\left[r,v\right]$ pair was projected on a coordinate plane by applying cosine and sine functions to $\left[r,v\right]$ as:

$$x\left(r,v\right)=v\times \cos \left(r\right)$$

(2)

$$y\left(r,v\right)=v\times \sin \left(r\right)$$

(3)

Several polar measures were calculated using the mean of $x$ and $y$, termed $\overline{x}$ and $\overline{y}$. These were calculated over $n$ samples within the period of interest, simply as:

$$\overline{x}={{\displaystyle \sum }}_{i=1}^n{x}_i/n$$

(4)

$$\overline{y}={{\displaystyle \sum }}_{i=1}^n{y}_i/n$$

(5)

Thus, $\left[\overline{x},\overline{y}\right]$ describes the coordinates of the average vector, which can be projected back into polar coordinates of angular displacement and distance, $\left[\overline{r},\overline{v}\right]$, using the relationships,

$$\overline{r}=\{\begin{array}{c}atan2\left(\overline{y},\overline{x}\right) atan2\left(\overline{y},\overline{x}\right)>0 \\ atan2\left(\overline{y},\overline{x}\right)+2\pi otherwise \end{array}$$

(6)

$$\overline{v}=\sqrt{{\overline{x}}^2+{\overline{y}}^2}$$

(7)

$\left[\overline{r},\overline{v}\right]$ is the resultant vector with a magnitude proportional to the normal (multi-annual) strength of seasonality, and direction indicating the central tendency of NDVI in the annual cycle. Phenological timing varies significantly across large landscapes due to latitudinal, elevational, hydrological and other environmental gradients affecting vegetation. For example, the phenology of deciduous locations across much of the southwest is six months out of phase with much of the Conterminous United States (Figure 2a). Polar coordinate transformation facilitates the calculation of standard metrics, normalized to the customized phenological year of any given pixel, which can clarify similarities in modes of phenological variability while preserving timing differences.

Following the polar transformations shown above, the next step identifies the unique starting date for the phenological year associated with each pixel. This corresponds to the period of least photosynthetic activity, when NDVI values are generally lowest (i.e., winter for most pixels), which rarely coincides exactly with the beginning of the calendar year. The period of least activity is simply the diametric opposite of the average vector $\left[\overline{r},\overline{v}\right]$, i.e., rotated exactly one-half year ($\pi$ radians) from $\overline{r}$. This new angle $\theta$ is a complement of $\overline{r}$ based on Equation (6) using the full multi-annual time series of NDVI values $\left(n=19\times 46=874\right)$. That is,

$$\theta =\{\begin{array}{c}\left(\overline{r}+\pi \right) \left(\overline{r}<\pi \right) \\ \left(\overline{r}-\pi \right) otherwise \end{array}$$

(8)

We recorded this parameter as the offset (in radians or back-transformed to days) between the beginning of the phenological year and that of the calendar year. The straight line formed by $\theta$ and $\overline{r}$ bisects the polar graph into two sections of equal area under the multi-annual NDVI profile; the sum of NDVI values clockwise from $\theta$ to $\overline{r}$ is equivalent to the sum clockwise from $\overline{r}$ to $\theta$.

Subsequent phenological metrics (Table 1) were calculated and are reported by their calendar date. The time series of coordinate pairs $\left[r,v\right]$ was first divided into phenological years, beginning with the first NDVI reading after $\theta$. The initial and final segments of values that occur before $\theta$ in the first year and after $\theta$ in the final year were omitted, thus resulting in $19–1$ phenological years with 46 dates each $\left(r\_1,1,\dots ,r\_18,46\right)$. Within each phenological year, NDVI values were converted to cumulative proportions of the year’s total NDVI accumulation, from zero to one. Proportional values were used to determine benchmarks that defined beginning, middle, and ending growing season dates. Any NDVI threshold, $t$ with range $\left(0,1\right)$, can be used to mark a phenological milestone, $J\left(t\right)$, corresponding to the earliest date $J$ when the cumulative proportion exceeds the threshold value $t$. We chose the phenological milestones $J\left(0.15\right)$, $J\left(0.5\right)$, and $J\left(0.8\right)$ to define the beginning (GSbegin), middle (GSmid), and end (GSend) of the growing season (Table 1). These correspond to 15%, 50%, and 80% completion thresholds of the phenological year; other threshold values can be chosen to suit a particular analysis. The length of the growing season (LOS) was simply the number of days (or radians) between $J\left(0.15\right)$ and $J\left(0.8\right)$ (Figure 1c).

Annual measures of NDVI’s central tendency include the mean (mean_NDVI_grw) and standard deviation (std_NDVI_grw) of observed NDVI within the growing season, indicating the average productivity and within-season variability, respectively. We also calculated an average vector for the growing season, AVgrw, $\left({\overline{r}}_{grw},{\overline{v}}_{grw}\right)$ using Equations (4)–(7) with appropriate n, whose vector length is also a measure of the degree of seasonality. For additional information about the seasonal pattern, we examined the early (GSbegin to GSmid) and late (GSmid to GSend) portions of the growing season separately. We calculated average vectors for both of these periods, indicated by $\left({\overline{r}}_{early},{\overline{v}}_{early}\right)$ and $\left({\overline{r}}_{late},{\overline{v}}_{late}\right)$, respectively. This yielded two additional vector directions (GSmid_early and GSmid_late) and two additional vector lengths (AVearly and AVlate). Collectively, these vectors and the differences among them provide insights about modality and symmetry throughout the growing season (Figure 1d). All phenological metrics are sensitive to different aspects of growing season variability, as discussed below. PCT analyses were performed in R [32]. Our repository of R code is accessible on GitHub at https://github.com/LandscapeDynamics/PolarMetrics; additional details of the PCT process are in [25].

2.3. Dimensional Reduction and Phenological Classification

Through the process described above, we produced 11 metrics to describe LSP (Table 1). These included six measures describing aspects of seasonality and greenness/productivity, and five timing milestones, i.e., growing season benchmarks recorded as a day of year ($r$, in radians). For further analysis as quantitative metrics, these circular timing milestones were transformed using the $sin\left(r\right)$ and $cos\left(r\right)$ functions such that dates near the beginning/end of the year had similar values. For example, 1 January and 31 December are distant when expressed in radians or days-of-year, but appropriately near in the circular sine-cosine coordinate space. This resulted in five timing variable pairs, e.g., $\left[sin\left(GSbegin\right), cos\left(GSbegin\right)\right]$. Therefore, the data set that we subjected to factor analysis consisted of 16 variables for each pixel-year combination.

We used factor analysis to reduce the 16 variables to a smaller set of linear combinations (factors) that explained the bulk of the variation across the data, and to help identify core dimensions of phenological variability. We performed a factor analysis using the factanal base stats function in R [32]. Pearson correlation coefficients among the 16 variables were used to calculate factor loadings. The factors having eigenvalues greater than one were retained and rotated using varimax rotation to improve interpretability. The rotated factor pattern (loadings) was then used to generate factor scores for each pixel-year combination.

Factors are continuous variables, and the combination of factor scores assigned to a particular pixel is potentially unique. However, the ecological significance of phenological differences measured by differences in one or more factor-scores is not known a priori. Pixels can be grouped by similarity to facilitate comparisons but implicit in such grouping is the presumption that inter-group differences are more meaningful than intra-group differences. Cluster analysis provided a statistically rigorous means of grouping pixel-years based on their factor scores, identified herein as phenological classes or phenoclasses. The term phenoregion is also used [18], but phenoclass is preferred when pixels can change group membership over time in response to environmental change and when pixels in each group may be widely disjunct geographically. Phenoclasses can be thought of as categorical types into which the multidimensional phenological space has been partitioned, and within which vegetation exhibits a common and characteristic phenological pattern.

We used non-hierarchical k-means clustering [33] where input variables were the factors at the pixel-year level. Objective clustering was used to minimize the error between each of the $k$ representative centroids and the cluster members (pixel-years) using the sum of squares criterion while maximizing the differences between centroids [34]. A k-means clustering approach was used because it is unsupervised and provides several options for seed selection [35,36,37,38,39].

Ensuring that the clustering process begins with high-quality initial seeds typically results in better-separated and more representative final clusters. For our purposes, the seed finding algorithm used in SAS’s “PROC FASTCLUS” [40] was suitable. Several of the more widely used seed finding algorithms [35,36,37,38] were tested and found to produce seeds whose distributions broadly mirrored those of the factor score values. That is, these methods placed the most seeds where most of the factor score values were located, with few seeds outside of that central domain. This, in turn, resulted in final cluster centroids with a similar relation to the original data distribution. Although desirable in many contexts, one consequence when measuring phenological change is that phenoclasses outside of the denser part of the data distribution exhibit both larger intra-cluster and inter-cluster differences, with more widely spaced centroids. For this reason, transitions (changes between phenoclasses) in the sparse data region would only be observed when phenological changes were unusually large, whereas in the denser part of the data distribution only comparatively small changes would be required to elicit phenoclass transitions. The result in ecological terms is that a typical landscape undergoing minor transitions between closely spaced phenoclasses can appear just as dynamic as a strongly disturbed landscape undergoing transitions between widely spaced, uncommon phenoclasses. Thus, a clustering process that tends to result in more uniformly distributed final cluster centroids, as employed here, has advantages in the context of monitoring phenological change.

The number of classes, $k,$ is negatively correlated with how large a step in the factor score space a phenoclass transition represents. Typically, $k$ is chosen a priori based on knowledge of or subjective preference for a given level of division, or it is iteratively determined based on a trade-off between within-cluster and between-cluster deviance. Deviance between groups (i.e., the between sum-of-squares, $BSS$) should be large relative to deviance within groups (i.e., the within sum-of-squares, $WSS$). This can be measured as $BSS/ \left(BSS+WSS\right)$ or $BSS / TSS$ or (total sum of squares), which approaches 1 as $k$ approaches the number of pixel-years. In the end we chose a 500-class analysis for demonstrative purposes based on (1) optimizing computation, and (2) ensuring that there was sufficient resolution to show relevant ecological divisions, without showing phenology changes due to minor climatic variation, which detracts from our focus on detecting dynamic change.

3. Results and Discussion

3.1. Phenology Metrics

Each PCT variable measures an intuitive aspect of phenology, e.g., the start of the growing season, the length of the season, and the average growing season greenness. Reducing the many original measurements in each year’s NDVI profile to a suite of PCT variables (Table 1) provides both a substantial reduction in data volume and intuitive ecological characterization. Phenology maps using PCT illustrate differences due to varying environmental drivers, including topography, climate, and land use (Figure 2).

Five of the 11 PCT variables in Table 1 are timing variables, which describe the juxtaposition of the growing season profile with the annual solar cycle. These variables quantify complex patterns in one of the most important components of phenological variability, particularly in the western part of the continent (Figure 2a). For example, the middle of the growing season (GSmid) is coincident with late summer across most of the continent but shows very different timing in the southwestern US, where vegetational phenology is more strongly asynchronous with the solar cycle. There, GSmid varies broadly from September–October in the east (e.g., the Chihuahuan Desert) to February–March in the west (e.g., the Mohave Desert). Along this broadly east–west timing gradient, there is an irregular band through southern Nevada, Arizona, and northwestern Mexico with growing season midpoints between December and January (near Julian days 365 and 1). This regional gradient illustrates the utility of transforming day-of-year phenology metrics to sine-cosine pairs whose values grade appropriately across the end-of-year transition. Moreover, the structure of spatial and temporal gradients in phenology timing is preserved in statistical products such as the factor scores, which take the sine-cosine variables as inputs.

Other PCT variables indicate strong contrast in seasonality and growing season length between much of the northern US (midwestern US through northeastern US and Mississippi River Valley) and the rest of the continent (Figure 2b,d). Landscapes across this broad northern belt are dominated by agricultural crops, prairies, and deciduous trees. The growing season length variable (LOS, Figure 2b) exhibits low values in this region, likely driven by the timing of agricultural activities as well as natural deciduousness. In comparison, more evergreen landscapes in the Pacific Northwest, Eastern Canada, and the southern and southwestern US show larger LOS values. Low seasonality in the latter regions is also indicated by seasonal magnitude (the average growing season vector length, AV_grw, Figure 2d). Although considerable landscape variability is present within regions, the agricultural landscapes of the midwestern US appear to be among the most seasonal systems on the continent, followed by regions dominated by deciduous vegetation such as the central Appalachians and the northeastern US.

Variables quantifying NDVI greenness were correlated with growing season productivity. For example, high growing season NDVI greenness (mean_NDVI_grw) values for the dense forests of the eastern US and Pacific Northwest contrasted with much of the western US (Figure 2c). NDVI greenness also aided in further discriminating regional and landscape phenologies. Consider, for example, that on the basis of length of season (LOS) alone, the Pacific Northwest and Mountain West were difficult to distinguish (Figure 2b). However, through mean_NDVI_grw (Figure 2c), the ecological distinction was clear. The high NDVI of productive, evergreen vegetation in the Pacific Northwest contrasted sharply with lower values in the intermountain and mountain west regions where sparser, more dry-adapted evergreen vegetation dominates.

3.2. Dimensional Reduction

Factor scores provide a small, orthogonal set of variables that describe phenological pattern, and are a concise means of quantifying spatial (among-pixels) and interannual (within-pixel) differences in phenology. While PCT reduced the data volume by two thirds (46 variables per year to 16), a further reduction in data volume and increased explanatory power per variable was achieved through factor analysis. There were four factors with eigenvalues greater than 1, which collectively explained over 95 percent of the variance among the 16 variables in the PCT data set (

Table 2. Factor loadings from factor analysis using varimax rotation. Coefficients smaller than $\left|0.2\right|$ are not shown. Variable definitions are given in Table 1.

		Factor 1	Factor 2	Factor 3	Factor 4
Timing Variables	GSbegin sin	−0.892	0.311	0.244
	GSbegin cos	0.283	0.911
	GSmid_early sin		0.959
	GSmid_eary cos	0.927		−0.245	0.202
	GSmid sin	0.675	0.702
	GSmid cos	0.672	−0.662		0.241
	GSmid_late sin	0.936		−0.231	0.215
	GSmid_late cos		−0.939		0.304
	GSend sin	0.568	−0.689		0.392
	GSend cos	−0.764	−0.579
Greenness & Seasonality Variables	LOS				0.898
	mean_NDVI_grw	−0.209		0.964
	std_NDVI_grw			0.321	−0.836
	AVearly			0.960
	AVgrw	−0.247		0.838	−0.457
	AVlate	−0.274		0.859	−0.367
		Factor 1	Factor 2	Factor 3	Factor 4
	Proportional Variance	0.294	0.282	0.231	0.145
	Cumulative Variance	0.294	0.576	0.807	0.953

).

Rotated factor loadings discriminate three main phenology characteristics: timing, productivity, and seasonality. Rotated factors 1 and 2 are heavily loaded on the timing variables (GSbegin_cos, GSbegin_sin, etc.), and together explain 58 percent of the variance in the data. This indicates that most of the observed differences in phenology are simply determined by when milestones occur. The complementary loadings of factors 1 and 2 (Table 2) indicate that two dimensions are required to represent circular annual timing as expressed through the sine and cosine variables (Figure 3). Factor 3 explains an additional 23% of the variance and was heavily loaded on variables associated with vegetative productivity as reflected in higher NDVI values, (i.e., mean_NDVI_grw) and the vector length variables (AVgrw, AVearly, AVlate). This is not surprising, as NDVI is one of two measures used to calculate the vector lengths. Factor 4 had heavy loadings for length of season (LOS) and variability in NDVI (std_NDVI_grw) and is referred to here as “seasonality.” Note that due to the signs of the loadings, factor 4 values related inversely to seasonality. Thus, high LOS values and low variation in NDVI result in high factor 4 values, whereas short growing seasons with high intra-season variability result in low factor 4 values. Our interpretation of factor 3 and 4 loadings is that while several of the original variables reflect mixtures of seasonality and productivity, these components are separable and result in two orthogonal factors.

The false-color composite continental map of phenological variability shown in Figure 4 was generated for an individual phenological year by assigning one of the timing factors (factor 1) in combination with the productivity and seasonality factors to blue, green, and red, respectively. This visualization illustrates that the combined factors provide a nuanced description of phenological variation, and collectively are a robust indicator of gradients in ecological and biophysical diversity at landscape, regional, and continental scales. For example, agriculture-dominated landscapes across the US Midwest corn belt and the Mississippi River Valley, where the timing of planting and harvest are synchronized and highly regular at regional scales, show internal similarity that distinguishes them from other regions.

The factor scores also quantify phenological change within pixels. We used a false-color composite mapping approach to illustrate magnitudes of inter-annual variability in the same three factors (Figure 5). For example, grassland and shrubland regions from the Dakotas to south Texas show high inter-annual variability across all factors, whereas forest and desert regions are more stable. Most regions show distinctive mixtures of variability in some phenological characteristics and stability in others.

3.3. Phenological Classification

Transitional changes of a pixel from one phenoclass to another between years is determined by the relative position of cluster centroids in the four-dimensional factor score space. Each centroid represents the “center of mass” of a cluster group (i.e., phenoclass) in that space. The phenoclass centroids exhibit different distributions along each of the four factor dimensions, and in three of those dimensions the values are concentrated about a central tendency (Figure 3). As discussed above in Section 2.3, we used a clustering approach with a seed finding algorithm that produced more uniformly spaced centroids, so that the centroid distributions are more uniform than the underlying factor score data. This was done to reduce bias in phenoclass change detection that would otherwise occur when counting transitions for phenoclasses that occurred near the central distribution of the data. This difference is especially evident in the factor 1 and 2 dimensions, wherein the circular joint distribution of timing factors is well-represented by fairly evenly distributed phenoclasses (Figure 3a). Utilizing a more uniform cluster centroid distribution allows for a more ecologically consistent measure of change from one phenoclass to another. Ultimately, classifying continuous factor scores into discrete phenoclasses provides a basis for measuring transitional jumps between phenological states, which can help in identifying meaningful ecological change.

Note that a phenoclass map will illustrate spatial and temporal patterns effectively identical to those evident in Figure 4 and Figure 5, which explicitly map three factor scores as red-green-blue (RGB) composites. This is because phenoclasses preserve factor score values in discretized form, with each phenoclass defined by its centroid’s four factor values. The greater the number of phenoclasses produced by a given analysis, the more closely a phenoclass map will approximate continuous, rather than discrete, factor score variability.

Our results demonstrate that temporal variation in phenology is geographically variable (Figure 5a). Moreover, phenoclass change through time can be used to quantify directional phenological change as reflected in the centroid values. For example, changes in the frequency of different phenoclasses at the continental scale suggest directional shifts in timing, productivity, and seasonality over the period of observation (Figure 5b). We can precisely quantify such changes, but we cannot say with certainty that the observed changes are ecologically important. It is possible to examine areas believed to have experienced substantive change, however, and to objectively assess whether land surface phenology provides a useful method for detection.

To illustrate, we examined factor scores and phenoclasses for an area of the south-central US, centered on Texas, known to have been affected by severe drought during October 2010 through September 2011 (Figure 6). Impacts of the drought on vegetation are captured by the phenological response across the immediate drought interval, in terms of both regional factor score variation and frequency of phenoclass changes. A strong difference between the 2011 growing season and the growing seasons of adjoining years 2010 and 2012 is obvious in the multitemporal factor score maps (Figure 6a,b). The extent and severity of short-term drought impacts, including reduced productivity and delayed growing season timing, are evident throughout eastern and central Texas and north across the High Plains. Anomalous factor values in 2011 also drove a more than 8% increase in pixel-level phenoclass transitions (Figure 6d).

Phenoclass transitions provide a wide variety of additional information about landscape condition and change not explored here. Information pertaining to dynamic landscape mosaics, directional ecological change, predictability, and stability over short and long terms can be gathered, for example, by examining multi-pixel landscapes in terms of their factor score and phenoclass composition and dynamics. Because pixels can change phenoclass assignment between years, it is possible to monitor phenological change in a framework of states and transitions. Preliminary analyses using this approach to address phenological evolution, disturbance, recovery, and land management applications look promising.

Note that phenoclasses are not explicitly land cover classes, although there may be congruence between many land-cover types and their signature phenology. Some land-cover types likely exhibit a mixture of phenological responses due to interannual meteorological differences or other factors. Alternatively, some different vegetation classes likely share phenological signatures, e.g., variation in species composition among pine-dominated ecosystems. Examining such relationships could readily be performed using the data products described here.

4. Conclusions

The three-step process of polar coordinate transformation, dimensional reduction using factor analysis, and statistical clustering into phenoclasses proved to be an effective and efficient means of using land surface phenology to characterize and monitor landscape dynamics. PCT allowed a ten-fold reduction in data volume while showing no loss in sensitivity to NDVI change. Multiple PCT-generated phenology metrics provide a rich characterization of the LSP cycle across large areas at high spatial and temporal resolution. Any of these can be used in focused analyses, based on their relevance for particular types of ecological investigations (e.g., climate drift, land degradation, ecological disturbance and recovery).

Developing factor score variables and phenoclasses from PCT variables resulted in a concise, integrative mapping and classification of LSP on an annual basis. Each phenoclass represents the gross phenological timing, productivity, and seasonality of the vegetation within a pixel, and variations in phenology are quantified across years. Movement along the factor dimensions and transitions between phenological states provide quantifiable evidence of ecological change. An advantage for landscape ecologists and others studying phenological variation and change is that landscapes are depicted using intuitive measures (PCT variables) rather than a complex time series of NDVI values.

Vegetation condition and dynamics mediate, in varying degrees, linkages between natural resources of conservation interest in terrestrial systems (e.g., endangered species, ecosystem services, forest products) and their stressors and drivers (e.g., climate change, land use change, wildland fire). Land surface phenology facilitates vegetation monitoring via remote sensing beyond thematically and temporally coarse land use/cover and towards ecologically nuanced gradients and dynamics. As such, the potential applications of large-scale LSP data sets for studying biodiversity and natural resource conservation and management are vast. These include studies of ecosystem resilience and vulnerability, species and resource distribution modeling, conservation planning, and other applications. The data presented here are available for such applications through an actively managed online database and viewer, namely the Landscape Dynamics Assessment Tool (LanDAT, https://landat.org).