Ha Long—Cam Pha Cities Evolution Analysis Utilizing Remote Sensing Data

: Socio-economic development has promoted the modiﬁcation of land cover patterns in the coastal area of Ha Long, Cam Pha cities since the 1990s. The urban growth, together with intensive coal mining activities, has improved the life quality of residents. However, it has also caused many environmental problems in this region. Change detection techniques based on post-classiﬁcation comparison were applied for monitoring the spatial and temporal evolution of land covers. The confusion matrix for 2001 and 2019 showed high overall accuracy (97.99%, 94.95%) and Kappa coefﬁcient (0.97, 0.92), respectively. Statistics from classiﬁed images have revealed that man-made features increased by about 15.32%, while natural features, mangrove jungles, and water bodies decreased 10.64%, 1.96%, 2.72%, respectively, and urban evolution presents various dynamics, soft in the ﬁrst period (1991–2001), but stronger in the second period (2001–2019) with different characteristics. The study also expresses the constraint of topographic and geologic resources, which have prevented the urban development in this coastal area. Such obtained results are very important for understanding interactions and relations between natural and human phenomena and they may help authorities by providing indicators and maps able to highlight necessary actions for sustainable development.


Introduction
The statistics of the United Nations showed that the Vietnamese urban population has accelerated from 14.7% in 1960 to 37% in 2020 [1], and the amount of urban centers has rapidly risen in recent years, reaching to 862 in 2020 [2] in comparison with 480 in 1986 [3]. Ha Long and Cam Pha cities have become the largest centers of the coal mining industry in the Northern East region of Vietnam owing to considerable coal reserves in the underground [4]. Over the past decades, rapid socio-economic development has promoted urbanization at Ha Long and Cam Pha cities. Specifically, the recent statistics show that their urban population comprised 76.5% of the total population in 2010, but it reached 88.0% in 2019 [5]. The urban development has interacted with the intensity of coal mining exploitation, causing significant pressures on the environment with many dynamics at the coastal area of Ha Long and Cam Pha cities, such as the expansion of man-made infrastructures [6], the modification of hydrologic network and topographic surface [7], the dissemination of environmental problems such as air pollution, water pollutions, soil degradation [8], and land subsidence hazard [9], etc. Thus, understanding and monitoring urban evolution over a long period in the coastal area of Ha Long and Cam Pha cities are crucial for various purposes such as land management, public service provision, spatial planning, territorial organization, and sustainable objectives achievements [10,11]. This might be an effective support to manage the urban development towards suitable land-use decisions in the future [12].
Urban growth can be highlighted through the environmental impacts caused by increased consumption of food, energy, water, land, etc. [13], which induces energy demands with the development of residential settlements [14], alters the local climate with the increase of urban heat flux [15], perturbs the natural water cycle with the extension of impervious surfaces [16], reduces biodiversity with the rise of habitat fragmentation [17], stimulates the expansion of the industrial zones with pollutant emission [18], converts land surface from non-urban to urban uses with the loss of arable soil, and intensive pressure on agricultural production [19]. However, the researchers noted that these issues could be described and traced notably via land use/land cover changes studies based on reflected and emitted measurements at various spectrum wavelengths [20,21]. The spectral variations between different image acquisitions due to the alteration of land use/land cover in urbanized areas are used for identifying the changes when the magnitude of spectral variations is more important than any noises [22].
As a matter of fact, monitoring urban evolution trends have started with the launches of Earth Observation satellites for spatially and temporarily determining land use/land cover pattern changes (location, intensity, period) through satellite data processing [23][24][25]. Studying the urban evolution over several decades requires a historical dataset with the repetitive acquisitions of multispectral images of the ground surface. The Landsat program can satisfy these requirements owing to the successive launches of some satellites from Landsat 1 to Landsat 8 (except Landsat 6 due to the failure after not reaching the velocity necessary) since the 1970s by the joint mission of the National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS) [26]. This Landsat program might be expected as a potential source of information for identifying land use/land cover changes for environmental management support at a regional scale [27] because it has successively provided identical image specifications with medium spatial and temporal resolutions [28]. Otherwise, Landsat images with enhanced spatial resolution at 30 m of multispectral bands could be considered as suitable and reliable datasets for analyzing urban evolution in a more timely and cost-effective manner [29]. Therefore, Landsat images have been exploited to investigate when, how, and where Earth's surface has been changing because of its historical archive availability over the last five decades. Although the number of studies using Landsat data has dramatically increased after the opening of entire Landsat collections for free access in 2008 [30], the studies on urban evolution have been confined in an inadequate amount in which the methods of change detection have been revolutionized depending on the number of handled images [31].
Several researchers noted that change detection is the process of identifying mutations in the state of an object or phenomenon by simultaneously processing two or more images of the same geographical region acquired at different times [32,33]. Meanwhile, Radke et al. [34] supposed that change detection was considered as a process to determine considerable differences in sequential pixel appearances due to the emergence, disappearance, movement, or shape alteration of objects on the ground surface [34]. For detecting changes in land use/land cover, many issues have been discussed and synthesized in literature reviews of several researchers' articles [31,[33][34][35][36][37][38][39][40][41]. Salah et al. [42] mentioned the different categorization schemes in previous reviews, which have usually focused on a few aspects or dimensions of change detection problems while ignoring others because of the complexity of this topic. However, Hemati et al. [31] and Zhu [41] have suggested six categorical change detection methods such as thresholding, differencing, segmentation, trajectory classification, statistical boundary, and regression, which have been applied in diverse applications using Landsat data to observe the dynamics caused by human activities including urban evolution analysis [43,44].
Land use/land cover changes can be detected by comparing satellite images acquired at two different moments-the bitemporal approach, which is a simple but less comprehen-sive understanding of the urban dynamics in comparison to the time-series approach [45]. Indeed, in this recent time-series approach, a set of satellite image scenes taken at many different times-Satellite Image Time Series (SITS) affords a large amount of information compared to a single image couple in the context of temporal tendencies of regional evolution [46]. Despite these benefits, it still raises specific challenges regarding: the irregular temporal phenological signature of different land cover types; the insufficient sampling used to train the supervised classification; the missing temporal data [42]; the network architectures or specific datasets shaping that need to be developed for exploiting the temporal information jointly with the spatial and spectral information of the data [47]. Thus, in a more classical way, other sets of approaches and methods can be used varying from manual change interpretation [48] to bi-temporal linear data transformation [49] or multi-temporal spectral mixture analysis [50] and deep learning [51].
Although in the coming years, together with Van Don Special Economic Zone, the two adjacent cities of Ha Long and Cam Pha will become the pole of development in the North-East region. However, there was not research on urban evolution in this dynamic region except at Ha Long city, where the relevant topic had been considered in the works of Brömme et al. [52], Hens et al. [6], Jalilov et al. [53]. Therefore, the present study examines the spatio-temporal land cover change utilizing Landsat multitemporal and multispectral images for a period of 30 years (i.e., 1991, 2001, and 2019). The aims are to (a) identify land cover changes using different sensors of Landsat satellites to characterize urban dynamics at the coastal area of Ha Long-Cam Pha cities over this period and highlight the different tendencies of urban evolution between these two cities; (b) verify the interest of topographic gradient issued from Shuttle Radar Topography Mission (SRTM) digital elevation model for analyzing the land suitability in urban development. It will then analyze the concern for urban planning policies in the region.
After the context of this study (Section 1), the presentation of the study area is provided (Section 2), the description of the data used, and the various processing steps are presented (Section 3). The results are followed with a validation procedure (Section 4). These parts are subsequent to a discussion (Section 5), a conclusion closes the paper (Section 6).

