Prediction of urban expansion by using land cover change detection approach.

Bangladesh has been experiencing rapid urban expansion over the last few decades, contributing much to the region's land cover transition into the urban area. The study aims to employ geospatial modeling techniques to investigate land cover scenarios in the Pabna municipality of Bangladesh. Therefore, the research examined Cellular Automata Markov and Multi-Layer Perceptron Markov models to detect land cover for 2023 and 2028. The study selected the Multi-Layer Perceptron Markov as the best fit model over Cellular Automata Markov based on the highest kappa value. The result reveals that urban area has increased from 3.39 to 8.79 km2 over 1998-2018. Urban expansion and its surrounding area are primarily occurring towards the northeast directions. However, the extent of urban build-up land will grow from 3.39 km2 in 1998 to 11.01 km2 in 2023 and 12.44 km2 in 2028. Moreover, the future land cover map delineated that the urban growth will expand in the northeast part of the study area. The scenario shown in this paper would assist urban planners in quantifying the urban growth under different land cover features and thus preparing proper strategic measures.

Introduction

Urbanization is a common trend worldwide (UN, 2018) that occurred by the concentration of population in urban areas, resulting in land transformation (Shao et al., 2021; Mosammam et al., 2017; Xu et al., 2019a; Cobbinah and Darkwak, 2016). UN (2014) reported that more than 50% world's population resides in urban areas, and this figure is forecasted to exceed 65% by 2050. This enormous population growth will add 2.5 billion more inhabitants in global urban areas, where 90% of urban expansion will occur in the less-developed territory (Fenta et al., 2017). The urbanization scenario is evident in Africa and Asia, and by 2050, these regions will have experienced about 90% urbanization (UN, 2016). The previous study reveals that the haphazard urban expansion is faster in developing countries than developed countries (Cohen, 2006; Grimm et al., 2008). The trend of urbanization is common in Africa, resulting from population agglomeration (Andreasen et al., 2017). UN-Habitat (2016) recorded the urban population growth rate of 3% annually in Africa, which accounted for 400 million in 2010 and projected 1 billion by 20140. Similarly, Asian countries follow rapid urbanization for faster population growth in city centers (UN, 2014). The Asian population is estimated to increase from 42.5% in 2010 to 64.6% in 2050 (UN, 2014). Bangladesh has also followed a similar pattern of urbanization and population growth in urban areas (Hasan and Nazem, 2015). The UN (2014) reported that the current population of Bangladesh in urban areas was 34%, which is forecasted to reach 56% by 2050. The previous survey has recorded that the city will be expanded to 10,712 km2, accounting for 7.25% of the country's land area (BBS, 2011). The high growth rate of the urban population accelerated faster urban expansion (Zhao, 2010). Shao et al. (2021) demonstrate that worldwide urban population growth is a prime factor responsible for the extreme rate of urban sprawl in global south cities. The previous study also indicates that the increasing urban population creates pressure for service facilities such as housing, transport, health, education, recreation, etc. to support their livelihood is contributing to urban sprawl (Shao et al., 2021; Fenta et al., 2017; Sumari et al., 2017; Tanveer et al., 2019; Ujoh et al., 2019).

Urbanization is considered a significant factor of land cover change (Ahmed et al., 2013; Kalnay and Cai, 2003; Chen et al., 2006). Hemani and Das (2016) examined that urbanization occurs in suburban areas by changing land cover features in Guwahati, India. Another report by Ahmed et al. (2013) found that urban areas will have almost half of the land of Dhaka city by 2019, where vegetation will be causing most among any other land cover classes. Hassan (2017) witnessed that build-up land of Khulna city corporation area (KCC), Bangladesh exceeded 2772 ha by 2014 while agriculture and vegetation cover reduced by 92% and 37% respectively during the study period of 1989–2000, and significantly contributed to residential development. A plethora of research has proved that the urban growth resulting from land use and land cover change has adverse effects on biodiversity, urban heat island, habitat fragmentation, and a threat of sustainable development (Pawe and Saiki, 2020; Griggs et al., 2014; Pawe and Saikia, 2017; Souza et al., 2016; Son and Thanh, 2017). Therefore, monitoring urban land cover contributes to assessing urban environmental conditions and managing natural and artificial resources (Li and Avissar, 1994). Many researchers have tried to identify the factors contributing to local and regional climate change caused by changing land cover features (Choudhury et al., 2019; Karakus, 2019; Ogunjobi et al., 2018). Research shows that the changing land cover influences urban climate (Jahan et al., 2021; Nagarajan and Basil, 2014; Grimm et al., 2008). Li and Wang (2019) reveals that the land cover change directly affects land surface temperature (LST) through changing land biophysical characteristics. Land cover changes (e.g., from vegetation to impenetrable covers, such as asphalt, rooftops, black-top) are the primary driver of LST change. This build-up area, known as urbanization, transforms natural land cover into artificial structures such as housing, commercial building, transport services, etc (Babalola and Akinsanola, 2016; Pu et al., 2006; Patra et al., 2018). These changes of land cover features affect air humidity correlated to atmospheric temperature (Ibrahim, 2017). Tewolde and Cabral (2011) mentioned that developing countries are more vulnerable to rapid urban expansion resulting from land cover changes than developed countries. People from developing countries are impacted most due to their limited capacity to cope with environmental and social consequences of urban expansion (Fenta et al., 2017; Cohen, 2006). Therefore, monitoring and detecting urban growth have become a vital concern for urban sustainability worldwide (Zhang et al., 2019a; Shao et al., 2021).

