Land use/land cover changes in the central part of the Chitwan Annapurna Landscape, Nepal

Study area

The CHAL in the central Nepal is drained by eight major rivers (Kali Gandaki, Seti, Madi, Marshyandi, Daraudi, Budi Gandaki, Trishuli and Rapti) and their tributaries. This landscape covers all or parts of six protected areas and 19 districts (WWF, 2013a). We have chosen the central part of CHAL that connects two biologically important protected areas, the Chitwan National Park (CNP) in the south and the Annapurna Conservation Area (ACA) in the north. This part of CHAL has given the highest priority corridor for landscape level connectivity (WWF, 2013b). The intensive study area covers Chitwan (around Barandabhar Corridor and surrounding areas), Tanahun (Seti River basin), Kaski and some parts of Syanja and Parbat districts (Panchase and part of Annapurna Conservation Area) with an area of 2,749.48 km² (Fig. 1). The elevation ranges from 150 to 3,300 m. The lowland part has tropical and subtropical types of climate, whereas mid-hills have the temperate type of climate and the upper mountain region has subalpine type has of climate.

This landscape is rich in biodiversity, including three Global 846 Ecoregions (Terai–duar Savanna and Grasslands, Himalayan Subtropical Broadleaf Forests, Himalayan Sub-tropical pine forest) (Dinerstein et al., 2017; Wikramanayake, Dinerstein & Loucks, 2002) and two Ramsar sites (Beeshazari and associated lakes, Chitwan and Lake Clusters of Pokhara valley, Kaski) (National Lake Conservation Development Committee of Nepal (NLCDC), 2020). This area is prime habitat for many important mammal species, birds, herpetofauna, fish and many other micros and macroinvertebrates (Bhuju et al., 2007; WWF, 2013b).

Data sources

