Time-Series analysis of MODIS NDVI data along with ancillary data for Land use / Land cover mapping of Uttarakhand

Land use and land cover plays an important role in biogeochemical cycles, global climate and seasonal changes. Mapping land use and land cover at various spatial and temporal scales is thus required. Reliable and up to date land use/land cover data is of prime importance for Uttarakhand, which houses twelve national parks and wildlife sanctuaries and also has a vast potential in tourism sector. The research is aimed at mapping the land use/land cover for Uttarakhand state of India using Moderate Resolution Imaging Spectroradiometer (MODIS) data for the year 2010. The study also incorporated smoothening of time-series plots using filtering techniques, which helped in identifying phenological characteristics of various land cover types. Multi temporal Normalized Difference Vegetation Index (NDVI) data for the year 2010 was used for mapping the Land use/land cover at 250m coarse resolution. A total of 23 images covering a single year were layer stacked and 150 clusters were generated using unsupervised classification (ISODATA) on the yearly composite. To identify different types of land cover classes, the temporal pattern (or) phenological information observed from the MODIS (MOD13Q1) NDVI, elevation data from Shuttle Radar Topography Mission (SRTM), MODIS water mask (MOD44W), Nighttime Lights Time Series data from Defense Meteorological Satellite Program (DMSP) and Indian Remote Sensing (IRS) Advanced Wide Field Sensor (AWiFS) data were used. Final map product is generated by adopting hybrid classification approach, which resulted in detailed and accurate land use and land cover map.


INTRODUCTION
Land Use and Land Cover (LULC) play a vital role in understanding the complex earth's ecosystem processes, biogeochemical cycles (Achard et al., 2004;Liu et al., 2008;Sellers et al., 1997), climate change studies (Liu et al., 2004) and deforestation studies (Hansen et al., 2008;Skole and Tucker, 1993;Steininger et al., 2002).Land use and land cover information is extensively used by various governmental and non-governmental organizations and thus placing a great demand on detailed, reliable and up to date LULC maps (Perera and Tsuchiya, 2009).International Geosphere Biosphere Program (IGBP 1990), Intergovernmental Panel on Climate Change (IPCC 1990), National Research Council (NRC 2005), and Global Terrestrial Observation systems (GTOS 2009) etc., have identified land use and land cover as an important driver for climate change.Land use/land cover also impacts the socioeconomic conditions of the people (Hietel et al., 2005;Lo and Yang, 2002;Ribeiro and Lovett, 2009;Swetnam et al., 2011) and hence detailed scientific study is required in this field for better environmental planning, management and conservation.
Before the era of satellite remote sensing, mapping of land resources started with manual interpretation of aerial photographs (Colwell, 1960).These techniques were tedious and time consuming.Global land cover data sets were first prepared using various published maps, atlases, national sources and ground level surveys (Matthews, 1983;Olson et al., 1983).The global land cover maps generated by this method suffered lack of consistency in classification schema used, differences in spatial resolutions and variables in measurement techniques (Running et al., 1994;Townshend et al., 1991).
Land use land cover mapping of larger areas using satellite remote sensing started with National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) data.(Nemani and Running, 1997).
Various Global land cover data sets were prepared using multi temporal AVHRR NDVI data sets (De Fries et al., 1998;IGBP DISCover, Loveland et al., 2000; UMD GLCC, Hansen et al., 2000).Major disadvantages of using AVHRR data is its poor geometric accuracy and limited radiometric calibration (Cihlar et al., 1997;Meyer et al., 1995).
Land use land cover mapping is also possible using multi resolution data obtained from various satellite platforms like Landsat (Bakr et al., 2010;Huang et al., 2009;Olthof and Fraser, 2007;Steininger et al., 2002), Indian Remote Sensing IRS-1C satellite based Wide Field Sensor (WiFS) data with 189 meter resolution (Joshi et al., 2001(Joshi et al., , 2006) ) and Indian Remote Sensing IRS-P6 satellite based Advanced Wide Field Sensor (AWiFS) data with 56 meter resolution (Kandrika and Roy, 2008).However, these data sets suffer limitations in terms of spatial coverage, temporal resolution and availability of cloud free data (Asner, 2001;Hansen et al., 2008;Song et al., 2001).
Present study has the advantage of (i) Using temporal plots of various land cover classes to understand the status of the vegetation, (ii) Smoothening the time-series curves to eliminate the erroneous interpretations, (iii) Use of ancillary data to Sandeep Kumar Patakamuri a* , Shefali Agrawal b , and M. Krishnaveni a a Centre for Water Resources, Anna University, Guindy, Chennai, India-600 025 b Indian Institute of Remote Sensing, 4-Kalidas Road, Dehradun, India-248001 support classification decision making and (iv) Use of land use land cover maps derived from AWiFS data alongside with high resolution Google earth imagery as reference data for classification accuracy assessment.The approach presented in this study is expected to have better accuracy and is adoptable for land use land cover mapping at different geographical locations.Statistics, 2010).Topography and climate plays a major role in ecology, socio economic status and culture of the people in this region (SoE, 2004).Agriculture is the major sector in this area and many industries are forest dependent.There are over fifteen important rivers and over a dozen glaciers in the state.Uttarakhand has over two hundred medium to large sized hydro-projects and these projects highly influence the environment and living conditions of the people (Rana et al., 2007;WSMD, 2009).