Earlier studies have investigated the urban growth and land cover change scenario for densely and fast-growing urban areas (Gazi et al., 2020; Kafy et al., 2020; Rahman et al., 2020; Roy et al., 2020; Trotter et al., 2017; Ahmed et al., 2013). Several studies have considered fixed land cover classes to detect land cover change (Hassan, 2017; Imran et al., 2021). Most of the previous research areas were experiencing rapid urbanization and almost homogenous land use characteristics within city boundaries (Pal and Ziaul, 2017; Imran et al., 2021; Fenta et al., 2017; Nzoiwu et al., 2017). This study focuses on assessing the urban expansion of the Pabna municipality due to mixed land use characteristics within the city boundary (Parvez and Islam, 2019) and the continuous changing trend of land cover classes in a suburban area (Abir and Saha, 2021). Pabna Municipality of Bangladesh, experiencing unprecedented urban growth and continuously changing land cover classes due to practical development (Parvez and Islam, 2019). However, the recent product of a public university and medical college with supporting structural setup adjacent to Pabna municipality acts as a pull factor of suburban development. However, the urban growth is prominent in the Pabna municipality area (Abir and Saha, 2021), the current trend of urbanization towards the suburbs of Pabna town in the last decades (Pabna Municipality, 2018). The chaotic population growth and urban expansion of Pabna municipality and its surrounding area give a new dimension to this research. This research has monitored the urban growth pattern of Pabna municipality in three historical periods for comparative land cover change analysis and stimulated the land cover maps based on selected land cover classes. Pabna municipality (Figure 1) has mixed land cover features like water body (1.50 km2), bare ground (2.15 km2), agriculture (4.57 km2), vegetation (4.46 km2), and urban area (3.07 km2) (Parvez and Islam, 2019).

Figure 1

Land use and land cover map of Pabna Municipality.

Abir and Saha (2021) witnessed that the built-up area of Pabna municipality has increased significantly during the study period of 1990–2020, where water body, bare land, and vegetation decreased by 52%, 52%, and 23%, respectively. None of these papers did examine the direction-wise urban growth pattern and predicted land cover model for two separate stimulating years. In addition, the previous study did not read the urban growth by considering agricultural land and bare surface together in the Pabna context to stimulate future land cover maps.

Different geospatial methods and techniques help to monitor urban growth and various land cover features (Huang et al., 2018). Remote Sensing (RS) and Geographic Information System (GIS) have become two commonly used geospatial tools to detect land cover change (Jahan et al., 2021; Bhuiyan et al., 2020b; Hegazy and Kaloop, 2015). Landsat satellite images are popular secondary data sources applied for change detection analysis (Nagne et al., 2018). Remote sensing and GIS techniques help measure land cover transformation through satellite images (Shao et al., 2021). Though the percentage of clouds in an image affects data interpretation results (Ahmed et al., 2013), previous studies mentioned that the land cover change detection through remote sensing is accurate and effective (Shao et al., 2021; Fu et al., 2013; Ligate et al., 2018; Nzunda and Midtgaard, 2019). Epstein et al. (2002) demonstrated that the quantitative analysis of urban expansion rate through remote sensing and GIS is a cost-effective and time-saving method due to the availability of good quality satellite data. Researchers have successfully used supervised and unsupervised classification of remote sensing methods to measure urban expansion at various spatiotemporal scales (Bhatta 2009; Mundia and Murayama 2010; Tewolde and Cabral 2011; Sharma et al., 2012; Schneider and Woodcock, 2008). This study used a supervised classification method under five land cover classes for analyzing urban expansion, similar to the previous research. In this regard, specific objectives of the research include (a) monitoring the rate of urban expansion between the years 1998–2018, (b) analyzing land cover categories for the year 1998–2018 (c) forecasting future changes in land cover for the projected years of 2023 and 2028.

Method and materials

Study area

This research is considered the Pabna municipality of Bangladesh as the subject region. It was founded in 1876 and raised to municipality type 'A' in 1989 (Parvez and Islam, 2019). Pabna municipality is a local administrative unit located in the Pabna district of the Rajshahi division, which lies between the latitudes of 23°59′15′'N and 24°02′15′'N, and the longitudes of 89°12′00′'E and 89°15′45″ E (Figure 2). It lies 161 km northwest of the capital city of Bangladesh and 110 km east of the Rajshahi division (Abir and Saha, 2021). The study area covers 27.20 km2, with 15 wards and 46 mahallas under the jurisdiction. The total population of Pabna municipality is 133403, with a growth rate of 3.15% annually (DPHE, 2021). According to the DPHE (2021) report, the population density of Pabna municipality is 30 people per acre, which is higher than the other municipalities. The study region has a tropical monsoon climate, which means summers are hot and humid, while winters are frigid (DPHE, 2021). According to BBS (2011), the annual average temperature of Pabna municipality ranges from 23.3 °C to 32.2 °C, with a yearly average rainfall of 1656 mm. The humidity level was 76.80% in April and about 85% in July (BBS, 2011). The study area's average elevation is 43 feet, and it's primarily flat plain (Parvez and Islam, 2019). Although the residents preferred a considerable portion of the northern region for residential purposes, the vast majority of the land area remains undeveloped (Abir and Saha, 2021). Citizens have employed the southern part mainly for residential and commercial purposes, and there are two industrial zones in the southwest and northeast outskirts of the site.

Figure 2

Study area map.

Data collection

