Using Hidden Markov Models for Land Surface Phenology : An Evaluation Across a Range of Land Cover Types in Southeast Spain

Land Surface Phenology (LSP) metrics are increasingly being used as indicators of climate change impacts in ecosystems. For this purpose, it is necessary to use methods that can be applied to large areas with different types of vegetation, including vulnerable semiarid ecosystems that exhibit high spatial variability and low signal-to-noise ratio in seasonality. In this work, we evaluated the use of hidden Markov models (HMM) to extract phenological parameters from Moderate Resolution Imaging Spectroradiometer (MODIS) derived Normalized Difference Vegetation Index (NDVI). We analyzed NDVI time-series data for the period 2000–2018 across a range of land cover types in Southeast Spain, including rice croplands, shrublands, mixed pine forests, and semiarid steppes. Start of Season (SOS) and End of Season (EOS) metrics derived from HMM were compared with those obtained using well-established smoothing methods. When a clear and consistent seasonal variation was present, as was the case in the rice croplands, and when adjusting average curves, the smoothing methods performed as well as expected, with HMM providing consistent results. When spatial variability was high and seasonality was less clearly defined, as in the semiarid shrublands and steppe, the performance of the smoothing methods degraded. In these cases, the results from HMM were also less consistent, yet they were able to provide pixel-wise estimations of the metrics even when comparison methods did not.