STUDY AREA AND DATA DESCRIPTION
The state is housing six national parks, six wild life sanctuaries and two conservation reserves.The state is one of the major contributors of the tourism sector in India, attracting several domestic and international tourists.The state also has a rich and ancient cultural heritage and is also called as Dev-Bhoomi which means land of Gods.Uttarakhand state is sensitive to natural calamities like flash floods; landslides and earth quakes and requires well planned environmental management system.In September 2010 Uttarakhand has faced floods due to heavy rain fall resulting in loss of lives and property (Sharma, 2012).This study can help planners and decision makers to better manage and protect the land resources of Uttarakhand.

METHODS
The methods adopted in this study can be broadly divided into five stages: (i) Data preparation, (ii) Image classification, (iii) Generating time-series plots and smoothening, (iv) Grouping the classes based on similarity and (v) Accuracy assessment.The flow chart depicted in Figure .2 represents the systematic approach adopted in the present study.

Data Preparation
A total of 23 tiles for H24V05 and 23 tiles of H24V06 covering the study area were downloaded for the year 2010 starting from January to December.These images were then mosaicked together to form 23 images, and later layer stacked to form annual composite image.

Image classification
In the present study ISODATA clustering technique was used to exploit its ability to further split and merge the clusters (Jensen, 1996;Richards and Jia, 1999).MODIS data was grouped to 150 clusters with 25 iterations and a convergence threshold of 0.95.The classification resulted in 105 useful clusters with non-zero values.

Generating time-series plots and smoothening
Once the classification process was finished, a set of four random samples were collected from each of the 105 clusters and multi-temporal NDVI plots were generated for each of the 420 sample locations.
There are several factors such as atmospheric conditions, sensor viewing geometry, cloud cover, aerosol concentrations, offnadir viewing and low sun zenith angles that causes noise in NDVI time series (Gutman, 1991;Holben and Fraser, 1984;Kobayashi and Dye, 2005;Li and Strahler, 1992).Savitzky-Golay filtering technique was used in this study to filter the noise in NDVI response.It is a simplified least squares fit convolution for smoothing and computing derivatives of a spectrum (Savitzky and Golay, 1964).

Grouping the classes based on similarity
The clusters with similar temporal characteristics were grouped together based on the theory that similar land cover classes exhibit similar temporal patterns.Time-series patterns of sample points were extracted from each cluster, elevation data, image interpretation of AWiFS false color composites and Google earth visual interpretation helped in proper identification of land use/land cover classes.Table .1 below shows the criteria followed in attributing suitable land use/land cover class to each cluster.

