Quantitative Analysis of Land Use and Land Cover Dynamics using Geoinformatics Techniques: A Case Study on Kolkata Metropolitan Development Authority (KMDA) in West Bengal, India

Abstract

One of the most valuable approaches in spatial analysis for a better understanding of the hydrological response of a region or a watershed is certainly the analysis of the well-known land use land cover (LULC) dynamicity. The present case study delves deeper into the analysis of LULC dynamicity by using digital Landsat TM and Landsat OLI data to classify the Kolkata Metropolitan Development Authority (KMDA) into seven classes with over 90% classification accuracy for decadal level assessments of 30 years (for the years 1989, 1999, 2009, and 2019). The change index, the Dematel method for analyzing the cause-effect relationship among the LULC classes, the Jaccard Similarity Index for measuring the nature of similarity among the LULC classes, and the Adherence Index for measuring the consistency of the LULC classes after the transition was used in this study to analyze the LULC transformation. In more detail, the present study considers how urban land use is altering at the expense of other land uses. Besides the shifting pattern of mean centers of the LULC classes through time, also gives a very significant insight into the LULC dynamics over 30 years of span. The current study of LULC dynamicity and transformation patterns over the 30 years of the KMDA area is expected to assist land and urban planners, engineers, and administrators in sustainable decisions and policies to ensure inclusive urbanization that accommodates population growth while minimizing the impact on potential natural resources within the whole study area.

1. Introduction

Since the beginning of civilization, man and nature have had an intimate relationship, and the rapid increase in population through the decades has led to anthropogenic impact on land and water resources management. Therefore, planning and management of land use pattern and their spatial distribution and change over time play a crucial role in the sustainable development of any region [1].

One of the most crucial measures of a region’s economic development is the dynamism of land use and land cover (LULC). Its development demonstrates how human behavior and needs change over time, as well as how they prioritize some land uses over others. Its dynamism is associated with several processes, including urbanization, industrialization, commercialization, and other related ones, each of which has a unique impact on the current environment. A crucial step in the analysis of the given spatial problem is choosing the proper spatial-temporal scale for the study of its variations.

Most of the earlier research focuses on using modeling methods to investigate the pattern of LULC alterations. Jiangle’s potential future land use patterns in China were predicted by Liping et al. [2] using a Cellular Automata (CA) Markov analysis. Similarly to this, Hiping et al. [3] have predicted LULC changes in Beijing (China) using the CLUE-S model in conjunction with a Markov model. Li et al. [4] have used a novel method consisting of an artificial neural network (ANN)-based CA analysis in a geographic information system (GIS) platform, to predict the land use pattern over the territory of southern China.

The conventional techniques typically employed for detecting changes in land use and land cover are expensive, imprecise, and only effective in a small region. Due to the multiple challenges that regularly occur in environmental research [5,6,7], new technologies such as satellite remote sensing and GIS are needed. These technologies enable researchers to investigate and track the dynamics of natural resources, providing information for environmental management. Applications of remote sensing analysis can assist planners or environmental managers in making decisions about sustainability by utilizing a quantitative and model-based approach, in addition to the object-oriented component of the result. Land use or land cover change is intrinsically related to the interaction of natural and human impacts on environmental change. Different types of land use and their persistence are what influence changes in the biosphere’s condition and bio-geochemical cycles [8,9,10].

Monitoring and modelling of land use/land cover patterns have become more consistent because of the introduction of high spatial resolution satellite imagery as well as improved image processing and GIS technology over the wide applicability based on spatio-temporal analysis. The multiple advantages of remote sensing lead to the chance of covering three decades of land use dynamics along with a higher resolution from satellite platforms [11,12,13,14,15,16,17]. In more detail, Landsat-TM images are freely available and provide valuable and continuous records of the earth’s surface, as well as a wealth of information for identifying and monitoring changes in manmade and physical environments over the last three decades [18,19]. In India, maps at a scale of 1:250,000 are created using the multitemporal Indian Remote Sensing (IRS) satellites for spatial accounting and trend monitoring of land usage (NRSA, 1989). By providing a fundamental understanding of the current land use pattern, satellite remote sensing can assist in determining future improvements and management. In the modern world, satellite remote sensing technology is very helpful for creating maps of land use and land cover and evaluating them [20]. It is easy to update an existing database for various land use planning and design purposes due to the availability of repetitive data. As a result, the compilation of LULC maps and information has become more effective with the application of geospatial technologies such as remote sensing, GIS, GPS, and computational approaches [21,22,23,24]. Geospatial technology ensures that real-time data and geospatial information are available and accessible rapidly for resource mapping. Several change detection algorithms based on remote sensing images have recently been introduced. A multitude of change detection strategies and algorithms have been developed and their advantages and disadvantages assessed. The most used techniques in classification analyses include unsupervised hard classification or clustering, supervised hard classification, principal components analysis (PCA), hybrid classification, and fuzzy classification [25,26,27,28,29,30].

It is essential to have information about existing land use and land cover over regional areas in making optimal use of land. It is also important to have the ability to monitor the changing dynamics of land use because of both increasing population demands and natural factors shaping the landscape. As a result of the variety of natural and man-made processes, the land is always transforming. Understanding the evolution of land use/cover systems and studying spatio-temporal patterns of intra- and inter-land use/cover systems are still important goals in land use studies [31]. The method of detecting changes in the process of an object or situation by observing it at different times is known as land use/land cover change detection [32,33,34]. According to Macleod and Congation [35], four aspects of change detection are relevant when monitoring natural resources. The authors entail, first and primarily, recognizing the changes that have occurred; second, defining the type of the change; third, evaluating the change’s area extent; and last, rating the change’s spatial pattern.

The objective of the study is to determine the changing trajectory, stability of each land use/cover element, identifying the land use/cover elements having alike changing characteristics and the interaction between the land use/cover elements. Therefore, considering all these above-mentioned aspects as a baseline, in this present study the status of land use/cover and its reconfiguration during the last 30 years (1989, 1999, 2009, and 2019) in the Kolkata Metropolitan Development Authority (KMDA) region of West Bengal has been analyzed in a quantitative manner using geoinformatics. The outcomes of this study represent a very useful tool for both hydrology and hydraulic researchers in the assessment of land use/cover impacts on the ecohydrological behavior of both urban and natural watersheds.