Introduction
Land Surface Phenology (LSP), the seasonal pattern of variation in vegetated land surfaces observed from remote sensing [1,2], has been proven in the last decades to be a basic tool to investigate climate change impacts in ecosystems (e.g., References [3][4][5]).The availability of long time series of vegetation indices derived from different space-borne sensors has allowed for the description of phenology metrics at global, regional, or local scales, and to analyze their spatio-temporal changes in relation to climate trends and human or natural disturbances (e.g., References [4][5][6][7][8][9][10][11][12][13][14]).
Moderate Resolution Imaging Spectroradiometer (MODIS) derived vegetation index (VI) time series have provided a convenient balance between the length of the available time series and temporal and spatial resolution, being available at 250 m pixel size, every 16 days since the year 2000 [7,15,16], which is suitable for describing and analyzing seasonal patterns and trends in regional-to-global LSP metrics.
Normalized Difference Vegetation Index (NDVI) is the most commonly used vegetation index to extract phenology metrics, although alternative vegetation indices can also be used and their possible advantages and drawbacks to derive phenology metrics in specific land cover types have been discussed (e.g., References [17][18][19]).Different methods have been proposed to extract phenological parameters from NDVI or other VI time series (e.g., References [6,7,[20][21][22][23][24]; see a review in Reference [25]), with some form of smoothing method, polynomial or model function fitting, and threshold or inflexion point estimation being the most commonly employed strategy, with active discussion on the performance of the different methodologies [26,27].
HMM were proposed more than two decades ago by Viovy and Saint [54] as a useful method to obtain LSP metrics.However, alternative methods of analysis were preferred at the time, and only recently, a very limited use of HMM for applications in LSP can be found.Cerqueira Leite et al. [47] and Siachalou et al. [49] used HMM modeling of crop phenology with the aim of classifying different types of crops from time series of remote sensed images.Shen et al. [55] used HMM modeling of corn phenological phases to estimate corn progress from NDVI and meteorological data.García et al. [56] discussed the application of HMM to extract LSP metrics from the time series of NDVI data.Given the successful experiences with applications of HMM based methods in different fields in the last decades, as well as the expansion in computing capabilities that allows for efficient implementation in any modern desktop computer, it is timely to reconsider the use of HMM for LSP and to test their performance in comparison with usual methodologies.There are two main expected potential advantages in the use of HMM to derive LSP metrics.Firstly, as discussed in Reference [56], HMM can model the different phenological states of the vegetation and the transitions between them, and so they could be tailored to include a priori information about the phenology of specific types of vegetation or land cover types.Secondly, should the estimation of LSP metrics using basic types of HMM provide satisfactory results, there would be wide opportunities for further improvements, as there are many variants and extensions of HMM that have already been employed in related signal analysis problems (see, e.g., References [29,42]).
The aim of this paper is to evaluate the potential of using HMM to derive LSP metrics from time series of vegetation indices.To this end, we analyzed NDVI time series data for the period 2000-2018 across a range of land cover types in Southeast Spain, including rice croplands, shrublands, mixed pine forests, and semiarid steppes, obtaining SOS and EOS metrics derived from HMM that were compared with those obtained using well-established smoothing methods.Our main objective was to check the feasibility of applying HMM to large areas with different types of vegetation, including semiarid ecosystems that exhibit high spatial variability and low signal-to-noise ratio in seasonality, where smoothing methods could possibly fall short of exhibiting the best performance.
From the results obtained in this work, we concluded that phenology metrics derived from the simple type of HMM, considered in this paper, were consistent with those obtained from commonly employed smoothing methods, that they were able to provide metrics estimations for all pixels even in some cases where smoothing methods partly failed, and that they showed correlations with climate variables in a similar way to the best smoothing methods, suggesting the validity of the HMM derived metrics.

Hidden Markov Models Basis
For the sake of clarity, we recall in this section the basic elements of HMM, in the context of their use in this paper.More general and detailed descriptions can be found in References [28,54].For technical aspects of inference and algorithms, we refer the reader to the specialized literature (see, e.g., Reference [57]).
An HMM consists of a sequence of observable states, {Y t }, a sequence of unobservable or hidden states, {S t }, and two probabilistic models determining how both sequences are generated (Figure 1a).The observable states, called "emissions" from the origins of HMM in the field of speech recognition, correspond in this paper to a time series of NDVI increments.The probability distributions of these observed states depend on the hidden states (Figure 1a), which could refer to the corresponding "phenological states" of the vegetation, i.e., growing, declining, or stationary.We assume that, for a given hidden state, the observed NDVI increments follow a Gaussian distribution whose mean and standard deviation depend solely on the hidden state.Thus, for instance, the expected mean values of NDVI changes would be positive or negative for the growing or declining hidden states, respectively.The probabilistic model for the sequence of hidden states is a Markov chain, defined by the transition probabilities between states.The basic type of HMM used in this paper is a first order HMM (Figure 1a), where the transition probabilities from a hidden state to any other depend only on the current state, and not on the history of the system.Thus, in Figure 1a, the absence of an arrow connecting two states, hidden or observed, implies their conditional independence, so that, given S t , both S t+1 and Y t are independent of previous hidden or observed states.
The so-called topology of the HMM is defined when the number of hidden states is given, as well as the possible existence of forbidden transitions between states (called blocking).Figure 1b presents the topology of the HMM used in this work, an HMM with four states, increasing (S + ), stationary high (S 0+ ), decreasing (S − ), and stationary low (S 0− ), with transitions between states only allowed for those connected with arrows.Since the sum of the transition probabilities from a given hidden state to any other is one, there are only four independent transition probabilities.In total, there are twelve independent parameters in the model, as for each of the four hidden states the Gaussian density for the emissions is defined by its mean and standard deviation.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 17 climate variables in a similar way to the best smoothing methods, suggesting the validity of the HMM derived metrics.

Hidden Markov Models Basis
For the sake of clarity, we recall in this section the basic elements of HMM, in the context of their use in this paper.More general and detailed descriptions can be found in References [28,54].
For technical aspects of inference and algorithms, we refer the reader to the specialized literature (see, e.g., Reference [57]).
An HMM consists of a sequence of observable states, {Yt}, a sequence of unobservable or hidden states, {St}, and two probabilistic models determining how both sequences are generated (Figure 1a).
The observable states, called "emissions" from the origins of HMM in the field of speech recognition, correspond in this paper to a time series of NDVI increments.The probability distributions of these observed states depend on the hidden states (Figure 1a), which could refer to the corresponding "phenological states" of the vegetation, i.e., growing, declining, or stationary.We assume that, for a given hidden state, the observed NDVI increments follow a Gaussian distribution whose mean and standard deviation depend solely on the hidden state.Thus, for instance, the expected mean values of NDVI changes would be positive or negative for the growing or declining hidden states, respectively.The probabilistic model for the sequence of hidden states is a Markov chain, defined by the transition probabilities between states.The basic type of HMM used in this paper is a first order HMM (Figure 1a), where the transition probabilities from a hidden state to any other depend only on the current state, and not on the history of the system.Thus, in Figure 1a, the absence of an arrow connecting two states, hidden or observed, implies their conditional independence, so that, given St, both St+1 and Yt are independent of previous hidden or observed states.
The so-called topology of the HMM is defined when the number of hidden states is given, as well as the possible existence of forbidden transitions between states (called blocking).Figure 1b presents the topology of the HMM used in this work, an HMM with four states, increasing (S+), stationary high (S0+), decreasing (S-), and stationary low (S0-), with transitions between states only allowed for those connected with arrows.Since the sum of the transition probabilities from a given hidden state to any other is one, there are only four independent transition probabilities.In total, there are twelve independent parameters in the model, as for each of the four hidden states the Gaussian density for the emissions is defined by its mean and standard deviation.Once the topology of the HMM is defined, as well as the type of probability distributions for the emissions, given a sequence of observed states the parameters of the model can be fitted using an expectation maximization algorithm [58,59] that provides the values of the parameters that maximize the likelihood of the observations.Then, having values for the parameters of the model, Once the topology of the HMM is defined, as well as the type of probability distributions for the emissions, given a sequence of observed states the parameters of the model can be fitted using an expectation maximization algorithm [58,59] that provides the values of the parameters that maximize the likelihood of the observations.Then, having values for the parameters of the model, the sequence of hidden states can be inferred using a dynamic programming algorithm known as the Viterbi algorithm [60].In this way, given a time series of NDVI changes, phenology metrics such as the Start of Season (SOS) and End of Season (EOS) can be defined from the inferred sequence of hidden states, as SOS and EOS would correspond to some moments where the hidden state of the system is S + or S − , respectively.

Study Areas
Five study areas in Southeast Spain, located in three different regions, were selected from general on-the-ground information regarding their predominant land cover types (Figure 2).In Valencia region, Albufera site comprises an area of traditional rice crops surrounding a coastal lagoon; Millares and Ayora sites are inland mountain areas, with altitudes ranging from 600 to 800 m in Millares, and 700 to 1000 m in Ayora, with vegetation cover dominated by shrublands.In Murcia region, Espuña site, located in Espuña range, with altitudes ranging from 650 to 1300 m, is a mixed pine forest derived from a reforestation project dating back to the end of the 19th century.In Andalucía region, the Níjar site is located close to the Southern Mediterranean coast, with altitudes ranging from 250 to 900 m with grass steppes as the main land cover.
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 17 the sequence of hidden states can be inferred using a dynamic programming algorithm known as the Viterbi algorithm [60].In this way, given a time series of NDVI changes, phenology metrics such as the Start of Season (SOS) and End of Season (EOS) can be defined from the inferred sequence of hidden states, as SOS and EOS would correspond to some moments where the hidden state of the system is S+ or S-, respectively.

Study Areas
Five study areas in Southeast Spain, located in three different regions, were selected from general on-the-ground information regarding their predominant land cover types (Figure 2).In  In each study area, around 500 MODIS pixels were visually selected using Google Earth TM images, which provided enough resolution to try to avoid areas with heterogeneous land cover (e.g., crops in areas of natural vegetation, very steep slopes in mountain areas, and civil constructions) without the need of any special software for image analysis.Approximate central coordinates and number of selected pixels for each site are shown in Table 1; kmz files showing the grid of pixels for each area are included as supplementary material (File S1).In each study area, around 500 MODIS pixels were visually selected using Google Earth TM images, which provided enough resolution to try to avoid areas with heterogeneous land cover (e.g., crops in areas of natural vegetation, very steep slopes in mountain areas, and civil constructions) without the need of any special software for image analysis.Approximate central coordinates and number of selected pixels for each site are shown in Table 1; kmz files showing the grid of pixels for each area are included as supplementary material (File S1).

Data Adquisition
MODIS data were downloaded from the Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov/)managed by the NASA Earth Science Data and Information System (ESDIS) project.We used a time series of NDVI data, from MODIS product MOD13Q1 Version 6 [16], corresponding to the best pixels in 16 day intervals (NDVI 16-days composite band) at 250 m resolution, from February 2000 to June 2018.This data is stored in scientific data sets (SDS) under the NASA Hierarchical Data Format Earth Observing System (HDF-EOS) format, which is the standard archive format for all EOS Data Information System (EOSDIS) products [16].We used the NASA Application for Extracting and Exploring Analysis Ready Samples (A EEARS) and R software to access and download the NDVI-MODIS data.Pixel-reliability (PR) and Day-of-the-year (DOY) layers were also downloaded and used in some analyses.PR provides a rank of the overall quality of the pixel, and DOY provides the actual date of acquisition of the reflectance's used to compute each pixel NDVI.

