1. Introduction
By 2050, the world’s urban population is expected to nearly double, from 3.5 billion in 2010 [
1]. Studies suggest that most of the growth will take place in developing countries, and cities may expand to as much as three times their sizes [
2,
3]. Urbanization brings forth the replacement of natural landscapes with built-up surfaces for accommodation [
4,
5,
6,
7]. Impervious surface knowledge is an important indicator of urbanization degree [
8] as well as environmental quality [
9]. The consequences of the surge in impervious cover include increases in the volume, duration, and intensity of surface runoff [
10,
11,
12]; decreases of groundwater recharge [
13,
14,
15,
16,
17,
18,
19,
20]; and the degradation of receiving water sources, as they directly impact the flow of nonpoint source pollutants [
21], which also affect the hydrologic cycle. Timely and accurate information on the changing land cover of developing cities will aid in decision making processes. This interests scientists, resources managers, and planners [
22]. Such information can also benefit the development of Sustainable Development Goals (SDG Goal) targeted for the developing world by the United Nations. SDG Goal 11 aims to make cities safe, resilient, and sustainable, while SDG goal 6 aims to provide proper water supply and sanitation for all by 2030, demanding comprehensive knowledge of growth rates and expansion patterns of cities. The authors of [
23] bring attention to the overwhelmed capacities of urban infrastructure such as water supply networks, solid waste management, and storm water drainage, due to lack of information on the expected amount of impervious surface and the corresponding surface runoff.
Both remotely sensed data and field measurements are used to quantify impervious surfaces. Ground truth acquisition, though reliable, is expensive and time consuming, especially for developing countries. Research in [
24] highlights the popularity of satellite datasets in environmental studies with impervious-surface applications. Land use and land cover classification has advanced rapidly with remote sensing data [
25,
26,
27,
28,
29,
30,
31].
Despite advances in space-borne sensors, limitations such as data not linking uniquely to land cover, and ambiguity with land uses persists. The authors’ data from [
32] identified fusion as a promising tool, due to its ability to combine information on land properties from sensors operating with fundamentally different physical principles. They found 28 studies supporting improved results of fusion relative to single data source in their review. For example, the spectral reflectance and various indices (e.g., normalized difference built-up index (NDBI)) of built-up surfaces from several bands of optical sensor depend on the material. Similarly, microwave energy, scattered by those surfaces from radar sensor, depends on orientation and dielectric properties. Image fusion techniques can be distinguished into three categories based on their integration level: pixel level, feature level, and decision level. The first category comprises geocoding and co-registration to stack images’ pixels [
33], while the second refers to each dataset’s feature extraction [
34,
35,
36,
37].
The European Space Agency (ESA) developed the Copernicus Mission, with a set of satellite families called Sentinel (S). S-1A and S-2A were launched on 3 April 2014, and 23 June 2015, respectively, followed by their twins (1B and 2B) on 25 April 2015 and 7 March 2017, respectively. Mapping of land surfaces is one of their applications. S-1 are active radar pairs, immune to weather conditions, while S-2 are passive optical satellites, with 13 bands [
38]. Both satellite-based radar and optical data were constrained by low spatial and temporal coverage of medium resolution data in the past, which can be overcome with high spatial-, temporal-, and spectral-resolution S-1 and S-2 data [
32]. The free and public data access policy of Sentinel data makes them more prospective for studying the additive value of data fusion.
While there are studies available for impervious surface mapping (ISM) of national, regional and provincial levels, studies focusing on spatiotemporal change at city scales are still scarce. A smaller spatial scale demands more attention to details as a city is full of dynamic land cover classes as there is a higher chance of mixed classes within a pixel. For comprehensive summary of impervious cover worldwide the following land cover dataset exists: MODIS Land Cover type product (annual 500 m resolution), GlobCover (GLC) 2009 (300 m resolution, 1 January 2009–1 January 2010), the China impervious fraction map (500 m resolution, 2000, 2005, 2010), and Copernicus Global Land Service (100 m resolution, annual from 2015–2019) (Copernicus Service Information). Such coarse temporal and spatial resolution render them inadequate for in depth analysis.
Conventional parametric statistical techniques, common with remote sensing, are not appropriate for classification [
39]. Thus, more general ensemble classifiers are gaining attention, such as random forest (RF). They are known to have improved classification accuracy significantly [
40]. Lighter computation and insensitivity to noise or overtraining with features to estimate the importance of each feature are some of its attractive characteristics.
Past studies have combined optical and radar data and utilized machine learning for land use land cover change assessment in cloud-prone regions [
41], winter vegetation cover extraction [
42], crop mapping [
43,
44], plastic-mulched landcover extraction [
45], monitoring of tropical rainforest [
46], monitoring wetlands [
47], urban mapping [
48], and change detection [
49]. Land cover classification for various purposes such as flood mapping [
50,
51], agricultural landscapes [
52], forest cover derivation [
53] also took advantage of the fusion of S-1 and S-2 data. The authors of [
54] developed various fusion strategies to combine different remote sensing data. They focused on improved discrimination power of remote sensing data due to additional information from SAR images, such as height of terrain. Such integration has proven improvement over optical data alone, with an increase in overall accuracy of up to 6% [
41,
42] and up to 13.5% [
52]. There are not many studies emphasizing the extraction of man-made built-up surfaces for urban monitoring except with reference to filling in the temporal gaps of cloudy optical radar imagery [
50]. For single-date multi-source data, the availability of S-2 data in the vicinity of S-1 data is a challenge. The present study thus assesses the contribution of radar data in the fused dataset in the form of complementary information, such as texture and backscatter, and emphasizes the effectiveness of the method in extracting built-up surfaces. Synergetic use of optical and radar sensors were popular in the past [
55,
56,
57], however, very little exploitation is seen in the fusion of S-1 and S-2 data in particular.
This study used pixel-based fusion of S-1 and S-2 datasets, with the RF land cover classification method, to map the IS of nine cities of Punjab, the most populous province of India. McNemar, a non-parametric test [
58] was used to confirm the better performance of the fused data in land cover classification than optical data alone. Some peculiar geographical cases were also assessed where the fused data can outperform optical. Each city’s impervious surface area and growth information can be used to update Punjab’s impervious surfaces inventory. Faisalabad, Gujranwala, Sialkot, Sheikhupura, and Rawalpindi are identified to have high urban growth [
59]. These cities cover the different sizes, climatic conditions, locations, and geography within the province. The results of ISM in 2016 and 2020 were used to compute the growth rate of each city.
4. Discussion
It is apparent from the results in
Table 7 that in all cities either in one or both time periods, the OA and kappa coefficients of the classified maps from fused datasets show improvement. There is a 2 to 9% increase in OA and a 3 to 12 % increase in kappa coefficient. Such instances emphasize the potentiality of S-1 data to add values to the cities’ impervious surfaces mapping relative to S-2 data alone.
The low resolution (100 m) Proba-V images used by the global dataset tended to overestimate the IS in most of the cities (refer
Table 8), as the pixels were likely to mix the coexisting heterogeneous covers. In [
26,
27,
86,
87,
88] the authors also argue on the inadequacy of coarse spatial resolution in representing spatial disparity within cities. From the HRI, many city locations were identified with a mixture of tree cover, grassland, water canal, buildings, and other open spaces. Different land covers in such proximity require high resolution images for clear discrimination. Thus, in comparison to global datasets, the results from fused data underestimated the values ranging from 0.5% to 15%.
Authors from [
89] used the McNemar test in their study to test the statistical significance in the difference of classification result from maximum likelihood classification (MLC) and Iterative Self-Organizing Data Analysis Technique and texture analysis (Iso-Tex) classification. Similar to this study, [
89] also emphasizes on extraction of built-up cover to map impervious surfaces. They have utilized SAR images and derived texture bands (energy, mean and variance) from Sentinel-1. According to [
90], heterogenous urban areas have distinct spectral signature on SAR images. [
89] concluded the potentiality of texture information in discriminating built-up areas from other land covers such as green areas, bare soil, and water bodies in urban areas. They found the most significant difference on variance texture band for VV based on the highest mean percentage difference on spectral signature between urban and non-urban areas. Significant contribution of variance texture band was also identified in this study through the spectrum plots.
From the spectral signature plots, valuable information added by the radar bands in the land cover classification is identified. In the sigma dB backscattering band and the variance texture derived from radar image, built up cover had the highest value, which could be important information for the RF classifiers. Such a distinct feature is helpful for the classifier in differentiating different land covers. This might have improved the overall accuracy and kappa coefficient of classified maps from the fused dataset. One of the constraints of optical images is limited spectral resolution, which causes different building roofs, roads, and parking lots to appear as different colors. Similarly, high spectral variation within the same land cover type and low spectral variation between different land cover types may occur in urban areas. Therefore, sole reliance on optical bands for land cover classification is not recommended. This makes the automatic extraction of impervious surfaces challenging [
91]. Therefore, [
92,
93,
94,
95,
96] suggest texture as one of the measures to reduce the impact of spectral variation within the same land cover (as cited in [
70]). With free access S-1 datasets, calculation of texture measures is possible with low computational cost, unlike in the past.
Figure 6 shows land cover maps of Rawalpindi–Islamabad from optical and fused data, respectively. They are the northernmost cities considered in this study. Unlike all other cities with flat topography, a part of the Siwalik Range surrounds Rawalpindi and neighboring Islamabad to the north. The encircled region in
Figure 6a depicts that the map from optical data misclassified vegetation into water. According to [
97,
98], topographic shadows of irregular mountains create a major obstacle in accurate land cover classification while using optical images. Such shadows are often misclassified as water as both have very low reflectance [
99]. Therefore, shadow removal techniques were applied during classification [
97,
98]. Shadows commonly occur in high resolution optical images in the form of dark features, reducing the spectral values of shaded objects affecting the land cover classification. One of the sources of shadow in urban areas are tall buildings [
100]. Thus, built up surfaces in such areas are likely to be misclassified as water. In their study [
100] performs a hierarchical classification method to first extract shadow and non-shadow areas and divide the latter into parking lots and non-parking lots to finally classify buildings and roads. They use band indices such as brightness, normalized difference water index (NDWI), and Ration G, through thresholding and vector data, to discriminate shadows. A study in [
97] used DEM data and various topographic correction methods, such as improved cosine correction, C-correction, Minnaert correction, statistical empirical corrections, and a variable empirical coefficient algorithm to reduce relief effects in forest cover classification of a complex mountainous topography. In contrast
Figure 6b shows that misclassified water pixels were significantly reduced by the fused dataset. This demonstrates the usefulness of radar bands to eliminate shadow effects without any correction algorithm. Unlike the optical image, radar images consider the geometry of the object backscattering the signal. Also surface roughness and internal structure of the surface can be retrieved. Therefore, inclusion of derived texture measures from the radar band can reduce the shadow effects without any complex topographic correction or complex shadow removal techniques.
Most of the cities considered in this study are surrounded by cultivation and other vegetation cover. In contrast, Bahawalpur has significant barren cover to its south and east that can also be seen in
Figure 7a. The 2020 optical data overestimated Bahawalpur’s impervious surface area by over 20%. Therefore, the city’s impervious surface map from fused data was overlaid on the impervious surface map from optical data for further examination. In
Figure 7b it is seen that considerable barren pixels were misclassified as built-up to the south and east, as well as on the riverbanks. Misclassification of barren surface into built up was also observed in other cities lying on the riverbank (Multan and Rawalpindi). The impervious surfaces overestimation from optical data alone in other cities ranged from 1.5% in Sahiwal to 7.1% in Multan, for 2020. A similar case was also observed in a study conducted by [
101] where river sand class was misclassified into built-up class. The built-up area extraction from spectral or spatial indexes or spectral alone data does not address the confusion between built-up and other land cover types as it emphasizes a single class [
102]. Due to heterogeneous spectral characteristics of urban areas, spectral confusion is created with other land cover classes. For instance, similar spectral signatures of barren land and asphalt concrete results in misclassification of barren into built-up and vice versa during land cover classification. Therefore, texture measures are widely used in built-up area extraction [
103,
104,
105,
106,
107,
108]. Each texture measure obtained from the radar band provides unique spatial information to discriminate between land cover types [
102] and can complement reflectance information from optical images.
IS expansion in these cities occurred at a cost of replacing natural surfaces, such as agricultural land, open spaces, and tree cover with built-up surfaces. Increasing such surfaces surges the wastewater and stormwater effluents. If such effluent is directed to streams without proper treatment, both humans and aquatic lives are at risk [
109,
110]. Authors in [
111] performed water quality modeling of the Ravi River in their study. Drainage of various urban areas of Sheikhupura, Faisalabad, and Sahiwal districts ultimately carry the untreated effluent into the river [
112,
113]. They identified a distinct difference between the water quality upstream of the river and at downstream locations in vicinity to urban settlements. The latter resulted in high concentrations of total dissolved solids (TDS), ammonia, and other toxic compounds [
111]. Like the Ravi, the Chenab River also faced adverse contamination with cadmium and lead, as identified by [
114]. Their results revealed rapid urbanization as one of the major sources of metal pollution. Urban heat island (UHI) effects are another consequence [
115,
116]. For instance, [
117] identified a 1.52 °C average warming effect in Islamabad from 1993 to 2018, due to the city’s increased IS. Thus, knowledge of IS spatial distribution in cities can aid in addressing the urban climate and environmental pollution issues.
6. Conclusions
In this study S-1 combined with S-2 were used to quantify the IS of nine Pakistani cities (Rawalpindi–Islamabad, Multan, Faisalabad, Gujranwala, Bahawalpur, Sahiwal, Sheikhupura, and Khanewal) in 2016 and 2020. It also estimated the rates of expansion in those cities. An RF classifier was used for land cover classification defined by four classes: barren, built-up, vegetation, and water. The maps were reclassified into binary maps to emphasize the IS. The results from optical and combined data were analyzed and values added by S-1 were interpreted with the support of land cover maps.
Results show that the IS expansion took place in all cities under consideration. Rawalpindi–Islamabad had the highest growth among the cities considered, with the rate of 2.5%, while Khanewal had the least growth, with a rate of 0.5%. The annual expansion rates of the remaining cities ranged from 1% to 2.25%. From CM there was a 2% to 9% increase in OA and a 3% to 12 % increase in kappa coefficient of the classified maps. Spectrum plots showed built up surfaces had the highest values in backscattering and variance texture bands. This could have contributed to improving map accuracies. Besides, in several optical bands, the signatures of barren soil were in the vicinity of the built-up surface. This likely resulted in the overestimation of IS from optical data. The maximum overestimation was seen in the city with significant barren cover on the periphery, i.e., Bahawalpur. The overestimation was more than 20%. However, S-1 was able to remove the misclassified topography shadows into water in the land cover maps of Rawalpindi. Such phenomena are common in places with irregular topography, according to past studies. Estimation of IS in cities surrounded with irregular topography and significant barren surfaces can be improved by using fused radar and optical dataset. Thus, this study was able to develop a cost and labor effective methodology to combine radar and optical data to improve impervious surface mapping at a city scale. This will aid in updating information on impervious surfaces, as well as present and future trends of a city’s expansion.