Accuracy assessment
Accuracy assessment is a key stage in thematic mapping using remotely sensed data.Classification accuracy assessment indicates the qualitative aspects of the land cover map.It determines whether the generated map serves the purpose of the users or not.Accuracy assessment is a complex process involving several problems of consideration (Foody, 2008(Foody, , 2002)).
In this study, accuracy assessment is carried out using error matrix.It accounts for commission errors, as well as omission errors and kappa statistics which reflects the difference between actual agreement and the agreement expected by chance (Lillesand and Kiefer, 2008).
Reference data for accuracy assessment is collected with the help of National land use and land cover maps derived by using multi-temporal AWiFS for the year 2009-2010 and Google Earth (GE) image interpretation.Stratified random sampling is used with a minimum of 25 samples per class, and a total of 300 sample points were generated.Grids were generated with 250 x 250 m size for each of the sample point to match the pixel size of the MODIS NDVI data and exported to Google Earth for visualization.

Plotting NDVI time-series and smoothening
Different land cover classes portray dissimilar temporal characteristics at various stages of the year.This helps in identification and classification of land use/land cover using multi-temporal remotely sensed satellite imagery.A matlab routine was coded to generate the plots and to smoothen their responses.Time series plots were generated for a set of four randomly selected points from each of the 105 clusters.Savitzky-Golay filtering was applied on the plots to obtain smoothened, noise free values to ease the interpretation process as this method is found to be suitable for smoothening the time series without much disturbing the seasonality parameters (Chen et al., 2004;Jonsson and Eklundh, 2002).Time-series plots of various land use/land cover classes and their corresponding smoothened plots were shown in the  Reference grids were created to match the spatial resolution of MODIS pixels.Many researchers have used google earth imagery for classification and validation of land use/land cover maps (Clark et al., 2010;Colditz et al., 2012).When a mixed class is observed inside the grid, the majority land cover type (>70%) present in the grid is allocated for the sample.Snow/Ice class was observed at higher altitudes and fresh snow appeared brighter than permanent snow.Barren lands support minimum to no vegetation cover and these samples were easily differentiated from other classes based on texture information.Shrub lands/degraded forest is a land cover class consisting of trees and shrubs.Crop lands were identified by rectilinear shapes and plow lines and slope agriculture was identified with the help of elevation information.Identification of grass lands was a difficult task experienced in the study.Grass lands were mixed with crop lands and shrub land classes, and differentiating these classes was found to be a challenging task.Built-up lands contained man made buildings and road networks.Evergreen forests were differentiated from deciduous frosts with the help of elevation information.Deciduous forests are located at the lower altitude zones compared to evergreen forests as deciduous forest could not adapt to the cold and dry weather at high altitudes.No shedding of leaves was observed in evergreen forests during autumn season.

Final Classification and map compilation
All the clusters belonging to a particular land use/land cover class were recoded at the final stage of classification Urban mask derived from DMSP-OLS Night time lights data and water mask derived from MOD44W product were appended to the final land use/land cover map of Uttarakhand for the year 2010.Map generated in the present study is presented on a reduced scale in Figure .4.
Snow/ice cover was evident in the higher altitudes and contributed about 12.39% of the study area.Ever green forests are evident throughout the state and constituted about 24.08% of the area.Deciduous forests covered about 13.07% of the area.Single crop lands, double/triple crop lands and fallow lands together are labeled as crop lands accounted for 20.22% of the entire area.Shrub/degraded forests covered 12.5% of the area.Water bodies contributed 0.43% of the area and built-up lands occupied about 0.16% of the area.Grass lands and barren lands accounted for 9.56% and 7.55% of Uttarakhand area respectively.

Evaluation of mapping results
Results obtained in the present study were evaluated in both qualitative and quantitative aspects.Reference map classes were aggregated to match the classification schema of the present study for quantitative validation.Area statistics for all the classes were calculated for the present map and compared with the area extent of reference map obtained from Bhuvan portal (Aggarwal et al., 2012) as shown in Qualitative aspects of the thematic mapping effort were validated using error matrix (Lillesand and Kiefer, 2008) generated with the help of three hundred sample points obtained by stratified random sampling method.NRSC-ISRO land use/land cover map for the year 2009-2010 was used as reference image for the purpose of validation.Rectangular reference grids were generated for all three hundred points with dimensions of 250m x 250m for the purpose of visualization on Google Earth.Spatial distortions occurred in the google earth reference will be negligible compared to the pixel size of the imagery used in the study (Clark et al., 2010).Table .3 represent the error matrix, producers accuracy, users accuracy and kappa statistics of the accuracy assessment process.
Error matrix obtained in the present study clearly highlights the complex nature of accuracy assessment.Present study incorporated careful examination of sample locations on multiple source like published reference maps, google earth as well as AWiFS imagery for assuring the quality of classification.The result showed an overall accuracy of 78.67% with over all kappa statistics of 0.7600.