Methods of Analysis
General statistical analyses, including descriptive statistics, and correlation and regression analyses, were carried out using the software R [61].

Pixel Homogeneity Check
We checked similarity in phenological behaviors of the selected pixels in each study area, by decomposing the time series of NDVI for each pixel (NDVI curve) as the sum of a trend or secular component (computed as the average in a moving window of size corresponding to one year) and an annual component, and computing correlations between the annual component for each pixel and the average for the whole area (mean annual curve).

Hidden Markov Models Fitting and Analysis
For each NDVI curve, data with NA values or with PR > 1 were excluded.Using DOY data, NDVI values were interpolated every four days, and a mild smoothing with kernels [0.0370, 0.1111, 0.2222, 0.2593, 0.2222, 0.1111, and 0.0370] was applied.From the resulting equally spaced NDVI time series, finite differences were computed.The set of finite differences series for all pixels in each study area was used to fit a four state HMM with blocking (Figure 1b) and Gaussian emissions, using the HMM Toolbox for Matlab TM [62].For each study area, once the fitted parameters were obtained, Viterbi algorithm was applied to obtain the most probable sequence of hidden states for each pixel.Finally, for each pixel and season, SOS (respectively, EOS) was defined to be the day of the year corresponding to a given percentile (25 or 50) in the series of hidden S + states (respectively, S − ).

