Applicability of Smoothing Techniques in Generation of Phenological Metrics of Tectona grandis L. Using NDVI Time Series Data

Abstract

Information on phenological metrics of individual plant species is meager. Phenological metrics generation for a specific plant species can prove beneficial if the species is ecologically or economically important. Teak, a dominating tree in most regions of the world has been focused on in the present study due to its multiple benefits. Forecasts on such species can attain a substantial improvement in their productivity. MODIS NDVI time series when subjected to statistical smoothing techniques exhibited good output with Tukey’s smoothing (TS) with a low RMSE of 0.042 compared to single exponential (SE) and double exponential (DE). Phenological metrics, namely, the start of the season (SOS), end of the season (EOS), maximum of the season (MAX), and length of the season (LOS) were generated using Tukey-smoothed MODIS NDVI data for the years 2003–2004 and 2013–2014. Post shifts in SOS and EOS by 14 and 37 days respectively with a preshift of 28 days in MAX were observed in the year 2013–2014. Preshift in MAX was accompanied by an increase in greenness exhibiting increased NDVI value.LOS increased by 24 days in the year 2013–2014, showing an increase in the duration of the season of teak. Dates of these satellite-retrieved phenological occurrences were validated with ground phenological data calculated using crown cover assessment. The present study demonstrated the potential of a spatial approach in the generation of phenometrics for an individual plant species, which is significant in determining productivity or a crucial trophic link for a given region.

1. Introduction

The occurrence dates of phenophases such as blooming, full leaf expansion, leaf coloration, or senescence are keys for the determination of phenological metrics, viz., start of the season (SOS), end of the season (EOS), and maximum of the season (MAX) of any tree [1]. These metrics can prove to be of crucial importance in the tree’s productivity assessment. They can also prove significant in modeling and monitoring climate change [2], explaining the seasonal changes, determining net terrestrial carbon dioxide flux [3,4,5,6], and assessing leaf cycle functioning.

Information is meager on phenological data of individual plant species in many parts of the world [7]. More focus is placed on generating phenological metrics of forest types [8,9,10,11] rather than on a particular natural species [12,13]. Hence, evaluating the capability of satellite data for effective and reliable phenological metrics generation for one individual forest species, namely, Tectona grandis L. (teak) will be of great importance. Teak, a deciduous tree of the Verbenaceae family, is considered a valuable timber and has been widely used in India for more than 200 years due to its elegance and durability. In addition, it is an economically important tree as it has a tight grain and high oil content and tensile strength and can also be grown under various climatic regimes. Moreover, this tree occupies major worldwide areas and is also a dominant tree of tropical dry deciduous forests (DDF) of the Narmada district. Due to a paucity of ground data on this system, there is a great need to develop remote sensing techniques to characterize its biology. Phenological metrics generation of teak growing in the DDF can provide useful information about both spatial and temporal patterns of its productivity. Such forecasts can attain a substantial improvement in its productivity.

Appraisal of phenophases of any plant species mainly involves two approaches; ground-based phenology measurements and satellite-based phenology monitoring [14]. Ground-based measurements entail direct traditional insitu visual recordings of phenological events [15,16,17,18,19], manual periodic photography, or fixed-position camera-based digital repeated photography [20,21,22]. These methods, though they provide good outputs, have their limitations concerning lack of synoptic coverage and time [23]. An alternative space approach overcomes such limitations both with respect to time and coverage. Recent studies, therefore, are based on using the potential of spatial data in understanding plant phenology [24,25,26,27,28,29,30,31]. Specifically, time series normalized difference vegetation index (NDVI) data derived from moderate resolution imaging spectroradiometer (MODIS) [32,33,34,35,36,37,38] and advanced very high resolution radiometer (AVHRR) [39,40,41,42,43,44] are largely used for the generation of phenological metrics. This satellite-based vegetation index is recognized as an effective tool in the quantification of vegetation greenness [45,46,47], canopy carbon [48,49,50,51,52], and water fluxes [53,54,55]. It has proved its utility in monitoring the phenological changes in vegetation across large spatial scales and over long time periods. Many researchers have used NDVI and enhanced vegetation index (EVI) time series in the identification of different phenological parameters at both single pixel scale and large spatial scales by developing various methodologies [25,56,57,58,59,60,61]. Studies are also available wherein Landsat time series data have been used for retrieving phenological metrics [62,63,64]. Some have built statistical relationships between NDVI or EVI and climatic parameters for monitoring the entire growing season [65,66,67,68]. The phenological parameters derived from vegetation indices may not correspond directly to conventional ground-based phenological events but do provide indications for understanding these phenophases which now have become pertinent in climate change analysis [69,70]. In addition, differences among biomes and their drivers of phenology (e.g., dry deciduous vs. cold deciduous) may necessitate different interpretations of image-derived phenophases. As such, detailed comparisons of these satellite measures and ground-based phenological events are needed across a range of biomes with different ecological drivers.