The research used Landsat satellite images (Figure 3) collected from the U. S. Geological Survey (USGS) public domain website. As a result, we acquired three-time satellite image data for 1998, 2008, and 2018 to examine land cover change detection and future prediction. The authors downloaded images with a pixel size of 30m from path code 138 and row code 43. Table 1 shows the comprehensive information of acquired data.

Figure 3

Methodological flow chart.

Table 1

Attributes of Landsat data.

Year	Date acquired	Sensor	Number of bands	Pixel Size	Path and Row	Cloud %
1998	27/02/1998	TM	8	30m	138 & 43	0
2008	17/02/2008	TM	8	30m	138 & 43	0
2018	28/02/2018	OLI/TIRS	11	30m	138 & 43	0

Map projection has been made through Universal Transverse Mercator (UTM) within zone 45-N Datum World Geodetic System-84. Satellite images from 1998, 2008, and 2018 indicate the band (4, 3, 2), Red-Green-Blue (RGB). The Landsat 8 OLI images are 11 bands, with any three bands from the same sensor forming a false-color composite image. False-color composite image created using any three bands of the OLI sensor (USGS, 2016). RGB stands for the primary three colors of red, green, and blue, and therefore we have been utilized to create a false-color composite band (4, 3, 2). This false-color composite shows the urban area as cyan blue; vegetation is red; water body as blue; and bare land as brown.

Radiometric calibration

Firstly, radiometric calibration is applied to obtain the reflection value. Then, the ENVI software FLAASH method employed atmospheric correction techniques to correct wavelengths. A copy of raw images was produced, and considered DOS corrections were to analyze the impacts of the DOS algorithm on Landsat 8 OLI and Landsat 5 TM (Chavez, 1988). DOS algorithm helps to detect dark objects in a satellite image, having zero to minimal reflectance numbers in the pixel (Song et al., 2001). The minimum digital number (DN) value indicates the subtraction of atmospheric scattering effects within a satellite image. We applied ENVI 5.3 software to perform DOS correction, and it automatically generates corrected multispectral images.

Band composition

Building the signature detect common areas of the classified image using band composition. Figure 4 represents the band composition for the land cover map of 2018. Landsat 8 OLI true color composite (4, 3, 2) has shown land cover as actual color. False-color composite (5, 6, 4) represents the urban area as light blue, vegetative area as reddish-brown, and agricultural land as orange. On the other hand, false-color composite (7, 5, 3) has delineated the urban area as magenta, forest as olive to light green, bare land as pink, coniferous forest as light green, and water as blue. Another band of false-color composite (5, 4, 3) for Landsat 8 OLI represents urban as cyan blue, vegetation as dark red, agricultural land as pink to red, and the water body dark blue to black.

Figure 4

Landsat 8 OLI false and true color band composition (2018).

Figures 5 and 6 delineate the false and true color band composition for the land cover map of 1998 and 2008. Landsat 5 TM true color composite (3, 2, 1) represents land cover as actual cover. RGB stands for red, green, and blue and creates false-color composites in bands 4,3,2. (FCC). False-color band composition (4, 5, 3) for Landsat 5 TM delineates urban areas as light blue, vegetative areas as reddish-brown, agricultural land as the orange patch, and water bodies as blue. In contrast, the false-color band composite (7, 4, 2) shows the urban area as magenta, forest as olive to the bright green, coniferous forest as light green, bare land as pink, and sparsely vegetation orange-brown, and water as blue.

Figure 5

Landsat 5 OLI false and true color band composition (2008).

Figure 6

Landsat 5 OLI false and true color band composition (1998).

Building structure

Digitized polygon has measured the signature value through band composition of images. The signature value of these images has a DN value. The reflectance value of the pixels is known as signature development. Table 2 depicts the selected land cover categories considering specific characteristics.

Table 2

Land covers classes.

Land cover types	Descriptions
Urban	All residential, commercial, industrial, transportation, and settlements
Agricultural land	Cropland, cultivated land, urban agricultural land, low lying area, and marshy land
Vegetation	Trees, shrubland, semi-natural vegetation, deciduous and coniferous forest, inter-city recreational areas, parks and playground, grassland, vegetation area, garden
Water body	River, lakes, pond, permanent open water, cannel
Bare ground	The landfilling site, constructions site, development land, bare and exposed soils, open space

Image classification

The signature file is raw data for supervised image classification. Parametric and non-parametric rules are the two common approaches to the supervised classification method. In this study, we considered the maximum likelihood parametric rule to classify images. The statistical parameters covariance matrix of the pixels are considered for the parametric approach.

Accuracy assessment

Accuracy assessment is the final stage of the image classification process. The ground truth data evaluate and compare the classified data through accuracy assessment. However, it is not typical to assess each pixel ground truth data from 1998 to 2018. The images produced several reference points with the specific geographic coordinate. There are 250 references points have been generated for each classified image to perform an accuracy assessment. Both user and producer accuracy measured the overall accuracy of a classified image. Finally, the agreement of classified map is determined by Kappa coefficient whether the images are suitable for analysis or not.

User accuracy is a technique to measure the performance of a classified map based on field data by category (NOAA, 2018). It is calculated by dividing a class's correct pixels by the total pixel assigned to that class (Canada Centre for Remote Sensing, 2010).

Producer accuracy measures an analyst's performance in the image classification by category (NOAA, 2018). It shows the percentage of the corrected ground class. It is measured by the ratio of correct pixels and the reference pixels of a classified image (Canada Centre for Remote Sensing, 2010).