Study Area
The study area is mainly located in the Quang Ninh province in the North of Vietnam, focused on the coastal area of Ha Long and Cam Pha cities, which are adjacent to each other. The Quang Ninh province is closed to the Chinese frontier and expanded along the Tonkin gulf. This area has remarkably diverse landscapes: hills and mountains predominating with elevations ranging from 0 to 1498 m and narrow coastal alluvial plains, where urban development can be extended toward shallow intertidal zone ( Figure 1). The remaining part is a vast gulf including hundreds of big and small limestone islands.
The territory of these two cities covers an area of 758 km 2 from the latitude 20 • 54 35 N to 21 • 14 05 N and from the longitude 106 • 54 08 E to 107 • 24 08 E. It is influenced by the tropical monsoon climate, in which the annual average temperature reaches 22.6-24.4 • C, and the annual average rainfall is approximately 2200-2700 mm [54]. The climate is separated into two distinctive seasons: hot in the summer (lasting from May to October) with heavy rainfall (accounting for 80% to 85% of the total annual precipitation) and cold in the winter (lasting from November to April) with low rainfall (15% to 20% of the total annual precipitation) [55]. The hydrographic system is mainly composed of small, short rivers originating from Northern mountainous districts characterized by a fast-rising flow and quick drainage to the sea due to the steep terrain and heavy rains in the rainy season.
Many coal mines are located at the central part of the study area, and the residential settlements were developed on both sides of the national road No18 with various man-made structures extending along the shoreline. However, the tropical rainforest still dominates in the Northern part. Besides, carbonate sediments with thick limestone layers have developed in the marine environment, generating unique landscapes with high geomorphological values. The Ha Long bay is recognized by the United Nations Educational, Scientific and Cultural Organization (UNESCO) as the World's Natural Heritage site since 1994 [56]. These resources provide comparative advantages for the Quang Ninh province in several developing economic sectors, including coal mining, aquaculture, tourism, etc. For example, the total number of tourist arrivals reached around 21.78 million, and the tourism sector contributed 8% to the Gross Domestic Product (GDP) of Quang Ninh province in the period of 2013-2015 [57]. Otherwise, in the coal mining sector, 80,000 employees have worked in 23 major coal mines yielding more than 38.5 million tons of raw coal in 2020 [58]. By 2020, the population of the two cities reached 518,300 people, prevailed in 21 wards, 12 communes (Ha Long city), and 13 wards, three communes (Cam Pha city) corresponding to an average density of 292 people km −2 and 570 people km −2 respectively [5]. In accompanying the dynamic socio-economic activities, local authorities have faced institutional challenges in environmental management [61], especially when these activities have intertwined with the rapid urban development on narrow piedmont situated along the coast (Figure 1).