The time series NDVI or EVI data when subjected to smoothing techniques remove the outliers [71] occurring due to cloud contamination, atmospheric perturbations, variable illumination, and viewing geometry. These include methods like moving average time series (MATS) [72], a weighted least-squares linear regression approach [73]. Fourier harmonics-based methods (classical and discrete) [74,75,76,77,78], threshold-based methods [79,80,81], curve fitting methods or fitting of polynomial functions [82,83], point of inflection methods [81,84], simple sliding windows [85], piecewise linear regression [86,87], single exponential [88], double exponential [89], Tukey smoothing [88], etc.

The objective of the present work is to: (1) investigate three statistical smoothing techniques, viz., single exponential (SE), double exponential (DE), and Tukey’s Smoothing (TS) for diminishing noise effects and removing outliers in MODIS NDVI time series; (2) use the best reconstructed time series NDVI MODIS data to generate different phenological metrics, namely, SOS, EOS, MAX, and LOS of teak growing in Narmada forests, Gujarat, India.

2. Materials and Methods

2.1. In Situ Phenological Data Collection

The study area selected for the present study is DDF of Narmada district, Gujarat, India (21.86 N, 73.56 E) (Figure 1). This area exhibits the dominance of teak trees. Occurrence dates of different phenophase events for the year 2003–2004 were obtained from Sardar Sarovar Narmada Limited (SSNL) report [90]. The data for the year 2013–2014 were generated by monitoring ten 30 m × 30 m (900 m²) homogenous teak patches at regular intervals of 20–25 days from May 2013 till April 2014. Analysis of within-population phenophase frequency at every sampling date was performed based on the modified qualitative method of Ohshan [91]. The population was randomly selected and labeled before the beginning of the sampling for ten healthy adult teak trees. During each field visit, the degree of incidence was determined for the different phenophases of teak, viz., leaf initiation, maximum greenness, and leaf fall, for which each tree was examined precisely. A frequency index was assigned for particular phenophases at the time of their presence on a tree, depending on the percentage of the crown where it occurred: 1 = presence less than 5% (leaf fall), 2 = presence between 5–25% (leaf initiation), and 3 = presence over 25% of the crown (maximum greenness).The average of the ten sampled plants’ frequency indicesrepresents the calculation of the incidence of a phenophase in the whole population for each date. These ground-based phenological observations were then compared to phenological metrics derived using time series satellite data for the validation purpose.

2.2. Satellite Data

Annual spatiotemporal MODIS vegetation indices-terra (MOD13A1)’s NDVI product of the year 2003–2004 and 2013–2014 with repetivity of 16 days, and 500 m spatial resolution was used to extract the phenological metrics of teak. Before the extraction of these metrics, NDVI time series data were subjected to different smoothing techniques.

2.2.1. Application of Smoothing Techniques

Different smoothing algorithms were applied to the time series NDVI MODIS data to diminish noise effects: single exponential (SE), double exponential (DE), and Tukey’s smoothing (TS). All three techniques are in ascending order in terms of sophistication.

Single Exponential Smoothing (SE)

SE is a simple and accessible tool for smoothing time series data. A simple average calculation is used to assign exponentially decreasing weights, starting with the most recent observations. New observations are weighted more in the average calculation than earlier observations. This method is for univariate data and does not include trend or seasonality. It uses only one parameter named alpha (α), even known as smoothing factor or smoothing coefficient. This smoothing technique is widely employed due to its simplicity and success. The notation for the same is as given below:

(1)$S_{t+1}=\alpha y_t+(1-\alpha )S_t$