Phenology Metrics Using Smoothing Methods
Two different tools, based on smoothing methods and commonly used to extract phenology metrics, were employed to compare HMM results: The R package greenbrown [63][64][65] and the stand-alone program TIMESAT [22,66].
Phenology analyses with the R package greenbrown were carried out, excluding NDVI values with PR > 1, in an automated way using the package function "Phenology" with default parameters.As described in Reference [64], this function performs three successive steps: 1) Filling gaps in the time series; 2) smoothing and interpolation; and 3) extraction of phenology metrics SOS and EOS.Three different methods of smoothing and interpolation were applied, using linear interpolation, splines, or double logistic functions; the resulting metrics from each method are denoted in what follows as gb_lin, gb_spl, and gb_dbl, respectively.For the SOS and EOS definitions, we used the approach "white" [67], based on 50% thresholds between minimum and maximum NDVI values.
For TIMESAT analyses, NDVI values were first weighted according to PR, such that values with a reliability of 0 ("good data") were weighted as 1, reliability of 1 ("marginal data") were weighted as 0.5, and all other reliability values were weighted as 0.1.Second, the data were smoothed by removing spikes from the time series.Spikes were identified as observations that deviated from the median of a moving window by more than one standard deviation.Once the time series was cleaned and smoothed, we used a double logistic fitting method to describe the seasonal curve, as described by Zhang et al. [7].SOS and EOS were established as 20% of the annual amplitude.