Landsat Images
In change detection analysis, seasonal variation due to solar angles and phenological changes is taken into consideration as the main source of noise, and it can be minimized by image selection [41]. In the present study, three aspects were necessary to consider: (a) the selection of cloud-free images, as in this cloudy region, at least around 10% was acceptable for the mainland part of the scene; (b) the selected images were captured in the period between Summer and Autumn when crop makes agricultural land in suburban areas highly distinguishing with the residence in urban areas; (c) the acquisition date must correspond to the moment of changing administrative structures at both cities.
Landsat images were chosen in moderate weather conditions corresponding to 3 acquisition dates. Cloud had partly covered some sections of the scenes, but this had not seriously impacted the study (Figure 4). The parameters from Landsat images are shown in Table 1.  In this study, Landsat images were downloaded from the website of the USGS: https://earthexplorer.usgs.gov/ (accessed on 15 August 2021). The images were preprocessed at the 2A level, providing surface reflectance products (Bottom Of Atmosphere reflectance-BOA), and they were geo-referenced in the UTM coordinate system (WGS84). Landsat 5 satellite provides 7 bands in a visible spectrum, infrared, and thermal parts with different spatial resolutions: 30 m and 120 m, respectively. Landsat 7 satellite has an additional panchromatic band with a spatial resolution of 15 m. Landsat 8 satellite provides 11 bands with identical characteristics, but it also has a panchromatic band ( Table 2). All spectral bands were exploited to identify the land covers in the study area.

High Spatial Resolution Imagery from Google Earth
Google Earth integrates several types of spatial data such as satellite images in natural color composite, aerial photos, and geographic information layers. Images can be displayed at the present time or at a specific time in the past. Currently, Google Earth images are provided with a super high spatial resolution of less than 1.0 m [63] and are regularly updated over time. Users can recognize features on the ground using Google Earth Pro applications. Google Earth images are suitable for interpreting land covers because all features on the images are clearly identified with prior knowledge via field surveys [64]. In this study, the images captured on 17 August 2001 and 8 November 2019 were used for selecting the ground truth region of interest (ROIs) by visual interpretation, while the image of 1991 was not available. Figure 1 presents the SRTM Digital Elevation Model Version 1 from NASA with a spatial resolution of 90 m. It was exploited for identifying the topographic gradient of the study area [65].