where S_i is the smoothed value of time series at time i, y_i is the actual value of time series at time i, and α is the smoothing constant and 0.0 < α < 1.0.

Double Exponential Smoothing (DE)

DE is used on data sets involving seasonality and for handling trend analysis. It is an extension to exponential smoothing, adding explicit support for trends in the univariate time series. It is used when there is a linear trend in the data. This involves an additional smoothing factor along with alpha (α) parameter. This is to control the decay of the influence of the change in a trend called beta (β).

For data exhibiting linear trend as:

(2)$y_t=b_0+b_{1}t+e_t$

where b₀ and b₁ can change with time at a slow pace. The basic equations named Holt’s method are as below:

(3)${\mu }_t=\alpha y_t+(1-\alpha )({\mu }_{t-1}+T_{t-1})$

(4)$T_t=\beta ({\mu }_t-{\mu }_{t-1})+(1-\beta )T_{t-1}$

where µ_t is the exponentially smoothed value of time series at time t, y_t is the actual observation of time series at time t, T_t is the trend estimate, α is the exponential smoothing constant for the data, and β is the smoothing constant for the trend.

Tukey’s Smoothing (TS)

TS uses running medians to provide flexible but straightforward curves and is a robust smoother. Median smoothing methods were introduced by Tukey in 1977 for extracting smooth patterns which tend to hide due to non-linear spikes in time series data [92]. Such filtering smooths any existing volatile behavior that occurs in trends or seasonal behavior. This method smooths out the data acquired from equally spaced, linearly ordered intervals such as every year, every month, every quarter, etc.

NDVI data reconstructed using all three types of smoothing techniques were compared based on their potential in removing outliers.

Performance assessment of the smoothing techniques was also carried out by calculating root mean squared error (RMSE). RMSE between observed raw NDVI time series data and predicted smoothed NDVI time series data was used for evaluating the performance of each smoothing technique. The best optimal technique that would reconstruct the best denoised data sets was considered to be the one generating the least RMSE. The best-reconstructed data was used further for extracting phenological metrics of teak by detecting the inflection point (i.e., date) when the NDVI time series begins to ascend or descend for the specific year. Occurrence dates obtained using these smoothed NDVI time series were compared to the ground-based phenological observation and these along with RMSE were considered for selecting the optimum smoothing technique.

2.2.2. Determination of Phenological Metrics from NDVI Time Series Data

The intra-annual variations in the NDVI time series were used as the base for determining SOS, MAX, EOS, and LOS, using the NDVI ratio. This is the derivative method where the maximum value of NDVI ratio corresponds to the greatest change of NDVI time series. Equation (5) is given as:

(5)$NDVIratio(t)=\frac{NDVI(t+1)-NDVI(t)}{NDVI(t)}$

where NDVI(t) is the NDVI value at time t, and NDVI ratio (t) is the calculated relative change at time t. SOS was determined as the time t or the day with the maximum NDVI ratio [93]. Likewise, EOS date was determined as the time t or the day having the minimum NDVI ratio. The time or duration between SOS and EOS with NDVI ratio closest to zero was identified as MAX date. Furthermore, the LOS was obtained by defining the period between SOS and EOS in each grid point.

3. Results and Discussion

RMSEs generated to check the potential of SE, DE, and TS are 0.072, 0.097, and 0.048, respectively. Comparison of reconstructed NDVI time series using SE, DE, and TS techniques showed that both SE and DE smoothed higher values but could not encompass all the outliers (Figure 2 and Figure 3). These methods are comparatively less effective in accounting for missing values or correcting outliers. TS smoothing approach is identified as a more effective and robust smoothing method, bringing out high-quality NDVI time series as it fairly handled the outliers (Figure 4). The scatterplot diagram distinctly highlights the presence of outliers that needed to be removed to generate precise output. Statistically, the smoothing techniques are supposed to remove these outliers. At this point, the selection of Tukey’s technique is performed in comparison to the other two smoothing techniques as it is considered to be a resistant smoothing technique which uses running medians. Despite the fact that it lacks mathematical generalization, the main purpose of the smoothing technique is to provide a general idea of relatively slow changes in values with little attention to the close matching of data values. Phenological attributes exhibit such characteristics. Further, running medians are thought to be fast exploratory tools to allow a quick view of data components. On comparing with ground data, TS again exhibited the greatest effectiveness over SE and DE (Figure 5, Table 1). The Tukey-smoothed NDVI data provided greater vegetation phenology information, i.e., even minor changes in phenological dates which cannot be correctly identified from raw NDVI data can be obtained from Tukey-smoothed NDVI data. TS retrieved phenological data and ground phenological data showed very close values, thereby validating the TS results. Thus, TS potentially proved to be the optimal technique to best reconstruct the NDVI time series data and hence TS reconstructed NDVI time series is used for delineating shifts in SOS, EOS, MAX, and LOS of teak.