2. Materials and Methods

2.1. Study Area

Kolkata Urban Agglomeration (UA) incorporates the administrative areas of the Kolkata Metropolitan Development Authority (KMDA) encompassing 3 municipal corporations (Howrah, Kolkata, and Chandan Nagar), 38 municipalities, 77 non-municipal urban towns, 16 outgrowths, and 445 rural villages [36]. KMDA functions as a statutory authority under the administrative control of the Urban Development Department of the Government of the State of West Bengal in India.

As shown in Figure 1, KMDA incorporates 6 districts i.e., Kolkata, Howrah, Hooghly, Nadia, North 24 Parganas, and South 24 Parganas under its jurisdiction. Kolkata UA is located between 22°00′19″ N and 23°00′01″ N latitudes and 88°00′04″ E and 88°00′33″ E longitudes with a surface coverage of 1851 km². The agglomeration is characterized by a concentrated population and dense settlements along both sides of river Hooghly. It has a total population of 14.72 million people, with an average density of 7950 persons/km² [32], and a annual population growth rate of 1.8% [33,34,35,36].

The population dynamics of the region witnessed a massive migration influx, which has caused many problems such as lack of land, overpopulation, pressures upon the existing resource potentials, etc., to name a few, and those problems are now leaving its indelible mark on the anthropogenic landscape over the years. These processes lead to massive LULC changes over the years, especially in such a metropolitan area as Kolkata [37,38].

2.2. Data Sources

