Assessment of values and trends in coarse spatial resolution NDVI datasets in Southeast Asia landscapes

ABSTRACT Normalized difference vegetation index (NDVI) has been widely applied for monitoring vegetation dynamics. However, NDVI values are known to be profoundly affected by various external factors. In this study, the variation of NDVI values and trends among the several long-term NDVI datasets with resolution of 1, 4 and 8 km were assessed to understand the differences between the available datasets. The assessment items were 1) Pearson’s correlation coefficient, 2) trend map and breakpoint spatial similarities and 3) comparison of NDVI from Landsat and Flux tower in 2007–2015. The comparison revealed a maximum correlation coefficient of 0.67 among NDVI datasets and average spatial similarity of 37.2% among the trend maps estimated from NDVI datasets. Furthermore, there was a possibility of having significantly opposite trends between two trend maps from different NDVI products. Comparisons with NDVI from vegetation pixel in Landsat 5 TM and 8 OLI resulted in the R2 between 0.06 and 0.68 and RMSE of 0.07–0.2, while comparison with NDVI from flux tower data yielded the RMSE of 0.04–0.41, although the R2 was relatively weak at 0–0.18. Our study highlights the possibility of differences among NDVI datasets, and suggests that these differences should be reconciled especially in time-series analysis.


Introduction
Monitoring and measuring the vegetation over large areas using satellite remote sensing has frequently been performed. To correlate values obtained by satellite imaging with actual vegetation characteristics, vegetation indices were developed as spectral transformations that mostly utilize the red and infrared bands, which are sensitive to leaf pigments. Among the different vegetation indices, one index that has been widely used for long-term datasets is the normalized difference vegetation index (NDVI). NDVI can be used to analyze various aspects of vegetation such as vegetation dynamics, productivity and the effects of climate on vegetation, allowing mapping of the spatial and temporal distributions of vertical and horizontal structures such as vegetation community and vegetation biomass (Pettorelli et al., 2005) and land cover mapping (Defries & Townshend, 1994).
NDVI has been viewed as a proxy for vegetation greenness and has been widely used for predicting biophysical properties of vegetation such as the leaf area index (LAI) (Turner, Cohen, Kennedy, Fassnacht, & Briggs, 1999;Wang, Adiku, Tenhunen, & Granier, 2005) and the fraction of photosynthesis active radiation (FAPAR) (Myneni & Williams, 1994) across different vegetation landscapes. Although NDVI has been successfully applied, the values are known to be sensitive to external factors such as soil and canopy background, atmospheric disturbances, and sun-sensor viewing angle (Baret & Guyot, 1991;Epiphanio & Huete, 1995;Myneni & Williams, 1994;Roujean, Leroy, & Deschamps, 1992). Moreover, application of NDVI in areas with dense vegetation results in signal saturation, particularly for LAI values above 2-3 (Carlson & Ripley, 1997), and thus makes it less likely to provide accurate information of biophysical features vegetation and phenology in this particular area.
While instantaneous NDVI values are also susceptible to the same problems, time series NDVI observations are sensitive to sensor changes, especially when data from several satellite missions are used to construct the time series observations, particularly in advanced very high-resolution radiometer (AVHRR)based NDVI data and to orbital drift, which affects the solar zenith angle (Nagol, Vermote, & Prince, 2014). Further, inconsistency in time series NDVI observations from different sensors also appear due to differences in the spectral response function (SRF) in the red and infrared bands, and this adds another level of complexity for datasets collected from multiple sensors (Teillet, Staenz, & William, 1997;Trishchenko, Cihlar, & Li, 2002). Thus, evaluation of these datasets for use with long-term NDVI applications is needed to better understand the differences and deviations between pairwise comparisons of NDVI datasets.
Assessment of consistency of NDVI application across sensors has been performed by several methods such as pairwise comparison of sensors (Detsch, Otte, Appelhans, & Nauss, 2016), comparison with in situ measurements (Fensholt, Sandholt, & Stisen, 2006), trend analysis (Alcaraz-Segura, Liras, Tabik, Paruelo, & Cabello, 2010;Detsch et al., 2016;Rasmus Fensholt et al., 2012;Rasmus Fensholt & Proud, 2012;Rasmus Fensholt, Rasmussen, Nielsen, & Mbow, 2009) and detecting structural breaks from sensor measurements which correspond to sensor change, orbital drift or sensor degradation (Tian et al., 2015). Analysis performed in previous studies showed an artefact in the data due to sensor change in AVHRR-based NDVI applications, and different NDVI responses and trends were estimated from different sensors. Nonconformity of trends in NDVI time series, with even competing greening and browning trends reported in some studies, has triggered many debates (Alcaraz-Segura et al., 2010;Baldi et al., 2008). Likewise the finding of greening in the Amazon forest by the enhanced vegetation index (EVI) from MODIS (Saleska, Didan, Huete, & Da Rocha, 2007) countered the finding of Morton et al. (2014) that greening was merely an artefact of sunsensor viewing angle and that greenness remains consistent throughout the dry season. Based on these examples, the results from different remote sensors need to be compared before drawing final conclusions about vegetation dynamics.
Studying the trends in vegetation is important for understanding the response of vegetation to climate change and identifying the relationships between greening and browning to climate variables (Cong et al., 2013). The Indonesia and Southeast Asia region is a hotspot for the distribution of tropical forest areas that are crucial for the global carbon cycle. However, the forest and vegetation ecosystems in this area are susceptible to climate anomalies (Poussart, Myneni, & Lanzirotti, 2006) such as El Nino and La Nina, which affect temperature and rainfall intensity. Further regional studies using satellite data in this area have been undertaken to identify the trends and areas that are vulnerable to climate anomalies (Erasmi, Propastin, Kappas, & Panferov, 2009;Gutman, Csiszar, & Romanov, 2000). Both studies were based on data from the AVHRR-based NDVI products. However, an exploration of the consistency of NDVI data using different sensors in this region is needed, particularly because NDVI values in tropical forest areas were categorized as having great uncertainty (Brown, Pinzón, Didan, Morisette, & Tucker, 2006;Rasmus Fensholt & Proud, 2012).
Several long-term NDVI datasets obtained with AVHRR campaigns have been operational and are currently providing the access for monitoring vegetation phenology; these include GIMMS3G  and STAR NESDIS (Kogan, 2001) which have the longest continuous observations and newer long-term systems such as Système Pour l'Observation de la Terre (SPOT) -VGT, METOP, MERIS, PROBA V, MODIS TERRA -MOD13A3 (Didan, Munoz, & Huete, 2015), MODIS AQUA -MYD13A3 (Didan et al., 2015) and STAR NESDIS VIIRS. In this article, we performed several statistical analyses aiming to compare and contrast differences in NDVI values and trends among moderate and coarse long-term NDVI applications in South East Asia (Table 1). We used several commonly available NDVI applications, such as GIMMS3G ver.0 (Pinzon & Tucker, 2014), SPOT Vegetation (VGT) and MODIS-Terra and AQUA (Didan et al., 2015), and other newer applications such as METOP, MERIS, PROBA V, NESDIS AVHRR and NESDIS VIIRS. The analyses were conducted in three parts: 1) intra-comparison among NDVI data by correlation analysis, 2) comparison of coarse NDVI dataset values with NDVI observations from the flux tower and 3) spatial similarities analysis using trends estimated from the Mann-Kendall test. These analyses aimed to comprehensively assess variation among NDVI applications and their performance for determining trends of the pixel values.