MODIS data analyzed for the monitoring of the phenological metrics of teak for the years 2003–2004 and 2013–2014 showed significant patterns despite their coarse resolution. Annual time series of NDVI data generated for teak enabled the differentiation of various phenological metrics such as the SOS (Figure 6 and Figure 7), EOS (Figure 8 and Figure 9), MAX (Figure 10 and Figure 11), and LOS (Figure 12 and Figure 13). Such data applications have proven to be useful in several studies on global environmental change [94,95,96,97,98].

Differences in phenological patterns are notable. SOS and EOS in DDF teak were delayed by 14 and 37 days, respectively, in the year 2013–2014 (SOS date—11 August 2013; EOS date—22 March 2014) when compared to 2003–2004 (SOS date—28 July 2003; EOS date—13 February 2004). The SOS in the year 2003–2004 was observed in late July which shifted to mid-August in 2013–2014. The EOS was observed in mid-February in 2003–2004 while in 2013–2014, the season ended in the end of March. Advancement of 28 days in MAX of teak was observed in the year 2013–2014 compared to 2003–2004. Greenness in teak in the district during the year 2003–2004 reaches its peak at the end of September, while in 2013–2014, teak reached its maximum greenness at the end of August. Results highlighted the fact that teak reached its MAX early in 2013–2014, indicating the shift in the phenology of the tree. LOS increased by 24 days in the year 2013–2014 (LOS—224 days) (Figure 14) compared to the year 2003–2004 (LOS—200 days) (Figure 13). Except for MAX, the year 2003–2004 showed earlier dates than the year 2013–2014 for the occurrence of all phenological metrics.

For the validation of the phenological metrics generated from MODIS TS NDVI time series, the dates of different phenological occurrences are compared to ground measured crown cover. In situ crown cover assessment of 2013–2014 provides the information on different phenophases (Table 1). A plot of crown cover versus sampling dates matches the NDVI curve of the period (Figure 14). Results of the crown cover assessment can be used for the validation of the results generated through MODIS data. Average crown cover between 5 and 25% that corresponds to leaf initiation on ground and SOS on satellite data was observed on 30 July 2013. Maximum greenness or MAX that corresponds to crown cover over 25% was noted on 27 August 2013 and leaf fall on the ground, i.e., less than 5% crown cover that corresponds to EOS was recorded on 6 March 2014. These ground-based phenological observations are comparable to phenological metrics derived using time series data, thereby validating our results. Results of in situ observations also show that teak trees considered for the present study shifted their phenophases between the years 2003–2004 and 2013–2014. Leaf initiation on the ground showed post shift of 8 days. In situ leaf fall recorded showed a delay of 24 days. Maximum greenness noted on the ground also compliments satellite-retrieved MAX results which showed a preshift of 15 days. The results highlight that phenological metrics from satellite data such as SOS, EOS, and MAX with temporal dynamics similar to those of ground phenology measurements such as leaf initiation, leaf fall, and maximum greenness through crown cover assessment can be generated.

Substantial interannual variability in the start and end of the growing season of the teak of DDF forests can be explained by studying many factors such as year-to-year variability in weather, especially climatic factors including rainfall, temperature, elevated CO₂, altered precipitation regimes, etc. Such types of variability may be responsible for the observed shifts in phenology which ultimately can be robust indicators of the impacts of climate change. This makes derivation of phenological metrics imperative as they can serve as important inputs for better understanding of the climatic drivers controlling phenology. These alterations or shifts in phenophases can influence the climate system at both global and regional levels through feedback mechanisms of surface albedo, CO₂ fluxes, and evaporation [99]. Improved inputs of phenological measurements into global biosphere models conducted via satellite data will also enhance the understanding of temporal and spatial global carbon dynamics. This will be more useful for the areas where there exists a lack of ground phenological observations due to inaccessibility [38].