The study area was covered by Landsat digital data from the TM, ETM+, and OLI/TIRS sensors for the years 1989, 1999, 2009, and 2019 on the USGS Earth Explorer website (http://earthexplorer.usgs.gov, accessed on 5 December 2022). The availability of the Landsat images and their medium to high spatial resolution led to their selection for this study. The detailed information regarding the data is provided in

Table 1. Detailed features (i.e., metadata) of the datasets analyzed in the present case study, with different sensors.

Year	Sensor	Path/Row	Band Count	Spatial Resolution	Radiometric Resolution
1989	TM	138/44, 138/45	7	Optical 30 m, Thermal 120 m	8 bits
1999	ETM+	138/44, 138/45	9	Optical 30 m, Thermal 60 m, Pan 15 m	8 bits
2009	TM	138/44, 138/45	9	Optical 30 m, Thermal 60 m, Pan 15 m	8 bits
2019	OLI, TIRS	138/44, 138/45	11	Optical 30 m, Thermal 100 m, Pan 15 m	16 bits

. These three Landsat digital data are referenced to the UTM map projection (Zone 45 N) using WGS84 geodetic datum and are Level-1 Terrain corrected (L1T).

2.2.1. Land Use/Cover Mapping Based on Digital Classification

The maximum likelihood algorithm and the hard classification approach were used to the FCCs with the best band combinations during the three decades to create the land use/cover classification. This method assumes that all classes’ probability densities are equal and is based on the probability density of pixels belonging to a given class. The maximum likelihood algorithm determines the spectral distance between the measurement vector for the candidate pixel and the mean vector for each signature within a pixel cluster, just like other supervised classification techniques. The study area analyzed in the present study case was classified by the means of this classifier algorithm into the following seven land use/cover classes: urban built-up area, homestead with the plantation, fallow land, bare land, waterbody, agricultural land, and vegetation by assigning per-pixel signatures from the corresponding satellite digital data on the basis of the specific pixel reflectance value.

Aiming at improving the accuracy of the classification and at reducing the presence of misclassification errors, post-classification refinement was therefore employed for guaranteeing the effectiveness and reliability to this method [39]. Moreover, the use of data characterized medium spatial resolution such as that of Landsat digital image pixels, is a very common problem still unsolved in the literature [40]; especially for urbanized surfaces such as the study area examined in the present study case where a heterogeneous mixture of features can be recognized [41].

Assessment of classification accuracy of analyzed and classified images of three decades years was carried out to determine the quality of information derived from the data by employing the error matrix and Kappa (K-hat) statistics. To use the classified outputs in detecting changes, it is essential to perform an accuracy assessment for individual classification [42]. For the accuracy assessment of land use/cover maps extracted from satellite images, the stratified random method was used to represent different land cover classes of the area. The accuracy assessment was carried out using 275 points, based on ground truth data. The comparison of reference data and classification results was carried out statistically using error matrices. In addition, a nonparametric Kappa test (K-hat) was also performed here to measure the extent of classification accuracy as it not only accounts for diagonal elements but for all the elements in the confusion matrix [43].

Kappa (k) is a measure of agreement between predefined producer ratings and user-assigned ratings. It is calculated as follows [1,2,3]:

$$k=\frac{\mathrm{Pr}\left(a\right)-\mathrm{Pr}\left(e\right)}{1-\mathrm{Pr}\left(e\right)},$$

(1)

where Pr(a) represents the observed agreement, while Pr(e) indicates the chance agreement.

2.2.2. Reconfiguration Detection of Land Use/Cover Categories through Time

Reconfiguration in this study is defined as the change from one land use/cover element to another (s). The post-classification outputs were integrated during the re-configuration study using a GIS platform. Reconfiguration is a local effect within an event. It involves the movement and change in event [33]. In the context of the temporal land use and land cover study, the term reconfiguration is pertinent. That would convey the arrangement of land use and land cover elements in different forms or combinations, actually.

Post-classification reconfiguration has been used successfully by many researchers in land use/cover dynamics assessments due to its effectiveness in detecting the location, nature, and rate of changes [44]. The quantitative determination of conversions from one land use/cover category to another on a “one-to-many” basis and their corresponding area across the assessed phase were determined using cross-tabulation analysis on a pixel-by-pixel basis. The three temporal classified maps, each of which had a unique “from–to” transition class, were consequently combined to produce a new thematic layer.

The change magnitude and direction of land use and land cover elements can be explained by temporal shifting of mean centers of each LULC element. Plotting of the weighted mean centers (X_w, Y_w) of each LULC element on a temporal basis has seemed useful for tracking changes in the distribution or comparing the distributions of each feature. The weighted mean centers have been calculated as [6,7,8]:

$$Xw=\frac{{{\displaystyle \sum }}_{i=1}^{n}wixi}{{{\displaystyle \sum }}_{i=1}^{n}wi}$$

(2)

$$Yw=\frac{{{\displaystyle \sum }}_{i=1}^{n}wiyi}{{{\displaystyle \sum }}_{i=1}^{n}wi}$$

(3)

where w_i is the weight of i^th element and x_i, y_i are the x and y coordinate of the i^th element.

Besides this, the rate of change of LULC elements is an important indicator that expresses the pace at which the LULC elements change either by expansion or as being encroached. The higher the positive value, the higher the pace to be expanded and the higher the negative value, the higher the areal loss as being encroached.

On the other hand, the value of it as “0” or within 0 to 1 is indicating the stability or near to stable state. In this present study, the rate of change in area for each land use/cover class is calculated using the following formula [8,9]:

$$C_i=\frac{\left(W_{ei}-W_{bi}\right)}{W_{bi}}\times \left(\frac{1}{t}\right)\times 100,$$

(4)

where C_i is the rate of change of impervious surface during the study phase, W_ei and W_bi are impervious surface area or percentage at the end and beginning of the study phase, respectively, t is the length of the study phase measured in years.

2.2.3. Calculation of Change Index

Such land use and land cover reconfiguration enumeration have provided effective inputs to calculate the transferred index T₁ and transformed index T₂. Transferred index T₁ is the quantitative explanation of fractional area transferred from a specific LULC element to other LULC elements within the total studied span retaining fractional unchanged area. Whereas the transformed index (T₂) is quite a straightforward quantitative temporal approach. It is a quantitative expression of fractional area transformed into a specific LULC element from the other LULC elements within the studied span. Both indices can be calculated using Equations (2) and (3) [4,5,6,7]:

$$T_1=\frac{{{\displaystyle \sum }}_{i=1}^n\left(A-a_i\right)}{{{\displaystyle \sum }}_{i=1}^n{a}_i}\frac{1}{a_j}\frac{1}{t},$$

(5)

where A is the total area, a_i is the fractional area transferred from the i^th LULC element in the beginning year to other LULC elements in the ending year, and j is the fractional unchanged LULC element in the ending year [4,5,6,7]:

$$T_2={{\displaystyle \sum }}_{i=1}^n\frac{a_e{}_i}{a_b{}_i}\frac{1}{t},$$

(6)

where the fractional area of different LULC elements $a_e{}_i$ has transformed into a specific LULC element $a_b{}_i$ from the beginning year to the ending year.

Finally, the so-called change index CI can be calculated as follows by the ratio between T₁ and T₂ [9,10,11]:

$$CI=\frac{T_1}{T_2}$$

(7)

2.2.4. Decision-Making Trial and Evaluation Laboratory (DEMATEL) Methodology

An extremely widespread methodology used for locating cause-and-effect chain elements in complex systems is represented by the so-called Decision-making trial and evaluation laboratory (hereinafter referred to as DEMATEL). In more detail, it deals with the assessment of interdependent interactions among components and identifying the significant ones using a visual structural model, and numerous studies on the use of DEMATEL have been undertaken in the past ten years.

This kind approach makes use of impact relation diagrams to identify the key components of a complex structural system as well as matrices to transform interdependency relationships into a cause-and-effect group. The formulating steps of DEMATEL are as follows:

i.

Finding the direct-relation (Average) matrix:

First, each respondent was asked to evaluate the direct influence between any two factors by an integer score from 0–4 representing 0 as no influence, 1 indicates low influence, 2 indicates medium influence, 3 indicates high influence, and 4 indicates very high influence.

The notation of x_ij indicates the degree to which the respondent believes factor i affects factor. For i = j, the diagonal elements are set to zero. For each respondent, a n × n non-negative matrix can be established as follows [14,15,16]:

$$x_k = x_{i,j}^k,$$

(8)

where k is the number of respondents with 1 ≤ k ≤ H, and n represents the number of factors. Thus, X¹, X², X³, …., X^H are the matrices from H respondents. To incorporate all opinions from H respondents, the average matrix (A) = [𝑎_𝑖𝑗] can be properly constructed as follows [14,15,16]:

$$\mathrm{A} = \left[a_{ij}\right] = \frac{1}{H}{{\displaystyle \sum }}_{k=1}^H{x}_{i,j}^k.$$

(9)

ii.

Calculation of the normalized initial direct-relation matrix:

Normalized initial direct-relation matrix D can be calculated by [14,15,16]:

$$D = mA,$$

(10)

where m = $\min \left(\frac{1}{ma{x}_i{{\displaystyle \sum }}_{j=1}^n{a}_{ij}} , \frac{1}{ma{x}_j{{\displaystyle \sum }}_{j=1}^n{a}_{ij}}\right)$ $\mathrm{with} i,j\in \left\{1,2,\dots ,n\right\}$, and element belonging to matrix D falls between zero and one.

iii.

Calculation of total relation matrix:

The total-influence matrix T is then obtained by utilizing Equation (11), in which, I is an n × n identity matrix. The $t_{ij}$ element summarizes the indirect effects that factor i had on factor j, and then the T matrix reflects the overall relationship between each pair of system factors [17,18,19,20]:

$$\begin{array}{c}T = \underset{m\to \infty }{\lim }\left(D+D2+\dots +D^m\right)\\ = {{\displaystyle \sum }}_{m=1}^{\infty }\mathrm{Di},\end{array}$$

(11)

where ${{\displaystyle \sum }}_{m=1}^{\infty }{D}^i$

= $D^1+D^2+\dots +D^m$

= $D\left(I+D^1+D^2+\dots +D^{m-1}\right)$

= $D{\left(I-D\right)}^{-1}\left(I-D\right)\left(I+D^1+D^2+\dots +D^{m-1}\right)$

= $D{\left(I-D\right)}^{-1}{\left(1-D\right)}^m$

$$T=D{\left(I-D\right)}^{-1}$$

iv.

Calculation of Influential relational map (IRM):

At this step, the vectors R and C, representing the sum of the rows and the sum of the columns from the total-influence matrix T, are defined by the following formulas [17,18,19,20]:

$$R={\left[r_i\right]}_{n\times 1}={\left[{{\displaystyle \sum }}_{j=1}^n{t}_{ij}\right]}_{n\times 1}$$

(12)

$$C={\left[c_j\right]}_{1\times n}={\left[{{\displaystyle \sum }}_{i=1}^n{t}_{ij}\right]}_{n\times 1}$$

(13)

Here ${\left[r_i\right]}_{n\times 1}$ demonstrates the total effects, both direct and indirect, given by criterion i to the other criteria j = 1, 2,…, n. Similarly,${\left[r_i\right]}_{n\times 1}$ represents total effects, direct and indirect, received by criterion j from the other criteria i = 1, 2,…, n. The sum (R + C) shows the total effects given and received by factor i. That is, (R + C) indicates the degree of importance that factor i plays in the entire system. On the contrary, the difference (R − C) depicts the net effect that factor i contributes to the system. Specifically, if (R − C) is positive, factor i is a net cause and if (R − C) is negative, factor i is a net effect.

In this present study case, DEMATEL has been used to enumerate the cause-and-effect relation between the LULC elements, which leads to the land use/cover reconfiguration.

2.2.5. Jaccard Similarity Index

The Jaccard index, also known as the Jaccard similarity coefficient, is a statistic that is used to determine how similar sample sets are (A and B). The size of the intersection divided by the size of the union of the sample sets is formally defined as the size of the intersection divided by the size of the union of the sample sets, by emphasizing the similarity between finite sample sets. The mathematical representation of the index is written as [21,22]:

$$J\left(A,B\right)= \frac{\mid A{\displaystyle \cap }B\mid }{\mid A{\displaystyle \cup }B\mid }= \frac{\mid A{\displaystyle \cap }B\mid }{\mid A\mid +\mid B\mid -\mid A{\displaystyle \cup }B\mid }.$$

(14)

The Jaccard index is essentially the number in both sets, divided by the number in either set, multiplied by 100. This will produce a percentage measurement of similarity between the two sample sets.

2.2.6. Adherence Index

The Adherence index is a metric that shows how much of each land use/cover category had not transitioned to another LULC category. The percent consistency of any land use/cover element after the transition is represented by this index. It is the ratio of the land use/cover area that had not changed $A_u{}_c$ to the mean of the area at the beginning A_b and ending year A_e associated with that specific land use/cover element [22,23]:

$$\mathrm{Adherence} \mathrm{Index} = \frac{A_u{}_c}{\frac{Ab+Ae}{2}}$$

(15)

3. Results and Discussion

3.1. Status of Land Use/Cover Classification

The classified LULC maps of the study area corresponding to 1989, 1999, 2009, and 2019 are illustrated in Figure 2.

The achieved classification accuracies were respectively 98%, 98%, 95%, and 95%, and the overall kappa statistics were 0.98, 0.96, 0.96, and 0.91 for the classification for 1989, 1999, 2009, and 2019 images, respectively. In the land use/cover classifications for the examined four years, user’s and producer’s accuracies of individual classes were extremely high, ranging between 81% and 99%, and 90% and 99% respectively (see

Table 2. Contingency matrix for land use and land cover, 1989.

Data	Urban Built-Up	Water Body/River	Vegetation	Bare Land	Fallow Land	Home Stead with Plantation	Agricultural Land	Row Total
Urban built-up	5314	35	0	0	415	0	0	5770
Water body/river	2	4102	0	0	0	0	0	4104
Vegetation	0	1	2883	0	0	60	2	2046
Bare land	0	0	0	2152	39	3	0	2194
Fallow land	684	13	0	236	2139	8	0	3080
Home stead with plantation	0	0	93	0	36	2233	1	2363
Agricultural land	0	0	16	10	0	0	576	602
Column Total	6000	4151	2992	2398	2635	2304	579	21,059

,

Table 3. Contingency matrix for land use and land cover, 1999.

Data	Urban Built-Up	Water Body/River	Vegetation	Bare Land	Fallow Land	Home Stead with Plantation	Agricultural Land	Row Total
Urban built-up	5429	40	0	0	415	0	0	5884
Water body/river	2	4602	0	0	0	0	0	4604
Vegetation	0	1	2083	0	0	60	2	2146
Bare land	0	0	0	2172	39	3	0	2194
Fallow land	684	13	0	236	2639	8	0	3580
Home stead with plantation	0	0	93	0	36	2243	1	2373
Agricultural land	0	0	16	10	0	0	676	702
Column Total	6000	4151	2992	2398	2635	2304	579	21,483

,

Table 4. Contingency matrix for land use and land cover, 2009.

Data	Urban Built-Up	Water Body/River	Vegetation	Bare Land	Fallow Land	Home Stead with Plantation	Agricultural Land	Row Total
Urban built-up	5731	40	0	0	415	0	0	6186
Water body/river	2	4611	0	0	0	0	0	4613
Vegetation	0	1	1183	0	0	60	2	1246
Bare land	0	0	0	2872	39	3	0	2914
Fallow land	684	13	0	236	2639	8	0	3580
Home stead with plantation	0	0	93	0	36	2143	1	2273
Agricultural land	0	0	16	10	0	0	684	710
Column Total	6000	4151	2992	2398	2635	2304	579	21,522

and

Table 5. Contingency matrix for land use and land cover, 2019.

Data	Urban Built-Up	Water Body/River	Vegetation	Bare Land	Fallow Land	Home Stead with Plantation	Agricultural Land	Row Total
Urban built-up	6531	49	0	0	410	0	0	6990
Water body/river	2	4611	0	0	0	0	0	4613
Vegetation	0	1	1083	0	0	60	2	1146
Bare land	0	0	0	2072	39	3	0	2114
Fallow land	684	13	0	236	2645	8	0	3586
Home stead with plantation	0	0	93	0	36	1943	1	2073
Agricultural land	0	0	16	10	0	0	884	910
Column Total	7217	4674	1192	2318	3130	2014	887	21,432

).

According to Li et al. [45] among others, an accurate classification should have an overall classification accuracy of 92% and kappa statistics above 0.9, which were successfully achieved in the present study.

3.2. Land Use/Cover Scenario

The classification area statistics are summarized in

Table 6. Temporal areal account on LULC elements.

Class Name	Area in km2	Area in Percentage	Area in km2	Area in Percentage	Area in km2	Area in Percentage	Area in km2	Area in Percentage
	1989 (km2)	1989 (%)	1999 (km2)	1999 (%)	2009 (km2)	2009 (%)	2019 (km2)	2019 (%)
Agricultural land	104.40	5.820	234.68	13.08	311.80	17.38	189.20	10.54
Bare land	388.00	21.62	194.70	10.85	137.76	7.68	206.66	11.52
Urban built-up area	362.71	20.21	473.36	26.38	493.85	27.53	539.10	30.05
Fallow land	62.88	3.50	101.54	5.66	76.45	4.26	124.74	6.95
Vegetation	182.66	10.18	187.50	10.45	224.30	12.50	174.66	9.73
Homestead with plantation	574.89	32.04	477.75	26.632	370.54	20.65	351.18	19.57
Water bodies	118.33	6.59	124.33	6.93	179.15	9.98	208.32	11.61

. The classified areas were measured by multiplying the count of pixels per class with a spatial resolution of remote sensing data (i.e., 30 m), in which the pixel counts were determined after applying post-classification analysis.

The data reported in Table 6 reveal that in 1989, about 5.82% (104.40 km²) area of the KMDA region was under agricultural land, 21.63% (388 km²) under bare land, 20.22% (362.71 km²) under built-up area, 3.51% (62.88 km²) under fallow land, 10.18% (182.66 km²) under vegetation, 32.05% (572.89 km²) under homestead with plantation, and 6.60% (118.33 km²) under water body. During 1999 the area under these land categories was found about 13.08% (234.68 km²) under agricultural land, 10.85% (194.70 km²) under bare land, 26.39% (473.36 km²) under built-up area, 5.66% (101.54 km²) under fallow land, 10.45% (187.50 km²) under vegetation, 26.63% (477.75 km²) under homestead with plantation, and 6.93% (124.33 km²) under water body. During 2009, the area under these land categories was found about 17.38% (311.80 km²) under agricultural land,7.68% (137.76 km²) under bare land, 27.53% (493.85 km²) under built-up area, 4.26% (76.45 km²) under fallow land, 12.50% (224.30 km²) under vegetation, 20.66% (370.54 km²) under homestead with plantation, and 9.99% (179.15 km²) under water body. Finally, during 2019 the area under these land categories was found about 10.55% (189.20 km²) under agricultural land, 11.52% (206.66 km²) under bare land, 30.05% (539.10 km²) under built-up area, 6.95% (124.74 km²) under fallow land, 9.74% (174.66 km²) under vegetation, 19.58% (351.18 km²) under homestead with plantation, and 11.61% (208.32 km²) under water body as shown in Figure 3.

3.3. Land Use-Land Cover Changes and Reconfiguration

3.3.1. Temporal Shifting of Mean Centers of Each LULC Element

The post-classification comparison method has been applied in this study, which is the most adopted approach to detect the spatial as well as the areal change of LULC elements [46,47,48,49,50,51,52,53].

The plot of the weighted mean centers of each LULC element on a temporal basis has seemed useful for tracking changes in the distribution or comparing the distributions of each feature, as shown in Figure 4.

Figure 4 shows that within 30 years of span the agricultural land, bare land, urban built-up, fallow land, vegetation, homestead with plantation and water bodies have shifted by 10.8 km north-westward, 6.0 km northward, 5.64 km southward, 13.47 km south-westward, 7.32 km northward, 1.51 km eastward and 3.22 km south-eastward, respectively, as shown in

Table 7. Shifting distances (expressed in km) and directions of each LULC element.

LULC Elements	1989 to 1999	1999 to 2009	2009 to 2019	1989 to 2019
Agricultural land	9.67 (NW)	1.48 (NE)	1.64 (W)	10.8 (NW)
Bare land	2.63 (E)	3.49 (NW)	3.55 (NE)	6.00 (N)
Urban built-up area	2.36 (SW)	2.41 (S)	1.07 (SW)	5.64 (S)
Fallow land	7.13 (S)	12.29 (SE)	9.44 (W)	13.47 (SW)
Vegetation	4.72 (NE)	3.41 (NW)	4.11 (NE)	7.32 (N)
Homestead with plantation	1.58 (S)	3.04 (NE)	1.57 (S)	1.51 (E)
Water bodies	0.23 (NE)	1.99 (N)	4.21 (S)	3.22 (SE)

.

Within this span, the fallow land has shifted by the longest distance (expressed in Table 4 in km) followed by agricultural land; minimum shifting distance is seen for the homestead with plantation. Apart from the overall scenario, if the total span is studied under three quarters, (1989 to 1999, 1999 to 2000, and 2009 to 2019) in the first phase, maximum transition has happened with agricultural land so that it has shifted by maximum distance (9.67 km northward) and minimum shifting distance is seen for water (0.23 km north-eastward direction). Whereas in the rest of the two spans, maximum shifting is seen for the fallow land (12.29 km south-eastward and 9.44 km westward in the second and last phase respectively), though in the last phase it is seen to be decreased and the shifting distance has seen to be minimized also. The analysis of LULC areal transition focuses on the fact that the mean centers shifting distance is seen to be minimized gradually, especially for the agricultural land and the fallow land, but the transition effects are seen on bare land, HSP, vegetation, and water in the span of 1999 to 2009 and vegetation and water in the span of 2009 to 2019 as their shifting distance has increased [54,55,56,57,58,59,60]. However, a noteworthy fact can be established from the explanation that within the 30 years of span, such mean enters of each LULC element have shifted within 13.5 km from their prior locations. Therefore, it can be concluded that maximum LULC dynamicity has been experienced by the 13.5 km area under KMDA. Among all these mean centers shifting about the land use /cover dynamicity, there are some similarities such as mean center shifting viz LULC dynamicity. As shown in

Table 8. Jaccard similarity index values per LULC element.

	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
Agricultural land	1.000	0.558	0.523	0.565	0.746	0.200	0.234
Bare land		1.000	0.874	0.483	0.928	0.537	0.721
Urban built-up area			1.000	0.345	0.777	0.537	0.578
Fallow land				1.000	0.585	0.207	0.276
Vegetation					1.000	0.413	0.556
Homestead with plantation						1.000	0.583
Water bodies							1.000

, implementing Jaccard Similarity Index similar nature in mean center shifting of respective LULC elements has been established. Here the similarity index value greater than or equal to 0.5 has been considered significant.

From Table 8 it is seen that the vegetation cover of the study area is the most dynamic of all and that is quite like the land use/cover such as agricultural land, built-up areas, bare land, and fallow land. The vegetation cover dynamics are the most alike the bare land (0.928) followed by built-up areas (0.77), agricultural land (0.74), and fallow land (0.58) respectively. On the other hand, strong similarities are seen between built-up areas and bare land (0.87) and between the water body and bare land (0.72).

Irrespective of aforesaid LULC elements, a near similarity in dynamicity can also be seen between the rest of the LULC elements. Agricultural land dynamics have a nearness with bare land (0.55), built-up (0.52), and fallow land (0.56) dynamicity. Bare land dynamics are found to be nearly similar to homestead with plantation (0.53). Built-up area dynamics have near similarity with the dynamics of homestead with plantation and water body by the similarity index value of 0.53 and 0.57 respectively.

The dynamics of both the vegetation area and homestead with plantation have the nearness to the water body dynamics by the index value of 0.55 and 0.58, respectively.

3.3.2. Study on the Annual Rate of Areal Change and Areal Loss and Gain Temporally

It has been revealed from

Table 9. Rate of changes (in percentage, %) in LULC areas over the time.

Class Name	1989 to 1999	1999 to 2009	2009 to 2019
Agricultural land	13.865	3.651	−4.369
Bare land	−5.535	−3.249	5.557
Urban built-up area	3.390	0.481	1.018
Fallow land	6.830	−2.745	7.017
Vegetation	0.295	2.181	−2.459
Homestead with plantation	−1.878	−2.493	−0.580
Water bodies	0.564	4.899	1.809

that within the overall span of 30 years, the maximum areal increase has been observed for the fallow land at the rate of 3.64% followed by agricultural land, water bodies, and built-up area at the rate of 3%, 2.81%, and 1.80% respectively. Whereas maximum areal encroachment has been observed in the case of bare land at the rate of 1.73%, followed by Homestead with plantation and vegetation at the rate of 1.41% and 0.16% respectively [61,62,63,64,65,66,67].

Though throughout the study built-up and water body has shown consistency with respect to the areal expansion, the rate of expansion is not alike. In the span of 1989 to 1999 built-up area has been seen to be expanded at a maximum rate (3.39%). But such expansion has a slowdown in the span of 1999 to 2009 at 0.48%. A little bit of expansion has taken place in this phase (20.49 km² or 1.14%), though that has been seen to be increased at the rate of 1.01%. On the other hand, water bodies have been observed to be consistent to some extent throughout the three phases. The rates of areal increase accounted for 0.56%, 4.89%, and 1.80% respectively. Another consistency is seen for homestead with plantation with respect to the decreasing areal extent.

Throughout the three spans they have been seen to be decreased by the rate of 1.87%, 2.49%, and 0.58% respectively. However apart from the aforesaid LULC elements, from 1989 to 1999 maximum rate of areal expansion is seen for the agricultural land (13.86%), followed by fallow land (6.83%). Whereas vegetation areas have almost no change (the rate of change is 0.29%). Apart from those LULC elements, for agricultural land and vegetation, a decreasing rate of areal expansion is evident. In the case of both, despite having areal expansion at the rate of 3.65% and 2.18% in the span of 1999 to 2009, those are seen to be decreased by 4.36% and 2.49% in the span of 2009 to 2019 respectively. Irrespective of all the elements, bare land area has decreased at decreasing rate of 5.53% and 3.24% in the span of 1989 to 1999 and 1999 to 2009 respectively though that has increased at the rate of 5.55% in the span of 2009 to 2019. On the other hand, fallow land area is seen to be increased in the spans of 1989 to 1999 and 1999 to 2019 at the rates of 6.83% and 7.01% respectively followed by a decrease in the area at the rate of 2.74% in the span of 1999 to 2009. LULC areal dynamics is the function of areal loss and gain temporally. The gain is the expression of expansion of any LULC element at the cost of the others; whereas loss is the encroachment of any LULC element by others. The post-classification comparisons between temporal areal databases can provide an apparent insight regarding the areal losses and gains of different LULC elements. Such areal gains and losses are an effective outcome of the interaction between each LULC element. A negative correlation between each can express the LULC areal growth and encroachment as well. As emerges from the detailed analysis of

Table 10. Correlation statistics between LULC temporal areas.

	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
Agricultural land	1
Bare land	−0.938	1
Urban built-up area	0.636	−0.856	1
Fallow land	0.147	−0.479	0.819	1
Vegetation	0.791	−0.544	0.114	−0.465	1
Homestead with plantation	−0.687	0.843	−0.948	−0.632	−0.312	1
Water bodies	0.394	−0.579	0.830	0.601	0.116	−0.927	1

, a high negative correlation can be observed between agricultural land and bare land (−0.93), built-up and bare land (−0.85), homestead with plantation and built-up (−0.98), and homestead with plantation and water bodies (−0.93). Besides the interaction between homestead with plantation and fallow land (−0.63), vegetation and bare land (−0.54) are worthy enough to be mentioned. For those cases, each LULC is experiencing a growth in the area at the cost of others.

Although it is evident from

Table 11. Account for areal gains (km²) per LULC elements over the time.

Gain	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
1989 to 1999	192.79	76.08	196.35	76.63	96.30	152.72	27.06
1999 to 2009	209.88	90.02	151.85	56.50	189.56	147.46	73.32
2009 to 2019	104.44	174.15	174.15	107.80	204.25	165.22	78.88
1989 to 2019	171.78	128.94	290.01	115.06	151.44	180.73	112.58

and

Table 12. Account for areal loss (km²) per LULC elements over the time.

Loss	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
1989 to 1999	62.51	269.38	85.70	37.98	91.46	249.87	21.06
1999 to 2009	132.76	146.96	131.35	81.59	152.76	254.66	18.49
2009 to 2019	227.05	96.32	128.90	59.52	150.46	223.61	49.72
1989 to 2019	86.98	310.27	113.61	53.20	159.44	404.44	22.59

that, within the 30 years of span, urban area has gained the maximum areas (290 km²) and homestead with plantation has experienced the maximum areal loss (404.44 km²), they did not remain same. In the period 1989 to 1999, maximum areal gain and loss were of built-up area (196.35 km²) and bare land (269.37 km²) respectively; from 1999 to 2009, maximum areal gain and loss were of agricultural land (209.88 km²) and homestead with plantation (254.66 km²) respectively; and within 2009 to 2019, the maximum areal gain and loss were that of the homestead with plantation (204.25 km²) and agricultural land (227.04 km²).

The transition matrix (

Table 13. Transition matrix (1989 to 2019).

1989 to 2019	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
Agricultural land	17.424	13.124	21.951	12.045	9.320	18.961	11.579
Bare land	47.467	77.730	90.308	37.664	33.309	78.743	22.775
Built-up	11.561	29.020	249.099	19.927	9.969	30.773	12.357
Fallow land	9.719	16.569	12.024	9.681	3.308	8.273	3.310
Vegetation	39.747	12.06	20.759	8.911	23.222	40.349	37.608
Homestead with plantation	56.908	55.384	139.823	34.369	93.008	170.454	24.949
Water bodies	6.371	2.777	5.141	2.141	2.525	3.630	95.742

,

Table 14. Transition matrix (1989 to 1999).

1989 to 1999	Agricultural Land	Bare Land	Urban Built up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
Agricultural land	41.893	15.138	8.19	10.279	15.433	12.492	0.979
Bare land	88.299	118.623	58.199	45.847	30.202	44.591	2.237
Built-up	11.292	15.798	277.014	9.920	5.121	37.337	6.225
Fallow land	5.553	19.476	5.643	24.909	3.969	2.730	0.603
Vegetation	19.921	5.631	4.937	2.601	91.197	49.822	8.548
Homestead with plantation	65.169	18.956	114.575	5.479	37.218	325.03	8.469
Water bodies	2.558	1.078	4.805	2.505	4.358	5.749	97.274

,

Table 15. Transition matrix (1999 to 2009).

1999 to 2009	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
Agricultural land	101.923	31.567	30.131	9.998	29.328	25.560	6.180
Bare land	47.952	47.741	40.276	29.166	8.069	17.122	4.377
Built-up	20.912	5.875	342.010	5.180	14.865	76.991	7.530
Fallow land	19.784	32.307	21.280	19.957	1.900	3.430	2.886
Vegetation	68.279	10.793	8.565	7.174	34.745	21.441	36.505
Homestead with plantation	49.468	9.144	45.532	4.212	130.472	223.088	15.837
Water bodies	3.487	0.338	6.063	0.770	4.926	2.912	105.843

,

Table 16. Transition matrix (2009 to 2019).

2009 to 2019	Agricultural Land	Bare Land	Urban Built-up Area	Fallow Land	Vegetation	Homestead with Plantation	Water Bodies
Agricultural land	84.757	39.714	37.151	24.535	34.960	69.594	21.093
Bare land	19.670	41.446	22.511	19.274	5.684	21.827	7.352
Built-up	10.938	40.183	364.962	39.007	3.207	22.556	13.005
Vegetation	27.451	10.385	15.675	5.425	73.847	75.120	16.403
Fallow land	4.921	21.417	19.229	16.937	1.333	4.971	7.648
Homestead with plantation	27.281	46.509	69.269	16.047	51.118	146.936	13.383
Water bodies	14.182	7.011	10.310	3.516	4.516	10.183	129.440

) is another useful method for expressing the areal gain and loss of LULC elements, from which areal losses and gains per LULC element may be derived from the corresponding rows and columns values.

3.3.3. Study on Change Index Per LULC Elements

The post-classification LULC transition matrix is a commonly adopted approach to represent the LULC reconfiguration. This is basically a cross-tabulation matrix that expresses detailed “from-to” change class information. From the transition matrix (

Table 17. Land use and land cover transferred index and transformed index.

Land Use/Cover Category	Transferred			Transformed
	1989 to 1999	1999 to 2009	2009 to 2019	1989 to 1999	1999 to 2009	2009 to 2019
Agricultural land	0.021	0.015	0.008	0.034	0.058	0.115
Bare land	0.006	0.006	0.018	0.200	0.069	0.049
Built-up	0.007	0.037	0.006	0.045	0.047	0.097
Fallow land	0.036	0.032	0.040	0.015	0.034	0.027
Forest	0.012	0.005	0.011	0.042	0.077	0.056
Homestead with plantation	0.004	0.009	0.006	0.135	0.122	0.108
Water bodies	0.034	0.018	0.016	0.011	0.007	0.021

), three types of information can be provided as the most dominant land use element which is effective, the most affected category, and the amount of persistent area. As per the areal account of the LULC elements registered in the transition matrix (Table 17) for 30 years, the urban built-up area is seen to be most dominant at about 290.01 km² (16%) area of the total area (1793.89 km²), reconfigured as urban built-up. Whereas homestead with plantation is the most affected category as out of its total area, 70% (404.44 km²) has transformed into an urban built-up area. However, apart from all these conversions out of the whole, 35.76% (643.35 km²) area has shown persistence within the span of the study period.

The LULC reconfiguration is a two-fold process, one element is “transferring” to others and many elements are “transforming” to one. Both processes in combination imply a fact as “change”. Therefore, in the LULC reconfiguration phenomenon, the “from-to” change can be renamed as “transfer” and “transform” respectively. From Table 9, it is seen that the maximum areas of existing LULC elements in 1989 have transformed to urban built-up in 2019 and hence the “transformed index” (Table 17) of urban built-up area has been accounted as 0.083 which is the highest. Whereas the vegetation area has been observed to have undergone maximum areal encroachment by the other LULC elements especially homestead with plantation (40.35 km²). The transfer index (Table 10) of the vegetation has been accounted as 0.023. However, such a pattern of class conversion viz., class reconfiguration on a “from–to” basis has changed the LULC nature from 1989 to 2019 and the vegetation cover has experienced the maximum change as indicated by the highest index value of 1.540 and agricultural land use has the lowest index value as 0.065 (

Table 18. Land use and land cover change index and adherence index.

Land Use/Cover Category	Change Index			Adherence Index
	1989 to 1999	1999 to 2009	2009 to 2019	1989 to 1999	1999 to 2009	2009 to 2019
Agricultural land	0.634	0.249	0.075	24.709	37.301	33.835
Bare land	0.032	0.087	0.376	40.715	28.719	24.066
Built-up	0.149	0.789	0.062	66.265	70.720	70.663
Fallow land	2.383	0.965	1.455	30.297	22.423	16.836
Vegetation	0.287	0.059	0.192	49.275	16.875	37.019
Homestead with plantation	0.029	0.077	0.056	61.754	52.597	40.718
Water bodies	2.944	2.642	0.756	80.171	69.749	66.811

).

However very specifically, from 1989 to 1999 water bodies had the maximum change index (2.94) and homestead with plantation had the minimum (0.029), within 1999 to 2009 the maximum and minimum change index were associated with water bodies (2.64) and vegetation (0.05) respectively and within 1999 to 2009 the maximum and minimum change index were for fallow land (1.45) and homestead with plantation (0.05) respectively (Table 18).

3.3.4. Identification of Cause-Effect Chain among LULC Elements

The LULC change associated with the study area for the 30 years of the span is an obvious outcome of the interaction between LULC elements as revealed by the transition matrix and the correlation table. From the overall analysis, it has been noted that a particular LULC element is changing in response to others. Some are expanding at the cost of others causing definite areal gains and losses.

Therefore, Dematel method had been ideally applied to deduce the cause-and-effect chain. As per this method, the highest R_i + C_j value (

Table 19. Cause and effect chain calculated from DEMATEL model.

LULC	ri	cj	(ri + cj)	(ri − cj)
Agricultural land	0.564	0.606	1.169	−0.042
Bare land	0.976	0.464	1.44	0.513
Urban built-up	0.694	0.3	0.994	0.394
Fallow land	0.348	0.856	1.205	−0.508
vegetation	0.12	1.32	1.44	−1.2
Homestead with plantation	1.512	0.133	1.645	1.378
Water	0.248	0.782	1.03	−0.534

) of bare land is indicative of its higher degree of relationship with the other. On the other hand, the positive values of R_i − C_j for built-up and HSP is indicating those as the causal elements, whereas the negative values of the rest of the elements are indicating those as the effect.

The plotting of R_i − C_j against R_i + C_j (Figure 5) depicts the LULC positions as per cause-and-effect assignments.

Therefore, apart from the inter land use/cover element interaction, it can be particularized that population growth and related urban expansion are the basic driver behind such land use/cover change viz., reconfiguration. Over time, KMDA has seen continual built-up and mixed built-up (that has been defined as homestead with plantation) spread at the expense of non-built-up land covering. However, there was a significant difference in such growth between the KMDA-urban and KMDA-non-urban settlement (that has been pointed as homestead with plantation). As a result of the peripheral expansion in KMDA-non-urban areas, the urban built-up and mixed built-up areas increased significantly. Urban built-up cover grew at a positive rate, whereas mixed built-up cover grew at a negative rate within the KMDA region. This can be explained by the growth of census towns within the area under KMDA. According to the Census of India data, the number of Census Towns had considerably increased from 113 in the year of 1991 to 449 in the year of 2011. Moreover, the global increase in urban population from 1991 to 2011 (Census of India) is a piece of evident information concerning urban growth as the driver behind land use/cover reconfiguration. As per Census of India data from 1991 to 2011, the total urban population of KMDA had increased from 10,481,539 to 12,969,556 out of which in 38 municipalities that had increased from 500,242 in 1991 to 7,228,920 in 2011 and in 3 Municipal corporations that had increased from 5,481,297 in 1991 to 5,740,636 in 2011.

4. Conclusions

In regard to the various observations and existing literatures related to LULC study, it can be generalized that changes in LULC pattern is inevitable. In any urbanized area with the continual urbanization process and associated essential infrastructural development, the changing LULC is very common. From the overall study, it is evident that the LULC in the KMDA area have changed significantly over the last 30 years. According to relevant quantitative assessments on temporal LULC scenario, the increase in urban built-up area was leading. The expansion is notably centered on the fringes of cities. The growth rate in southern and northeastern Kolkata has been very high in recent decades. As a result, a substantial amount of previously undeveloped land is now used for residential purposes. Cities are generally expanding toward semi-urban areas due to a lack of open space in the core areas. As the consequences of this, vegetative areas, agricultural land, and bare land have reduced. However, such reductions in area may not be the effect of areal increase of any single LULC element. LULC change is an outcome of interactive response among existing LULC elements. Such interaction may include an impetus in strategic planning and management of land and environment as well.

Therefore, the quantitative measures as adopted in this present study to assess the LULC dynamics as the function of interaction among LULC elements could help in abstracting the LULC scenario of a region and in tracking the effective and effected elements which could extend a scope of decision-making to the researchers, planners, and managers to ensure sustainability.