Study areas
This study was conducted in a region of South East Asia centred on Indonesia and bounded by 94°to 154°E and 13°S to 13°N (Figure 1). In this region, there is one Flux Tower site located at Pasoh Forest Malaysia precisely at 102°E and 3°N that collects the daily NDVI observations from October 2008 to present.

Data
Coarse and moderate resolution NDVI data processed with different sampling time intervals and post-processing algorithms from the applications listed in Table 1 were retrieved from the respective online repositories. We aggregated weekly, trimonthly and bi-monthly data into monthly data using a maximum value composites (MVC) analysis (Holben, 1986). In the case of data sampled weekly, weeks with data in 2 months were treated as follows: NDVI values at the start and end of the month were linearly interpolated using the last weekly values, and assuming that the weekly data represented the end of the week value, data were composited into monthly data by taking the maximum value. Further, post-processing was performed on all datasets (except for the AVHRR NESDIS and VIIRS NESDIS applications) using the adaptive Savitzky-Golay filter in Timesat (Jönsson & Eklundh, 2004) to minimize possible effects of atmospheric disturbances such as clouds and haze. A filter of adjacent data points was applied to smooth time series data to avoid massive change of the original values. NDVI data from different satellite missions have different starting and ending dates, and also different temporal archives were selected to have the greatest number of datasets overlapping and the longest period of continuous data and statistical analysis of pairwise comparisons and trend analysis were performed. For pairwise comparison and trend analysis, the data period was divided into two sub-periods to permit analysis of newly released datasets such as Proba-V and VIIRS NESDIS and currently operating datasets such as AVHRR NESDIS, METOP, MODIS TERRA, MODIS AQUA, VIIRS NESDIS and PROBA V as follows: 1) April 2007 to March 2012 (60 months), and 2) January 2014 to December 2015 (24 months). In addition, pairwise assessments were performed for the period of November 2013 to May 2014 (7 months) to analyze the compatibility of SPOT VGT and PROBA-V.
There are three analyses conducted to compare and evaluate the similarity and performance of coarse NDVI products, such as 1) Pearson's correlation, 2) trend and breakpoint analysis and 3) comparison with vegetation indices data from moderate spatial resolution data and flux tower data.