Kappa is a final accuracy measure between the classified and ground reference maps. The following formula measures the Kappa coefficient (Hagen, 2002).Here, TS = Total number of sample

TCS = Total correctly number of sample

Land cover change detection analysis

Monitoring land cover through remote sensing geospatial tools helps assess land cover change scenarios in the different periods. Thus, the researcher could understand the pattern of urbanization in the study region. Both pre-classification and post-classification methods measured land cover change. Post-classification is the commonly used method for detecting land cover change (Ahmed et al., 2013). Therefore, this research applied the post-classification method.

Land cover suitability analysis for future urban expansion

Suitability analysis helps stimulate land cover based on existing growth patterns and driving factors (Hopkins, 1977). Suitability analysis is the most appropriate way of land cover prediction according to purpose. The criterion of suitability varies according to land cover types (Malczewski, 2004). This research aims to predict future urban expansion based on current growth patterns. Different areas have different growth patterns with specific trends and directions (Malczewski, 2004). Therefore, the study considered past and present urban growth trends for suitability assessment.

Multi-Criteria Evaluation

Multi-Criteria Evaluation is a popular method of suitability analysis (Drobne and Lisec, 2009). To delineate suitable areas for different land cover requires defining factors (Figure 3). The variables are combined with importance weighting to indicate a suitability map (Drobne and Lisec, 2009). MCE must be considered constraints and factors for suitability assessment.

Constraints: Constraints are those variables that refer to the restricted region for development. Therefore, the value is either zero (not suitable) or one (suitable) (Eastman, 2012).

Factors: Factors are those variables that have a significant role in urban growth. These factors have different values means that diverse suitability. So, factors have continuous suitability values (Eastman, 2012).

Eastman (2012) mentioned that the most common techniques of MCE include Boolean Intersection, Weighted Linear Combination (WLC), and Ordered Weighted Averaging (OWA) (Eastman, 2012). In this research, we utilized WLC techniques for the feasibility of the study.

Weighted linear combination

WLC is a technique for the weighted overlay operation of criteria. The weight and criterion value generated the suitable area in the process of overlay operation. The following formula calculated the WLC (Eastman, 2012):Where, S = Suitability, W = Weight of factor i, X = Criterion score of factor i, C = Criterion score of constraint J, = Production.

After that, the factors were standardized within a numeric range and combined according to the weights. The combined value is overlaid with constraints for the final suitability map (Eastman, 2012). The standardization formula is:Where, R = Raw score

The fuzzy set membership approach is implemented for standardizing the factors. The fuzzy truth shows the membership in an inarticulately defined set, and the defined group does not have a defined boundary (Jiang and Eastman 2000). Fuzzy set functions vary from Sigmoidal, J-shaped, and linear (Eastman, 2009), and this research considered the linear fuzzy set function.

Model calibration and validation

This research has chosen Pattern-based models over Agent-based models (Agarwal et al., 2002), which is the aggregation of the CA Markov method and the multitude of MLP Markov methods (Balzter, 2000; Zalta, 2012; Basharin et al., 2004). The more justified method for the simulation of 2023 and 2028 is verified and considered according to findings. Here, the principal variables are the land cover data and Markov transition areas discovered by the Markov model. The transition region created suitability images, and the CA filter is the primary variable here (Purves and Pacala, 2008).

Model validation includes comparing the simulated and base maps, which is a crucial step as it generates numerical value to compare two maps in the case of predictive change modeling (Eastman, 2009; Vliet, 2009). Though it is challenging to anticipate changes in land cover using model validation (Ahmed, 2011), the best technique to justify a model is to match the prediction map (Pontious and Chen, 2006). Kappa statistics value validated the predicted map.

Data analysis and results

Image classification and land cover change detection

In 1998, the most considerable amount of land covered by vegetation and agriculture was 8.10 km2 (29.41%) and 9.39 km2 (34.10%), respectively, which reduced to 7.71 km2 (28%) and 8.62 km2 (31.31%) in 2008 (Table 3). Between 1998 and 2008, the urban area rose by 8.42% from 3.39 km2 in 1998 to 5.71 km2 in 2008. The classified map also revealed that between 1998 and 2008, the quantity of bare ground in the study region dropped by roughly 3.37%. Table 4 demonstrates that the vegetative land cover had declined considerably between 2008 and 2018. Nearly 11.21 % of the land was converted to the build-up area, resulting in a considerable increase in the size of the urban area.

Table 3

Land cover transfer matrix (1998–2008).

Land Cover Types	1998		2008		Change in the area (%)
	Area = (km2)	(%)	Area = (km2)	(%)
Urban	3.39	12.30	5.71	20.72	+8.42
Agricultural land	9.39	34.10	8.62	31.31	-2.79
Vegetation	8.10	29.41	7.71	28.00	-1.41
Water body	5.55	20.16	5.31	19.10	-1.06
Bare ground	1.11	4.04	0.18	0.67	-3.37
Total	27.54	100	27.54	100	0

Table 4

Land cover transfer matrix (2008–2018).

Land Cover Types	2008		2018		Change in the area (%)
	Area = (km2)	(%)	Area = (km2)	(%)
Urban	5.71	20.72	8.79	31.93	+11.21
Agricultural land	8.62	31.31	8.30	30.15	-1.16
Vegetation	7.71	28.00	5.85	21.25	-6.75
Water body	5.31	19.10	3.96	14.37	-4.73
Bare ground	0.18	0.67	0.63	2.29	+1.62
Total	27.54	100	27.54	100

