Impervious Surfaces Mapping at City Scale by Fusion of Radar and Optical Data through a Random Forest Classiﬁer

: Urbanization increases the amount of impervious surfaces, making accurate information on spatial and temporal expansion trends essential; the challenge is to develop a cost- and labor-effective technique that is compatible with the assessment of multiple geographical locations in developing countries. Several studies have identiﬁed the potential of remote sensing and multiple source information in impervious surface quantiﬁcation. Therefore, this study aims to fuse datasets from the Sentinel 1 and 2 Satellites to map the impervious surfaces of nine Pakistani cities and estimate their growth rates from 2016 to 2020 utilizing the random forest algorithm. All bands in the optical and radar images were resampled to 10 m resolution, projected to same coordinate system and geometrically aligned to stack into a single product. The models were then trained, and classiﬁcations were validated with land cover samples from Google Earth’s high-resolution images. Overall accuracies of classiﬁed maps ranged from 85% to 98% with the resultant quantities showing a strong linear relationship (R-squared value of 0.998) with the Copernicus Global Land Services data. There was up to 9% increase in accuracy and up to 12 % increase in kappa coefﬁcient from the fused data with respect to optical alone. A McNemar test conﬁrmed the superiority of fused data. Finally, the cities had growth rates ranging from 0.5% to 2.5%, with an average of 1.8%. The information obtained can alert urban planners and environmentalists to assess impervious surface impacts in the cities.


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 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.

Study Area
The study area consists of nine growing cities of Punjab province, given in Figure 1. There was a total of 194 cities in the province as of 2015. With an overview of all those cities, the Urban Unit conducted a thorough study of 50 of them. These cities represent the dynamics of the province in terms of "urban extent, area, expansion, urban population, population growth, land consumption, densities and economic regions" (Punjab Cities Growth Atlas (2015), Section 3). Nine cities were chosen for this study based on varying urban extent, which represent the overall spatial cover of the province for this study. Since Rawalpindi and Islamabad are connected at their physical boundaries, they are considered a single entity, Rawalpindi-Islamabad, during the assessment. It represents the 21st largest city in the group. The studied cities lie from the northernmost region of the province down to the south. Multan and Bahawalpur lie on the left bank of the Chenab and Sutlej rivers, respectively. Both cities, as well as Khanewal, are in the southern part of the province, which has an arid climate. Cities in the middle reach, such as Sahiwal, Faisalabad, Gujranwala, and Sheikhupura, have semi-arid climates. Rawalpindi-Islamabad has a humid subtropical climate. They lie on the banks of the Haro River, in the northern part of Punjab. Rawalpindi-Islamabad, Multan, Faisalabad, Gujranwala, and Bahawalpur rank 2nd, 3rd, 4th, 5th, and 6th, respectively, among the 50 cities in terms of area, while Sahiwal is the 11th, Sheikhupura is the 13th, and Khanewal is the 21st largest city on the list.

Dataset
This section describes the datasets used for this research. The details of sources and the purposes of the datasets used are given in Table 1. The coordinates of the cities used to define the image extents during preprocessing are listed in Table 2. Land cover classification scheme in random forest classifier is given in Table 3. The acquisition dates of radar and optical images of the nine cities for 2016 and 2020, respectively, are listed in Table 4.

SAR Data
Two S-1 images of each city were obtained from the Copernicus Open Access Hub. They were acquired for February or March 2016 and 2020. S-1 carries a C-band (5.405 GHz) imaging radar, called synthetic aperture radar (C-SAR), which is an active sensor. For land monitoring, ground range detected (GRD) products obtained in interferometric wide (IW) swath mode with 240 km swath are recommended by the ESA S-1 observation strategy.
This mode supports dual-polarization modes (VV and VH) bands. However, for 2016 images, only VV mode was present. The spatial resolution is 10 m, temporal resolution is six days, and swath width is 400 km. All the S-1 data can be freely accessed, according to the free data policy of ESA. The raw S-1 images were in jp2 format. These images were combined with S-2 images to classify the land cover map of each city.