Pairwise comparison between NDVI products
The pairwise comparison was used to estimate the degree of differences between NDVI observations for overlapping observations. This analysis was performed using Pearson correlation at several spatial resolutions. We divided the analysis into resolution categories of 8, 4 and 1 km and resampled data collected at 1 and 4 km resolution using the bilinear method to obtain data in the coarsest resolutions of 4 and 8 km resolutions categories. Further zonal analysis using the land cover map using the type 3 classification scheme excluding the water category ( Table 2) of MODIS land cover MCD12Q1 (Friedl et al., 2010) in 2013 was performed with an original spatial resolution of 500 m in order to understand the general agreement of correlation for different land cover categories.

Trend and breakpoint analysis
Aside from performing a correlation analysis, we also performed a trend analysis to estimate the directional changes in NDVI time series to estimate whether different sensors produced the same trend over the same observation period. Trend analysis was performed using the Mann-Kendall test in Rstudio utilizing the Kendall package (McLeod, 2005) after deseasonalizing NDVI time series data by subtracting the seasonal monthly means from the raw values over a same time periods (60 months from April 2007 to March 2012 and 24 months from Jan 2014 to December 2015). We used the Mann-Kendall test to estimate trends in time series data, including NDVI time series, because of its ability for detecting significant trends (Li et al., 2014), which are useful for comparing trend maps among different datasets.
Using the time series data, additional analysis to detect breakpoint or abrupt change from vegetation time series data can be performed. The breakpoints detected from vegetation indices including NDVI were identified to relate to the cover change (i.e. deforestation), and were also used for consistency assessment (Tian et al., 2015). In this study, the biggest point resembling the breakpoint with the biggest magnitude of change was detected by using bfast package in R (Verbesselt, Hyndman, Newnham, & Culvenor, 2010). Within the time periods, the month of breakpoints during 2007-2012 (5 years), and 2014-2015 (2 years) identified from different NDVI coarse resolution products were classified annually and compared one to another.
Besides comparing the annual breakpoint maps, we also compared trend maps and the annual breakpoints from different NDVI time series by reclassifying the type of trend based on slopes and p-values as shown in Table 3. From these categories, we calculated the percentage of area that was labelled with the same trend on different maps by graphing the trend distribution into a matrix as illustrated in Table 4.
From the matrix in the Table 3, we calculated the percentage of trend map by dividing the pixels that fall into the same categories in both maps divided by the total pixels as showed in the following equation 1: where A, F, K and P denote the diagonal cells in Table 3 that hold the number of pixels that possess the same trend categories in both maps.