Landsat images from 2000, 2010 and 2020 were used to detect the LULCC within the 10-year time interval. The Landsat 7-ETM (Enhanced Thematic Mapper) for 2000, Landsat 5-TM (Thematic Mapper) for 2010, and Landsat 8-OLI (Operational Land Imager) for 2020, images with same 30 m spatial resolution were downloaded from the United States Geological Survey (USGS) (https://glovis.usgs.gov/app) geoportal. A total of six scenes of satellite images of two from each year were downloaded (Table 1). The entire Landsat images consist of around 3–10% of cloud cover, but this was less than 1% in our study area. We also used the topographic maps with 1:25,000 and 1:50,000 scales developed by the Department of Survey, Government of Nepal.

In addition, the Google Earth and a classified map of 2010 developed by the International Centre for Integrated Mountain Development (ICIMOD) (http://rds.icimod.org) used as a reference for verification. The reference field data were collected using a Global Positioning System (GPS) during the field study and used as ground-truthing points during classification of images and accuracy assessments (Table 1).

Image pre-processing

Each band of Landsat image was checked using metadata and georeferenced to the WGS_84 datum and Universal Transverse Mercator (UTM) Zone 44 or 45 North coordinate system. The details of bands and resolution are mentioned in Table S1. We only used the bands with 30 m resolution for further analysis. Landsat 5 TM images have seven bands, Landsat 7 ETM images have eight spectral bands, and Landsat ETM has 8 bands in which 1 to 7 bands have 30 m resolution (Barsi et al., 2014). Similarly, Landsat 8 OLI images have 11 bands in which eight bands 1 to 7 and 9 have 30 m resolution. (https://www.usgs.gov) (Table S1). For the natural color composite of Landsat 8 OLI images, band 4 (red), 3 (green) and blue (2) were combined for the natural color, whereas bands 7, 6, 4 were used for false color (urban). Similarly, bands 5, 4, 3 for vegetation composition, bands 6, 5, 2 for agriculture, 5, 6, 4 for land and water (Vermote et al., 2016; Barsi et al., 2014).

The images were processed in ERDAS IMAGINE 9.2. Bands of each satellite image (2000, 2010 and 2020) were stacked within Raster main icon with layer stack function as a single layer. In this study, we selected band 1 to 5 and 7 (blue, green, red, near infrared (NIR), shortwave infrared I (SWIR1) and shortwave infrared II (SWIR2)) for Landsat 5 TM and Landsat 7 ETM; band 1 to 7 (coastal, blue, green, red, NIR, SWIR1 and SWIR2) for Landsat 8 OLI in land use and land cover classification. Band 8 to 11 of Landsat image 2020 are less used in LULCC (Yu et al., 2019). The images of each scene were masked using the Area of Interest (AOI) of the study area using mask function (Fig. 2).

Ground-truthing points

For the data collection (ground truthing coordinates) in the central part of Chitwan-Annapurna Landscape, we obtained permission from the Department of National Parks and Wildlife Conservation (Permission letter number 3372), Chitwan National Park (Permission letter number 2723), Division Forest Offices of Chitwan (Permission letter number 2723), Tanahun (Permission letter number 749), Kaski (Permission letter number 200) districts and Annapurna Conservation Area Project (Permission letter number 66). The field survey which was carried out from 2018 to 2020 provided a clear idea about the field, forest types and land cover types. For ground-truthing, geographic coordinates were collected during the sign survey of large mammals, including leopard and their prey using GPS (Garmin eTrex 10). These geographic coordinates represented all land cover types along the landscape. Each coordinate was taken from the central point of the land cover patches which was more than 30 m radius. A total of 1,350 coordinates were collected (259 from Sal dominated forest, 125 from riverine forest, 299 from mixed forest, 125 from grasslands, 88 from barren areas, 135 from developed areas, 92 from water bodies and 229 from cropland). Out of the total sampling coordinates, half of the coordinates (667) were used for supervised classification and the remaining coordinates (683) were used for accuracy assessment. In addition to this, we also used printed versions of topographic maps to locate the different land cover types including changes over there through participatory GIS (pGIS) techniques. pGIS studies consider that the local people are familiar and experience with change to their surroundings and provide the greater spatial information about the area (Aynekulu et al., 2006; Zolkafli, Brown & Liu, 2017). For this purpose, focus group discussions were performed with members of community forests and elderly people who inhabited for a long time in that area and easily felt the changes in their surroundings. Twenty group discussions were arranged in different locations of the landscape (five discussions on Barandabhar and associated area, ten on the Seti River basin of Tanahun, five on Panchase and lower part of the ACA).

Image classification

The consistency of the land cover classes at national, regional and international level is not same (Chettri et al., 2013; Nepal Ministry of Land Reform and Management (MoLRM), 2015; Uddin et al., 2015b; Wang et al., 2020). In the present study, land cover classification was established with the help of published literature and maps (Khanal et al., 2020; Nepal Ministry of Forests and Environment (MoFE), 2019; Nepal Ministry of Land Reform and Management (MoLRM), 2015; Thapa, 2011; Uddin et al., 2015b; Zomer, Ustin & Carpenter, 2001). We classified the land cover of the central part of CHAL into eight major classes, based on the dominant plant species, human settlements, landscape and agriculture. We categorized the forest types as Sal dominated forest, riverine forest and mixed forest (Table 2). The dominant plant species composition in the mid-hills is of mixed type and difficult to separate into other subcategory, hence, we classified such forest as mixed forest.

Unsupervised classification

In the beginning, the unsupervised classification of the multi-temporal Landsat images of 2000, 2010, 2020 was performed. This classification is based on the automatic identification and assignment of image pixels to spectral grouping. It starts with a spectral plot of the whole image and group the pixels with similar features. Two common algorithms are used for the creation of the clusters in unsupervised classification (Duda & Canty, 2002). They are k-means clustering and Iterative Self-Organizing Data Analysis Technique (ISODATA) (Ragettli, Herberz & Siegfried, 2018). In this classification, we used k-means algorithm. The nearest likelihood with 10 iterations were used to group the pixels having similar features. The images were classified into 40 classes with a convergence threshold 0.90. Then, the similar classes were merged into eight different classes using recoding of classes (Table 2 and Fig. 2). The unsupervised classification of images was used for the planning of field data collection that provided the basic field knowledge. The unsupervised classes were revised after the collection of ground-truthing points.

Supervised classification

The supervised classification was performed using the widely used parametric classification algorithm namely Maximum Likelihood Classification (MLC) (Chamling & Bera, 2020; Rai et al., 2018). The signature classes or training sets were prepared from ground-truthing points for 2020 and Google Earth map for 2000 and 2010 were used to prepare signature classes for supervised classification. Two separately classified Landsat images were mosaicked to make a single image. Finally, the images were filtered fixing the pixels 3 × 3 for smoothing the image and avoid the errors of misclassification. The images were again recoded based on field knowledge to minimize the errors of misclassification. We selected five sites, two from low land (Barandabhar and associate area), two from mid hill (Seti River basin and Panchase area) and one from an upland area (lower part of the ACA) for the separate analysis where the land use/land cover was changed drastically within the land 20 years.

Accuracy assessment

Accuracy assessment increases the quality of the remotely sensed data on classified thematic maps. It compares the classified image with ground truthing points (Congalton, 2001; Rai et al., 2018; Siddique et al., 2020; Song et al., 2001; Thapa, 2011). Another common method to assess the accuracy of the classified map is to generate stratified random points as the classified class. These random points compared with the Google Earth and topographic maps as reference for verification (Crowley & Cardille, 2020). The topographic maps of Nepal were used as reference of settlements or developed area, water resources, croplands and forest area. In this study, ground-truthing points (n = 683) were used as reference for the accuracy assessment of classified images of 2020 (Fig. 3). For Landsat images of 2000 and 2010, 500 stratified random points were generated and compared them with references such as Google Earth, topographic maps of Nepal (for water bodies, settlements, urban or developed area and forest) and the classified maps of ICIMOD (for the classification of forest and grassland). The evaluation was performed computing confusion matrix or error matrix and Kappa Coefficient (Congalton, 2001; Foody, 2002). The user’s accuracy, producer’s accuracy, overall accuracy was obtained from the error matrix. The user’s accuracy provides the reliability that the classified pixels of the map match with the ground-truthing points (Eq. (2)). Similarly, the producer’s accuracy determines the probability of correctly classified reference pixels (Eq. (3)). The overall accuracy was calculated by dividing the correctly classified pixels by the total number of reference points (Eq. (1)) (Congalton, 2001; Foody, 2002). Kappa Coefficient ( $\hat{K}$) is used to measure the agreements between model prediction and reality (Congalton, 2001). It is the multivariate analysis technique to evaluate the accuracy of the classified map statistically. The Kappa Coefficient ( $\hat{K})$ ranges from 0 to 1. If the value of $\hat{K}$ is 0, this reflects there is no agreements, 0–0.2 signifies as slight, 0.21–0.40 as fair, 0.41–0.60 as moderate, 0.61–0.80 as satisfactory or good and 0.81 to 1 as almost perfect agreements (Maingi, Kepner & Edmonds, 2002). Statistically, the $\hat{K}$ was calculated using Eq. (4).

(1) $${\displaystyle \mathrm{O}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{Totalnumberofcorrectlyclassifiedpixels}{Totalnumberofreferencepixels}\times 100}$$

(2) $${\displaystyle \mathrm{U}\mathrm{s}\mathrm{e}\mathrm{r}{}^{\prime }\mathrm{s}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{Numberofcorrectlyclassifiedpixelsineachcategory}{Totalnumberofclassifiedpixelsthatcategory\left(rowtotal\right)}\times 100}$$

(3) $${\displaystyle \mathrm{P}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{Numberofcorrectlyclassifiedpixelsineachcategory}{Totalnumberofclassifiedpixelsthatcategory\left(columntotal\right)}\times 100}$$

(4) $${\displaystyle \mathrm{K}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}(\hat{\mathrm{K}})=\frac{N\left({\sum }_{i=1}^r{X}_{ii}\right)-{\sum }_{i=1}^r\left(X_{1+}\times X_{+i}\right)}{N^2-{\sum }_{i=1}^r\left(X_{i+}\times X_{+i}\right)}}$$where, r = Number of rows in the error matrix

Xii = number of observations in row i and column i (on the major diagonals)

Xi+ = Total number of observations in rows i

X+i = Total number of observations in column i

N = Total number of observations included in matrix

Land use/land cover classes and change

Out of eight land cover classes of 2020, mixed forest was the most dominant (44.95%) followed by croplands (28.3%), Sal dominated forest (15.4%) and developed area (4.13%) (Table 3, Fig. 4).

The results of LULCC from 2000 to 2010 indicated that there was a decrease in water bodies, barren land, grassland, riverine forest and croplands by 0.9%, 7.7%, 6.2%, 13% and 16%, respectively; build-up or developed area, Sal dominated forest and the mixed forest were increased by 19.1%, 4.62% and 18.2%, respectively. Similarly, from 2010 to 2020, water bodies, riverine forest, Sal dominated forest, developed area and mixed forest were increased by 2.54%, 2.09%, 3%, 10.3% and 16.3%, respectively. Barren area, cropland and grasslands were decreased by 13.3%, 16.3% and 46.4% respectively (Table 4, Figs. 5 and 6). Overall, from 2000 to 2020, the areas of grassland, riverine forest, cropland and barren area were drastically decreased, whereas developed area, mixed forest and Sal dominated forest were increased (Table 4).

The separate analysis of LULCC between 2000 to 2020 in old Padampur and associated areas (low land) clearly showed that more than 93% of the total cultivated land was changed into the grassland and forest. Similarly, the barren area (flood plain of Rapti River) was reduced by 74.67%. However, grassland, riverine forest and mixed forest in the old Padampur and associated areas were increased by 94.45%, 91.26% and 62.5%, respectively (Figs. 7A1–7A3, 8A, Table S2). The trend of land cover change from 2000 to 2020 in new Padampur and associated areas (low land) indicated that the riverine forest, Sal dominated forest and grassland were drastically reduced by 61.21%, 54.14% and 64.88%, respectively, whereas the cropland and developed areas were increased by 88.17% and 1433.33%, respectively (Figs. 7B1–7B3, 8B, Table S2). Land cover change from 2000 to 2020 in Byas municipality of Tanahun district and surrounding areas showed a significant reduction in the cropland by 40.86%, whereas there was a significant increase in developed areas and mixed forest by 86.55% and 62.14%, respectively. The trend of land cover change in Byas and surrounding areas was more between 2010 to 2020 than 2000 to 2010 (Figs. 7C1–7C3, 8C, Table S2). The results of land cover change analysis of Panchase Protected Forest and associate areas between 2000 to 2020 showed a reduction in cropland by 51.92% and grassland by 43.22%, whereas an increase in mixed forest and Sal dominated forest by 68.1% and 23.29%, respectively (Figs. 7D1–7D3, 8D, Table S2).

The results of land cover change analysis of a part of the ACA between 2000 to 2020 clearly showed an increase in mixed forest and developed area by 14.93% and 166.66%, respectively, whereas a decrease in cropland, barren area and grassland were decreased by 40.97%, 24.09% and 19.94%, respectively (Figs. 7E1–7E3, 8E, Table S2).

Accuracy assessment

The overall accuracy of classified images of 2000, 2010 and 2020 was 81%, 81.6% and 84.77%, respectively. The user’s accuracy ranged from 73.33% to 87.09% in 2000, 73.68% to 83.33% in 2010 and 80.26% to 90.69% in 2020. The low range of user’s accuracy in barren area in 2000 (73.33%), in a developed area in 2010 (73.68%) indicated confusion during land cover classification (Tables 5). Riverine forest in 2000, mixed forest in 2010 and Sal dominated forest in 2020 were more reliable with user accuracy of 87.09%, 83.77% and 90.69%, respectively (Tables 5, Tables S3–S5). The Kappa coefficient for the years 2000, 2010 and 2020 were 0.76, 0.79 and 0.82, respectively.