Regression Analysis of Normalized Difference Vegetation Index (NDVI) to Compare Seasonal Patterns and 15 Year Trend of Vegetation from East to West of Nepal

Understanding the changing patterns and trend of vegetation is essential for its socio-environmental values. Normalized difference vegetation index (NDVI), a satellite based data obtained from Moderate Resolutions Imaging Spectro-radiometer (MODIS) were analysed. The data have a characteristic resolution of 250 × 250 m2 and a 16-day composite period. They were ordered separately for each sample plot from east, centre and west of Nepal, for 15 years period, 2000 to 2015. MODIS, Land Surface Temperature (LST) data were used to identify unreliable NDVI values and hence eliminated. Also, the unusually fluctuating NDVI values during the rainy season were removed. A cubic spline function (for seasonal patterns), linear regression model (for NDVI trend) and generalized estimating equations (GEE for comparison of the changing trends) were the models used. The results showed a patterned annual seasonal vegetation and significant trends during the 15 years. Seasonal growth showed a peak in rainy season and trough in the winter season, with slight temporal variation among the areas with a characteristic shift of seasonal greening (start of greening) and browning (end of greening) from east to west of Nepal. The NDVI trend was significantly rising in eastern and western suburban areas while the central urban city had a significant decline. The temporal shift of start and end of the season from east to west can be of value to agriculturalists. INTRODUCTION The global climate has been changing and so is the trend of vegetation throughout the world. On one hand, vegetation is one of the prone factors to climate change while on the other hand, the vegetation helps to curb the effects of climate change effects. Moreover, changes in vegetation constantly help to change the climate factors at local or regional level (Meng et al. 2019, Zhu et al. 2016). Both increasing (Li et al. 2011, Liu et al. 2019) or decreasing patterns of vegetation have been observed in Asia. Apart from being geo giants of Asia, China and India represent as ‘greening leaders’ (Dunne 2019). Nepal, a small landlocked country between China and India, owns agriculture and forest as a major national economic resource and essentially needs to be aware of the changing trends of vegetation in the country. The climate-changing factors, for instance, temperature and rainfall are mostly correlated with vegetation (Liu 2017, Anbazhagan & Paramasivam 2016, Kaufmann et al 2003). The changes in terrestrial vegetation can modify regional and global climate at diurnal, seasonal and long term (Bounoua et al. 2000). The Australian Government Bureau of Meteorology (BoM 2014) suggests that the ability to link events in the natural world to a cycle that predicts seasonal changes is essential in the successful development of any community. These natural barometers are not uniform across the land but instead use the reaction of plants and animals to gauge what is happening in the environment. The cropping intensity (yield, the timing of plantation and harvest etc.) response to climate can help the understanding of proper growth, development and net production, which is practically influenced by previous year’s knowledge in most of the developing countries (Cooper et al. 1997). To understand the climate and seasonal effect to vegetation and crops, its start and end are required to appropriately fix the calendar of vegetation or other crops cycles. This start or end of crop seasons would change by time and space, the knowledge of which can eventually yield better. The spatio-temporal characteristics of remote sensing data are considered to be the primary advantage in environmental studies. MODIS is a sensor, fitted aboard the Terra and Aqua satellites by the National Aeronautics and Space Administration (NASA). It monitors environmental changes due to vegetation, temperature, droughts and flood on Earth 2021 pp. 267-273 Vol. 20 p-ISSN: 0972-6268 (Print copies up to 2016) No. 1 Nature Environment and Pollution Technology An International Quarterly Scientific Journal Original Research Paper e-ISSN: 2395-3454 Open Access Journal Nat. Env. & Poll. Tech. Website: www.neptjournal.com Received: 26-11-2019 Revised: 02-01-2020 Accepted: 01-03-2020