Regression analysis with Landsat data and in situ measurements
The ability of coarse resolution NDVI products when presenting NDVI values was compared with the two data sets: 1) Landsat systems (Landsat 5 TM and 8 OLI) and 2) flux tower measurements by constructing the ordinary least square (OLS) regression. This regression analysis was conducted as an alternative due to the lack of field measurements of vegetation biophysical features during the time periods in Indonesia needed to estimate the performance of Table 2. Type 3 classification scheme of MOD12Q1 and the associated notation for each classes in this study.
coarse resolution NDVI products when estimating the biophysical properties. In addition, the comparison with the available global raster data of vegetation biophysical properties may create bias since those data were originated and constructed from the sensors of coarse NDVI products, for instance, LAI and fPAR estimated from MODIS, SPOT VGT and Proba V, and tree cover percentage from MODIS and AVHRR data (Baret et al., 2007;Defries, Hansen, Townshend, Janetos, & Loveland, 2000;DiMiceli et al., 2011;Myneni, Knyazikhin, & Park, 2015;Sánchez, Camacho, Lacaze, & Smets, 2015;Yang et al., 2006). The assessment with moderate and detailed spatial resolution data were performed with the assumption that Landsat systems and flux tower data have greater quality when used to predict biophysical properties with their detailed spatial resolution of Landsat system and flux tower data, and complete spectral bands of Landsat systems. Landsat data were employed along with other detailed resolution data as one of the input data to model fraction cover estimated from MODIS and also often used to create detailed resolution LAI map from ground sampling data to validate the coarse data estimation of vegetation biophysical properties (Chen et al., 2002;Defries et al., 2000;Fernandes, Butson, Leblanc, & Latifovic, 2003;Kang et al., 2016;Wang et al., 2004;Yang et al., 2007).
Landsat data and coarse resolution satellite imageries have a different temporal resolution or revisit time. The pixel values from Landsat system were recorded in every 16 days, while NDVI from coarse resolution products were constructed from daily, weekly, 10-day or bi-monthly data. To minimize the bias due to the difference in the revisit time, the analysis was conducted by comparing the average values from the minimum three continuous cloudmasked Landsat data with the average values of NDVI from coarse products at the same periods to generate the OLS regression line. We collected 46 scenes of time series data originated from 255 individual scenes of Landsat 5 TM during 2007-2012, and 73 scenes of time series data originated from 472 scenes of Landsat 8 OLI during 2014-2015 ( Figure 2). Both Landsat systems data obtained at Level-2 correction, which pixel values represent atsurface reflectance and were resampled to 8-km.
Besides using Landsat as comparison, we also measured the RMSE of NDVI calculated from satellite observations with measurements taken at the flux tower located in the Pasoh Tropical Rainforest, Malaysia (Figure 3). At the flux tower, NDVI measurements have been taken every 30 min on days with no rain since October 2008 along with monitoring of other environmental parameters (details on homepage: http://www.bluemoon.kais.kyoto-u.ac.jp/pasoh/ main.html) (Nakaji et al., 2014). The data enabled us to compare NDVI values collected at the flux tower with the values collected in satellite data over the period of overlapping data collection. The monthly flux tower data were analyzed by taking the average monthly values and the maximum monthly values ( Figure 3) and comparing them to the MVC data of NDVI measurement. Statistical analysis was performed by deriving the root mean square error (RMSE) and R 2 values for the comparison between observed values from satellites observation and actual values from the flux tower measurement.
In summary, the OLS regression was performed by using NDVI data as listed in Table 5. The regression analysis was performed by using 1) all pixels in the imageries, 2) vegetation pixels only, and, 3) major vegetation biomes such as a) evergreen broadleaf forests (EBF), b) savanna (SAV), and c) cropland (CRO) to assess the performance of the VI values of coarse resolution products in different vegetation cover types.