Field Campaign
The selection of training samples for a supervised classification must be independently applied to three Landsat images. However, the recognition of land cover categories was difficult due to the ambiguity of features on two ancient images (1991,2001). A field campaign needs to be taken place for referencing land cover categories. Each land cover category was identified with five sampling sites, which supported the selection of the training samples on Landsat images. Hence, the field campaign was carried out with five team members on 12-13 November 2019. Photos of some sites are shown in Figure 5.

Method
Zhu [41] noted that whether change detection approaches were time series or bitemporal, the differencing was most commonly adopted among change detection methods by using images at different acquisition dates. Changes in land cover will be identified by differencing classification results, also known as a post-classification comparison [35,39,66]. This principle was applied to urban change detection from the early 1970s, and it required the comparison of classified images corresponding to spectral variations over time due to the urban land use/land cover changes [67]. In recent years, many advanced classifiers have been taken into account to process satellite images, and image classification approaches can be distinguished into some categorical classifiers: supervised and unsupervised, or parametric and non-parametric, or hard and soft, or per-pixel, sub-pixel, and per-field. For detailed information about each category, readers can refer to the research of Lu and Weng [68]. In general, different image classifiers have their own advantages and disadvantages, and the selection of suitable classification methods depends on the sufficient training samples, the distribution of features on the ground, the spatial resolution of satellite images, the availability of classification software [69].

Pre-Processing
Landsat images were atmospherically corrected for the surface reflectance at 2A level, geo-referenced to the UTM coordinate system (WGS84) by the USGS. The multispectral bands of these images with the spatial resolution of 30 m were pre-processed as illustrated in Figure 6: 1.
Co-registration with respect to each image ensures that the images become spatially aligned thus that all features in one image overlap as well as possible with its footprint in other images. Accurate image co-registration is a prerequisite to the accurate extraction of features, and it may subsequently provide correct land cover mapping results.

2.
Cropping image aims to discard the unwanted portion outer areas from Landsat's scene for preserving the important part-the study area (Figure 1).

3.
Stacking layers address to combine several channels of Landsat image with the identical frame of reference, thus that multi-channels can be processed later in the image processing software. The reflectance values in integer format were converted into a floating format according to Formula (1) for classifying at each acquisition date.
Here, Ref n is the reflectance value of band n in a decimal number, B n is the apparent value of band n in an integer number [70].

Image Classification
Image classification is a favorite method to identify the land cover categories because it provides the ability to generate a series of land cover maps at different moments. Lu and Weng [68] described many classifiers, which have been successfully applied into image classification, such as artificial neural network [71], fuzzy classification [72], or object-based classification [73], maximum likelihood [74], etc. In the present research, we applied the Support Vector Machines (SVMs) classifier, an algorithm of Supervised Machine Learning developed by Cortes and Vapnik [75]. It was initially applied to machine vision domains such as handwriting digit and text recognition based on Statistical Learning Theory [76]. Subsequently, it was experienced in satellite image classification [77][78][79]. Although this classifier was developed for about 3 decades, the results of recent researches proved the prominence of SVMs in comparison with other usual classifiers due to its ability to generalize desirable events despite limited training sample sets [80][81][82][83]. Otherwise, Paneque-Gaslvez et al. proved that SVMs outperformed all the rest for establishing an efficient classification approach to accurate map all broad land cover classes in a large, heterogeneous tropical area [84], thus that Support Vector Machines (SVMs) were chosen for classifying land use/land cover categories in urban studies and this classifier yielded better performance regarding accuracy and generalization [85,86].
In order to proceed with the supervised classification, three sets of training samples were, respectively, determined here by on-screen digitization based on color composite images of Landsat 1991, 2001, and 2019. Each sample set corresponds to 10 samples with an average of 200 pixels/sample for each land cover category (Figure 7). A code corresponding to each category of level 2 in Table 2 was assigned to each class, of which a distinctive radiometric curve was calculated by the average values of pixels within the window sample for all spectral bands (Figure 8). The SVMs classifier was executed to separate them from each other using this training sample set.

