Monitoring deforestation, forest health, and environmental criticality in a protected area periphery using Geospatial Techniques

Image pre processing

Pre processing of data involved several key steps to ensure that the data is clean, accurate, and ready for analysis. Firstly, radiometric calibration was performed by converting raw digital numbers (DN) to radiance or reflectance values. This process corrected sensor-specific biases. Secondly, we removed the effects of the atmosphere on the reflectance values by employing Dark Object Subtraction (DOS). To maintain spatial consistency, we performed geometric correction by projecting all data layers into the WGS 1984, UTM zone 44N projection system. In the image subsetting phase, we clipped the images to the areas of interest (AOI) using the DSD boundary. Then, we adjusted the spatial resolution of the Landsat images using the nearest neighbor resampling technique while maintaining a 30 m spatial resolution.

Image classification

The features in the images were retrieved to predetermine the forest cover changes in Musali DSD. During classification, five land use classes were identified: forest, agriculture, water bodies, built-up/settlements, and others, as shown in Table S2. There are various statistically driven supervised classification algorithms, such as maximum likelihood, parallelepiped, and Mahalanobis. The chosen classification method for this study is maximum likelihood classification (MLC). MLC, a widely used algorithm for LULC classification, relies on statistical sampling using probability density functions to detect predefined sets of LULC classes (Alawamy et al., 2020). Moreover, this method is widely used because fine-resolution satellite remote sensing data are comparatively inexpensive sources for LULC mapping (Weng, 2002; Dissanayake et al., 2019). It relies upon the likelihood that a pixel belongs to a specific class. One well-known method, MLC, based on the Bayesian equation, calculates the likelihood D of unknown measurement vector X, belongs to Eq. (1) (Mohajane et al., 2018): (1)${\displaystyle \mathrm{D}=\mathrm{In}\left(a_c\right)-\left[0.5\mathrm{In}\left(\left|Co{v}_c\right|\right)\right]-\left[0.5\left\{X-M_c\right)\mathrm{T}\left(Co{v}_c-1\right)\left(X{M}_c\right)\right]}$

In the equation the weighted distance represent by D; particular class is c; the measurement vector of the candidate pixel is represented by X. The mean vector of the sample of class c represented by M_c. The percent probability of the candidate pixel being a member of class c represented by a_c. The covariance matrix of the pixels in the sample of class c is presented by Cov_c. $\left|Co{v}_c\right|$ is the determinant of Cov_c (matrix algebra). Cov_c − 1 is is the inverse of Cov_c (matrix algebra). ln = natural logarithm function; and T = transposition function (matrix algebra) (Mohajane et al., 2018).

The majority filtering method has been employed to rectify the errors that may arise during the classification process. This approach has been utilized by previous researchers as well to bolster the reliability of their results (Weng, 2002; Dissanayake et al., 2019). Figure S1 illustrates the LULC information provided in Table S2. These images were collected from the Google Earth Pro. August 25, 2022, was the date for selecting the images as it is the closest available date to offer an illustrative representation of the LULC information.

Accuracy assessment

Securing the quality of classified images is paramount in every LULC change detection process. To achieve this, a stratified random sampling technique was employed in the procedure, ensuring coverage of all LULC classes, with 600 Ground Control Points (GCP) created each year (Wu & Murray, 2003). Subsequently, Google Earth Pro historical imagery served as reference data for accuracy assessment. Four widely used accuracy metrics were computed: Overall Accuracy (OA), Producer’s Accuracy (PA), User’s Accuracy (UA), and Kappa Coefficients (K).

OA represents the overall percentage of correctly classified LULC classes, calculated by dividing the number of accurately classified land cover pixels by the total number of pixels in the datasets using Eq. (2) (Yuh et al., 2023; Olofsson et al., 2014; Weng, 2002; Dissanayake et al., 2019; Estoque & Murayama, 2017). (2)${\displaystyle OA=\frac{1}{N}\sum _{ii=1}^n{P}_{ii}}$

where, OA is overall accuracy; N is total samples number; n is total categories number; and Pii is correct classifications number of ith sample in confusion matrix.

PA measures the percentage accuracy of individual LULC classes within a map, determined by dividing the number of correctly classified pixels in a specific land cover class by the total number of pixels belonging to that class in the reference data. Misclassified pixels within this metric are known as errors of omission and was calculated using as Eq. (3): (3)${\displaystyle PA=\frac{\text{Correctly classified number pixel in each category}}{\text{Correctly classified total number pixels in that category (column total)}}}$