MSI Data
Two S-2 images of each city were also obtained from the same portal as the radar images. The S-2 data were acquired close to that of the S-1 data, with a maximum gap of 20 days. S-2 carries multispectral instrument (MSI) sensors that acquire optical imagery in 13 bands. There are three visible, four red edge, and three infrared, with one each of aerosol, vapor, and cirrus bands in the list. All the visible and near infrared bands are 10 m resolution, while red edge and short wave infrared are 20 m resolution, and the others are 60 m resolution as given in Table 5. The temporal resolution is five days, and the swath width is 290 km. S-2 images from 2020 consist of projected and processed to BOA (bottom of atmospheric) level-2 products, so there was no need to convert them into reflectance value as well. However, for images from 2016, only TOA (top of atmospheric) level-1 products were accessible, requiring further processing. S-2 images are also freely available from ESA. S-2 images alone were also used to classify the land cover maps of the city.

DEM Data
The latest DEM data available for the cities were downloaded from USGS Earth Explorer browser. They were 1 arc-sec shuttle radar topography mission (SRTM) models and were later used for geometrical correction of the S-1 images. The spatial resolution was 30 m. For some cities, two or more DEM data needed to be merged to cover the whole extent.

Google Earth's High-Resolution Images
The training and validation samples for land cover classification of optical and radar images were obtained from Google Earth. It allowed us to acquire the samples for 2016 and 2020 independently. Since the land cover goes through a dynamic process, it was important to extract samples from the nearest time the remote sensing images were obtained. The training and validation samples were not overlapped for the robust accuracy analysis of land cover maps later. Four land cover classes were distinguished prior to extracting impervious surfaces from the cities: water, vegetation, bare soil, and built-up. Later the water, vegetation, and bare soil classes were merged into non-built-up surfaces. Each polygon delineating the training and validation samples were roughly larger than the pixel size of the images.

Global Land Cover Maps
The Copernicus Global Land Service produced 100 m resolution annual global land cover maps from 2015-2019. They utilized images from Proba-V sensors. They were trained and validated using 168 K and 28 K points, respectively, from Geo-WIKI's crowdsourcing. The overall accuracy of their map for 2015 was 80.6%. These datasets were used to compare and correlate the quantity of each city's impervious surface cover, obtained from the land cover maps from the current study. Total number of training and validation samples extracted for each city are listed in Table 6.

Methodology
This section details the algorithms and processing used to combine the S-1 and S-2 images. It also covers the steps to map the impervious surfaces and quantify the growth rates of the nine Punjab cities for 2016 and 2020. There are five subsections of the method: image preprocessing, land cover classification, map validation, impervious surfaces quantification, and correlation analysis. Figure 2 shows a detailed flowchart of the research methodology and Figure 3 shows the steps for IS area quantification.