The classified images of different years' land cover maps imply that the urban area of Pabna municipality increased over time, from 12.30% in 1998 to 20.72% in 2008 and 31.93% in 2018. However, there has been a significant decrease in vegetative land cover. The urban area has dynamically increased in the study area due to successive reductions of vegetative coverage. Agricultural land has dropped, from 34.10% in 1998 to 31.31% in 2008, and then exhibited a modest decrease of roughly 30.15% in 2018. The water body has reduced from 20.16% in 1998 to 19.10% in 2008 and has continued to degrade with 14.37% in 2018. There was a high percentage reduction of bare ground from 1.11 km2 to 0.63 km2 over 20 years (Table 5).

Table 5

Land cover transfer matrix (1998–2018).

Land Cover Types	1998		2018		Change in the area (%)
	Area = (km2)	(%)	Area = (km2)	(%)
Urban	3.39	12.30	8.79	31.93	+19.63
Agricultural land	9.39	34.10	8.30	30.15	-3.95
Vegetation	8.10	29.41	5.85	21.25	-8.16
Water body	5.55	20.16	3.96	14.37	-5.79
Bare ground	1.11	4.04	0.63	2.29	-1.75
Total	27.54	100	27.54	100	0

Image classification accuracy should be at least 85% (Anderson et al., 1976); however, more than 80% accuracy is acceptable (Carletta, 1996). The higher value of accuracy assessment means that the base map is most suitable for future expansion prediction. The overall accuracy measured 81. 68%, 82.12%, and 88.21% for 1998, 2008, and 2018 land cover maps, respectively. The cappa coefficient value calculated 88.21%, 83.21% and 81.68% for land cover maps of 2018, 2008 and 1998, respectively. Here, the Kappa coefficient value for each land cover category was also more than 80%. The accuracy level of this study was greater than the acceptable range, and the finding from the classified map is applicable for urban expansion prediction.

Direction of urban expansion

Different forms of land cover turned into urban areas over time. Pabna municipality has divided into eight spatial directions from the center of the Central Business District (CBD) area. The spatial direction approach follows geographical trends, namely NW-Northwest, N-North, NE-Northeast, E-East, SE-Southeast, S-South, SW-Southwest, and W-West. Figure 7 illustrates that the historical trend of urban expansion occurred in the northeast direction, indicating the fast-growing zone of Pabna municipality. Urban growth in the West, Northwest, Northeast, and East portion of the study area is comparatively lower than in the Northeast direction. In contrast, the Southeast, South, and Southwest portion of Pabna municipality has no historical trend of urban growth.

Figure 7

A radar chart of urban areas at different times and directions.

Model selection

Suitability map

Suitability analysis defines all factors and constraints. Constraints were delineated regions that are not suitable for future urban expansion. The study considered the factors based on the study area profile (Figure 3), with a different scale. Therefore, the standardization of the factors was the same for suitability analysis. This research used a fuzzy membership approach to standardize the factors and create distance images of each land cover class. The distance images used a simple Euclidean distance (meter) function, which measures the distance between each cell from the featured images. After that, the weight given to standardized images according to expert opinion. The aggregate MCE resulted in images using a weighted linear combination (WLC) method.

Urban area suitability map (2008)

Constraints map for urban suitability (Figure 8) shows the regions which are not suitable for future development. The water body constraint map has revealed the restricted area for further development. The urban built-up constraint map indicates the existing occupied area. So that further urban expansion is not possible because the models were not considering vertical growth.

Figure 8

Constraints map for urban suitability.

Figure 9 shows the factors fuzzy set distance images generated according to standardized values. The urban expansion was mainly in the periphery of the existing urban area. Distance from the current urban area was on a continuous scale but not in 0–255. So the standardization has been done by the monotonically decreasing liner function. That implies that the urban suitability value decreases when we go far from the urban area, and after 255 m, it is close to 0. Agriculture, vegetation, and bare ground distance images used monotonically increasing linear function from 0 to 255 m. Distance images of water bodies also followed monotonically increasing linear function. The area within 30 m is considered suitable, and then suitability value increase with increasing distance up to 255 m. The commercial site's monotonically increasing linear function indicates that the land is not ideal for urban expansion between 0 m and 500 m, but >500 m is the highest suitability for urban development. The distance (0 m–2000 m) for urban center implies the land is not suitable for urban expansion, whereas the value of more than 2000 m highlighted the high suitability for urban growth. The fuzzy set membership symmetric linear function has been applied in this research to determine the suitability of proximity to the highway. In this case, the distance between 0 – 50 m indicated the land is not suitable for urban expansion. The suitability decreases continuously from 50 m to 2000 m, and after 2 km, it is close to 0. Distance image of the secondary road followed symmetric linear function, where 0–10 m distance area is not suitable for urban growth. The suitability for the secondary road is constantly decreasing from 10m to 48000 m, and after 4.8 km, it is close to 0. Distance image of tertiary used using monotonically decreasing linear function, where the suitability decreases with increasing distance from urban areas (0–255 m). The distance image of the river developed following a monotonically increasing function. Suitability for the intermittent river within 30 m is limited. Beyond that, suitability increases with distance up to 1000 m, then their effect terminates. The distance elevation image developed through monotonically increasing function and slope distance image by the decreasing process.

Figure 9

Fuzzy set distance images of the urban area.

The experts have assigned the weight of the factors based on their relative importance. Here, thirteen factors are considered and weighted to generate an urban suitability map of 2018. The factor weight should range from 0 to 1, with 1 being the overall weight. The weight of the distant image of elevation (0.14), tertiary road (0.12), water body (0.11), secondary road (0.11), and highway (0.10) is prominent due to their appropriateness of urban suitability. The WLC approach combines the weights of the factors, resulting in the urban suitability map (Figure 10).