CONCLUSION
During the land use and land cover mapping of Uttarakhand, barren lands, dry river beds, fallow lands and portions of the built-up land exhibited similar temporal pattern at some of the locations.In order to better delineate Built-up lands and Water bodies, the study made use of urban mask derived from DMSP-OLS nighttime lights time series data and water mask derived from MOD44W product.Identification of grass lands and shrub lands/mixed forest classes is observed to be a difficult task as the temporal patterns were difficult to interpret and differentiate in these classes.It is thereby suggested to consider for more number of sample locations at these clusters.Elevation information obtained from SRTM DEM data helped in differentiating various vegetation classes.Prior knowledge about the study area and suitable reference datasets also plays a major role in improving the quality of classification.It is of prime importance to have a clear definition and classification criteria for each land use/land cover type used in a particular mapping exercise.
Time series plots generated in the study may contain disturbances due to sensor geometry, atmospheric errors and sensor calibrations.Savitzky-Golay filtering was adopted in the present study to smoothen the time-series curves and there by avoid erroneous interpretations.The study incorporated the visual interpretation of google earth imagery in classification and accuracy assessment phases.With recent advancement in web based GIS it would be a great opportunity to deploy crowd sourcing in which experts from the field, research collaborators and volunteers can collectively work at various stages of land cover classification and validation.
Land use/land cover map generated in the present study was evaluated in both qualitative and quantitative aspects.Over all accuracy of the present map is found to be 78.67% with over all kappa coefficient of 0.7600.
study area is covered by two MODIS tiles numbered H24V05 and H24V06."Data required for the present study was obtained through the online Data Pool at the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota (https://lpdaac.usgs.gov/get_data)"with the help of United States Geological Survey (USGS) Earth Explorer (EE) tool (http://earthexplorer.usgs.gov/).

Figure 4 .
Figure 4. Land use/land cover map of Uttarakhand for the year 2010 on a reduced scale Land use land cover maps at 1:250,000 scale derived using AWiFS, obtained from National Remote Sensing Centre, Department of Space, Government of India, Hyderabad.
downloaded from LP DAAC, Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS) nighttime lights time series data downloaded from National Geophysical Data Center (http://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html),Advanced Wide Field Sensor (AWiFS) imagery downloaded from Bhuvan portal developed by National Remote Sensing Centre Open Earth Observation Data Archive (NOEDA) to visualize and disseminate earth observation data (http://bhuvannoeda.nrsc.gov.in/download/download/download.php), and

Table 1 .
Criteria followed in identification and attribution of land use/ land cover class names to each cluster

Table 2 .
Table. 2. Results were also compared with Forest Survey of India publication (FSI) statistics published in India State of Forest report 2011.Area statistics of land use/land cover map derived in the present study compared with reference map It is observed that the results obtained are in close agreement with the area statistics of the reference data.Minor over estimation in water bodies and built-up lands are attributed to the use of urban mask and water mask in delineating these land cover classes.Differences in the remaining land cover classes may be attributed to the moderate resolution of the imagery used in the study.According to Forest Survey of India report, the total forest area in Uttarakhand is 24,766 sq.km and nonforest area is 28,717 sq.km (Forest Survey of India, 2011).When aggregated, the present study results showed a total forest area of 26,563 sq.km and non-forest area of 26,920 sq.km.
The methods and techniques used in the study demonstrated the ability of multi-temporal NDVI data in mapping land use/land cover.Various governmental, non-governmental and research organizations can make use of the results (available athttp://doi.pangaea.de/10.1594/PANGAEA.816654)obtained in the present study for sustainable ecological, economical and environmental development planning and management./I: Snow/Ice, GL: Grass lands, BL: Barren lands, EF: Evergreen forest, SH/DF: Shrub/degraded forest, CL: crop lands, DF: Deciduous forest, WB: Water bodies and BU: Built-up lands. S

Table 3 .
Error matrix for the land use/land cover map of Uttarakhand