Image Preprocessing
The following section describes the preprocessing steps for radar and optical data, respectively. This study utilized the Sentinel Application Platform (SNAP) and ArcGIS pro for the process. Firstly, the coordinates given in Table 2 were used to subset the S-1 images for each city. The orbit files were then updated for accurate information on the position and velocity of the satellite. Then a calibration algorithm was applied to convert the raw digital number of each pixel to radiometrically calibrated SAR backscatter values called sigma naught. This was followed by the application of a range Doppler terrain correction algorithm to correct the distortion in the image due to the side looking geometry of the satellites. Finally, a logarithmic transformation of sigma naught values converted them to sigma dB values for end band composition. S-1 images are characterized by grainy structure called speckles which can be removed by speckle filtering in SNAP. However, [60] does not recommend the filtering when the interest is in the identification of small spatial structures, such as isolated houses or image textures, as the process may blur out important information.
For additional information from the radar images, grey level co-occurrence matrices (GLCM) were used to compute texture measures. According to [61,62], such measures capture the spatial relationships of pixels by identifying the pattern based on the neighborhood size provided. The resultant matrices store the occurrence frequency of pixel pairs with specific grey level (G) or pixel brightness values [61]. Such measures are common in image classification [63]. However, the texture measures must not be auto-correlated with each other according to [64]. Therefore, only contrast and variance were chosen, similar to [65]. For a neighborhood, a 5 × 5 moving widow was chosen, and single pixel displacement was used [65,66]. For a consistent classification result, VH polarization was eliminated from the 2020 images during the image sub-setting process. Then there were two additional bands in the form of texture measures in each radar image. Finally, all the SAR images were reprojected to UTM zone 42/43 to align with the extent and coordinate reference system (CRS) of their MSI pairs, using the bilinear interpolation resampling method in ArcGIS.  To convert the level 1C TOA products from 2016, sen2cor algorithms were applied in SNAP. They use the reflective properties of scene and cloud screening to create accurate atmospheric and surface parameters in the level 1 images. The corrections included cirrus correction, scene classification, aerosol optical thickness retrieval, and water vapor retrieval. The corrected image disregards band B10, as it only stores information on cirrus clouds [67]. A bilinear interpolation-based technique was used to resample 20 m and 60 m resolution optical bands to 10 m spatial resolution. Particularly bands B5, B6, B7, B8a, B11, and B12 to align with bands B2, B3, B4 and B8, and radar image ( Table 5). The same city coordinates from Table 2 were used for the optical image sub setting. Bands B1 and B9 were eliminated because of their coarse spatial resolutions (originally 60 m) and aerosol and cloud sensitivity (ESA, 2015). There were ten bands in the final product (B2, B3, B4, B5, B6, B7, B8, B8a, B11, and B12.
Apart from same spatial resolution, the radar and optical images must be in same coordinate system and geometric extent to stack them into a single product. Therefore, while reprojecting radar image, a bilinear resampling method was used to align the radar image geometrically with optical image.

Land Cover Classification
After obtaining the processed images, the 10 optical bands, one band index, and texture and backscattering bands from radar data were stacked, and land cover maps were classified in R. A band index was developed using Equation (1) to highlight the built-up surfaces called normalized built-up index (NDBI).
Built-up surfaces have higher reflectance in band B11 according to [68]. Each image was loaded in R as a raster file, with bands as individual layers. Then using the "stack" command in R, ten bands and an index from the optical image were put together with a VV-polarization band and two texture measures from the radar image. This resulted in 14 bands, or layers, in a final raster stack image. Except for reflectance bands, sigma db, texture bands and band index were normalized before feeding them to the RF classifier. The resolution was preserved to 10 m in the final product. Such layer stacking is a form of pixel-based fusion technique for multisensory images. The following 14 layers were contained in the stacked image: Sentinel-1A VV amplitude image (1); Sentinel-1A VV variance texture image (1); Sentinel-1A VV contrast texture image (1); Sentinel-2A BOA reflectance bands 2 to 8a and 11 to 12 (10); and Sentinel-2A NDBI image (1).
Four land covers, vegetation, barren, water, and built up classifications were defined for each city; these were the dominant types of land cover in the cities. Built up cover included parking lots, road networks, buildings, and pavements, while vegetation included all tree cover, agricultural fields, fallow lands, and grasslands. Optical images alone, as well as stacked images, were then used to classify each city using the RF algorithm in R v.4.0.2.
RF uses bagging to create an ensemble of classification and regression tree (CART). It uses a random subset of variables each time to split at tree nodes to keep their intercorrelation minimal. They output the mode of class identified by each classifier tree [69,70]. Like a supervised classification approach, the RF model is trained with the land covers of interest across the remote sensing image scene [71][72][73].
Four RF models were prepared. The task included feeding land cover samples to train and test the model, defining RF parameters, and classifying the whole image. RF algorithms were developed by [74] and were implemented on the number of trees (ntree) and number of features in each split (mtry). The features are the 11 optical image bands and 14 stacked image bands. Generally, the square root of the total features present is used as the value of mtry, i.e., 4 in our case. A study conducted by [75] concluded that 500 is the optimum number of trees required for an RF classifier; thus 500 trees were used in this study. The land cover was assigned to each pixel based on the nearness of band values in the training set. Then each of 500 trees decide a particular cover for the pixel, and the most voted cover is the result.

Correlation Analysis
The resultant land cover maps were then reclassified into built-up and non-built-up surfaces, including vegetation, barren, and water cover, to emphasize the impervious surfaces of the cities. Then the resultant IS areas of each city from the fused data were correlated with the area of built-up surface in the global land cover maps for 2016 to assess their closeness. A graph was then plotted to compute the R squared value between the areas.

Test of Statistical Significance of Two Models
Two RF models were trained for each study year for each city. McNemar test was used to test statistical significance in extraction of built-up class for impervious surface mapping of the cities. It is applied to 2 × 2 contingency table with matched pairs of subject, correctly and incorrect classified built-up pixels by fused and optical data in this study. Then Chi-square values were computed for each city according to [76] as given in Equation (2). It compares the sensitivity and specificity of two datasets on land cover classification by finding out where the built-up pixels in classified maps disagree with each other [77].
where n f,o are number of samples correctly classified by RF model trained with fused data but incorrectly classified by RF model trained with optical alone data. Then the Chi-square value was checked for a threshold of 3.84 [76][77][78][79][80].

Impervious Surface Area and Expansion Rate Computation
Then pixel level area computation was performed to quantify the IS areas in each city using the fused data. The number of pixels identified as built-up class in the maps were multiplied by the spatial resolution of the images, i.e., 10 m, and converted into square km. Finally, the rate of urban expansion was computed for the past four years using Equation (3). The binary maps from 2015 were then overlaid on the maps from 2020 for further spatial analysis. Equation (3) is given by where, R is the annual growth rate, n is the temporal gap, and I(2016) and I(2020) are the impervious surface area in 2015 and 2020 respectively for each city.

Results
In this section results from the current study are explained. The processed S-2 and combined S-1 and S-2 data were used to prepare land cover maps of each city with four classes: barren, vegetation, built-up, and water. Therefore, this section begins with RF models' statistics and validation results of land cover maps followed by spectrum plot analysis. After establishing a premise for fused data results, the cities' IS areas and expansion rates are presented. Finally, a correlation plot between the results from current study and Global Land Cover Maps are discussed.

RF Model Statistics and Land Cover Map Validation
Each of the four RF models trained for each city showed more than 95% overall accuracy (OA) and a kappa coefficient rate with the test datasets. This result is produced by the RF model based on test set prediction after adapting to the training sets. A confusion matrix was produced each time to compute the OA and kappa coefficients. The results support the claim that RF classifiers are known for their classification accuracy in the RS community [81].
After the RF models predicted each city's whole image, the resultant land cover map was assessed with unique validation polygons. Validation samples were acquired for each land cover class. Automatic extraction or classification techniques need to be evaluated to develop confidence in the maps. Table 7 presents the accuracy assessment results for both optical and fused data. The accuracy evaluations were performed by computing the OA and kappa coefficients from the CM. OA is the ratio of the sum of diagonal elements and total elements in the CM. The kappa coefficient is the ratio of correctly predicted elements and total elements, considering the misclassified elements as well, from the CM. According to [82] a kappa coefficient gives a more robust accuracy assessment, as it considers both predicted classification and referenced classification data. A kappa value ranging from 0.8 to 1 is considered a highly accurate match with the reference values [83].  Table 8 lists the IS areas of the cities in 2016 from fused and global datasets. For further assessment of land cover maps, quantitative validation with CLS was performed. In Figure 4, the IS area from CLS is given on the y-axis, and the result from the current study is plotted on the x-axis. The minimum impervious area estimated was 14.8 km 2 in the smallest city Khanewal and the maximum estimated impervious area was 324.9 km 2 in the largest city Rawalpindi-Islamabad. The correlation of the estimated impervious area by the RF models and the CLS data showed a strong linear relationship with an R-squared value of 0.99. The RMSD calculated between the estimation and the CLS data was 12.8 km 2 which is 11.8% of the average of the estimated areas. From the plots, an underestimation of the IS was seen compared to the global data, as the slope of the point data was 1.0789, with an intercept of 0.36 km 2 . The average underestimation was 9.05 km 2 . The Pearson correlation test was conducted with a 95% confidence limit. Thus, based on the accuracies of the classified maps, strong linearity to independent datasets, and an analysis from the spectrum plots, IS area was quantified and urban expansion rates were computed from fused datasets only.

Test of Statistical Significance
Although, from Table 7, the improvement in OA and kappa coefficients were not significant in many cities, McNemar test results highlight that the classification of builtup pixels by two different RF models: trained with fused data and trained with optical data are statistically significant. The Chi square value obtained from the 2 × 2 confusion matrices of Rawalpindi-Islamabad, Multan, Faisalabad, Gujranwala, Bahawalpur, Sahiwal, Sheikhupura and Khanewal are given in Table 9. The Chi-square values are larger than or equal to 3.84 and p-values are less than 0.05 for all cities at 95% confidence limit. The result confirms the better performance of fused dataset over optical alone.  Table 10 shows the IS areas of all nine cities from the fused data of 2016 and 2020 and corresponding cumulative growth of the cities in the past four years. From the Table, it is seen that the urbanized surface of Rawalpindi-Islamabad grew from 324.9 km 2 in 2016 to 358.1 km 2 in 2020. The total growth of the city is 10.2%, with an annual rate of 2.5%. It has the highest growth in terms of km 2 and the 2nd highest growth rate in terms of percentage among the cities considered. RI is also the largest in the list in terms of area as well as population. Since Rawalpindi and Islamabad share the same border, they are jointly known as a twin city. They have a combined population of more than 3 million as of the 2017 census. The smallest city in the list, Khanewal, had urbanized surface growth from 14.6 km 2 to 14.9 km 2 with the lowest expansion rate. Additionally, the cumulative growth of Multan was 6.3%, Faisalabad was 9.7% Gujranwala was 8.1%, Bahawalpur was 11.8%, Sahiwal was 8.3%, and Sheikhupura was 9.3%. Overall, IS expansion was seen in all nine cities studied.   (Figure 5f) also had significant growth. Since the temporal gap between the analyses of urban expansion was just four years in our study, a thorough review of the growth pattern is challenging. In summary the most predominant type of urban growths were infill growth and edge expansion. In the former development of new built-up surface is surrounded by preexisting built up surface [84]. And in the latter newly developed built-up surface spreads out from the edge of an existing built-up surface [84] also known as edge development [85]. Infill is common to places with preexisting urban services such as water, road, and sewers according to [84].

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.

Methodology Transfer and Limitation
Fusion of S-1 and S-2 data as a single, stacked product requires both images to have same spatial resolution, coordinate reference system, and geometric alignment. However, to replicate the process in tropical regions would be challenging for impervious surfaces mapping, with S-2 being data constantly contaminated with cloud covers [45,64]. The temporal gap between images must not be large to avoid biases in the classification. Besides, availability of Sentinel data, starting in 2014-2015, prevents the user from long term assessment of impervious surfaces expansion. For a larger study area, the method can be computationally time consuming. The technique would mainly benefit urban areas with open spaces within the settlement due to higher spatial resolution of Sentinel data [26,27,[84][85][86] despite the urban density.

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.
Author Contributions: All three authors contributed equally to conceptualization of the research. B.S. prepared the methodology which included curating data, training the classifier and validating the classified maps. Substantial technical help was provided by H.S. for software operation while both H.S. and S.A. contributed equally to formal analysis and investigation of the result. Original draft was prepared by B.S. while the writing was reviewed and edited by H.S. and S.A. Manuscript was finalized under the supervision of S.A. and H.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Data Availability Statement:
The data that support the finding of this study are available from the corresponding author upon the reasonable request. Due to privacy concern the link to the drive where the data are deposited is not provided publicly.