4. Conclusions

The present study tested the potential of three statistical techniques, viz., single exponential, double exponential, and Tukey smoothing. Comparison of these three techniques demonstrates the effectiveness of the Tukey smoothing technique in denoising the NDVI time series and removing the outliers. Statistical Tukey smoothing enhances the visibility of unique patterns present in the data and thereby aids in identifying distinct decadal shifts in the phenological metrics. Significant shifts in the growing seasons of DDF teak of Narmada district from the years 2003–2004 to 2013–2014 are identified on the basis of changes in SOS, EOS, LOS, and MAX delineated using Tukey-smoothed MODIS NDVI time series data. These shifts are validated by ground phenology measurements that are calculated using crown cover assessment. Earlier studies using remote sensing more commonly focused on phenological studies of forest types but the present study highlights the fact that satellite data has great capability and can be used for generating phenological metrics of individual plant species as well.

Author Contributions

Conceptualization, G.S.K. and R.K.M.M.; methodology, R.K.M.M., M.N.S. and P.A.T.; software, N.V.M., V.H.B. and R.K.M.M.; validation, N.V.M., V.H.B., R.K.M.M., S.M.; N.V.M., V.H.B. and R.K.M.M.; investigation, R.K.M.M., N.V.M., V.H.B. and G.S.K.; data curation, G.S.K.; writing—original draft preparation, R.K.M.M. and G.S.K.; writing—review and editing, G.S.K., P.A.T. and S.M.; visualization, S.M.; supervision, C.P.S. and B.K.B.; project administration, G.S.K.; funding acquisition, G.S.K. and M.N.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by ISRO RESPOND, project id OGP134” and we give thanks to MoEFCC for their added support.

Institutional Review Board Statement

Ethical review and approval were waived for this study.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

In this section, I acknowledge the funds received from ISRO RESPOND for carrying out the present work. I also extend my thanks to MoEF&CC for data and field work support.

Conflicts of Interest

Authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. A location map showing the study area, Narmada district, Gujarat, India; image used in the study area map is a Landsat 11 RGB composite.

Figure 2. NDVI time series for year 2013–2014 reconstructed using single exponential smoothing.

Figure 3. NDVI time series for year 2013–2014 reconstructed using double exponential smoothing.

Figure 4. NDVI time series for the year 2013–2014 reconstructed using Tukey’s smoothing.

Figure 5. Different smoothing techniques’ utility in derivation of maximum of the season.

Figure 6. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2003–2004 for derivation of start of the season in Tectona grandis indicated by green arrow.

Figure 7. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2013–2014 for derivation of start of the season in Tectona grandis indicated by green arrow.

Figure 8. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2003–2004 for derivation of end of the season in Tectona grandis indicated by orange arrow.

Figure 9. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2013–2014 for derivation ofend of the season in Tectona grandis indicated by orange arrow.

Figure 10. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2003–2004 for derivation of maximum of the season in Tectona grandis indicated by orange arrow.

Figure 11. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2013–2014 for derivation of maximum of the season in Tectona grandis indicated by orange arrow.

Figure 12. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2003–2004 for derivation of length of the season in Tectona grandis indicated between green and orange arrows.

Figure 12. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2003–2004 for derivation of length of the season in Tectona grandis indicated between green and orange arrows.

Figure 13. Raw MODIS NDVI and TS-smoothed reconstructed NDVI time series of year 2013–2014 for derivation oflength of the season in Tectona grandis indicated between green and orange arrows.

Figure 14. Crown cover in percentage during in situ sampling period.

Table 1

Phenological metrics derived using different data at different time periods.

Phenophases		SOS	EOS	MAX
Ground data	2003–2004	22 July	10 February	15 September
	2013–2014	30 July	6 March	31 August
	2015–2016	1 August	23 March	27 August
Tukey-smoothed NDVI	2003–2004	28 July	13 February	26 September
	2013–2014	11 August	22 March	29 August