UA evaluates the reliability of a given land cover map concerning its agreement with ground observations. It is calculated by dividing the number of correctly classified pixels in a specific land cover class by the total number of pixels classified within that class. Similar to Producer’s Accuracy, misclassified pixels in this metric are referred to as errors of omission. Equation (4) is as below: (4)${\displaystyle UA=\frac{\text{Correctly classified number pixel in each category}}{\text{Correctly classified total number pixels in that category(row total)}}}$

K indicates the level of agreement between test and validation data in a generated land cover map. It is based on the probability of the test data closely matching the validation data during the land cover mapping process and is highly correlated with overall accuracy. K coefficient was utilized as a metric to gauge the agreement between classified results and actual conditions, thereby determining the levels of accuracy. A threshold of 0.75 or higher was set, indicating sufficient agreement for actionable insights based on the image. The Eq. (5) was used to calculate the K (Alqurashi & Kumar, 2014). The results of the accuracy assessments are presented in Table S3. (5)${\displaystyle K=\frac{N\sum _{i=j}^k{x}_{ii}-\sum _{i=1}^k\left(x_{i+}\times x_{+i}\right)}{N^2-\sum _{i=1}^k\left(x_{i+}\times x_{+i}\right)}}$

where k is the number of rows in the matrix; x_ii is the number of observations in row i and column i; x_i+ and x+_i are the marginal totals of row k and column i. N is the number of observations (Alqurashi & Kumar, 2014).

Change detection

A change detection method was used to cross-tabulate LULC data between different periods after image classification: (i) 1988 vs. 1996, (ii) 1996 vs. 2009, (iii) 2009 vs. 2022, and (iv) 1988 vs. 2022.

Retrieval of LST

LST is commonly defined as the temperature at the interface between the Earth’s surface and its atmosphere and this serves as a critical indicator in all physical processes concerning surface energy and water balance at both micro and macro spatial scale (Malik, Shukla & Mishra, 2019; Sobrino, Jiménez-Muñoz & Paolini, 2004). This plays a pivotal role in surface processes and boasts a wide array of applications, including vegetation stress monitoring, climate change studies, environmental assessment, and urban climate analysis. With the availability of large-scale remote sensing data, near-surface air temperature measurements can be effectively monitored and utilized to retrieve the temperature of various LULC surfaces (Malik, Shukla & Mishra, 2019). In our study, Landsat 5 and Landsat 8 data were employed to analyze LST and investigate the impact of forest degradation. The LST was obtained by Thermal Infrared (TIR) band 6 of Landsat 5 and TIR band 10 and 11 of Landsat 8 by converting DN values into radiance values (Jaafar et al., 2020; Ranagalage, Estoque & Murayama, 2017; Dissanayake et al., 2019). Here we used thremal bands holding brightness temperatures which are represented in Kelvin. Before retrieval, the LST land surface emissivity values was derived using Eq. (6): (6)${\displaystyle \varepsilon =m{P}_{\mathrm{V}}+n}$

where ɛ represents land surface emissivity; m represents (ɛ_v − ɛ_s) − (1 − ɛ_s) ɛ_v; P_v represents the amount of vegetation; n represents ɛ_s + (1 − ɛ_s) ɛ_v ; ɛ_s is soil emissivity; ɛ_v is the vegetation emissivity; and F is a shape factor (Ranagalage, Estoque & Murayama, 2017; Dissanayake et al., 2019). Here we used m = 0.004 and n = 0.986 as in previous research (Jaafar et al., 2020; Ranagalage, Estoque & Murayama, 2017; Dissanayake et al., 2019). The proportion of vegetation (PV) is derived from Eq. (7). (7)${\displaystyle \text{PV}={\left(\left(\text{NDVI}-\text{NDVImin}\right)/\left(\text{NDVImax}-\text{NDVImin}\right)\right)}^2}$

where NDVI is the normalized difference vegetation index derived from Eq. (9) as in sub section below. The NDVImin and NDVImax are the minimum and maximum values of the NDVI, respectively. Then emissivity corrected LST were retrieved using Eq. (8). (8)${\displaystyle LST\left(°C\right)=\frac{TB}{1}+\left(\lambda \times TB/p\right)In\varepsilon }$

where TB = Landsat TM Band 6 at-satellite brightness temperature; λ = wavelength of emitted radiance (λ = 11.5 µm for Landsat TM Band 6, λ = 10.8 µm for Landsat TIRS Band 10) (Ranagalage, Estoque & Murayama, 2017; Dissanayake et al., 2019); p = h × c/σ(1. 438 × 10 –2 mK), σ = Boltzmann constant (1. 38 × 10 − 23 J/K), h = Planck’s constant (6. 626 × 10 − 34 Js),c = velocity of light (2. 998 × 108 m/s), ɛ is the land surface emissivity. Last, the LST values of Kelvin were converted into degrees Celsius (°C).