Pairwise correlation between NDVI datasets
Pairwise analysis of NDVI datasets collected by different satellites and prepared for comparison in the 7, 24 and 60 months periods in the spatial resolution categories of 1, 4 and 8 km had Pearson correlation coefficients that were positive but ranged from weak to moderate (Table 6). For 1 km spatial resolution, the strongest coefficient (0.52) was obtained for the comparison of METOP and Proba during the observation period 2014-2015, followed by the pairwise comparison of METOP and SPOT with a coefficient of 0.50. On the other hand, the weakest coefficient (0.17) was obtained for the comparison of Terra and AQUA in 2014-2015, and this value was smaller than in the first analysis in 2007-2012 when it was 0.35. PROBA-V and SPOT showed a relatively strong correlation coefficient (0.44) for datasets at the original spatial resolution. At a spatial resolution of 4 km, AVHRR NESDIS and VIIRS NESDIS showed a strong correlation with 0.67 while the comparison of PROBA-V and METOP was the second strongest with 0.54. The weakest correlation coefficient (0.27) at this spatial resolution was for the comparison of MERIS and AQUA. This combination indicates that extended measurement to expand NDVI observation is feasibly achieved by combining AVHRR NESDIS with VIIRS NESDIS. At 8 km spatial resolution, METOP and SPOT combination followed by SPOT and TERRA combination became the two strongest correlations with 0.63 and 0.58, respectively. In addition, the longest NDVI observation from different sensors can be generated by combining the AVHRR-based NDVI products from GIMMS3G data and NESDIS due to the strong correlation between those datasets with 0.42 (correlation of coefficient).
The average correlation of coefficient map, as shown in the Figure 4, was estimated by resampling  the 61 correlation of coefficient combination in the previous analyses into 8 km and compiling the average and standard deviation values from that pairwise combination. The map indicated the distribution of strong and weak values within specific categories. To investigate whether that distribution corresponds to specific land cover classes, we performed the zonal analysis using the land cover categories.

Variation of correlation coefficients across land cover types
Mapping pairwise comparison correlation coefficients by land cover type in a zonal analysis at a spatial resolution of 1, 4 and 8 km revealed that some land cover types showed consistently strong or weak correlation coefficients at different spatial resolutions as showed in the Figure 5. On the correlation coefficient maps for spatial resolution of 1 km, there was a strong agreement for grass/cereal crops, broadleaf crops and SAV with average correlation coefficients of 0.67, 0.66 and 0.63, respectively, whereas correlation coefficients showed weak agreement for EBF, evergreen needleleaf forest and non-vegetated classes with 0.32, 0.26 and 0.34, respectively. At 4 km spatial resolution, great agreement (strong average correlation coefficient) was found in the classes of broadleaf crops and SAV with 0.84 and 0.76, respectively, and a moderate agreement was found for grassland/cereal crops (0.68). The land cover types EBF, evergreen needleleaf forest and non-vegetated class showed weak correlation coefficients at the spatial resolution of 4 km and had average values of 0.42, 0.37 and 0.41, respectively, at 8 km spatial resolution. Likewise,  broadleaf crops, SAV and grass/cereal crops showed strong correlation coefficients of 0.81, 0.73 and 0.63, respectively, at 8 km spatial resolution. Based on correlation coefficient results, strong correlation values for SAV and grassland land cover types which have frequent vegetation cover change are mainly influenced by climatic factors when the distinct pattern of precipitation and temperatures appeared in the dry and wet season. Agreement between the phenology and climate parameters is strong for these land cover types due to the dependence of the growth of those land cover types to the water availability which differs at different seasons (dry and wet), while vegetation cover types such as evergreen broad/needleleaf forest are less influenced by climate variables, and the values are stable but different, thus producing weak average correlation coefficients. The weak agreement in those land cover types corroborates the research of Rasmus Fensholt and Proud (2012), which found great differences in trends in land cover types of needleleaf forest and broadleaf forest.