Validation
The accuracy of change detection highly depends on the classification accuracy, and errors in each classified image are present in the final change products. The obtained results are thus vulnerable to classification errors [41]. Moreover, classification accuracy can be affected by many reasons, such as the classification methods, the algorithms used, and the type of satellite images [68]. In our study, Landsat's SVMs classification accuracy was visually adjudged by knowledge of the study area. In accordance with this interest, the selection of ground truth samples was realized from Google Earth images in order to quantitatively assess the Landsat's SVMs classification accuracy using confusion matrix, on the one hand [88]. In addition, the outline of features issued from the classified images was visually explicit on Google Earth images to identify the coincidence between these two data sets, on the other hand [89].

Keys for Visual Interpretation
Land cover categories can be apparently recognized on Google Earth images through visual interpretation due to high spatial resolution in comparison to Landsat images. Figure 9 below presents the elements of visual interpretation from satellite images. However, some interpretation keys are used to choose ground truth samples on Google Earth images captured on 17 August 2001 and 8 November 2019 as follows: Site: The mangrove is often located near the natural water bodies such as the riverbanks or estuaries (Figure 10e). However, many mangrove jungles have been partly replaced by aquaculture ponds, and the rest of the mangrove jungle became adjacent ponds where water impoundment is constructed and managed for farming fish, shrimp, mollusks, etc. Structure and shape: It is easy to distinguish the structure of mangrove with other vegetation types because mangrove's foliage leaf often has a finer structure than forest but coarser structure than a paddy field (Figure 10b,e,f). Tone/color: Coalfields are usually black (Figure 10a), while mangrove is dark green, logging land is brown, rock dump is grey on the natural color composite. It helps to clearly distinguish mangrove from other vegetation such as terrestrial forest, paddy field, logging land, etc. (Figure 10b,c,e,f).
Association: It is an indirect sign allowing to recognize the presence of objects or phenomena that cannot be directly perceived. However, we can identify their existence via traces, which have been left on the ground. For example, there are often many different types of equipment such as trucks, excavators, conveyors in coal fields (Figure 10a).

Accuracy Assessment
The accuracy of established land cover maps is often assessed via the confusion matrix [91] and the Kappa coefficient [92]. In fact, the land cover maps issued from Landsat's SVMs classifications are compared with ground truth samples, which were independently selected by visual interpretation on high spatial resolution satellite of Google Earth images. Due to the unavailability of Google Earth images in 1991, only two sets of random samples were respectively determined on Google Earth images (17 August 2001 and 8 November 2019) corresponding to eight classes of land covers mentioned in Figure 7. Each set consists of 10 samples with an average of 200 pixels/sample for each land cover category. Various indices might be pointed out from the confusion matrix: the overall accuracy, producer's accuracy, user's accuracy, Kappa coefficient, which are summarized in Table 3. The statistical parameters for 2001 and 2019 show that the overall accuracy alternately reached 97.99%, and 94.95%, while the Kappa coefficients reached 0.97 and 0.92. The user's and producer's accuracies of individual classes are consistently high, ranging from 72% to 100%.

Comparison with High Spatial Resolution Image
In order to validate SVMs classification results, the outline of the land cover map issued from Landsat image (23 September 2019) was investigated by overlaying on Google Earth image (8 November 2019). The five random specific frames are closed-up in Figure 11, showing the properly geometric and thematic matching between two datasets in the validation context.

Filtering Noise
After SVMs classification, the patches of pixels representing detached small and discrete land covers are quite common in the raster dataset and considered as noise. This unwanted noise can cause disturbance in image quality assessment and affect the reader's ability to grasp the information. The filtering approach is introduced into the postclassification stage as an appropriate solution [93][94][95]. Thus, the filters help to eliminate discrete pixels in the image but still preserve the accuracy of results. In this study, the Median filter with a dimension of 3 × 3 window is applied to the SVMs classification products for all acquisition dates. The results show that noises were reduced in the final land cover maps. However, the general shape of objects in the classified images was still preserved (Figure 12).

