Exploring Long Term Spatial Vegetation Trends in Taiwan from AVHRR NDVI 3 g Dataset Using RDA and HCA Analyses

Due to 4000 m elevation variation with temperature differences equivalent to 50 degrees of latitudinal gradient, exploring Taiwan’s spatial vegetation trends is valuable in terms of diverse ecosystems and climatic types covering a relatively small island with an area of 36,000 km2. This study analyzed Taiwan’s spatial vegetation trends with controlling environmental variables through redundancy (RDA) and hierarchical cluster (HCA) analyses over three decades (1982–2012) of monthly normalized difference vegetation index (NDVI) derived from the Advanced Very High Resolution Radiometer (AVHRR) NDVI3g data for 19 selected weather stations over the island. Results showed two spatially distinct vegetation response groups. Group 1 comprises weather stations which remained relatively natural showing a slight increasing NDVI tendency accompanied with rising temperature, whereas Group 2 comprises stations with high level of human development showing a slight decreasing NDVI tendency associated with increasing temperature-induced moisture stress. Statistically significant controlling variables include climatic factors (temperature and precipitation), orographic factors (mean slope and aspects), and anthropogenic factor (population density). Given the potential trajectories for future warming, variable precipitation, and population pressure, challenges, such as land-cover and water-induced vegetation stress, need to be considered simultaneously for establishing adequate adaptation strategies to combat climate change challenges in Taiwan.