INTRODUCTION
The global climate has been changing and so is the trend of vegetation throughout the world. On one hand, vegetation is one of the prone factors to climate change while on the other hand, the vegetation helps to curb the effects of climate change effects. Moreover, changes in vegetation constantly help to change the climate factors at local or regional level (Meng et al. 2019. Both increasing (Li et al. 2011, Liu et al. 2019 or decreasing patterns of vegetation have been observed in Asia. Apart from being geo giants of Asia, China and India represent as 'greening leaders' (Dunne 2019). Nepal, a small landlocked country between China and India, owns agriculture and forest as a major national economic resource and essentially needs to be aware of the changing trends of vegetation in the country.
The climate-changing factors, for instance, temperature and rainfall are mostly correlated with vegetation (Liu 2017, Anbazhagan & Paramasivam 2016, Kaufmann et al 2003. The changes in terrestrial vegetation can modify regional and global climate at diurnal, seasonal and long term (Bounoua et al. 2000). The Australian Government Bureau of Meteorology (BoM 2014) suggests that the ability to link events in the natural world to a cycle that predicts seasonal changes is essential in the successful development of any community. These natural barometers are not uniform across the land but instead use the reaction of plants and animals to gauge what is happening in the environment. The cropping intensity (yield, the timing of plantation and harvest etc.) response to climate can help the understanding of proper growth, development and net production, which is practically influenced by previous year's knowledge in most of the developing countries (Cooper et al. 1997). To understand the climate and seasonal effect to vegetation and crops, its start and end are required to appropriately fix the calendar of vegetation or other crops cycles. This start or end of crop seasons would change by time and space, the knowledge of which can eventually yield better. (NASA 2015). The Normalized Difference Vegetation Index (NDVI) is a data product from MODIS which is much common (NASA 2015, Eckert et al. 2014, Suepa et al. 2016 and more reliable (Yin et al. 2012) than other data types and much useful for local (Zhang et al. 2013), regional (Piao et al. 2011) or global scale (Liu 2017) study purpose. NDVI is a type of vegetation index which is dependent on the spectral behaviour of vegetation that characteristically absorbs in the red and blue wavelengths, reflects in the green wavelength, strongly reflects in the near-infrared (NIR) wavelength. Based on this principle, the data for vegetation index is calculated by the given equation (Ozyabuz et al. 2015).
Various previous studies, showing positive and negative trends of vegetation change, have commonly applied linear regression model (Kaufmann et al. 2003, Karnieli et al. 2010, Zhang et al. 2013. Moreover, a lot of variation does exist among their works regarding the methods, data types and management (Eckert et al. 2014, Piao et al. 2011, Chandola et al. 2010, Mishra & Chaudhuri 2015. Nonetheless, the analysis of the bulk data in smaller temporal and spatial sections for understanding the characteristics of vegetation (such as phenology, the start of the season (greening) and end of the season (browning)) in a particular locality in detail is entirely lacking. Also, the knowledge of intra-annual pattern together with the annual trend of vegetation in Nepal is still rare. Therefore, this study aims to assess the seasonal pattern and the trend of vegetation in Nepal from 2000 to 2015.

MATERIALS AND METHODS
This study was carried out in three purposively selected districts of Nepal, Dhankuta (East), Kathmandu (Center) valley and Surkhet (West). All the three are from the same geographical area, the hilly region, so that the other influencing factors would be least in effect, east -Dhankuta (27.15°N, 87.35°E), centre -Kathmandu (27.59°N, 85.39°E) and west -Surkhet (28.62°N, 81.88°E) of Nepal (Fig. 1). Department of Forest Resource and Survey (DFRS 2015), Department of Hydrology and Meteorology (DHM 2015) and National Population and Housing Survey 2011 (NPHS 2012) of the government of Nepal have reported that Dhankuta is an eastern suburban area of 892 km 2 with a population density of 183/ km 2 . It has an annual temperature ranging from 14.6°C (January) to 24.9°C (April), annual rainfall 1,121 mm and has an average altitude of 1,192 m. Kathmandu is an urban area located in the central region, with the area of 899 km 2 and the population density of 4,416/ km 2 . Its annual temperature ranges from 6.6°C (January) to 16.6°C (May), annual rainfall 1,667 mm and has an average altitude of 1,337 m. Surkhet is a suburban area located in the western region, with the area of 2,451 km 2 and the population density 143/ km 2 . The annual temperature of this city ranges from 15.1°C (January) to 27.0°C (April), rainfall 1,392 mm (July) and has an average altitude of 875 m. The conventional seasons in the country were classified as summer, rainy and winter which falls from March to May, June to August and December to February respectively.
The NDVI data were downloaded from MODIS's website (ORNL DAAC 2015). Data were specified for a period from 2000, the starting point of MODIS service, to 2015. The NDVI data were observed in every 16-day interval period, there were around 23 observations every year, that made a total of 345 observations for 15 years. In each of the three selected areas, NDVI observation occurred in 6561 grids (81×81). The data for all three regions were retrieved by the same process.
The raw data of each NDVI grid were divided by 10000 to obtain the values ranging from -1 to +1. The negative values up to 0 correspond to water. The values from 0 to 0.1 mean soil, rocks or concrete, the snow land and barren land. Low positive values (0.2 to 0.4) represent shrubs and grassland. The value close to 1 (0.6 and above) means the forest (Weier & Herring 2000). A higher NDVI value represented the denser vegetation in the area.
The seasonal pattern of vegetation in an area does not vary much among its nearby grids unless the geography and climate features of the area is hugely changed. That is why, only the central position grid, out of 6561, is chosen for analysing the seasonal pattern of vegetation. However, even though the trend of NDVI (that means the change of vegetation by time) can vary in the shorter distance too, all the grids are not considered for analysis to avoid the special correlation of the data. Hence, 49 locations were systematically selected to represent the whole study site. Then, from each location, four surrounding grids were chosen for The data for all three regions were retrieved by the same process.
The raw data of each NDVI grid were divided by 10000 to obtain the values ra from -1 to +1. The negative values up to 0 correspond to water. The values from 0.1 mean soil, rocks or concrete, the snow land and barren land. Low positive v (0.2 to 0.4) represent shrubs and grassland. The value close to 1 (0.6 and above) m the forest (Weier & Herring 2000). A higher NDVI value represented the d vegetation in the area.
The seasonal pattern of vegetation in an area does not vary much amo nearby grids unless the geography and climate features of the area is hugely cha That is why, only the central position grid, out of 6561, is chosen for analysin seasonal pattern of vegetation. However, even though the trend of NDVI (that m the change of vegetation by time) can vary in the shorter distance too, all the gri not considered for analysis to avoid the special correlation of the data. Henc analysis and that made 196 grids altogether. The process was repeated for all three sites one by one.

Data Management
First of all, NDVI was plotted with respect to observation day as shown in Fig. 2 (a, b and c). Here, n is the total number of observations and every dot on each vertical grid line represents one observation value on the same recording period. Therefore, every vertical line on the x-axis displays 15 dots of NDVI values corresponding to a particular day in each of the 15 years. This plot starts from Julian day 1 and ends on day 365. The data structure needs to be further managed before going into analysis due to big data gaps (big NDVI gap between consecutive 8 day's series is practically impossible in case of vegetation growth) and unreliable values (all data of a particular satellite for a particular site would be unusual when temperature data is deemed faulty). Those data are believed to reduce the value outcome of the result and hence eliminated before analysis. From Fig. 2, cross marked points are the sparse NDVI data and blue dots are indicated to be unreliable values, all are eliminated before analysis.

The Statistical Methods
First of all, the NDVI data for the central grid were analysed to identify the annual seasonal pattern by using a cubic spline function. The form of the function is, 6 ts one observation value on the same recording period. Therefore, every on the x-axis displays 15 dots of NDVI values corresponding to a particular f the 15 years. This plot starts from Julian day 1 and ends on day 365. The e needs to be further managed before going into analysis due to big data VI gap between consecutive 8 day's series is practically impossible in case n growth) and unreliable values (all data of a particular satellite for a e would be unusual when temperature data is deemed faulty). Those data to reduce the value outcome of the result and hence eliminated before m Fig. 2, cross marked points are the sparse NDVI data and blue dots are be unreliable values, all are eliminated before analysis.

cal Methods
the NDVI data for the central grid were analysed to identify the annual tern by using a cubic spline function. The form of the function is, the spline function, α, b and c k are the parameters in the model. k is the not. t denotes time in Julian day, that is, specified day of the year. t1 < t2 < ecified knots and (ttk) + means that (ttk) is positive for (t > tk) and zero

…(1)
Where, S t is the spline function, α, b and c k are the parameters in the model. k is the location of knot. t denotes time in Julian day, that is, specified day of the year. t 1 < t 2 < ... < t p are specified knots and (t -t k ) + means that (t -t k ) is positive for (t > t k ) and zero otherwise.
The data were, then, seasonally adjusted to stabilize the mean of the data. The seasonal adjustment was computed by, then, seasonally adjusted to stabilize the mean of the data. The as computed by, lly adjusted NDVI, x is NDVI after removal of unreliable and he fitted values from cubic spline function and x is the mean of sted NDVI data were further used for analysing the time series d.
NDVI trends of selected 196 grids (out of n=6561) were ear models to seasonally adjusted data separately for each of 196 of the linear model is, Here, y is the seasonally adjusted NDVI, x is NDVI after removal of unreliable and doubtful values, S t is the fitted values from cubic spline function and x is the mean of x. The seasonally adjusted NDVI data were further used for analysing the time series trend in 15 years period.
Thereafter, the NDVI trends of selected 196 grids (out of n=6561) were identified by fitting linear models to seasonally adjusted data separately for each of 196 sample grids. The form of the linear model is, e, then, seasonally adjusted to stabilize the mean of the data. The was computed by, nally adjusted NDVI, x is NDVI after removal of unreliable and the fitted values from cubic spline function and x is the mean of justed NDVI data were further used for analysing the time series iod.
he NDVI trends of selected 196 grids (out of n=6561) were inear models to seasonally adjusted data separately for each of 196 rm of the linear model is, ally adjusted NDVI,  is intercept and  is the coefficient of t vation day and  is the error term.

…(3)
Here, y is the seasonally adjusted NDVI, α is intercept and b is the coefficient of t that represents observation day and e is the error term.
However, the linear trend in each grid had a slight spatial correlation with the nearby grids. To tackle this problem, the Generalised Estimation Equations (GEE) were finally applied in this study. GEE is an extension of a linear model that is specially designed for correlated data (Liang & Zeger 1986). Furthermore, from the fitted models explaining the NDVI changes for all 15 years, it was difficult to distinguish the details of trends within this period. Therefore, the data were divided into three periods of 5 years each (2000-2004, 2005-2009 and 2010-2015) before fitting the model. At last, 95% confidence intervals were calculated for each sub-period and finally plotted to show the NDVI changes in three subsections of 15 years' period in each of the three study regions, separately. The GEE (Dormann et al. 2007) equation is illustrated as follows, Here, Y is a vector of seasonally adjusted NDVI, E(Y) or  is an expected v of Y, 1  g is an inverse link function, T is a matrix of observation days and  is a v of regression coefficients.
All data analysis and graphical displays were carried out using R Statis Programming version 3.2.1(R Core Team 2015).

Seasonal Pattern from Cubic Spline Function
Based on trial and error, eight knots (blue plus sign at the bottom of Fig. 2 a,  58, (before) to 86, 50 and 82 after the data management process.
The NDVI in the east region extended from 4-8.5 (Fig. 2 a), meaning, it consis shrubs as well as the densely wooded forest. The central region showed a na variation, from 0.65 to 0.90, indicated that vegetation is dominantly the big trees Here, Y is a vector of seasonally adjusted NDVI, E(Y) or m is an expected value of Y, g -1 is an inverse link function, T is a matrix of observation days and g is a vector of regression coefficients.
All data analysis and graphical displays were carried out using R Statistical Programming version 3.2.1(R Core Team 2015).

Seasonal Pattern from Cubic Spline Function
Based on trial and error, eight knots (blue plus sign at the bottom of Fig. 2 a, b and c) were fixed at the position of 15,40,70,120,150,200,230 and 350 days of the year, to fit cubic spline function. The plots of fitted values from cubic spline function for east, centre and west areas showed that r squared increased from, respectively, 11, 15 and 58, (before) to 86, 50 and 82 after the data management process.
The NDVI in the east region extended from 4-8.5 (Fig. 2 a), meaning, it consists of shrubs as well as the densely wooded forest. The central region showed a narrow variation, from 0.65 to 0.90, indicated that vegetation is dominantly the big trees. The western region was identified with a range from 0.4 to 0.8 meaning the area equally had both shrub and trees.
Besides, a lot of variations in vegetation response towards the seasonal effects were observed. The three seasons in Nepal are broadly classified as, summer: March, day 61 to May, day 150, rainy: June, day 151 to August, day 240 and winter: December, day 360 to February, day 60. When the data were plotted from day 1 to day 365, the NDVI in each selected region [ Fig. 2 (a-c)] showed the growth pattern with a peak in the rainy season and gradually declined to trough in winter. The beginning of the season (greening) in east, centre and west took place in days 81, 97, 129 while the end of the season (browning) occurred on 257, 273, 273 days respectively. Here, both greening and browning are earlier in the east than the other areas. The result of seasonal patterns of vegetation in Kathmandu valley and Surkhet are similar with some difference in eastern Dhankuta. In former two areas, NDVI peaks in the rainy season and declines to trough in the winter season [ Fig. 2(b) and 2(c)] where the greening is remarkably late by 60 to 70 days as compared to the eastern region. These graphs clear that start of the season (greening of vegetation) had a trend to move from east to the westward area and same was the case of browning.

Time Series Trends from Linear Model
The results from the linear regression models of 196 selected grids, were plotted separately for 15 years. In Fig. 3, only two grid plots were selected to represent increasing and decreasing trends from each of the three locations. The green horizontal line across the graph explains the NDVI trend for 15 years period the annual seasonal fluctuation cycle of NDVI, derived from the spline function was added back to the plot and shown in the red line. In Fig. 3 to Fig. 5, black dots are data points plotted year wise. The cross dots are the data eliminated at the time of cubic spline fit. The increasing or decreasing trend (Inc/dec) per decade and respective p-values from linear regression show how much vegetation has been changed from 2000 to 2015. Here, n represents the number of observations in each plot.
In the eastern region [ Fig. 3 (a) and (b)], the total increasing trend was in 85.20% area and out of that 69.5% was found significant. The decreasing trend was in 14% with 47% being significant in the area. Hence, we can say that majority of vegetation was tending to rise in this area. The central region [ Fig. 4 (a) and (b)] has a mixed form of result. The total rise was seen in 56.50% area and only 40.70% was significant. The decreasing trend was in 40.10% out of which 60.30% were significant. The result is almost balanced between increasing and decreasing trends. Therefore, the net gain or trend seemed to be balanced here. Finally, in the western gradually declined to trough in winter. While in the centre, Kathmandu NDVI ranged from 7 to 9 [ Fig. 2 (b)] indicating a dominant wooded forest. In the west, Surkhet, the

GEE and Confidence Interval Plot
All the 196 grids' data for 15 years, were staked into one, separately for every area. The data were then divided into three time sections of five years each, they are, 2000-2004, 2005-2009 and 2010-2015. Finally, to take account of the autocorrelation of the data, GEE were fitted to each of five years of data in each site. A 95% CI plots were drawn from the mean of predicted values from the model to show vegetation trends in three time sections (Fig. 6). The red horizontal line is the level of no change of NDVI and the area above it indicates the increasing pattern and below it means the vegetation is decreasing state. The colourful vertical lines are CI plots of coefficients from the GEE model.
The results from CI plots are much interesting. The changes in vegetation at different time frames were compared among the groups using CI plots and the overall temporal change were compared among the three study areas.

East
The eastern area showed significant increasing trends of vegetation with the highest rate (0.022 per decade, p-value < 0.05) in the first period from 2000-2004 and the lowest rate (0.019 per decade, p-value < 0.05) in the latest phase, 2010-2015 (Fig. 6, green CI lines). Hence, the eastern area has an overall increasing trend but the rate of change is gradually declining. The recent trend showed that the mean increase in NDVI is 0.019 per decade.

Centre
The red CI lines (Fig. 6) all lying at or below 0 change line indicate that the net change is negative here. In the central area, the negative values of NDVI showed that vegetation is declining with a faster rate in the latest period (-0.006 per decade, p-value < 0.05).  part [ Fig. 5 (a) and (b)], the increase was seen in 81.63% of the area and 75% was significant. Also, out of the total decreasing trend, 80.76% was significant.

model.
The results from CI plots are much interesting. The changes in vegetation a time frames were compared among the groups using CI plots and the overal change were compared among the three study areas.

East
The eastern area showed significant increasing trends of vegetation with th rate (0.022 per decade, p-value<0.05) in the first period from 2000-2004 and rate (0.019 per decade, p-value<0.05) in the latest phase, 2010-2015 (Fig. 6 lines). Hence, the eastern area has an overall increasing trend but the rate of gradually declining. The recent trend showed that the mean increase in NDV per decade.

Centre
The red CI lines (Fig. 6) all lying at or below 0 change line indicate that the n

West
In the western area, the vegetation is significantly increasing in overall study period similar to in eastern area (Fig. 6, pink CI lines). However, this area showed a recent faster rate of increment as compared to the earlier periods. The recent five years showed that the mean increase rate is 0.014 per decade (p-value < 0.05).

DISCUSSION
The seasonal pattern of vegetation, in this study, showed that it grew the highest in the rainy season and gradually declined to the lowest in the winter season, probably due to defoliation and physiological dormancy. Minimum and maximum NDVI help us in identifying the intra-annual vegetation change (Liu 2017). This result is consistent with other studies in tropical or temperate vegetation (Suepa et al. 2016, Zang et al. 2016, Chen et al. 2014, Evrendilek & Gulbeyaz 2008, Yin et al. 2016. Also, it is scientifically proved that a relatively cooler air temperature reduces the plant physiology including its growth (Fitter & Hay 2002). The seasonal fluctuation of vegetation may occur due to several climatic factors like the temperature, rainfall, humidity and the seasons which this study does not explain. The east area is quicker to be affected by the seasonal change (both start of greening and browning). It might be due to its proximity of the area to the Bay of Bengal, where Nepal draws the monsoon from. That is why it receives earlier and more precipitation during rainy season resulting in a shift of seasonal changes from east to westward region. The phonological shift was studied (Suepa et al. 2016) in Southeast Asia but on the annual (temporal) basis. While in our study, the seasonal shift of vegetation pattern is intra-annual and spatial type. This result will particularly be helpful in the agricultural sector to understand the annual climate response of vegetation or crops, especially to understand the existing phenological characteristics of the area.
Actually, the vegetation trends have nothing to do with the monsoon and seasons. It reflects the direction of change over a period of time. This study identified that the annual vegetation trend is considerably rising in east and west areas of Nepal while the central Kathmandu showed a declining trend. The Global NDVI trend studied during 1982-2012 showed an increasing trend in many parts of the world including India and Southeast China (Liu et al. 2015). Nepal lies in between these two blocks and possibly have a similar increasing trend overall. But Kathmandu being a growing urban area and densely populated city, might have been affected by several local environmental factors to show a decline of vegetation and hence, detecting the local varia-tion of vegetation change within a country is unique in this study.
The method of data management had contributed to obtaining much improved result from cubic spline to obtain the seasonal patterns and also the method can be applicable for analysing other noisy data. Here, the spline function along with linear regression and GEE models have been successfully used to analyse the seasonal pattern and the time series trends of NDVI.
However, this study is limited to local sites and study does not answer the overall factors contributing to these particular trends and patterns. This can be further worked out.

CONCLUSION
To sum up, this study shows the changes in NDVI in two ways. First, the seasonal pattern at a grid level, that clears the local level intra-annual change and second, the trend of vegetation in 15 years period. Moreover, a comparison of rates of changes within specified time segments can be made clear by CI plots. Eastern part of Nepal shows a wider range of vegetation with faster annual greening and browning than the rest of the two areas. The vegetation in Nepal has been increasing from 2000 to 2015 in sub-urban areas (Dhankuta and Surkhet) while it decreases in the urban region (Kathmandu). The study has academic implication in assessing the seasonal changes and time-series trends. It can also provide important information for urban planners, agriculturalists and the environmental experts. The study underpins the focus on analysing the potential reason behind this kind of NDVI changes in Nepal.