Figure 10

Urban suitability map (2008).

Agricultural land suitability map (2008)

Figure 11 shows the water body and urban constraints map for agricultural suitability. The red color denotes a favorable area for agricultural land, whereas the black indicates a restricted region for agricultural suitability. It is important to note that a considerable amount of land indicates suitability for further expansion.

Figure 11

Constraints of agricultural suitability.

To develop the fuzzy set distance images (Figure 12), the researchers considered three factors: agricultural land, bare ground, and vegetation. The images followed monotonically increasing linear function from 0 to 255 m. The appropriateness increases with distance, reaching a maximum of 255 m.

Figure 12

Fuzzy set distance images of agricultural land.

A weight of 75% is assigned to the distant image of agricultural land since it is the most ideal for agricultural suitability. The importance of the fuzzy set distance images of vegetation and the bare ground were 0.15 and 0.10, respectively. Figure 13 implies that the suburban region is most suitable for the agricultural land cover than city centers.

Figure 13

Agricultural land suitability map (2008).

Vegetation suitability map (2008)

Figure 14 depicts the water body and urban constraints map for vegetation suitability. The black color represents the inappropriate region for vegetation suitability. In contrast, the red color delineates the most suitable area for vegetative land cover.

Figure 14

Constraints of vegetation suitability.

The study has generated the fuzzy set distance images of agricultural land, bare ground, and vegetation using the monotonically increasing liner function (Figure 15). That means that as we get further away from the city, the urban appropriateness rating increases, and close to 255 m, it is maximum.

Figure 15

Fuzzy set distance images of bare ground.

The fuzzy set distance image of vegetation obtained an 80% weighting, whereas agricultural land and barren ground obtained just 0.10 each. Figure 16 shows the vegetation suitability map after integrating the weight of the factors. The vegetation suitability map implies that the east portion of the study region is most suitable for vegetative land cover.

Figure 16

Vegetation suitability map (2008).

Water body suitability map (2008)

The study considered the non-water body class to develop the water body suitability map of 2008. The water body constraints map (Figure 17) displays the existing restricted land cover class that will not consider when assessing suitability. The suitability map (Figure 18) shows that the west portion of the study region is suitable for water bodies. The northeastern part of the city, on the other hand, is unsuitable for water body land cover.

Figure 17

Constraints of water body suitability.

Figure 18

Water body suitability map (2008).

Bare ground suitability map (2008)

The researchers have displayed the restricted region in the water body and urban restrictions map (Figure 19). The majority of the land in urban centers is occupied by build-up land. As a result, the sub-urban area suggests appropriateness.

Figure 19

Constraints of bare ground suitability.

The authors used the land cover map of agricultural, bare ground and vegetative land to create fuzzy set distance images of bare ground (Figure 20). A monotonically increasing liner function is applied to standardize the data. The rate of appropriateness rises with distance, peaking at 255 m.

Figure 20

Fuzzy Set distance images of bare ground.

The expert has weighted the factors to create the bare ground suitability map. The weight of the distance image of agricultural land, vegetation and the bare ground was 0.23, 0.22 and 0.55, respectively. Figure 21 delineates the scattered areas of suitability for bare ground.

Figure 21

Bare ground suitability map (2008).

CA Markov simulated final land cover image (2018)

At first, five land cover suitability maps merged into a single file through IDRISI Selva 17.02 software package. The software tools considered suitability maps of 2008 to stimulate land cover image. The decision support wizard of software tools created Markov transition area from 1998 and 2008 transition suitability images. The CA Markov process consists of a set of iterations for either six months or one year. This study utilized one year period with ten iterations. This research used the 5 × 5 CA contiguity filter and base map of Pabna municipality (2008) to land cover prediction map of 2018. Figure 22 shows the CA Markov projected land cover map of 2018 after integrating the five land cover classes.

Figure 22

CA Markov projected land cover map (2018).

MLP Markov potential land cover map (2018)

The weights of the transition in the probability matrix indicate the area for future urban expansion. The MLP potential land cover map of 2018 (Figure 23) has been stimulated based on the Markov chain analysis.

Figure 23

MLP Markov projected land cover map (2018).

Comparison of models

Two alternative models were used to investigate the trends of land cover. Firstly, it is needed to select the best appropriate model for the research between CA Markov and MLP Markov. The maximum kappa value implies the best fit model according to the assumption. The Kappa value for CA Markov and MLP Markov were 0.84 and 0.83, respectively. The MLP Markov model has a greater numerical value of Kappa coefficients between the two models. As a result, this research applied MLP Markov to monitor land cover change in the projected years.

Prediction of urban expansion

It is feasible to determine the transition weights in the set of possibilities (Table 6 for 2023 and Table 7 for 2028) of the Markov chain for prediction through MLP Markov analysis. The MLP Markov generated the final prediction map of 2023 and 2028.

Table 6

Transition probabilities grid for Markov chain (2023).

	Urban	Agricultural Land	Vegetation	Water Body	Bare Ground
Urban	0.80	0.10	0.07	0.01	0.02
Agricultural Land	0.18	0.55	0.17	0.07	0.03
Vegetation	0.29	0.24	0.40	0.06	0.01
Water Body	0.06	0.20	0.18	0.53	0.03
Bare Ground	0.22	0.38	0.01	0.00	0.39

Table 7

