Assessment of Eco-Environment Quality of Darjeeling and Kalimpong Districts, West Bengal, India using Remote Sensing based Ecological Index

Study Area

Darjeeling, located between 26°31´N-27°13´N latitude and 87°59´E-88°53´E, longitudes having 3149 sq.km area (as per 2011 census) is the northern most district of West Bengal. It consists of four sub-divisions i.e., Darjeeling, Kalimpong, Kurseong, and Siliguri. In the year of 2017, Kalimpong sub- division was declared as a separate district. Therefore, the new sub-divisions after 2017 are Darjeeling Sadar, Kurseong, Siliguri, and Mirik. The central attraction, principal town and administrative headquarter of the region is Darjeeling town, which has an area of nearly 10.75 sq.km. Geographically the district has two major divisions, hills in the north and plains in the south. The hills varying in elevation from 400 to 3,000 metres rise abruptly from the plains with a north-west to south-east alignment (Figure 1). They are part of lower or Shivalik Himalayas characterised by ridges, bold spurs, and narrow deep valleys. The plains with an average elevation of

80 to 300 meters lying in the south are marshy low-lying area called Terai. The district is drained by Teesta, Mahananda, Jaldhaka, Balason, Mechi, Rangeet, and Rammam rivers [21]. It has temperate subtropical highland climate (Köppen climate classification: Cwb). The maximum temperature reaches up to 24°C-27°C during summer and 14°C-17°C in winter, while minimum temperature ranges between 16°C-19°C and 5°C-6°C in summer and winter respectively. It receives an average annual rainfall of 220 cm per year, ranging from 1 inch in December-January to 34 inches in July-August.

Most of the ridges and valleys in Darjeeling hills are dominated by forest cover, and in many cases, they are altered into tea estates and other cultivated lands. For being a famous tourist spot in West Bengal, Darjeeling Himalaya has relatively high population density. Due to limitations in extending agricultural land, the forested and protected areas are getting encroached to meet the increasing population pressure on hills. Conversion of forest and grass lands into agriculture and settlements are mostly common in Darjeeling. Thus, changes in land use are leaving imprints on the ecology of the area. Hence, the assessment on regional ecological change can be a good example for understanding the ecological status of the region.

Click to enlarge

Database

This study is based on Landsat images (Table 1) downloaded freely from USGS (United States Geological Survey) website (https://earthexplorer.usgs.gov/) using Earth Explorer. These images were geometrically corrected and reprojected into WGS 1984 UTM Zone 45 in GIS environment.

Generation of Remote Sensing based Ecological Index (RSEI)

Water, vegetation and temperature are the major factors that affect the environment. Therefore, dryness, greenness, wetness and heat are the major parameters to study the status of environment. Dryness refers to build-induced land- area is evaluated using RSEI model and the indicators of RSEI are explained within the Pressure-State-Response (PSR) framework (Figure 2).

$$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \mathrm {E} _ {1} + \mathrm {E} _ {2} + \dots + \mathrm {E} _ {n} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$

$$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$

$$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$

$$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$

$$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} + \frac {1}{2} \mathrm {C} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \mathrm {E} _ {1} + \mathrm {E} _ {2} + \dots + \mathrm {E} _ {n} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \mathrm {E} _ {1} + \mathrm {E} _ {2} + \dots + \mathrm {E} _ {n} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} + \frac {1}{2} \mathrm {C} ^ {2} $$ $$ \mathrm {E} = \mathrm {E} _ {1} + \mathrm {E} _ {2} + \dots + \mathrm {E} _ {n} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \mathrm {E} _ {1} + \mathrm {E} _ {2} + \dots + \mathrm {E} _ {n} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \mathrm {B} ^ {2} $$ $$ \mathrm {E} = \frac {1}{2} \mathrm {A} ^ {2} + \frac {1}{2} \mathrm {B} ^ {2} $$ surface desiccation, greenness to the state of vegetation cover, wetness represents moisture condition and heat to changing temperature and climatic conditions. RSEI combines these four indicators to detect the ecological condition of a region. These indicators are first calculated individually and then integrated to obtain the synthetic value of RSEI. This is done using the ecological indices derived from Landsat data (Table

1) at the pixel level as discussed below:

**Calculation of Indicators**

Dryness

It is assessed using Normalised Difference Built-up and Bare Soil Index (NDBSI) using Equations 1-3. It is a proxy of "Dryness" which combines IBI (Index-based Built-up Index) and SI (Soil Index) to better represent built-up areas and bare soil lands [22].

Generally, NDBI (Normalized Difference Built-up Index) and IBI (Index-based Built-up Index) are used to extract built-up areas, but NDBI often mixes with plant noise, which needs further filtration using NDVI. On the other hand, IBI has higher accuracy than NDBI and SI (Soil Index) is useful to extract bare soil areas.

$$NDBSI = \frac{(IBI + SI)}{2}$$

(1)

$$IBI = \left\{ \begin{array}{l} \frac{2 \cdot \rho SWIR1}{\rho NIR + \rho NIR} - \frac{\rho NIR}{\rho NIR + \rho RED} + \frac{\rho GREEN}{\rho GREEN + \rho SWIR1} \\ \frac{2 \cdot \rho SWIR1}{\rho NIR + \rho NIR} + \frac{\rho NIR}{\rho NIR + \rho RED} + \frac{\rho GREEN}{\rho GREEN + \rho SWIR1} \end{array} \right$$

(2)

$$SI = \left( \frac{\rho SWIR1 + \rho RED}{\rho BLUE + \rho NIR} \right) / \left( \frac{\rho SWIR1 + \rho RED}{\rho BLUE + \rho NIR} \right)$$

(3)

Here, $\rho SWIR1$, $\rho NIR$, $\rho GREEN$, $\rho RED$, $\rho BLUE$ are the spectral reflectance values of the respective bands in both Landsat 5 TM and Landsat 9 OLI/TIRS images.

Greenness

Greenness is assessed through Normalised Difference Vegetation Index (NDVI) using Equation 4.

$$NDVI = \frac{\rho NIR - \rho RED}{\rho NIR + \rho RED}$$

(4)

Here, $\rho NIR$ and $\rho RED$ are the spectral reflectance values of NIR and RED bands in both Landsat 5 TM and Landsat 9 OLI/TIRS images.

Wetness

"Wetness" is described using Land Surface Moisture (LSM) which represents moisture content in soil, water bodies, and plants.

It is derived from Tasselled Cap Transformation (TCT) given by Kauth, et al. In this study the coefficients from the TCT for both TM and OLI sensors are used to get the Land Surface Moisture (LSM) index over the study area [23, 24]; using Equations 5 & 6.

$$LSM (TM) = 0.0315 \cdot \rho BLUE + 0.2021 \cdot \rho GREEN + 0.3012 \cdot \rho RED + 0.1594 \cdot \rho NIR + (-0.6806 \cdot \rho SWIR1) + (-0.6109 \cdot \rho SWIR2)$$

(5)

$$LSM (OLI) = 0.1511 \cdot \rho BLUE + 0.1973 \cdot \rho GREEN + 0.3283 \cdot \rho RED + 0.3407 \cdot \rho NIR + (-0.7117 \cdot \rho SWIR1) + (-0.4559 \cdot \rho SWIR2$$

(6)

Here, $\rho BLUE$, $\rho GREEN$, $\rho RED$, $\rho NIR$, $\rho SWIR1$, and $\rho SWIR2$ are spectral reflectance values of the respective bands for both TM and OLI/TIRS sensors of Landsat 5 and Landsat 9 respectively.

Heat

Land Surface Temperature (LST) is used as a proxy of "Heat". It is derived from the thermal bands of the Landsat Data of the study area (Table 2) using Radiative Transfer Equation (RTE) algorithm [25] using Equations 7-10.

$$L_{\lambda} = \text{gain} \cdot DN + \text{bias}$$

(For Landsat 5 TM) (7)

$$L_{\lambda} = M_L \cdot Q_{\text{cal}} + A_L$$

(For Landsat 9 OLI/TIRS) (8)

$$BT = [K_2 / \ln \left( \frac{K_1}{LA} + 1 \right) + (-273.15)]$$

(9)

$$LST = \frac{BT}{1 + \left( \lambda \cdot \frac{BT}{\rho} \right) \ln \varepsilon}$$

(10)

where $L_{\lambda}$ is the spectral radiance value (Watts/s/(m$^2$*sr*$\mu$m)). Gain is thermal infrared gain value and bias offset values; $M_L$ and $A_L$ are the band-specific Multi band and Add band values (Table 2), and D N and $Q_{\text{cal}}$ are the pixel values of band 6 and band 10 (thermal bands) of Landsat 5 TM and Landsat 9 images respectively.

BT is the Top of Atmosphere (TOA) brightness temperature (in Kelvin). $K_1$ and $K_2$ are band-specific coefficients and thermal constants in the metadata file (Table 2).

To convert the Kelvin temperature to Celsius, an absolute zero is added to the equation, i.e., -273.15.

LST is Land Surface Temperature (in °C), $\lambda$ is the average wavelength of thermal bands. $\rho$ is the constant value derived from $(h^*c/s) = 14387.685 \mu$m K; where $h = \text{Plank's constant} = 6.626 \times 10^{-34} \text{Js}$, $c = \text{Velocity of light} = 2.998 \times 10^{10} \text{m/s}$, $s = \text{Boltzmann constant} = 1.38 \times 10^{-23} \text{J/K}$, $\varepsilon$ is land surface emissivity calculated using Equation 11.

$$\varepsilon = m^* P_c + n$$

(11)

where $m$ is soil emissivity (value 0.004), $n$ is vegetation emissivity (value 0.986) and $P_c$ is the proportion of vegetation calculated using Equation 12.

$$P_c = \left[ (NDVI - NDVI_{\text{min}}) / (NDVI_{\text{max}} - NDVI_{\text{min}}) \right]^2$$

(12)

where $NDVI_{\text{max}}$ and $NDVI_{\text{min}}$ are the maximum and minimum NDVI values.

Synthesis of RSEI

The values for different factors obtained using the equations above were the inputs to derive the synthetic index (RSEI). For this, Principal Component Analysis (PCA) was performed to integrate the values for different factors. This is because it is more effective as it is able to find out the relative importance of each variable to the main component of RSEI [3, 5]. It removes the impact of co-linearity, avoids the disadvantages of subjective weight-putting methods (RSEI evaluation) and allocates weight to each factor on the basis of the factor’s contribution to the principal component [26]. Prior to performing PCA the values of all the indicators have been normalized between 0 to 1 to eliminate the unit differences among each other and get unitless values for synthesizing them into one indicator (RSEI) [27] using Equation 13.

$$ N I _ {i} = \left(I _ {i} - I _ {\min }\right) / \left(I _ {\max } - I _ {\min }\right) \tag {13} $$ Here, NIi is the Normalized value for each index. Ii denotes DN value of each image element i. Imax and Imin are the maximum and minimum values of each image element respectively.

After obtaining the PCA values RSEI values are calculated using Equation 14. However, since RSEI0 (PC1) has low values for good ecological conditions and high values representing poor ones, the RSEI0 is subtracted from one (equation 14.1) to reverse the condition. The RSEI values obtained are again normalised between 0 and 1. Here, 0 indicates very poor and 1 represents very good ecological quality. Normalized RSEI values are then classified into five categories, i.e., very poor, poor acceptable, good, and very good (Table 3).

( ) 0 RSEI PCA

NDVI, LSM, NDBSI, LST

$$ = \mathrm {P C A} \left[ f \left(\mathrm {N D V I}, \mathrm {L S M}, \mathrm {N D B S I}, \mathrm {L S T}\right) \right] \tag {14} $$ $$ \mathrm {R S E I} = 1 - \mathrm {R S E I} _ {0} \tag {14.1} $$

Where, RSEI denotes Remote Sensing based Ecological

Index, NDVI (greenness), LSM (wetness), NDBSI (dryness),

and LST (heat).

Pressure-State-Response (PSR) Framework

PSR framework is a widely accepted research tool used for spatial assessment modelling. It is built on the basis of different indicators under three major categories i.e., indicators of anthropogenic pressure, environmental state, and societal responses. Here, “Pressure” refers to the anthropogenic activities that exert pressure and stress on the natural resources or ecosystem; “State” represents the quality and health of the local environment, and “Response” indicates the ecological responses to environmental changes due to pressure and state index.

The merit of this model is that while on one hand it takes into consideration vibrant factors that critically affect the quality of environment on the other hand it links them with human activities that can provide more realistic results. Further the various indicators can be easily derived from and calculated using remote sensing data.

As the RSEI is carried out within the PSR framework, the indices of RSEI are restructured through the PSR framework by introducing NDBSI under “Pressure” category, NDVI under

Pressure State Response

Table 4 exhibits the descriptive statistics of RSEI indicators. From the table it can be observed that NDBSI assumed as a “Pressure” indicator in our study, exhibits an increase in values increasing trend. It means values increased from -0.0898 to 0.0242 in 2022 as compared to1999. Likewise, its minimum and maximum values also increased.

This indicates increasing dryness viz. expansion of human induced built-up area and bare soil during the study period. In case of NDVI, which represents “State” of the ecosystem, the average value has increased notably from 0.3782 to 0.5135 during this period. However, its maximum and minimum values have declined, indicating partial deterioration of vegetation cover which is the result of increase in built-up areas in lower elevations, expansion of tea plantation along hill slopes, and degradation of shrublands and scrubs [28].

State Forest Report of the year 2001, 2011 and 2021 also shows an increase in both open and dense forest cover during the last two decades (State Forest Report 2001, 2011 and 2021). On the other hand, the mean value of both the “Response” indicators viz. LSM and LST have increased from -0.1071 to -0.0329 and 23.63° to 24.92° respectively from 1999 to 2022. Though LSM has increased but its negative values indicate the lack of soil moisture over the region.

LST which was in the range of 4.5-39.2°C in 1999 has shown an upward trend 6.9-41.3°C in 2022. Its minimum and maximum values have increased by 2.4°C and 2.1°C respectively.

Kumar, et al. mentioned an increase of 0.5°C to 1.5°C in mean temperature in the entire Darjeeling region in last 15 years, especially in the block in the northern, western to the southern part of Darjeeling region.

Thus, the increase of pressure index (human induced LULC change) has contributed to reduction of vegetation cover (forests) in some places, but the value of greenness is high because densely vegetated areas have converted more into croplands and tea plantation as compared to built-up and bare soil areas.

Therefore, an increase in pressure (NDBSI) and state (NDVI) index have influenced response indices (LSM and LST); which ultimately have a significant and comprehensive contribution to RSEI (Figure 3).

Click to enlarge

Assessment of Ecological Environmental Quality (EEQ)

PCA Results The RSEIs of Darjeeling district for 1999 and 2022 was computed through Principal Component Analysis of the four indicators and their contributions. Table 5 exhibits the results of PCA. PC1 with highest eigen values (78.82% in 1999) and (81.6% in 2022) contains the most variable information of the four metrics, compared with PC2, PC3, and PC4 and is represented as RSEI. PC1 has both positive and negative signs, which indicates the positive and negative roles of the indicators in the ecological environment. Here NDVI and LSM have negative signs whereas NDBSI and LST have positive signs which contradicts the reality. In the real world NDVI and LSM are assumed to be positive indicators while NDBSI and LST are assumed to be negative indicators for an ecological environment. Hence PC1 is rectified using equation 14.1.

Note: PCA loads of NDVI and LSM have negative signs, causing RSEI0 (PC1) to show lower values for good ecological conditions and higher values for poor ecological conditions,

Click to enlarge

as mentioned earlier. Reverse process of subtraction of RSEI0 (PC1) from 1 will make the conversion (Equation 14.1), which will make NDVI and LSM values positive in PC1 or RSEI.

Spatio-Temporal Distribution and Change in Ecological Environmental Quality (EEQ) using RSEI

The mean RSEI value of the study area was 0.49 in 1999 and 0.58 in 2022 (Table 4) indicating an acceptable (Level 3) Ecological Environmental Quality of the study area (Table 3). Spatially higher Ecological Environmental Quality (EEQ) is observed in the north which is the part of Darjeeling Himalayas. Less human intervention owing to extremely rugged terrain and the presence of the two biosphere reserves viz. Singalila Forest in the north-west and Neora Valley National Park in the north-east has preserved the natural environment here. On the other hand, low EEQ is found in the southern part called the Terai region (Figure 4). It is acceptable to poor, and very poor in some cases. It is a plain area characterised by intense human intervention and bustling anthropogenic activities. There is very high concentration of rural and urban settlements, Siliguri the largest urban centre of the study area is located here. Besides, it has vast expanse of agricultural land, low vegetation cover and dense network of transportation which has significantly altered the natural ecosystem.

Area statistics of RSEI classes (Figure 5) shows that more than 60% area was under Level 3 (acceptable) and Level 4 (good) in 1999, which increased to 78% in 2022. Area under acceptable ecological conditions was 37.46% in 1999 which reduced to 33.3% during 2022. On the other hand, area under very good and good ecological conditions increased from 4.55% to 7.14% and 23.42% to 44.7% during 1999 and 2022 respectively. Likewise, area under poor, and very poor ecological conditions reduced from 31.49% to 13.94% and 3.08% to 0.91% respectively during the study period (Figure 5).

Click to enlarge

Table 6 exhibits the type of change that has occurred in the Ecological Environmental Quality (EEQ) of the study area in 2022 as compared to 1999. From the table it can be observed that 47.41% area had stable conditions with no change, there was degradation in 3.99 % area and 48.60% area had improvement in EEQ. Improvement in EEQ is observed in all parts of the study area. On the other hand, ecologically degraded areas are seen in the pockets of north- eastern, north-western and southern parts of the study area whereas ecologically stable area is concentrated in the northern part (Figure 6).

Click to enlarge

Relationship of RSEI with the Ecological Indicators

*Mean values are calculated by averaging the other 4 values, for example, Mean RSEI for 1999 = (0.94+0.91+0.86+0.71)/4 = 0.855 Table 7: Correlation matrix of RSEI and four ecological indicators in the year 1999 and 2022.

We have also assessed the relationship among RSEI and the four individual indicators in this study (Table 7 & Figure 7), RSEI is positively correlated with NDVI and LSM, while it has strong negative correlation with NDBSI and LST. The value is more than 0.80 in both the cases. Individually NDVI has positive correlation with LSM (r = 0.49 in 1999 and r = 0.63 in 2022), but is negatively correlated with NDBSI (r = -0.64 in 1999 and r = -0.89 in 2022) and LST (r = -0.52 in 1999, r = -0.56 in 2022). Likewise, LSM has high negative correlation with NDBSI. On the other hand, NDBSI has a positive correlation with LST (r = 0.75 and r = 0.70 for the 1999 and 2022 respectively), as built-up areas and bare soil are responsible for the trapping temperature and thus allowing the LST to rise.

Click to enlarge

Detecting a region’s ecological status and its spatio- temporal changes poses a significant challenge due to difficulty in producing time series of ecological images showing ecological status. The RSEI technique is distinct from other approaches because its results not only produce a single numeric value but also a comprehensive ecological- status image that depicts differences in regional ecological conditions. This study on the ecological quality of Darjeeling and Kalimpong Districts reveals an improvement in the ecological quality of the study area. This is because area under forests has increased by 154 km2 in the last two decades (State Forest Report 2001, 2011, 2021). Although the built- up area has also increased but it is less than the increase in area under forests. Moreover, the contribution of vegetation to the PC1-based RSEI measurement is much higher than that of built-up land. The study also found a significant spatial difference in the ecological quality, the northern parts of the study area which are a part of the Darjeeling Himalayas have high EEQ as compared to the southern parts. This is because while the northern parts are covered with thick forest cover and biosphere reserves, the southern part, the Tarai region is bustling with anthropogenic activities. Much of the forest cover has been transformed into human settlements, plantations or agricultural fields (Figures 8 & 9). It is for this very reason that the NDVI and LSM values are higher in the northern part whereas NDBSI and LST values are high in the southern part for both 1999 and 2022 (Figure 3).

Click to enlarge

Click to enlarge