Homogeneity Analysis and Seasonality Behavior in Different Land Cover Types
Homogeneity analysis showed generally high correlation values between NDVI seasonal curves for each pixel and the mean curve in each study area, as illustrated, for example, in Figure 3 for Níjar site; similar figures for the other study areas are included as supplementary material (Figure S1).As seen in Figure 3, there is certain variability in correlation values, probably related in this site with differences in slope aspect; in other areas (Figure S1) there are also a few pixels having lower correlations with the mean seasonal curve; however, no areas exhibiting markedly different phenology behaviors were detected in any of the sites.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 17 follows as gb_lin, gb_spl, and gb_dbl, respectively.For the SOS and EOS definitions, we used the approach "white" [67], based on 50% thresholds between minimum and maximum NDVI values.
For TIMESAT analyses, NDVI values were first weighted according to PR, such that values with a reliability of 0 ("good data") were weighted as 1, reliability of 1 ("marginal data") were weighted as 0.5, and all other reliability values were weighted as 0.1.Second, the data were smoothed by removing spikes from the time series.Spikes were identified as observations that deviated from the median of a moving window by more than one standard deviation.Once the time series was cleaned and smoothed, we used a double logistic fitting method to describe the seasonal curve, as described by Zhang et al. [7].SOS and EOS were established as 20% of the annual amplitude.

Homogeneity Analysis and Seasonality Behavior in Different Land Cover Types
Homogeneity analysis showed generally high correlation values between NDVI seasonal curves for each pixel and the mean curve in each study area, as illustrated, for example, in Figure 3 for Níjar site; similar figures for the other study areas are included as supplementary material (Figure S1).As seen in Figure 3, there is certain variability in correlation values, probably related in this site with differences in slope aspect; in other areas (Figure S1) there are also a few pixels having lower correlations with the mean seasonal curve; however, no areas exhibiting markedly different phenology behaviors were detected in any of the sites.Mean NDVI curves, showing different seasonal behaviors in croplands and in natural vegetation are presented in Figure 4.While in rice crops (Figure 4a) a clear and consistent intra-year seasonal variation was present, in the areas with natural vegetation (Figure 4b-e) seasonality was essentially defined by the extreme dry conditions in summer, with very low precipitation and high temperatures.
Thus, in natural areas, seasons extended through two consecutive years, and a high variability between different seasons was present.

Estimated Hidden Markov Models Parameters and Inferred Hidden States
Parameters for HMM estimated for each study area are presented in Table 2.As shown in Table 2, HMM parameters reflect the differences in NDVI seasonality between croplands and natural

Estimated Hidden Markov Models Parameters and Inferred Hidden States
Parameters for HMM estimated for each study area are presented in Table 2.As shown in Table 2, HMM parameters reflect the differences in NDVI seasonality between croplands and natural vegetation, with higher values in Albufera, the site with rice crops, for transition probabilities between identical states, implying more well-defined sequences of states, and also much higher absolute values for mean NDVI changes (emissions) in the states corresponding to NDVI increasing (S + ) and decreasing (S − ) values.For each study area, inferred hidden states for each pixel in the area are presented in Figure 5.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 17 vegetation, with higher values in Albufera, the site with rice crops, for transition probabilities between identical states, implying more well-defined sequences of states, and also much higher absolute values for mean NDVI changes (emissions) in the states corresponding to NDVI increasing (S+) and decreasing (S-) values.For each study area, inferred hidden states for each pixel in the area are presented in Figure 5.
(a) (b) As clearly shown in Figure 5, inferred states in croplands exhibit short and well-defined compact growing seasons (Figure 5a, red), consistent across pixels and years, which also applies with a bit more variability to declining states (Figure 5a, blue).In natural land cover types (Figure 5b-e) seasons are much less defined, extend wider in time, and present a much higher variability both between pixels and between different years.

Methods
In rice crops, where growing and declining seasons are marked and short, we defined HMM-derived SOS and EOS as the median (50 th percentile) of DOY with corresponding inferred hidden states (S+ and S-, respectively).As clearly shown in Figure 5, inferred states in croplands exhibit short and well-defined compact growing seasons (Figure 5a, red), consistent across pixels and years, which also applies with a bit more variability to declining states (Figure 5a, blue).In natural land cover types (Figure 5b-e) seasons are much less defined, extend wider in time, and present a much higher variability both between pixels and between different years.