Transition probabilities grid for Markov chain (2028).

	Urban	Agricultural Land	Vegetation	Water Body	Bare Ground
Urban	0.87	0.06	0.05	0.01	0.01
Agricultural Land	0.10	0.67	0.15	0.05	0.03
Vegetation	0.22	0.21	0.52	0.05	0.00
Water Body	0.00	0.14	0.16	0.68	0.02
Bare Ground	0.13	0.32	0.00	0.00	0.55

Analysis of prediction map

According to a quantitative examination of the projection map, the urban area will cover 12.44 km2 (almost 45%) of the entire area in 2028 (Table 8). Gains in urban land cover types are prominent while other land cover types are continuously decreasing. Almost half of the vegetative area will reduce, from 8.10 km2 in 1998 to 4.13 km2 in 2028. Moreover, bare ground, water bodies, and agricultural lands will follow a similar decreasing trend, contributing to urban expansion.

Table 8

Change in area (km2) (1998–2028).

Year
Types	Urban (km2)	Agricultural Land (km2)	Vegetation (km2)	Water Body (km2)	Bare Ground (km2)
1998	3.39	9.39	8.10	5.55	1.11
2008	5.71	8.62	7.71	5.31	0.18
2018	8.79	8.30	5.85	3.96	0.63
2023	11.01	7.44	4.58	3.95	0.55
2028	12.44	6.78	4.13	3.70	0.49

The urban area map (Figures 24 and 25) illustrates the Pabna municipality's predicted urban expansion (2023 and 2028).

Figure 24

Predicted land cover map of 2023.

Figure 25

Predicted land cover map of 2028.

Direction of future urban expansion

Pabna municipality and its surrounding area had seen rapid urbanization, resulting in considerable changes in land cover over the 30 years. The research period's spatial trajectory of urban expansion suggests that the North, Northeast, and East were the key routes of urban growth (Figure 26). The extension of urban areas in the West and Northwest was limited. The low-lying agricultural terrain controlled urban development in the South zone. The thematic map measured the rate of urban expansion over the last 30 years in the paths of 12.70 km2 in the Northeast and 8.09 km2 in the North (Table 9).

Figure 26

The spider diagram showing spatial patterns of urban expansion (2).

Table 9

Urban area based on direction (km2).

	N	NE	E	SE	S	SW	W	NW
2028	2.47	3.72	1.68	0.12	0.03	0.13	2.22	2.21
2023	2.16	3.37	1.51	0.11	0.03	0.12	1.89	1.87
2018	1.68	2.84	1.30	0.08	0.01	0.06	1.26	1.46
2008	1.11	1.77	0.78	0.05	0.02	0.04	0.68	1.04
1998	0.67	1.00	0.37	0.01	0.01	0.01	0.50	0.51
Total	8.09	12.70	5.64	0.37	0.10	0.36	6.55	7.09

NW-Northwest; N-North; NE-Northeast; E-East; SE-Southeast; S-South; SW-Southwest; W-West.

Discussion

Rapid urbanization is a significant cause of changing land cover features (Rasel and Parvez, 2021). The trend of urbanization enhances urban growth by transforming vegetative and agricultural land into urban build-up areas (Pal and Ziaul, 2017). Nzoiwu et al. (2017) found that the build-up area increased over time in Awka town of Nigeria, rising from 9.55 km2 in 1986 to 21.45 km2 in 2015. Consequently, vegetation decreased from 33.69 km2 in 1986 to 21.41 km2 in 2015, losing around 12.29 km2 (Nzoiwu et al., 2017). Hassan (2017) delineated that the cumulative built-up area in the five major cities of Bangladesh increased by 468% from 2,356 ha in 1973 to 13,435 ha in 2014, while vegetation cover and agricultural land diminished by 27.77% and 61.91%, correspondingly, during the same time. Another research by Ahmed et al. (2013) found that approximately 50% of the land cover of Dhaka city will convert into urban areas, where vegetation is affected most than other land cover classes. This research found that there have been noticeable changes in urban and vegetative land cover features, which is similar to previous studies (Hegazy and Kaloop, 2015; Bose and Chowdhury, 2020).

In Pabna Municipality, the urban area was 3.39 km2 in 1998, and it had increased 8.79 km2 in 2018. The vegetation land cover decreased from 8.10 km2 in 1998 to 5.85 km2 in 2018. A noticeable amount of water bodies has been transformed from 5.55 km2 in 1998 to 3.96 km2 in 2018. The northeastern part of the Pabna Municipality measured the central path of urban growth. All other land use categories, on the other hand, dropped, with the vegetation being the most affected. The change in land cover indicates that vegetation, agriculture, bare land, and water bodies have continued to shrink in size, contributing to urban growth. The predicted land cover map suggests that 11.01 km2 of land in 2023 and 12.44 km2 of land in 2028 will convert into urban land. Agricultural and vegetation land cover will reduce by 6.78 km2 and 4.13 km2 in 2028, respectively. Similarly, water bodies and bare land cover will reduce in size in the predicted year of 2023 and 2018. Vegetation, agriculture, bare ground, and water bodies have shrunk due to land use and land cover changes, leading to urban development. Evidence proved that migration from rural to urban areas forced suburban development in providing service facilities (e.g., housing, accessibility, etc.) for the migrated people (Hua and Ping, 2018). It is measured from the data that the urban areas of 2018 were nearly 2.5 times larger than in 1998, with vegetation losing the most. This substantial urban growth is linked with rising livelihood demand from migrated populations (Hassan, 2017). Imran et al. (2021) delineated that urban expansion of Dhaka city was uneven and occurring in the suburbs due to the availability of agricultural and vegetative land. This research observed similar findings, where urban growth occurs mainly in the North, Northeast, and Northwest suburbs. The suburb of Pabna municipality has agricultural, vegetative, and bare land covers. As a result, there is a need for structural development to accommodate rising population demand by converting various land use features into urban areas.