Derivation of NDVI

The results of spectral analysis are typically summarized into vegetation indices, which establish relationships between reflectance across two or more wavelength intervals or bands in satellite images. Measurement of such indices serves as an especially valuable tool for precision agriculture and vegetation analysis (Fernández-Alonso, Hernández & Torres-Costa, 2023). Various sensors and multi spectral devices are extensively utilized nowadays to assess the condition of vegetation, including chlorophyll meters, canopy reflectance sensors, and Plant-O-Meter (Padilla et al., 2018; Kitić et al., 2019). NDVI is the most common vegetation index, frequently employed as an indicator of chlorophyll content and overall vegetation health (Fernández-Alonso, Hernández & Torres-Costa, 2023). The NDVI is a vital remote sensing metric, acting as a key indicator for assessing vegetation coverage and growth. This quantifies vegetation through the spectral reflectance and derived from the ratio of the difference to the sum of values in the red and near-infrared bands, with values typically ranging from −1.0 to 1.0 (Yin et al., 2024). Near infrared (NIR) (band 4 in TM and band 5 in OLI) and red bands (band 3 in TM and band 4 in OLI) of Landsat are needed to retrieve the NDVI. Vegetation density and chlorophyll activity variations were measured using the Landsat images which were taken during the dry period of the study area using Eq. (9) (Jaafar et al., 2020; Ranagalage, Estoque & Murayama, 2017; Dissanayake et al., 2019). (9)${\displaystyle \text{NDVI}=\left(\text{NIR}-\text{RED}\right)/\left(\text{NIR}+\text{RED}\right)}$

The retrieved maps were classified upon their values.We used the thresholding method to classify NDVI into four categories: non vegetation (less than 0) low-density vegetation (0–0.2), moderate-density vegetation (values between 0.2 and 0.5), and high-density vegetation (values above 0.5) following the similar work conducted by Mohajane et al. (2018).

Derivation of VCI

VCI is a widely used drought index which developed by Kogan (1990) for monitoring vegetation drought stress and estimating drought trends. Kogan (1990) developed this index to enhance the weather-related components in the NDVI value (Ha, Uereyen & Kuenzer, 2023). It is constructed by normalizing long-term satellite-based NDVI data. The VCI offers several advantages, such as eliminating envelope signal variations in NDVI and accounting for regional climate diversity. It is also easy to construct using only NDVI data. Consequently, it has various applications, including drought investigation, crop yield estimation, and vegetation dynamics assessment (Yin et al., 2024; Ha, Uereyen & Kuenzer, 2023). This helps to compare current NDVI values with past year NDVI values of the same season. Equation (10) was used to calculate VCI for the selected years (Yin et al., 2024; Ha, Uereyen & Kuenzer, 2023; Kogan, 1990; Dutta et al., 2015). (10)${\displaystyle \text{VCI}=\left(\text{NDVI}-\text{NDVImin}\right)/\left(\text{NDVImax}-\text{NDVImin}\right)\ast 100}$

where NDVI indicates the value of a specific pixel in that particular month, NDVImax and NDVImin show the multi-layer highest and lowest NDVI values in the same period. The VCI values represented as % and 0 is for extreme drought conditions and 100 is optimal vegetation health. Following (Ha, Uereyen & Kuenzer, 2023; Kogan, 1990; Dutta et al., 2015), we classified VCI values into five categories upon drought severity: non drought (100%–50%), mild drought (50%–30%), moderate drought (30%–20%), severe drought (20%–10%), and extreme drought (less than 10%).

Derivation of ECI

Due to rising LST and declining NDVI, the environment is in a critical state, as measured by the ECI. Increases in LST have been shown to directly correlate with ECI, while decreases in NDVI have been reported to have an inverse correlation with ECI (Senanayake, Welivitiya & Nadeeka, 2013; Saputra, Jamadi & Sari, 2023). Based on the ratio between LST and NDVI the ECI is calculated to identify environmentally critical areas (Senanayake, Welivitiya & Nadeeka, 2013; Ranagalage, Estoque & Murayama, 2017)). The LST and NDVI layers are used to derive the ECI using Eq. (11). The retrieved NDVI and LST layers were first normalized using the histogram equalization method, which ranges pixel values between 1–255 (Ranagalage, Estoque & Murayama, 2017). The higher the ECI value, the more environmentally critical. The spatial variation of ECI over the study area was interpreted as high, medium, and low. (11)${\displaystyle \text{ECI}=\text{LST}\left(\text{Stretched}1-255\right)/\text{NDVI}\left(\text{Stretched}1-255\right).}$