Comparison of Phenology Metrics Derived from Hidden Markov Models and Smoothing Methods
In rice crops, where growing and declining seasons are marked and short, we defined HMM-derived SOS and EOS as the median (50th percentile) of DOY with corresponding inferred hidden states (S + and S − , respectively).
Table 3 presents Pearson correlation coefficients between mean SOS and EOS values (averages of SOS or EOS for all pixels) derived from different methods in Albufera site.Correlations vary in magnitude and significance between different pair comparisons, with those corresponding to HMM derived SOS showing, in general, intermediate values, significant in most cases, and lower values for EOS, but still significant when compared to greenbrown with linear interpolation.For the study areas with natural land cover types, with more fuzzy and extended seasons, we defined SOS and EOS from HMM as the first quartile (25th percentile) of DOY with hidden states S + and S − , respectively.Correlations between HMM derived SOS in these areas, in comparison with those derived from other methods, are presented in Table 4. Correlation coefficients for phenology metrics from HMM and TIMESAT showed high values, which were highly significant in all cases.Comparisons between either HMM or TIMESAT and the other three methods showed rather similar correlations in all four land cover types.

Performance of the Methods on Phenology Metrics Definition and Variability
We analyzed, for each method, site and year, the number of pixels where the method was not able to derive a defined metric, resulting in a NA (not available) metric.Table 5 shows the proportion of NA pixels for SOS metric, as the average over all the years.SOS metrics derived from HMM are always defined, as they only depend on the presence of inferred growing states.Except for TIMESAT in Albufera, there are always a certain proportion of pixels with NA metrics for the different smoothing methods, usually higher in natural land cover types.Mean yearly values of coefficients of variation (COV) among pixels, for the different study areas and methods, are presented in Table 6.While a higher variability among pixels in different areas could reflect a real higher spatial variability, higher values of COV could also be the result of a lower performance in the precision of the estimation method.As shown in Table 6, HMM and TS metrics exhibited, in most cases, lower values of COV, with either of them being the lowest depending on the study area.

Relation of Phenology Metrics with Climate Variables
Lacking on-the-ground phenology measurements to compare with those from LSP, the presence of significant relations between LSP derived metrics and climate variables, as well as their magnitudes, could be considered as indicative of their validity.Figure 6 presents the relations between HMM derived SOS metrics and selected climate variables in Níjar site, showing significant linear dependencies, with previous cumulative precipitation being related to earlier SOS, while differences in temperature between September and spring were related to higher values.

Relation of Phenology Metrics with Climate Variables
Lacking on-the-ground phenology measurements to compare with those from LSP, the presence of significant relations between LSP derived metrics and climate variables, as well as their magnitudes, could be considered as indicative of their validity.Figure 6 presents the relations between HMM derived SOS metrics and selected climate variables in Níjar site, showing significant linear dependencies, with previous cumulative precipitation being related to earlier SOS, while differences in temperature between September and spring were related to higher values.
(a) (b) Coefficients of determination (r 2 values) for similar relations between SOS metrics derived from different methods and climate variables are presented in Table 7. Apart from significant relations for HMM derived SOS with both precipitation and temperature, of the different smoothing methods only TS derived SOS was significantly related with precipitation, and in all cases their relations exhibited lower r 2 values than those from HMM.  Coefficients of determination (r 2 values) for similar relations between SOS metrics derived from different methods and climate variables are presented in Table 7. Apart from significant relations for HMM derived SOS with both precipitation and temperature, of the different smoothing methods only TS derived SOS was significantly related with precipitation, and in all cases their relations exhibited lower r 2 values than those from HMM.