The approach of monitoring urban growth helps to monitor the impacts of land cover changes on urban environmental components by transforming vegetative, agricultural, and water bodies into build-up land. The study reveals that the Pabna municipality has increased over time. Shao et al. (2021) mention that the urban expansion caused by population growth has adversely affected the ecosystem due to loss of vegetation, water body, and land degradation. Likewise, urban growth exerts pressure for new construction and service facilities (Shao et al., 2021) that significantly raise the urban temperature by absorbing and reflecting more heat (Fahad et al., 2018; Rahman et al., 2017). Alamgir et al. (2020) demonstrate that urban growth contributed to the rise of LST by reducing water bodies and vegetative land surfaces. The rising trend in LST impacts land use and land cover classes' thermal capacity, contributing to the urban heat island (UHI) effects (Bonafoni et al., 2017). In addition, the increasing rate of LST has climatic effects such as droughts and seasonal temperature changes (Hossain et al., 2019). A study has proved that urban temperature contributes to global warming (Roy et al., 2014). The land cover changing pattern of Pabna municipality is noticeable and needs to monitor for environmental sustainability. The enormous amount of vegetative, water body and agricultural land has transformed into the urban area, which has adverse effects on the natural ecosystem and biodiversity of the study region. The prominent features of land cover changes in Pabna municipality have direct impacts on LST and other climatic effects. So, monitoring environmental components associated with land cover changes is essential for urban sustainability, which is the scope of this research.

The researcher tried to show the connections between urbanization and land cover changes of Pabna municipality. There is a negative association between urban growth and land cover changes, where urban areas are increasing in size and other land covers are continuously decreasing in size. If the current rate of urbanization persists, approximately 45% of the land will be converted to an urban area by 2028. Patarkalashvili (2017) indicates that vegetative land cover is a natural filter and retaining noise to improve the urban microclimate and natural resources. Another study highlighted that the vegetative areas help reduce and regulate the high temperature in urban areas (Nzoiwu et al., 2017). Nzoiwu et al. (2017) also suggest that conservation of vegetative cover through reforestation, tree planting, parks, and botanical garden could minimize the adverse effects of urbanization and land cover changes. Thus, as evidenced by previous studies, we recommended conserving and vegetative cover and more planting to manage the arising problems of urban expansion. Another study suggests vertical urban growth for optimal utilization of land (Fenta et al., 2017). Fenta et al. (2017) also mentioned proper land use planning to ensure social, economic, and environmentally sustainable development through stakeholders, policymakers, and practitioners. This research also suggests considering vertical expansion rather than horizontal development due to the compact urban nature of Pabna municipality. As a result, appropriate land use and zoning guidelines are essential for urban sustainability, both controlling urban growth and conserving other land use categories.

Conclusion

This paper quantitatively assesses the urban expansion of Pabna municipality in Bangladesh. The result indicates that the Pabna municipality has experienced rapid urban growth from 1998 to 2018 because of rapid urbanization and practical development. The urban built-up area is increased by 19.63%, from 3.39 km2 in 1998 to 31.93 km2 in 2018. In contrast, all other land cover categories are decreased in size, contributing to urban expansion. Vegetation and water body have reduced by 8.16% and 5.79% respectively from 1998 to 2008. The urban growth is mainly occurring in the northeast periphery of the city due to the availability of vegetative and agricultural land. The urban development is noticeable in the West, Northwest, Northeast, and East parts of the city. The predicted land cover map of 2023 and 2028 suggests that the current trend of urban expansion will continue in the future. The area of urban build-up land will cover 11.01 km2 in 2023 and 12.44 km2 in 2028, which was 3.39 km2 in 1998. Bare ground and vegetative land cover have accounted for 0.49 km2 and 4.13 km2 in 2028, respectively, a gradual reduction from 1.11 km2 and 8.10 km2 in 1998. Agricultural land measured at 9.39 km2 in 1998, and in 2023 and 2028, it is forecasted to be 7.44 km2 and 6.78 km2, respectively. A similar decreasing trend is observed for a water body, from 5.55 km2 in 1998 to 3.70 km2 in 2028. Notably, the North, Northeast, Northwest, West, and East part of the city will continue to share the land of urban expansion, where the Northeast region will contribute most. This study demonstrates that satellite remote sensing data helps detect land cover changes and spatial patterns of urban growth. However, the images were three historical periods' data, which impacted interpretation results. The more reliable result is found for the short interval time-series satellite data. As a result, a robust and credible dataset should consider for more accuracy. The findings from the results highlighted that the changes in land cover affect the urban environmental condition. Therefore, the urban planner and city planning authority should monitor the rate and dimension of urban expansion to ensure sustainable development. The policymaker could consider urban growth factors under different land cover scenarios to alleviate the negative consequence caused by unplanned and haphazard urban development. Also, there is a need for further study to monitor the impacts of urban expansion on urban climate to take effective mitigation measures.

Declarations

Author contribution statement

Subroto Sarker & Md.

Sohel Rana: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

The Landsat satellite images analyzed during the study are collected from U. S. Geological Survey. These datasets were derived from the following public domain resources: https://www.usgs.gov/.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.