Introduction
The terrestrial ecosystem is subject to the unprecedented speed of global climate change and must be carefully monitored because of its central role in maintaining a healthy and sustainable planet.Because of the close relationship between vegetation phenology and the terrestrial ecosystem [1][2][3], vegetation is commonly used to determine the ecological response to environmental changes [4][5][6][7], such as climate change [8][9][10][11], land cover and land use change [12][13][14][15][16][17][18][19], drought [20,21], and changes in net primary productivity [22].Satellite remote sensing vegetation indices offer opportunities to monitor vegetation phenology and environmental changes in a repeatable manner [22].The normalized difference vegetation index (NDVI) is one of the most widely accepted vegetation indices based on a mechanism that chlorophyll in vegetation absorbs more strongly in the red light spectrum (RED) whereas the cell structures of leaves reflect and scatter light in the near-infrared spectrum (NIR).The NDVI value is calculated by (NIR ´RED)/(NIR + RED) and expressed on a scale from ´1 to +1, in which positive values correspond to the presence of vegetation and negative values correspond to the absence of vegetation [23].
Time series of NDVI data have proven to be an effective tool in assessing the phase of climate variability, and has been applied to worldwide geographical locations [3,[24][25][26][27][28][29].NDVI data derived from the Advanced Very High Resolution Radiometer (AVHRR) onboard the National Oceanic and Atmospheric Administration's (NOAA) satellite series are widely used and have been considered as the best dataset available for long-term analysis of vegetation dynamics at reasonable spatial and temporal scales [30].A number of studies investigating changes in AVHRR NDVI data found that regional climate variability differentially regulated vegetation phenology across different ecosystems at various spatial scales [31][32][33][34][35].For instance, Yu et al. (2003) [36] revealed that the response of vegetation photosynthetic activities to climate variations depends on season while other studies showed that climate warming is asymmetric among different seasons and continents [37,38].
In order to maintain a livable planet, the United Nations Development Programme Global Environment Facility developed the Adaptation Policy Framework (APF), which aimed to reduce negative impacts and enhance any beneficial consequences of a changing climate [39].Additionally, an Ecosystem-based Adaptation (EbA) approach has emerged in the international climate change arena, which has proved to provide flexible, cost effective and broadly applicable alternatives for reducing the impact of climate change in many countries, such the United States [40].For the major concept in APF and EbA, an understanding of the target natural terrestrial ecosystem pattern is required for a successful implement.
Located in East Asia, Taiwan is a unique island with relatively small land area of 36,000 km 2 and a variable topography ranging from sea level to nearly 4000 m.Due to the 4000 m elevation variation with temperature differences of 24 ˝C that could occur equivalent to 50 degrees of latitudinal gradient, Taiwan has high levels of species richness and biodiversity across various ecosystems and climatic types [41,42].As such, the vegetation trends in Taiwan could represent the typical plant communities across various ecosystems equivalent to 50 degrees of latitudinal gradient.Additionally, Taiwan regularly encounters many natural challenges, such as earthquakes, typhoons, mudslides, and flash floods.Based on numerous evidences around the world, those natural challenges would be highly possible exacerbated by global climate change [43].Besides natural challenges, Taiwan's landscape has been transformed extensively by human developments and the pressure from an increasing population density continues to grow.The compacted population coupling with exacerbated natural challenges, the ecosystems and the safety of people's life in Taiwan are vulnerable with great risk.Therefore, the need to establish adequate adaptation strategies based on a comprehensive understanding of Taiwan's natural ecosystems has never been greater [44][45][46].
Taiwan's government adapted the concepts of APF, and a series of comprehensive strategies was established since 2012.With considering and understanding the natural terrestrial ecosystems of Taiwan, the effectiveness of the current adaptation strategies would be improved in terms of their ability for handling emerging large-scale climate change challenges.NDVI data have been applied in Taiwan on many short-term land cover change detection cases [18], forest coverage evaluation [47], and crop yield evaluation [3].However, a long-term assessment of NDVI patterns remains unexplored.With the longest available NDVI series data, the present study intends to assess the long-term vegetation patterns in Taiwan by analyzing NDVI as a vegetation surrogate with associated environmental variables.The main objectives of this study are to resolve: (1) the spatial characteristics of the NDVI in Taiwan, (2) the environmental variables that control the spatial characteristics of NDVI variation; and (3) the current vulnerabilities and future climate risks in Taiwan.Consequently, proper adaptation strategies can now be established to improve Taiwan's adaptive capacity and avoid compounding the costs and damages resulting from climate change challenges [43].

Study Area
Taiwan straddles the Tropic of Cancer (22 ˝-25 ˝N, 120 ˝-122 ˝E, Figure 1) and is an island located on the northwestern edge of the Pacific Ocean.The annual average temperature is approximately 22-25 ˝C [48] and the annual precipitation is approximately 2500 mm.Taiwan is mainly dominated by the north-south trending Central Mountain Range (CMR) with the highest point at Yu Mountain (3952 m) and an average height of 3000 m [49].The climate of Taiwan is primarily influenced by lofty mountains and alternating monsoons [50,51].In addition, the low pressure system, which causes summer typhoons and the winter monsoon [49], has a large influence on the island's climate during the course of the year.Topographic factors such as the land-sea distribution further amplify the complexity of Taiwan's precipitation mechanism [52].For this study, 19 weather stations were selected on the basis of two criteria: complete records for NDVI and environmental variables from 1982-2012 and spatially even distribution across the island to ensure a general representation.Of the 19 weather stations, there are five, six, three, and five stations located in the northern, central, southern, and eastern counties, respectively.According to a detailed national vegetation community survey [53], nearly 60% of the land area in Taiwan is covered by forest and the majority vegetation formation for each station is provided in Figure 1.

Study Area
Taiwan straddles the Tropic of Cancer (22°-25°N, 120°-122°E, Figure 1) and is an island located on the northwestern edge of the Pacific Ocean.The annual average temperature is approximately 22-25 °C [48] and the annual precipitation is approximately 2500 mm.Taiwan is mainly dominated by the north-south trending Central Mountain Range (CMR) with the highest point at Yu Mountain (3952 m) and an average height of 3000 m [49].The climate of Taiwan is primarily influenced by lofty mountains and alternating monsoons [50,51].In addition, the low pressure system, which causes summer typhoons and the winter monsoon [49], has a large influence on the island's climate during the course of the year.Topographic factors such as the land-sea distribution further amplify the complexity of Taiwan's precipitation mechanism [52].For this study, 19 weather stations were selected on the basis of two criteria: complete records for NDVI and environmental variables from 1982-2012 and spatially even distribution across the island to ensure a general representation.Of the 19 weather stations, there are five, six, three, and five stations located in the northern, central, southern, and eastern counties, respectively.According to a detailed national vegetation community survey [53], nearly 60% of the land area in Taiwan is covered by forest and the majority vegetation formation for each station is provided in Figure 1.

NDVI
This study utilizes the third generation of the AVHRR NDVI3g dataset developed by the Global Inventory Modeling and Mapping Studies group from the National Aeronautics and Space Administration Goddard Space Flight Center.The long temporal coverage, spanning from July 1981 to December 2012, makes this dataset valuable for assessing long-term vegetation activities at a reasonable spatial resolution (approximately 8 km).The AVHRR NDVI3g dataset is calibrated using the atmospheric Rayleigh-scattering-over-oceans approach by Vermote and Kaufman (1995) [54].

NDVI
This study utilizes the third generation of the AVHRR NDVI3g dataset developed by the Global Inventory Modeling and Mapping Studies group from the National Aeronautics and Space Administration Goddard Space Flight Center.The long temporal coverage, spanning from July 1981 to December 2012, makes this dataset valuable for assessing long-term vegetation activities at a reasonable spatial resolution (approximately 8 km).The AVHRR NDVI3g dataset is calibrated using the atmospheric Rayleigh-scattering-over-oceans approach by Vermote and Kaufman (1995) [54].Volcanic stratospheric aerosol periods (1982-1984 and 1991-1994) were established and subjected to atmospheric corrections.Empirical decomposition reconstruction methods were applied to correct the satellite orbital drift effect [55].Further improved cloud masking and sea-viewing wide field-of-view sensor NDVI data were used to calibrate the AVHRR NDVI3g data using Bayesian methods [56] to reduce the factors unrelated to changes in vegetation greenness.
To minimize the corruption of vegetation signals [57], the AVHRR NDVI3g data obtained two maximum NDVI values.The first value was the maximum NDVI value obtained over the first 15 days, and the second was the maximum NDVI value obtained for the rest of the month.Thus, each month contained two NDVI images.For this study, the image with the higher NDVI value for each month was chosen to represent the monthly NDVI.These NDVI values were extracted for 19 selected cells where weather stations present for further analysis (Figure 2).Volcanic stratospheric aerosol periods (1982-1984 and 1991-1994) were established and subjected to atmospheric corrections.Empirical decomposition reconstruction methods were applied to correct the satellite orbital drift effect [55].Further improved cloud masking and sea-viewing wide field-ofview sensor NDVI data were used to calibrate the AVHRR NDVI3g data using Bayesian methods [56] to reduce the factors unrelated to changes in vegetation greenness.
To minimize the corruption of vegetation signals [57], the AVHRR NDVI3g data obtained two maximum NDVI values.The first value was the maximum NDVI value obtained over the first 15 days, and the second was the maximum NDVI value obtained for the rest of the month.Thus, each month contained two NDVI images.For this study, the image with the higher NDVI value for each month was chosen to represent the monthly NDVI.These NDVI values were extracted for 19 selected cells where weather stations present for further analysis (Figure 2).

Environmental Variables
Environmental variables are characterized as natural and anthropogenic variables.Natural variables can be categorized into climatic and orographic variables.Climatic variables include temperature, precipitation, cloud amount, and sunshine duration, whereas orographic variables include mean slope, aspects, and elevation.Population density is the only anthropogenic variable considered in the present study.
Monthly observational temperature and precipitation data between 1982 and 2012 were derived from the Taiwan Climate Change Projection and Information Platform (TCCIP).The spatial resolution of 5 km × 5 km for these data sets was generated by the inverse distance-weighted interpolation and weighted-average methods suggested by Watson (1992) [58].Altitude, latitude/longitude, distance, and azimuth have been involved into the interpolation process by the TCCIP.Cloud amounts and sunshine duration were collected from first-order weather stations operated by Taiwan's Central Weather Bureau.The mean slope, aspect, and elevation data (spatial resolution of 40 km 2 ) were calculated and extracted from digital elevation model data from the Department of Land Administration, Ministry of the Interior.Due to the actual population density information for all NDVI cells are unavailable, the population density data were adopted on a township basis depending on where each station is located by deriving the latest census records.
For this study, 19 weather stations were selected on the basis of two criteria: they had to have complete records from 1982-2012, and they needed to be relatively distributed across the island to ensure a general representation of Taiwan.All environmental variables for each station were extracted and their basic information can be found in Figure 3.

Environmental Variables
Environmental variables are characterized as natural and anthropogenic variables.Natural variables can be categorized into climatic and orographic variables.Climatic variables include temperature, precipitation, cloud amount, and sunshine duration, whereas orographic variables include mean slope, aspects, and elevation.Population density is the only anthropogenic variable considered in the present study.
Monthly observational temperature and precipitation data between 1982 and 2012 were derived from the Taiwan Climate Change Projection and Information Platform (TCCIP).The spatial resolution of 5 km ˆ5 km for these data sets was generated by the inverse distance-weighted interpolation and weighted-average methods suggested by Watson (1992) [58].Altitude, latitude/longitude, distance, and azimuth have been involved into the interpolation process by the TCCIP.Cloud amounts and sunshine duration were collected from first-order weather stations operated by Taiwan's Central Weather Bureau.The mean slope, aspect, and elevation data (spatial resolution of 40 km 2 ) were calculated and extracted from digital elevation model data from the Department of Land Administration, Ministry of the Interior.Due to the actual population density information for all NDVI cells are unavailable, the population density data were adopted on a township basis depending on where each station is located by deriving the latest census records.
For this study, 19 weather stations were selected on the basis of two criteria: they had to have complete records from 1982-2012, and they needed to be relatively distributed across the island to ensure a general representation of Taiwan.All environmental variables for each station were extracted and their basic information can be found in Figure 3.

Hierarchical Cluster Analysis (HCA)
To discover the spatial characteristics of NDVI in Taiwan, the extracted NDVI values from the AVHRR NDVI3g dataset for each of the 19 weather stations were assessed through clustering methods.Rather than segmentation approaches [59,60], Hierarchical cluster analysis (HCA) is suitable for clustering targets with an unknown number of groups a prior.An agglomerative HCA was applied to group the mean monthly values to reveal natural groupings [61,62].In the present study, the Ward linkage method was used with the squared Euclidean distance as the similarity measurement [63].Detailed calculation processes can be found in a classical paper conducted by Punj

Hierarchical Cluster Analysis (HCA)
To discover the spatial characteristics of NDVI in Taiwan, the extracted NDVI values from the AVHRR NDVI3g dataset for each of the 19 weather stations were assessed through clustering methods.Rather than segmentation approaches [59,60], Hierarchical cluster analysis (HCA) is suitable for clustering targets with an unknown number of groups a prior.An agglomerative HCA was applied to group the mean monthly values to reveal natural groupings [61,62].In the present study, the Ward linkage method was used with the squared Euclidean distance as the similarity measurement [63].Detailed calculation processes can be found in a classical paper conducted by Punj and Stewart (1983) [64].Initially, the agglomerative HCA starts with individual weather stations.The most similar stations are first grouped, and these initial groups are merged according to their similarity measurements.Groups (or clusters) are formed by merging stations one stage at a time, until all stations are merged in one big cluster.The number of stages is bounded by the number of stations in the initial partition, which is one less than the total number of stations [65].The total number of stations is 19 in this study, so the number of stages is 18.At each stage, the corresponding HCA coefficient indicates the distance between the two clusters.A large difference between the coefficients of two consecutive stages indicates that the clusters being merged possess increasing heterogeneity so to terminate the merging process before the clusters become too different [66].As the merging process being terminated, the final number of clusters can be obtained by deducting the value of the terminated stage from the total number of stations.

Redundancy Analysis (RDA)
To identify the environmental variables controlling the spatial characteristics of NDVI variation, Redundancy analysis (RDA) was utilized to explore the relationship between NDVI and environmental variables.The RDA was first proposed by Rao (1964) [67] and later rediscovered by Van Den Wollenberg (1977) [68].Redundancy is often synonymous with explained variance [69].The RDA is a direct extension of multiple linear regression by allowing the regression of multiple response variables (RVs) on multiple explanatory variables (EVs) [70,71].In addition, an RDA may also be seen as an extension of a principal component analysis because the canonical ordination vectors are linear combinations of the RVs.In general, RDA is a method based on canonical multivariate analyses by assuming a linear response between variables [72][73][74] to extract and summarize the variation in a set of RVs that can be explained by a set of EVs [71].The assumed linear response between variables can be expressed by the Eigen analysis equation as follows, and be decomposed using a standard eigenvalue-eigenvector algorithm.
where S YX is the covariance matrix among RVs and EVs, S ´1 XX represents the inverse covariance matrix among standardized EVs, I denotes a unit matrix, λ k represents the eigenvalues of the corresponding axis k, and u k denotes normalized canonical eigenvectors (details can be found in Legendre and Legendre 2012 [70]).In this research, the suitability of an RDA was confirmed by a detrended correspondence analysis with a measured length of turnover units smaller than three.An RDA can assess how much of the variance in the mean monthly NDVI values (RVs) can be explained by the environmental variables (EVs).Arcsine transformation was applied for the percentage data, including the cloud amount and all aspects data.Moreover, the score scaling type of this study was set to focus on RV correlations, and RV scores were divided by standard deviations.The explanatory power of each environmental variable was evaluated by their simple term and conditional term effects.The relative relationships between the weather stations, NDVI, and environmental variables were demonstrated with biplot and triplot diagrams.
Basic interpretation rules and examples for RDA biplots and triplots are listed as follows and can be found on page 186-195 in Smilauer and Leps 2014 [75].The arrows represent RVs and EVs, and the symbols represent cases.All interpretations of ordination diagrams focus on the concept of relative relationships, such as the relative distances of symbols, the relative directions of arrows, or the relative ordering of projection points.In this study, the NDVI values from each month stood for the RVs, the environmental variables were the EVs, and the 19 weather stations were the cases.

‚
The arrows point in the direction of maximum increase in the value of the variable across the diagram, and its length is proportional to this maximum rate of change.

‚
An approximate ordering of the value of one RV or EV across cases is obtained by projecting the case points perpendicular to the RV or EV arrow.

‚
A case point projecting onto the origin of the coordinate system (perpendicular to an RV or EV arrow) is predicted to have an average value of the corresponding variable.The cases projecting further from zero in the direction of the arrow are predicted to have above-average values, and the case points projecting in the opposite direction are predicted to have below-average values.

‚
The relative directions of arrows approximate the linear correlation coefficients among the variables.If an EV arrow points in a similar direction to an RV arrow, the values of that RV are predicted to be positively correlated with the EV values.

‚
The cosine of the angles between any two arrows indicates their respective relationship.If the arrows meet nearly at a right angle, these two variables are predicted to have a low (near to zero) correlation.

Results and Discussion
To evaluate the relationship between the selected environmental variables, Pearson productmoment correlation analysis was adapted [76].Table 1 represents the correlation matrix of all environmental variables except the aspect variables.The positive correlation (p < 0.01) between rain and clouds suggested that more clouds would result in more rainfall.In addition, the positively significant correlation (p < 0.05) between elevation and the mean slope reflected the mountainous terrain of Taiwan.Moreover, the evidence of orographic lifting effects on precipitation observed in Taiwan with a positively significant relationship (p < 0.05) was found between rain and mean slopes.
The negative and significant correlation (p < 0.01) between elevation and temperature directly reflected the lapse rate concept, which is defined as the decrease of atmospheric temperature caused by an increase in altitude.The degree of the mean slope also increased as altitude increased in Taiwan, resulting in the significantly negative correlation between the mean slope and temperature (p < 0.01).The sunshine duration was found to be negatively correlated with rain and cloud amounts (p < 0.01), when more clouds produced more rain, the length of sunshine duration was reduced.In addition, the negative correlations between rain and temperature (p < 0.05) can be explained by physical science: dry conditions favor more sunshine and less evaporative cooling, whereas wet summers are cool [77].

HCA
HCA coefficients indicating the distance between two clusters joined at each stage can objectively determine the final number of clusters.A large increase of the HCA coefficient values indicates joining clusters being less homogeneous so to terminate the merging process.From the result of HCA in Figure 4, the largest increase in the coefficients is observed between stages 17 and 18, therefore, stage 17 would be the most appropriate stage for terminating the merging process.By deducing the value 17 from the total number of stations, 19, the final cluster number is two.As a result, NDVI patterns in Taiwan can be categorized into two groups spatially (Figure 5).Group 1 comprises weather stations located in the northern, central, and eastern part of Taiwan, whereas Group 2 comprises stations located in the western, southern, and one station from the northeastern Taiwan (Yilan).In terms of land cover types (Figure 5), Group 1 comprises 7 stations with forest types, 3 stations with urban types, and 2 stations with cultivated land types.Group 2 constitutes 4 stations with urban types and 2 stations with cultivated land types.The different composition of land cover types within each group provides possible explanation to the revealed NDVI values.Group 1 stations have higher average NDVI values than Group 2 stations (Table 2, 0.70 and 0.46, respectively).The stations located in the CMR area in Group 1 (Yushan, Alishan, and Sun Moon Lake) reveal high average NDVI values   In terms of land cover types (Figure 5), Group 1 comprises 7 stations with forest types, 3 stations with urban types, and 2 stations with cultivated land types.Group 2 constitutes 4 stations with urban types and 2 stations with cultivated land types.The different composition of land cover types within each group provides possible explanation to the revealed NDVI values.Group 1 stations have higher average NDVI values than Group 2 stations (Table 2, 0.70 and 0.46, respectively).The stations located in the CMR area in Group 1 (Yushan, Alishan, and Sun Moon Lake) reveal high average NDVI values In terms of land cover types (Figure 5), Group 1 comprises 7 stations with forest types, 3 stations with urban types, and 2 stations with cultivated land types.Group 2 constitutes 4 stations with urban types and 2 stations with cultivated land types.The different composition of land cover types within each group provides possible explanation to the revealed NDVI values.Group 1 stations have higher average NDVI values than Group 2 stations (Table 2, 0.70 and 0.46, respectively).The stations located in the CMR area in Group 1 (Yushan, Alishan, and Sun Moon Lake) reveal high average NDVI values ranging from 0.73 to 0.78 (Figure 3).Land development in the CMR area is restricted by law in Taiwan so that vegetation in this area maintains its natural condition.In addition, the relatively high altitude, steep slopes, and low temperatures make these stations inaccessible for large populations, and thus result in low human intervention so to yield high NDVI.The stations located in the eastern part of Taiwan in Group 1 (Dawu, Taitung, Hualien, and Su'ao) also reveal high average NDVI values ranging from 0.67 to 0.84 (Figure 3).Due to blocking effects of the CMR and limited lowland areas, human development in the eastern part of Taiwan is slower than the western part so to result in less human disturbance, less infrastructure development, and less land use changes, which would exert less impact on vegetation and may yield high NDVI.Furthermore, this study adopted the latest census on a township/district basis depending on where each station is located due to unavailable population density for all NDVI cells.Group 1 constitutes one densely populated station in Zhongzheng district (21,478 people/km 2 ) in Taipei city and several sparsely populated sites in mountainous area, such as Yushan (12 people/km 2 ) and Alishan (13 people/km 2 ).Taipei weather station locates in a restricted district within 300 m nearby the President Hall of Taiwan so that the NDVI cell containing Taipei weather station has a great reason to maintain a relatively high NDVI than other urban areas within Zhongzheng district.Since the population density of Zhongzheng district cannot represent the population density of Taipei weather station neighborhood, the averaged population density of Group 1 was obtained as 1796.98 people/km 2 by excluding Taipei station.According to Taiwan's history, human immigrated into lowlands and forests were gradually cleared for human settlements [78].Most of Group 2 stations locate in the western and southern lowlands of Taiwan (average elevation 32 m) where are easily accessible for human settlements since the 17th century and nowadays are occupied by approximately 95% of the island's population [79].Furthermore, the average population density of Group 2 stations (9243 people/km 2 ) is much higher than Group 1 (1796.98people/km 2 ).It is reasonable to assume that human population with accompanying human activities, such as land clearing, agricultural practices, industrialization, and urbanization processes, result in a massive landscape modification and a negative effect on vegetation so to yield low NDVI values.Overall, the result of HCA implies that the level of anthropogenic activity plays an important role in controlling the variation of NDVI in Taiwan.
From the annual average NDVI time series (Figure 6a), Group 1 shows a slight increasing tendency whereas Group 2 expresses a slight decreasing tendency.In addition, both groups illustrate increasing temperature trends over the period 1982-2012 (Figure 6b).The total precipitation for both groups represents fluctuated patterns with no discernable trends (Figure 6c).Although these tendencies were not very strong in terms of low R 2 values, the trajectories are extremely critical while facing future climate change.In a recent study of tropical mountain regions, Krishnaswamy et al. (2014) [25] found mild greening trends followed by browning trends around the mid-1990s.They also observed strong associations between vegetation with increasing temperature and temperature-induced moisture stress from their study period of 1982-2006.In the present study, Group 1 comprising stations with high altitude possesses a slight increasing tendency from 1982 to early 1990s (Figure 6d) and a noticeable decreasing tendency from 1991 to 1994 and 1994 to 2006 (Figure 6e,f) that agrees well with Krishnaswamy et al. (2014) [25].Therefore, it is reasonable to assume that the observed rising temperature plays an important role in vegetation trends for Group 1. Furthermore, based on a century of precipitation data, Chen et al. (2009) [80] found the incidence of meteorological drought has increased in the central and southern Taiwan.Because water is naturally associated with vegetation state and cover, the increasing incidence of meteorological droughts would negatively correlate with NDVI [81].Therefore, the decreasing NDVI tendency found in Group 2 may associate with the increasing incidence of meteorological droughts.In short, both groups illustrate significant influence of climate variability.
The scatterplots of NDVI, temperature, and precipitation for each group are shown in Figure 7a-d.None of the scatterplots show significant correlations.However, both temperature and precipitation illustrate positive associations with NDVI for Group 1 and negative associations for Group 2. This may imply that these two groups respond to temperature and precipitation differently.
From the annual average NDVI time series (Figure 6a), Group 1 shows a slight increasing tendency whereas Group 2 expresses a slight decreasing tendency.In addition, both groups illustrate increasing temperature trends over the period 1982-2012 (Figure 6b).The total precipitation for both groups represents fluctuated patterns with no discernable trends (Figure 6c).Although these tendencies were not very strong in terms of low R 2 values, the trajectories are extremely critical while facing future climate change.In a recent study of tropical mountain regions, Krishnaswamy et al. (2014) [25] found mild greening trends followed by browning trends around the mid-1990s.They also observed strong associations between vegetation with increasing temperature and temperatureinduced moisture stress from their study period of 1982-2006.In the present study, Group 1 comprising stations with high altitude possesses a slight increasing tendency from 1982 to early 1990s (Figure 6d) and a noticeable decreasing tendency from 1991 to 1994 and 1994 to 2006 (Figure 6e,f) that agrees well with Krishnaswamy et al. (2014) [25].Therefore, it is reasonable to assume that the observed rising temperature plays an important role in vegetation trends for Group 1. Furthermore, based on a century of precipitation data, Chen et al. (2009) [80] found the incidence of meteorological drought has increased in the central and southern Taiwan.Because water is naturally associated with vegetation state and cover, the increasing incidence of meteorological droughts would negatively correlate with NDVI [81].Therefore, the decreasing NDVI tendency found in Group 2 may associate with the increasing incidence of meteorological droughts.In short, both groups illustrate significant influence of climate variability.
The scatterplots of NDVI, temperature, and precipitation for each group are shown in Figure 7a-d.None of the scatterplots show significant correlations.However, both temperature and precipitation illustrate positive associations with NDVI for Group 1 and negative associations for Group 2. This may imply that these two groups respond to temperature and precipitation differently.(e) (f)

RDA
In addition to climatic factors represented by temperature and precipitation, other environmental variables, including cloud amounts, sunshine duration, elevation, mean slope, aspects, and population density, were explored to explain the NDVI variations in the last three decades.The variance in NDVI values explained by the environmental variables was fairly high, as (e) (f)

RDA
In addition to climatic factors represented by temperature and precipitation, other environmental variables, including cloud amounts, sunshine duration, elevation, mean slope, aspects, and population density, were explored to explain the NDVI variations in the last three decades.The variance in NDVI values explained by the environmental variables was fairly high, as

RDA
In addition to climatic factors represented by temperature and precipitation, other environmental variables, including cloud amounts, sunshine duration, elevation, mean slope, aspects, and population density, were explored to explain the NDVI variations in the last three decades.The variance in NDVI values explained by the environmental variables was fairly high, as indicated by the eigenvalues obtained from the first two canonical axes in the RDA ordinations (Figure 8).In a direct ordination with the RDA, the first axis (RDA Axis 1) and the second axis (RDA Axis 2) explained 90.06% of the variation in NDVI values.For each month's NDVI variations (Table 3), Axis 1 explained over 80% of the NDVI variation.With all axes together, more than 96% of the NDVI variations could be explained.
indicated by the eigenvalues obtained from the first two canonical axes the RDA ordinations (Figure 8).In a direct ordination with the RDA, the first axis (RDA Axis 1) and the second axis (RDA Axis 2) explained 90.06% of the variation in NDVI values.For each month's NDVI variations (Table 3), Axis 1 explained over 80% of the NDVI variation.With all axes together, more than 96% of the NDVI variations could be explained.In Figure 8a, the length of May's NDVI was the shortest, indicating that May's NDVI had the smallest standard deviation.Furthermore, the angles between each month's NDVIs were less than 90 degrees, which indicated that each month's NDVIs were positively correlated.An approximate NDVI value for each weather station could be inferred by projecting the station points perpendicularly onto the NDVI arrow.By comparing the perpendicular projections of weather stations 1 and 3 onto the line overlaying the arrow for December's NDVI, we deduced that the NDVI value in Yushan (0.86) was higher than that of Sun Moon Lake (0.81) in December.Furthermore, weather stations 4 and 5 (Anbu and Zhusihu) located closer to the coordinate origin indicates that their December's NDVI values were closer to the average of NDVI value in December of all 19 stations.The angle between In Figure 8a, the length of May's NDVI was the shortest, indicating that May's NDVI had the smallest standard deviation.Furthermore, the angles between each month's NDVIs were less than 90 degrees, which indicated that each month's NDVIs were positively correlated.An approximate NDVI value for each weather station could be inferred by projecting the station points perpendicularly onto the NDVI arrow.By comparing the perpendicular projections of weather stations 1 and 3 onto the line overlaying the arrow for December's NDVI, we deduced that the NDVI value in Yushan (0.86) was higher than that of Sun Moon Lake (0.81) in December.Furthermore, weather stations 4 and 5 (Anbu and Zhusihu) located closer to the coordinate origin indicates that their December's NDVI values were closer to the average of NDVI value in December of all 19 stations.The angle between rain and sunshine duration was larger than 90 degrees and indicated a negative correlation between these two variables (Figure 8b).Likewise, a positive correlation between the mean slope and elevation could be discerned from the in-between small angle (Figure 8b). Figure 8c is a triplot presenting the overall relationship among weather stations, NDVIs, and environmental variables.The correlation between each month's NDVI and environmental variables can be depicted from the cosine of the angle between them.For instance, population density and sunshine duration showed a negative correlation with NDVI (the cosine of the angle being larger than 90 degrees), whereas rainfall, cloud amount, the mean slope, and elevation displayed positive relationships with NDVI (the cosine of the angles being smaller than 90 degrees).The effects of each environmental variable on weather station clusters can be inferred from Figure 8d.The NDVI values are generally higher for Group 1 than for Group 2 stations.The values for mean slope, elevation, rain, and cloud amounts are higher in Group 1 stations.For Group 2 stations, the population density, sunshine duration, and temperature are higher.Focusing on the characteristics of aspect, Group 1 stations have higher percentages in north-and northeast-facing slopes, whereas Group 2 stations show higher percentages in west-facing, southwest-facing, and flat slopes.
Table 4 summarizes the explanatory power of all variables to the NDVI variations, which were evaluated by their simple and conditional term effects.The simple term effect is the amount of variability in the response data that could be explained by a constrained ordination model using that variable as the only EV.Hence, the conditional term effect of the variable tested is dependent on the variables already being selected and the size and significance is influenced by previously selected variables.The explanatory power of each variable needs to be interpreted with caution and discussed as follows.
The mean slope had the highest explanatory power (34.4%, p < 0.01) in both simple and conditional term effects, which may be explained by the inverse relationship between mean slopes and human disturbances (Table 4).As the lowland areas of Taiwan have been developed and urbanized, most of the pristine forests have been eliminated, whereas the steeper slope areas may have maintained their original habitats.Thus, the NDVI increases with the slope.Moreover, the population and flat aspect represent a certain level of human disturbance.The explanatory power of the flat aspect and the population density to the NDVI was 30.2% (p < 0.01) and 17.7% (p < 0.10), respectively, thus, the level of human disturbance played a key role in determining the NDVI (see Table 4).In addition, southwest-, west-, and southeast-facing slopes were significant contributors in terms of the explanatory power of southwest-and west with 25.1% and 21.7% for their simple term effect (p < 0.05), respectively, and 10.1% of southeast for its conditional term effect (p < 0.05).The effect from these aspects can be explained by the received amount of solar radiation and climate phenomena, such as monsoon systems with the effects on local microclimate, microenvironment, and corresponding local phenology [82,83].In the present study, temperature had an 18% explanatory power (p < 0.10) in its simple term effect (Table 4).Under the global warming trend [84], Hsu and Chen (2002) [85] estimated that the mean temperatures near Taiwan will likely increase between 0.9 and 2.7 ˝C relative to the 1961-1990 average.Additionally, Guan et al. (2009) [86] predicted the monthly mean temperatures for the mountain regions in the central and southern of Taiwan could be more vulnerable to future winter temperature increases than the northern region because of the chilling requirements of plant species.In this study, Group 1 stations located in the CMR area were expected to experience a more pronounced change because the vegetation status of those stations was relatively undisturbed by human activities and the response of vegetation would reveal relatively pure influences from climate variability.Biologically, some plant species need a minimum period of cold weather to ensure that vernalization occurs, a strong warming of winter temperatures could slow the fulfillment of chilling requirements and thus postpone spring phenology [87].Because the mountain systems in Taiwan harbor high levels of biodiversity and are often rich in endemic species, the changes in temperature are expected to represent a significant threat to the plant diversity of Taiwan [41].Group 1 stations in Northern Taiwan would need further investigation because the evidence shows that Northern Taiwan exhibited an increase in annual rainfall [88] and a decrease in meteorological droughts [80].This increasing trend of precipitation can be associated with other weather phenomena such as monsoon systems [89], and could be exaggerated by the topographic effect of Taiwan's lofty terrain [90].Group 2 stations, located in the western, central, and southern parts of Taiwan, may reflect part of the vegetation responses of common paddy fields.Yu et al. (2002) [91] revealed that evapotranspiration will be greatly influenced by temperature and relative humidity and would pose additional stress on agricultural water demand and associated primary productivity in Taiwan.However, further information on the detailed plant species in each group is needed to distinguish highly temperature-responsive vegetation, this is out of the scope of this research.
Besides temperature, precipitation had a 16.3% explanatory power (p < 0.10) in its simple term effect (Table 4).Because several studies delineated a drier future of the earth [92,93], many scholars declared that Taiwan is experiencing a more variable precipitation pattern.Shiu et al. (2009) [94], using 45 years (1961-2005) of hourly meteorological data in Taiwan, found a reduction in light precipitation and an increase in heavy precipitation.Another study conducted by Tu and Chou (2013) [95] also found a pattern of decreased lighter rain corresponding with an increase in heavier rain.Moreover, Chu et al. (2014) [96] observed an upward trend in precipitation intensity during the typhoon season in Taiwan contributed by the increased precipitation induced by typhoon and monsoon systems.Chu et al. (2014) [96] also noticed longer drought durations in Southern Taiwan, in agreement with Yu et al. (2002) [91].In the present study, stations in Group 2 received relatively lower amounts of precipitation, whereas stations in Group 1 received more rain.With an increasing variability in precipitation, the influence on the NDVI is expected to be varied depending on direct effects such as typhoons, extreme drought, and extreme rainfall events, and indirect effects such as rainfall-induced landslides.For instance, the severe drought in 2002-2004 directly caused an island-wide water shortage that resulted in tremendous impacts on agriculture and primary productivity in Central and Southern Taiwan (Group 2), where major agricultural activities take place [78,97].As such, establishing adaptation strategy based on the characteristics of these two groups with consideration of climate variability effects is crucial to effectively tackle climate change challenges.

Conclusions
In this study, we derived the NDVI value and relevant environmental variables for 19 weather stations in Taiwan for the period of 1982-2012, to address three questions.
What is the spatial characteristic of the NDVI value in Taiwan?As shown in the results of HCA, the spatial characteristics of the NDVI in Taiwan can be delineated into two groups.Spatially, Group 1 comprises weather stations located in the northern, central, and eastern part of Taiwan, whereas Group 2 comprises stations located in the western, southern, and one station from the northeastern Taiwan (Yilan).Group 1 stations have higher NDVI values than Group 2. Additionally, stations in Group 1 have higher values in elevation and mean slope, and lower values in temperature.Stations in Group 2 were characterized by relatively high values in temperature, sunshine duration, and population density, and relatively low values in cloud amounts, rain, elevation, and mean slope.Generally, Group 1 stands out as stations remaining their natural status with relatively less human disturbances, whereas Group 2 appears high level of human development.These two groups clearly represent two distinct vegetation characteristics of Taiwan, which suggests that the level of anthropogenic activity plays an important role in controlling the variation of NDVI in Taiwan.
What environmental variables control the spatial characteristics of NDVI values?Over 96% of the NDVI variations could be efficiently explained by the selected environmental variables based on the RDA results.Statistically significant controlling variables include climatic factors-temperature and precipitation, orographic factors-mean slope and aspects, and anthropogenic factor-population density.Temperature and precipitation would be obviously direct link to the nature of vegetation growth.The aspects of southwest-, west-, and southeast-facing slopes could be related to effects of the amount of solar radiation and climate phenomena such as monsoon systems on local microclimate, microenvironment, and corresponding local phenology.Population density, mean slope, and flat slope may associate with human disturbances and routine agricultural activities.
What are the current vulnerabilities and future climate risks in Taiwan?Group 1 shows a slight increasing NDVI tendency whereas Group 2 expresses a slight decreasing NDVI tendency.Both groups illustrate increasing temperature trends over the period 1982-2012.The NDVI tendencies and warming trends correspond well to relevant studies [25,80], suggesting significant climate variability influence can be found in both groups.With a projected rising temperature and variable precipitation [48], the estimated larger amplitudes in temperature and precipitation are expected to threaten the ecosystems and water resources in Taiwan.Weather stations in the CMR region of Group 1 may experience a more pronounced alteration because their relatively undisturbed natural environments would reveal purer responses to natural climate change.Group 1 stations in Northern Taiwan are expected to receive more rain in the future, thus, the interactions among precipitation, topographic effects, and other weather phenomena such as monsoon systems must be considered simultaneously.Group 2 stations corresponding to the major agricultural areas in Central and Southern Taiwan are expected to face challenges of water shortage due to rising temperature and variable precipitation.Because of water shortage are naturally associated with vegetation state and cover, declined water availability would impact healthy vegetation and yield low NDVI with associated primary productivity.
In summary, two groups were identified in terms of their different long term NDVI patterns, spatial characteristics, significant controlling environmental variables, and possible vulnerabilities with future climate risks.Overall, the present study contributed a valuable detailed assessment on the vegetation trends of Taiwan.On the basis of this assessment, the current vulnerability and future climate risks for Taiwan can now be estimated, evaluated, and considered in policy making and spatial planning processes by national, regional, and local agencies.Therefore, proper, effective, and efficient adaptation strategies can be implemented to tackle the emerging large-scale threats that climate change poses to people's lives and livelihoods in Taiwan.

Figure 1 .
Figure 1.Study area map of Taiwan with the 19 selected weather stations and corresponding vegetation formation.

Figure 1 .
Figure 1.Study area map of Taiwan with the 19 selected weather stations and corresponding vegetation formation.

Figure 3 .
Figure 3. Temporal variations of variables for each weather station: (a) NDVI values; (b) sunshine duration; (c) temperature; (d) cloud amount and (e) precipitation.Variables for each weather station: (f) elevation; (g) mean slope and (h) population density.

Figure 3 .
Figure 3. Temporal variations of variables for each weather station: (a) NDVI values; (b) sunshine duration; (c) temperature; (d) cloud amount and (e) precipitation.Variables for each weather station: (f) elevation; (g) mean slope and (h) population density.

Figure 4 .
Figure 4. HCA coefficient and cluster number variations at different stages.The gray box shows the observed largest gap between coefficients and star symbol indicates the stage for determining cluster number.

Figure 5 .
Figure 5. HCA clusters spatial distribution with land cover types.

Figure 4 .
Figure 4. HCA coefficient and cluster number variations at different stages.The gray box shows the observed largest gap between coefficients and star symbol indicates the stage for determining cluster number.

Figure 4 .
Figure 4. HCA coefficient and cluster number variations at different stages.The gray box shows the observed largest gap between coefficients and star symbol indicates the stage for determining cluster number.

Figure 5 .
Figure 5. HCA clusters spatial distribution with land cover types.

Figure 5 .
Figure 5. HCA clusters spatial distribution with land cover types.

Figure 7 .
Figure 7. NDVI, temperature, and precipitation scatterplots: (a) NDVI with temperature for Group 1; (b) NDVI with temperature for Group 2; (c) NDVI with precipitation for Group 1; and (d) NDVI with precipitation for Group 2.

Figure 7 .
Figure 7. NDVI, temperature, and precipitation scatterplots: (a) NDVI with temperature for Group 1; (b) NDVI with temperature for Group 2; (c) NDVI with precipitation for Group 1; and (d) NDVI with precipitation for Group 2.

Figure 7 .
Figure 7. NDVI, temperature, and precipitation scatterplots: (a) NDVI with temperature for Group 1; (b) NDVI with temperature for Group 2; (c) NDVI with precipitation for Group 1; and (d) NDVI with precipitation for Group 2.

Figure 8 .
Figure 8. RDA results: (a) RDA biplot showing NDVI and weather stations; (b) RDA biplot showing environmental variables and weather stations; (c) RDA triplot showing NDVI, environmental variables and weather stations; and (d) HCA clusters incorporate with RDA triplot (red oval represents Group 1 stations and green oval denotes Group 2 stations).

Figure 8 .
Figure 8. RDA results: (a) RDA biplot showing NDVI and weather stations; (b) RDA biplot showing environmental variables and weather stations; (c) RDA triplot showing NDVI, environmental variables and weather stations; and (d) HCA clusters incorporate with RDA triplot (red oval represents Group 1 stations and green oval denotes Group 2 stations).

Table 1 .
Correlation matrix of environmental variables.

Table 2 .
HCA clusters with corresponding average value and standard deviation (provided in parenthesis) for each variable.
* Taipei weather station was excluded from the calculation of population density for Group 1.

Table 3 .
Cumulative fraction of variation in individual monthly NDVI explained by the first, the first two, and all RDA

Table 4 .
Explanatory power of each environmental variable to the NDVI variations.