Discussion
HMM derived phenology metrics, in different land cover types, were consistent with those derived from usual smoothing methods, providing estimations for all pixels in all of the study areas considered in this work, even when comparison methods did not.
Usual curve fitting methods can be applied either to NDVI mean curves for a whole area or to individual pixel curves.In the first case, no information on the spatial distribution of phenology metrics is obtained, while in the second, the fitting of individual curves might result in high variances and the lack of metrics for a number of pixels, especially with noisy individual NDVI curves.Conversely, there are two scales involved in the HMM fitting process.Given an area with homogeneous NDVI behavior, the whole set of pixels in the area can be used to obtain estimations of the HMM parameters, which makes the estimation process and the resulting model very robust.This model is then used to infer the sequence of hidden states for each pixel in the area, which is the base to obtain the LSP metrics at the pixel level.
The spatial distribution of the resulting HMM derived metrics is expected to be well-behaved, assuming that NDVI curves in the area are sufficiently homogeneous.In this work, we wanted to evaluate the use of HMM for LSP in different land cover types, and so we selected a priori homogeneous areas, and also checked, a posteriori, their homogeneity.In general, as suggested and carried out in References [54,56], HMM could be applied in large heterogeneous areas by using a previous automatic clustering process of the pixels based on NDVI similarities.
Despite the HMM being proposed very early as a tool to investigate vegetation dynamics from remote sensing observations [54], and although there have been a few examples of applications of HMM in phenology in recent years [47,49,55,56], to our knowledge there has been no previous published evaluation of the performance of HMM to define phenology metrics in comparison with other methods.In this work we have employed a very simple model of HMM and tested its ability to provide phenology metrics in comparison to ready applications of different smoothing methods, without trying to optimize any of the methods or adjusting them to particular study areas.Thus, our results should not be considered as an evaluation of the smoothing methods used as reference, but as indication that HMM are able to provide consistent results.
Since NDVI and vegetation phenology are expected to be related to climate variables (e.g., References [68,69]), the relations of SOS derived from HMM and other methods with climate variables presented in Figure 6 and Table 7 are indications of the validity of HMM derived metrics.Notwithstanding, further work including comparisons with ground-based data would be required to confirm these indications.The results presented in this paper strongly support the potential of HMM as an alternative or complementary method for deriving LSP metrics, thus opening new lines of research and possible improvements by using many of the variants and extensions of the basic type of HMM used in this work.For instance, we have used the same topology of the HMM in all land cover types, but they could be adapted to particular areas with different vegetation types, either in the number of states or the allowed transitions between states.Additionally, we have employed homogeneous HMM, i.e., with transition probabilities between states constant in time, but they can be allowed to be time-dependent [55].More complex models, as higher order or semi-Markov models [70], could also be used to model the duration of the different phases of the season.

Conclusions
Our results show that HMM derived phenology metrics can be obtained in large areas with different land cover types, producing consistent values with those obtained from commonly employed smoothing methods.HMM were able to provide metrics estimations for all pixels, while smoothing methods partly failed in some cases.Spatial variability of HMM metrics were similar to the best behaving smoothing method, and lower than other methods in many cases.Correlations of HMM derived SOS metrics with climate variables were similar or higher than those with the best smoothing methods, suggesting the validity of the HMM derived metrics.As a general conclusion, our results indicate that HMM can be very useful for the derivation of LSP metrics, providing pixel-wise estimations in large areas, although further work is needed to assess their performance in relation with on-the-ground observations.Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/5/507/s1,

Figure 1 .
Figure 1.(a) Structure of a first order HMM (Hidden Markov Model), as a temporal sequence of hidden (St) and observed states (Yt); and (b) topology of the HMM used in this work to model NDVI changes with four states: increasing (S+), stationary high (S0+), decreasing (S-), and stationary low (S0-).

Figure 1 .
Figure 1.(a) Structure of a first order HMM (Hidden Markov Model), as a temporal sequence of hidden (S t ) and observed states (Y t ); and (b) topology of the HMM used in this work to model NDVI changes with four states: increasing (S + ), stationary high (S 0+ ), decreasing (S − ), and stationary low (S 0− ).
Valencia region, Albufera site comprises an area of traditional rice crops surrounding a coastal lagoon; Millares and Ayora sites are inland mountain areas, with altitudes ranging from 600 to 800 m in Millares, and 700 to 1000 m in Ayora, with vegetation cover dominated by shrublands.In Murcia region, Espuña site, located in Espuña range, with altitudes ranging from 650 to 1300 m, is a mixed pine forest derived from a reforestation project dating back to the end of the 19 th century.In Andalucía region, the Níjar site is located close to the Southern Mediterranean coast, with altitudes ranging from 250 to 900 m with grass steppes as the main land cover.

Figure 2 .
Figure 2. Location of the five study areas in Southeast Spain.Borders of Spain's autonomous regions are displayed.

Figure 2 .
Figure 2. Location of the five study areas in Southeast Spain.Borders of Spain's autonomous regions are displayed.