Trend and breakpoint analysis
Trend analysis using the Mann-Kendall test generated Kendall's slope map with corresponding p-values. Mapping the associated p-values produces a spatial distribution showing areas with greening and browning trends that are statistically significant (p-value < 0.05) and areas with no significant trend for each NDVI application for the 60-month dataset ( Figure 6). Pairwise comparison of trend results from different NDVI applications in terms of spatial agreement expressed as a percentage of the same or similar class producing maps that show agreement are shown in Table 7. The average spatial agreement between trend maps of different NDVI applications from different sensors and spatial resolutions was weak at 37.2 ± 2.5%.
In 8 km NDVI data, agreement varied from 33.13% (GIMMS vs. SPOT VGT) to 46.44% (Metop vs. SPOT VGT). The spatial agreement with GIMMS3G data on the trend map was relatively weak with only 35.2% on average, with the strongest agreement for the comparison of MODIS Aqua with GIMMS data (37.45%). In 4 km spatial resolution analysis, a strong agreement was observed between SPOT VGT and Metop (44.12%) and between Aqua and Terra (42.88%). Data assessed over the period from April 2007 to March 2012 revealed that the strongest agreement by using AVHRR NESDIS was shown with Aqua (37.04%), although other datasets showed more or less similar level of agreement. In the assessment for the second period (January 2013 to December 2014), which involved the VIIRS datasets and Proba-V, AVHRR NESDIS had the greatest agreement with VIIRS NESDIS (38.23%), while the second strong agreement was obtained from the pairwise comparison of Metop and Proba-V (40.25%).
The trend map assessed for the original spatial resolution scale of 1 km revealed the strongest agreement between trend maps with an average of 37.32% (Apr 2007 to Mar 2012) and 35.98% (Jan 2013 to Dec 2014). The strongest agreement was obtained from METOP and SPOT VGT with 40.94% of similarity percentage. In the period of January 2014 to December 2015, the trend map from Proba-V showed slightly smaller agreement with another trend map from other datasets in this assessment, the strongest agreement was observed with the Metop data (37.98%).
This analysis also revealed the possibility of producing contradictory results in trend maps estimated from different satellite observation, such as achieving significantly negative trends in one trend map but significantly positive trends in another trend map. By taking the average of the matrix of similarity tables and measuring the distribution of contradictory category results from 16, 30 and 21 combinations of trend maps at 1, 4 and 8 km spatial resolution scales, respectively, we found a small percentage of significant contradictory trends from SN to SP (around 2.5-5.39%) as shown in Table 8. The percentage of contradictory trends seems small if we look at the proportion from the total; however, this percentage of contradictory trends may account for a 10-23% portion of detected significant greening or browning (SP and SN). In conclusion, there was approximately 0.20-0.25 (20-25%) possibility that a trend detected as being significant from one NDVI observation might be classified as a significant opposite trend if another NDVI observation was used.
The similar discrepancy among the coarse resolution NDVI products also appeared in the trend map analysis. In the first analysis to detect annual breakpoints during 2007 -2012 (Table 9a), the similarity percentages were fallen between 29.12 and 40.51%. The greatest level of similarity percentage when identifying breakpoint was achieved by using SPOT VGT and MODIS Terra with Table 7. Percentage of spatial similarities in trend maps in vegetation areas estimated from different NDVI time series applications at spatial resolution of 1, 4 and 8 km.  Figure 6. Greening and browning trends over 60 months mapped using different NDVI time series applications at different spatial resolutions.
40.51% of similarity, while the smallest came from comparing STAR NESDIS with AVHRR GIMMS. In the period of 2014-2015 (Table 9b), the percentages were much bigger considering that there were only two categories of 2014 and 2015 for the detected breakpoints. The similar percentages fell within the range of 51.1-58.42%. The greatest similarity percentage in this period was produced by comparing STAR VIIRS and NESDIS data with 58.42%, while the weakest came from comparing STAR NESDIS with MODIS TERRA data.