Combining Classes
Mapping land cover products based on Landsat's SVMs supervised classification identified eight land cover categories for each year (1991,2001,2019). However, change detection in urban areas focused on man-made features analysis because they are the most important artifacts reflecting the interactions between human activities with the natural environment. Therefore, eight original classes were selectively aggregated for responding to the study objectives, and new land cover maps were established with four derived classes: natural features, man-made features, mangrove, water (Figure 13a,b). Specifically, man-made features are composed of artificial land with coal fields, settlements, and rock dumps. In contrast, natural features consist of woodland (forest), bare land (logging land), and cropland (rice field). Meanwhile, mangroves and water were considered as two independent classes with natural features strongly affected by human activities in the study area ( Figure 14).

Change Detection
The change detection matrix of land covers is deduced from 3 classified images. They are synthesized in Tables 4 and 5. The land cover categories can be determined when and what has been changed to what during an interval of time. In the present study, two periods are considered I: from 1991 to 2001 and II: from 2001 to 2019, corresponding to the administrative adjustments of Ha Long, Cam Pha cities. For change detection analysis, the post-classification approach was applied to the comparison of independent classified images for each period. The main benefit of the post-classification approach is the capability to provide the matrix of change information and to minimize external impact from atmospheric and environmental differences between the multi-temporal images [38]. The statistic numbers show that land covers had been slightly changed over the first period, but their mutation in the second period was considerably faster. For instance, the absolute values indicate that the areas of man-made features increased 2.08% for 10 years of the first period, but they have reached 13.24%, corresponding to 6 times for 18 years of the second period. However, the average area of man-made features increased from 248 ha/year in the first period to 828 ha/year in the second period corresponding to 3.3 times, while the area of natural features, mangrove, and water bodies decreased due to the coal exploitation or sea encroachment from 148 ha, 74 ha, and 25 ha/year in the first period to 620 ha, 88 ha, and 166 ha/year in the second period, respectively. That manifests that the temporal interval of the second period is almost twice as long as the first period, the change of land covers does not occur linearly.

Analysis of Urban Evolution
After integrating two classified maps of each period, the situation of each land cover category over two periods is expressed in Table 5.
Ha Long city was established in the early 1990s based on Hon Gai town and then followed by Cam Pha city in the early 21st century based on Cam Pha bourg. These major urban centers provide labor sources for coal mines in the region. After the national opening economy plan in the middle 1980s, the urban evolution at Ha Long-Cam Pha cities emerged with slight changes during the first period and stronger dynamics after 2001. The changes in Ha Long city were bigger than in Cam Pha city ( Figure 15).  (Figure 16a,b). Otherwise, many small patches of man-made features were expanded by pouring waste rock into the sea with the replacement of water bodies (473.94 ha) at Cam Trung and Cam Thinh wards (Figure 16d). The disappearance of man-made features at two cities is explained by environmental reclamation when coal mines have been exhausted (Figure 16b). In the second period, the population growth and economic development promoted the rapid increase of urban change. Specifically, in 2010 the population of Ha Long and Cam Pha cities reached 223,700 people and 177,400 people, respectively. In 2020, the population reached 327,000 and 191,300 people [5]. Moreover, Table 6 illustrates the differentiation of economic structures between the two cities: industrial sectors have been concentred at Cam Pha city, where coal exploitation became dominant. Meanwhile, service sectors have been concentred in Ha Long city of which the potential is tourism with recreation activities. However, the changes are identical between the two cities ( Figure 15).  (Figure 17d), where the demand of real estate has risen with the installation of infrastructures for recreation areas, residential settlements and luxury hotels due to the close proximity to world heritage site [99]. This has constituted a driving force to stimulate the widening of man-made features toward the intertidal zone, which will be described in next Section 5.3 (Figure 17b). Although provincial authorities have perceived the ecological role of mangrove jungles and they desired to preserve them, mangroves have lost 1142.55 ha in comparison with the previous period (Figure 17c). In Cam Pha city, man-made features have continued to expand toward the intertidal zone at the Cam Trung, Cam Thinh, Quang Hanh wards, occupying 744.37 ha of mangrove jungles ( Figure 17c) and 1747.45 ha of water bodies; making the National Route No18 no longer adjacent to the sea, but located in the interior with the farthest distance about 1 km (Figure 17d). During this period, open-cast coal mining flourished in the territory of both cities with 23 major coal mines, which caused the expansion of man-made features for encroaching upon the natural features (Figure 17a,b). These areas are mainly located in the mountains in the central and northern part of the study area.