Figure 3 .
Figure 3. Correlations between annual component of NDVI (Normalized Difference Vegetation Index) time series for each pixel and global mean annual curve in Níjar site (grass steppe).Pixel position in UTM (Universal Transverse Mercator) coordinates.

Figure 3 .
Figure 3. Correlations between annual component of NDVI (Normalized Difference Vegetation Index) time series for each pixel and global mean annual curve in Níjar site (grass steppe).Pixel position in UTM (Universal Transverse Mercator) coordinates.

Figure 6 .
Figure 6.Relations between HMM derived mean SOS metrics and climate variables in Níjar site; (a) Relation with precipitation (cumulative August-October); and (b) relation with temperature (difference between mean temperature in September and average of mean temperatures in February-April).

Figure 6 .
Figure 6.Relations between HMM derived mean SOS metrics and climate variables in Níjar site; (a) Relation with precipitation (cumulative August-October); and (b) relation with temperature (difference between mean temperature in September and average of mean temperatures in February-April).

Figure S1 :
Correlations between annual component of NDVI time series for each pixel and global mean annual curve in Albufera, Millares, Ayora, and Espuña sites, File S1: Compressed file including kmz files for the five study areas.Author Contributions: Conceptualization, M.A.G, F.R. and S.B.; Methodology and Formal Analysis, M.A.G., H.M., G.M.C. and F.R.; Writing-Original Draft Preparation, F.R.; Writing-Review & Editing, all the authors; Visualization, M.A.G. and H.M.; Funding Acquisition, F.R. and S.B.

Table 1 .
Basic data on the study areas: Name of the site, after geographical reference (Albufera lagoon, Espuña range) or nearby town (Millares, Ayora, Níjar); central approximate coordinates; number of selected pixels; and type of land cover.

Table 1 .
Basic data on the study areas: Name of the site, after geographical reference (Albufera lagoon, Espuña range) or nearby town (Millares, Ayora, Níjar); central approximate coordinates; number of selected pixels; and type of land cover.

Table 2 .
Estimated parameters for HMM fitted in the five study areas.Transition probabilities between hidden states (from state in column to state in row), and mean and standard deviation (SD) of Gaussian emissions for each hidden state.

Table 2 .
Estimated parameters for HMM fitted in the five study areas.Transition probabilities between hidden states (from state in column to state in row), and mean and standard deviation (SD) of Gaussian emissions for each hidden state.

Table 3 presents
Pearson correlation coefficients between mean SOS and EOS values (averages of SOS or EOS for all pixels) derived from different methods in Albufera site.Correlations vary in magnitude and significance between different pair comparisons, with those corresponding to HMM derived SOS showing, in general, intermediate values, significant in most cases, and lower values for EOS, but still significant when compared to greenbrown with linear interpolation.

Table 3 .
Correlations between phenological metrics derived from different methods in the site with rice crops (Albufera).Correlation coefficients for mean values of SOS (upper pair comparisons, grey shaded) and EOS (lower pair comparisons), and statistical significance.

Table 3 .
Correlations between phenological metrics derived from different methods in the site with rice crops (Albufera).Correlation coefficients for mean values of SOS (upper pair comparisons, grey shaded) and EOS (lower pair comparisons), and statistical significance.

Table 4 .
Correlation coefficients between mean values of SOS (Start of Season) metrics derived from different methods in the four sites with natural land cover types.

Table 5 .
Proportion of pixels (%; mean yearly values) with NA SOS metrics derived from different methods in the five study areas.

Table 6 .
Coefficients of variation among pixels (%; mean yearly values) for SOS metrics derived from different methods in the five study areas.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 17 metrics exhibited, in most cases, lower values of COV, with either of them being the lowest depending on the study area.

Table 6 .
Coefficients of variation among pixels (%; mean yearly values) for SOS metrics derived from different methods in the five study areas.

Table 7 .
Linear regression of mean SOS metrics derived from different methods on climate variables in Níjar site.Coefficients of determination and statistical significance of Simple Linear Regression.

Table 7 .
Linear regression of mean SOS metrics derived from different methods on climate variables in Níjar site.Coefficients of determination and statistical significance of Simple Linear Regression.