Comparisons with Landsat and flux tower measurements
The comparison between NDVI from coarse resolution products and Landsat data revealed different level R 2 from 0.19 (VIIRS) to 0.64 (SPOT VGT and PROBA-V). Further removal of non-vegetation pixel data from the regression analysis still put SPOT VGT and PROBA-V with the data that have strong relationships with Landsat with the R 2 of 0.68 (SPOT VGT) and 0.67 (PROBA-V), respectively (Figure 7). In addition, MODIS-Terra and METOP were the third and fourth strongest with 0.54 and 0.53 when all pixels were considered, while MODIS-Terra and Aqua were the third and fourth strongest with 0.58 and 0.53 when only vegetated pixels were included in the regression analysis. The complete regression models from all coarse products can be seen in the Supplementary Materials 1 and 2 (S1). This analysis concluded that the newer products such as SPOT VGT, Proba-V, MODIS-Terra, Metop and Aqua   matched strongly with the Landsat data compared than the AVHRR-based NDVI index (i.e. GIMMS or NESDIS). The regression in specific vegetation cover types revealed that generally MODIS Terra/Aqua, SPOT VGT, Proba-V and METOP generated stronger model fit in different vegetation cover types. Generally, weakest linear fit can be found when modelling the EBF class with the maximum 34% of variance in EBF can be explained by using SPOT VGT data. In modelling SAV class, overall greater linear fits were displayed from NDVI coarse products with the range of 31% (GIMMS3g) to 76% (SPOT VGT) of the explained variance. Overall strong fit also can be found in modelling CRO class with the exception from GIMMS3g v0 data that only able to explain 5% of the variance. The regression results and the corresponding RMSE values can be discovered at Table 10.
Further analysis was conducted between NDVI data collected at the flux tower in Pasoh Forest, Malaysia with the pixel value of NDVI from satellite data in the same coordinate where the Flux Tower is located for monthly averages, and the comparison was evaluated in terms of RMSE and R 2 values (Table 11). GIMMS3G ver. 0, SPOT VGT, METOP, MODIS TERRA, MODIS AQUA and PROBA V showed comparable average RMSE with a range of 0.06 to 0.09 with SPOT VGT being the greatest and Proba V being the smallest, and the maximum RMSE for these data was 0.04 to 0.11 with the same pattern of results. On the other hand, a significant difference from flux tower measurements was observed for AVHRR NESDIS and VIIRS NESDIS, which showed similar average RMSE (0.41 and 0.16, respectively) and maximum RMSE (0.44 and 0.2, respectively), which differed from comparison results for the other NDVI applications. Weak R 2 values in this analysis indicate a weak correlation between data collected at the flux tower and satellite data used in this study. Although a weak R 2 value might be due to the density of the canopy, NDVI applications were most likely saturated in the area of the flux tower, which also results in smaller R 2 values. The density of canopy cover in this area was reflected in the comparison of NDVI values from the flux tower datasets and satellites with the RMSE of the average and maximum values being not markedly different, indicating constant and stable NDVI values in this area.
Comparing pixel and point measurements from the flux tower observations was apparently not a proper comparison since satellite values were influenced by the surrounding environment and were more likely to be exposed to atmospheric disturbances than the flux tower measurements. In addition, observations from more than one site in different areas should be used to estimate the RMSE. However, for this analysis, most of the error deviation for all NDVI applications appears to be relatively similar, with the exception in AVHRR NESDIS and VIIRS NESDIS, which somehow have much bigger RMSE.