Analysis of the Urban Development
Urban development requires a suitable space to respond to the specific needs of construction. However, suitable spaces have been limited at Ha Long and Cam Pha cities due to the topographic constraint with high and steep mountains. These limited spaces are also occupied by coal mining fields in the central part of the study area. These settings make both cities face the sea and lean against the cliffs (Figure 18a). Lofrano et al. [100] have assumed that urban development is initially compelled by underlying geologic resources and overlying geomorphologic forms. In addition, Yang, He, and Xu [101] also recognized the role of terrain slope as one of the most important factors for evaluating land suitability for urban development in the mountainous area, and the perfectly suitable steepness ranges from 0 • to 6 • . For analyzing the relationship of urban evolution and urban development at Ha Long and Cam Pha cities, the SRTM digital elevotion model was used to calculate slope map. The values of slope have ranged from 0 • to 71.5 • (Figure 18a) and two gradient classes of suitability for urban development were derived as illustrated in Figure 18b together with the distribution of coal mines. Rivkin [102] mentioned that urban evolution was the adaptation of organisms to heavily populated urban areas. Anthropogenic environmental and geographical changes brought to habitats when urbanization takes place have a significant evolutionary impact on organisms inhabiting these city areas. Whereby, people initially exploited the flat ground located in narrow ravines and along the coast to establish urban settlements, which are connected by the inadequate road network. The national route No18 is a unique backbone for serving the communication of two cities for a long time. The analysis of urban evolution in the previous section expressed that urban development in the next stages occupied suitable places in terms of steepness by filling mining waste rock into water bodies or mangrove jungles at intertidal zones. This tendency is clearly observed in Ha Long city, where man-made features extended across the Cua Luc estuary and at the Northern part of Ha Long bay, while man-made features have principally spread in the Northern part of Bai Tu Long bay in Cam Pha city. Figure 18b presents the disappearance of water bodies and mangroves jungles corresponding to the emergence of man-made features across the period from 1991 till 2019.
In fact, this phenomenon matched with the policy of provincial authority relating to the general land use planning from 2000 to 2020 of Ha Long city [103] and Cam Pha city [104]. By 2019, the flat terrain for urban development in both cities was depleted. Therefore, Ha Long city has been merged with Hoang Bo district to be able to exploit flat terrain at the Northern part of Cua Luc bay and the areas located along with highway No. 9 connecting Ha Long city with Van Don airport. However, the open space for urban development of Cam Pha city has still rested a critical problem. However, over the last decades, on the one hand, the opencast coal mining has emptied coal deposits, and on the other hand, the provincial authority would like to protect the environment by fostering tourism at the World's Natural Heritage site-Ha Long Bay. In consequence, all opencast coal mines in this coastal area have been gradually closed in 2025 to minimize environmental pollution [105]. The reclamation then returns the disturbed land from mining activities into natural habitats. A vast surface near cities's old centers can be exploited to decrease the pressure on open space, although it does not belong to suitable places for urban development.

Conclusions
This study has proved the ability of Landsat datasets to produce accurate land cover maps and to analyze urban evolution at the coastal area of Ha Long-Cam Pha cities using change detection techniques with post-classification comparison approach. The application of SVMs classifier from spectral bands of Landsat images allows classifying land-cover categories at different periods of time. The accuracy assessment based on visual interpretation from Google Earth image captured at the identical period and in-situ observation provides suitable confidence with an overall accuracy of 97.99% and 94.95%, while Kappa coefficients reached 0.97 and 0.92. The study also expresses the constraint of topographic and geologic resources to urban development at the coastal area of Ha Long-Cam Pha cities. The steep slope mountains have prevented the enlargement of these cities to the North direction, while the bays and estuaries characterized by flat terrain are a lifeline for urban expansion to the South direction. Although the spatial resolution of Landsat images up to 30 m has been applied to distinguish the land cover categories, the resolution is not fine enough for features recognition because many small features are not identified in comparison with Google Earth images. However, Landsat images are the single remote sensing dataset providing freely historical image series with long enough duration. The results could help the local authorities to propose appropriate solutions and policies for conserving natural resources and developing socio-economic in a more harmonized manner for this coastal area.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviation are used in this paper.