Discussion
The comparison analysis in this study, which comprises of three sub-analyses of correlation coefficients, trend and breakpoint analysis using the Mann-Kendall test, and bfast package, and comparison of satellite data with NDVI from Landsat and flux tower, showed the variability in NDVI values and their consistencies in trends, and performance over South East Asia landscape. The correlation coefficients from the measurement of different NDVI time series products showed weak to medium strength with AVHRR GIMMS, and monthly STAR NESDIS had an overall smaller correlation with other datasets. In order to ensure seamless integration of multiple datasets including inter-sensor NDVI time series observations, especially in areas covered by broadleaf or needleleaf forest which were common to be found in this area, appropriate pre-processing algorithms to address factors that affect NDVI values are needed. It is also worth noting the possibility of performing calibrations between different NDVI time series products, particularly in this region of South East Asia, and these should be performed using values in climate-dependent vegetation types such as SAV, broadleaf crops and grassland areas.
The variability of NDVI values between different NDVI products as shown in the pairwise comparisons with weak to medium correlation coefficients also gets reflected in the trend and breakpoint map similarity analysis, which exhibits a similar weak to medium percentage of a similar trend with AVHRR GIMMS and STAR NESDIS data became the data with smaller overall similarity. The possibility of smaller paired similarity between STAR NESDIS and STAR VIIRS can be caused by the aggregation process from 7-day data into monthly data.
The comparison of coarse resolution NDVI products and Landsat data depicted that some sensors such as SPOT VGT, PROBA V and MODIS Terra/ Aqua had an overall strong relationship when linked with Landsat data with overall strong linear in representing the values of different vegetation cover types. On the contrary, weaker relationships were showed by GIMMS, NESDIS and VIIRS data. It should be noted that at the time of conducting this study, the newer version of AVHRR GIMMS ver.1 was launched which may have an improvement over the GIMMS ver.0 data that was used in this study.
The performance of coarse resolution NDVI data in the dense canopy forest as represented by the flux Table 10. Regression results from NDVI coarse resolution products versus Landsat systems. tower data showed that different NDVI time series observations have similar RMSE when compared to NDVI datasets at the flux tower in the forest area. The correlation coefficient was quite weak as NDVI tended to become saturated in the area where the comparison was undertaken. Furthermore, this assessment comparing flux tower measurements is not sufficient to make a conclusion on the overall quality of NDVI time series products since data from only a single site of flux tower was used. In addition, linking the time series values of NDVI with the actual field datasets is also essential to assess the performance of this index when used to model vegetation biophysical properties and to understand what kind of properties in vegetation being represented by this vegetation index.

Conclusion
This study emphasized the discrepancies among the available and existing NDVI products and stressed the importance of fully understanding the limitation and possible differences in outputs that might be present, especially when time series analysis is performed using a single dataset. Other researchers aiming to use the coarse data to estimate the trends from NDVI are suggested to employ multiple satellite time series observations to confirm conclusions and to identify methods for reconciling discordant findings by various NDVI applications. Furthermore, a vegetation index that possesses the smaller sensitivity to the atmospheric gasses compositions and background variations is beneficial to generate stable and comparable vegetation index products (Huete et al., 2002). Therefore, another effort can be made by contrasting time series analysis using NDVI with other proxy parameters that are less impacted by atmospheric disturbances; such as vegetation optical depth (VOD) or Enhanced Vegetation Index (EVI), which does not saturate in areas of dense vegetation. In future studies, assessing the quality of remote-sensing-based vegetation proxy by linking the values with the ground datasets is essential to further understand the performance of vegetation indices when used to model the vegetation biophysical properties.