Built-up area extraction using Landsat 8 OLI imagery

Abstract

The normalized difference built-up index (NDBI) has been useful for mapping urban built-up areas using Landsat Thematic Mapper (TM) data. The applicability of this index to the newer Landsat-8 Operational Land Imager (OLI) data was examined during this study, and a new method for built-up area extraction has been proposed. OLI imagery of urban areas of Lahore, Pakistan, was used to extract built-up areas through a modified NDBI approach and the proposed built-up area extraction method (BAEM). Instead of using individual bands, BAEM employed principal component analysis images of the highly correlated bands pertinent to NDBI computation. Through integration of temperature data, normalized difference vegetation index (NDVI) and modified normalized difference water index (MNDWI), BAEM was able to improve the overall accuracy of built-up area extraction by 11.84% compared to the modified NDBI approach. Rather than employing the binary NDBI, NDVI and MNDWI images, continuous images of these indices were used, and the final output was recoded by determining the threshold value through a double-window flexible pace search (DFPS) method. Results indicate that BAEM was more accurate at mapping urban built-up areas when applied to OLI imagery as compared to the modified NDBI approach; omission and commission errors were reduced by 75.96% and 33.36%, respectively. Moreover, the use of DFPS improved robustness of the proposed approach by enhancing user control over the segmentation of the output.

Keywords:

• built-up area extraction method
• normalized difference built-up index
• thermal infrared
• double-window flexible pace search
• Landsat 8 OLI

1. Introduction

Urban areas constitute just a small subset of the world’s overall land mass; however, population density and resource utilization intensity are very high in these areas as compared to their surroundings, compelling improved resource management practices (Cohen 2006; Lambin and Geist 2001). Understanding the spatial distribution and growth of urban areas is essential for urban planning and resource management, and one of the basic activities required for this purpose is mapping the built-up areas (Bertrand-Krajewski, Barraud, and Chocat 2000). Such mapping activity requires a lot of resources if carried out through conventional means like ground surveying and aerial photography (Richards 2013). The often rapid urban expansion also makes the task of a timely and accurate mapping of urban built-up areas quite difficult (Small 2003; Perepechko et al. 2005; Lein 2006).

Remote sensing (RS) data, especially from satellite remote sensing (SRS) systems, is an invaluable resource for mapping built-up areas for several reasons. It provides a synoptic and comprehensive view, which, in the case of large urban areas, is not possible through ground surveys (Richards 2013). Another practical advantage of utilizing RS data for urban studies is the availability of historical archives that can help in mapping and understanding urban sprawl over time (Guindon, Zhang, and Dillabaugh 2004; Maktav, Erbek, and Jürgens 2005; Xu 2008; Griffiths et al. 2010).

SRS data has been increasingly used for automated and semi-automated mapping of vegetation, snow, water, and other land-cover features (Lillesand, Kiefer, and Chipman 2004; Joseph 2005; Jensen 2006). A number of techniques have been formulated for this purpose, which can be grouped into two generic categories. The first one is based on the classification of input imagery that broadly encompasses pixel- and object-based classification methods (Guindon, Zhang, and Dillabaugh 2004; Cleve et al. 2008; Gao 2008), whereas the second type involves the direct segmentation of the imagery through indices (Zha, Gao, and Ni 2003; Zhang et al. 2005; Knight et al. 2006). Each method has its own set of advantages and limitations; however, indices have a certain edge over other classification methods in terms of time needed to generate results. A variety of indices have been developed for the extraction of features of interest from satellite imagery. The normalized difference vegetation index (NDVI) is the most commonly used for the extraction of vegetation. Other indices include the normalized difference snow index (NDSI), normalized difference water index (NDWI) and normalized difference built-up index (NDBI) for snow, water and built-up area extraction, respectively (Bannari et al. 1995; Hall, Riggs, and Salomonson 1995; Mcfeeters 1996; Zha, Gao, and Ni 2003). The calculations of these indices are based on the specific properties of the features of interest in terms of strong absorption or reflection in different spectral bands of multispectral imagery (Jensen 2006). The examination of the spectral signature of the object under observation helps in developing an index for its identification and extraction from an image (Gitelson and Merzlyak 1996; Huete and Jackson 1987; Dozier 1989).

The extent of urban land area is often linked to the impervious surface area that can be used to quantify and map urbanization (Bauer et al. 2004; Civco et al. 2002; Rosa and Wiesmann 2013). However, such mapping is still challenging due to spatial, spectral and temporal variability in built-up areas (Powell et al. 2007). To assist urban built-up area mapping, Zha, Gao, and Ni (2003) developed the NDBI on the footprints of the NDVI that used Landsat Thematic Mapper (TM) near-infrared (NIR) band 4 (low reflectance in built-up area) and mid-infrared (MIR) band 5 (high reflectance in built-up area). The output NDBI was further refined by removing vegetation noise using NDVI. Instead of using continuous raster images of NDBI and NDVI for final computation, this approach recoded them into binary images with the presumption that positive values of NDBI and NDVI represented built-up areas and vegetation, respectively. Although this recoding simplified the interpretation of the final output (positive values indicating built-up areas), Zha, Gao, and Ni (2003) mentioned that the usage of binary images restricted any refinement of results and due to this constraint, their approach was unable to separate urban from bare areas.

The shortcomings of the NDBI approach of Zha, Gao, and Ni (2003) in terms of refinement of the output were addressed by He et al. (2010) as they proposed a modified version of the conventional NDBI method. The core methodology remained the same, however, instead of using the binary images of NDBI and NDVI, this modified approach employed continuous images of both indices. The output was a continuous raster in which the pixels with higher values indicated a higher probability of them to represent built-up areas. The final segmentation into built-up and non-built-up classes was carried out using the double-window flexible pace search (DFPS) technique of Chen et al. (2003). Although not as straightforward as the customary NDBI method, this approach was comparatively better in terms of overall accuracy (He et al. 2010).

Since the basic foundation of indices discussed above is the spectral response of built-up areas in different bands of Landsat TM data, the output of these approaches is supposed to be different if applied to the latest Landsat-8 Operational Land Imager (OLI) data. The basic reason is that OLI uses narrower wavelength ranges to acquire data for NIR and MIR bands as compared to those in TM data (USGS 2013a). Considering the sensitivity of these indices to the wavelengths being used (Teillet, Staenz, and William 1997; Gitelson and Kaufman 1998; Maynard et al. 2007), it becomes important to examine the spectral response of built-up areas in different bands of OLI data before applying any TM-data derived method to it.

Spatial resolution is another important factor in extracting and mapping built-up areas (Lwin and Murayama 2013). A decent spatial resolution is required for these purposes (Welch 1982; Thomas, Hendrix, and Congalton 2003); however, limited coverage area and high cost preclude the utilization of high-resolution satellite data. Enhancing the spatial resolution of multispectral satellite data through a resolution merging process is a common and cost-effective way to increase the level of detail in satellite images with low spatial resolution. While enhancing the spatial characteristics, resolution merging methods should preserve the equally useful spectral properties of original multispectral data. Studies indicate that high-spatial resolution images produced using a high-pass filter (HPF) resolution merging method have the least spectral distortion (Showengerdt 1980; Wang et al. 2005; Ahmad and Singh 2002). Chavez, Sides, and Anderson (1991) compared different resolution merging methods (hue-intensity-saturation, principal component analysis (PCA) and HPF) through statistical, visual and graphical analyses, and found that the distortions in spectral properties of the output of the HPF method were the lowest.

When using indices for the extraction of desired information by processing the spectral bands of an image, it is important not to distort its spectral characteristics beforehand. In a study conducted by Jawak and Luis (2013), it was found that the NDVI values of high-spatial resolution imagery produced through resolution merging highly correlated with the ones in the original high-spectral resolution imagery, which indicated that resolution merging did not affect the spectral properties (less spectral distortion). In another study, it has been reported that the spectral characteristics were not altered when NDVI was applied to resolution-merged high-spatial resolution data (Acerbi-Junior et al. 2005). This study concluded that the NDVI profile of a resolution-merged image was very similar to the one derived from the original multispectral image. The successful application of NDVI to resolution-merged imagery has also been demonstrated in a study conducted by Walker et al. (2012). Another study points out that the quality of NDVI computed through resolution-merged imagery depended on the resolution ratio of the merged bands (Vrabel 2000). For instance, a few distortions were observed for smaller or linear features in NDVI image produced through resolution-merged image that was sharpened from 30 m to 1 m (30:1), but were not found in sharpening combinations such as 10:1 or 6:1. These examples indicate that spectral indices can be safely applied to resolution-merged imagery if the resolution ratio between the merged bands is small.

Another important characteristic of urban built-up areas is an increased land surface temperature which is directly related to urban intensity (Oke 1973; Weng 2001; Yuan and Bauer 2007; Lo 2004) and thus can be used as an additional indicator of urban areas for mapping. The response of a built-up area is somewhat different from other land-cover types in the thermal range of the electromagnetic spectrum (Zhangyan, Yunhao, and Jing 2006; Tran et al. 2006; Mallick, Kant, and Bharath 2008). Obviously, the urban heat island effect is basically responsible for higher temperature values of built-up areas as compared to their surroundings (Jusuf et al. 2007; Weng, Lu, and Schubring 2004; Kim 1992; Gartland 2010). This characteristic suggests that temperature can be used as an indicator to separate built-up from non-built-up areas.

The objectives of this study comprise: (1) examining the spectral response of built-up areas in different bands of OLI data to evaluate the applicability of previously developed built-up area extraction approaches; and (2) developing an improved approach that can be applied to Landsat-8 OLI data for urban built-up area extraction.

2. Study area and data

Lahore, the capital city of Pakistan’s Punjab province, was selected as the study area to develop this built-up area extraction method (BAEM). Located along the bank of the Ravi River towards the west and near the Indian–Pakistani border towards the east, the city has a total population of around 9 million (estimated in 2012) with about 82% residing in urban and the rest in peri-urban and rural areas (Bureau of Statistics 2012). Covering an area of around 1772 km² with an average surface elevation of about 217 m above mean sea level and characterized by a flat terrain, the city is administratively divided into nine towns and a cantonment (Figure 1(A)). For the purpose of this research, five urban towns, namely, Data Gunj Baksh, Gulberg, Ravi, Samanabad and Shalimar, covering an area of around 167 km², were selected as they comprised areas having built-up densities ranging from very high to very low (Figure 1(B)).

Figure 1. Study area: Lahore, Pakistan (A) towns in Lahore city and (B) towns selected for this study in Landsat-8 OLI false-colour composite imagery of 18 May 2013.

This research used Landsat-8 OLI imagery of Lahore (path: 149, row: 38) acquired on 18 May 2013. Spectral band 1 (visible, 0.43–0.45 µm) and band 9 (cirrus, 1.36–1.38 µm), which were not pertinent to this study, were excluded from further processing. The OLI optical bands 2–7 had a spatial resolution of 30 m, whereas that of the thermal bands 10 and 11 was 100 m which had already been downscaled to 30 m by the data vendor, Earth Resources Observation and Science Center (EROS), United States Geological Survey (USGS). Panchromatic band 8 had a spatial resolution of 15 m. Universal Transverse Mercator (UTM) projection with zone 43 North and WGS84 datum were used for the data-sets in this study.

3. Methodology

The proposed method for extraction of built-up areas using Landsat-8 OLI imagery comprised four major steps: preprocessing and examination of satellite data, image enhancement through resolution merging, development of the BAEM and accuracy assessment.

3.1. Preprocessing of Landsat-8 OLI imagery

The file format and data type of the Landsat-8 OLI imagery were tagged image file format (TIFF) and integer, respectively. All bands of the OLI imagery of Lahore (bands 2–8 and 10–11) were clipped to the study area boundary. Since the imagery used was cloud-free, similar to the situation mentioned by Deng and Wu (2013) in their study, no atmospheric corrections were performed. Digital numbers (DNs) of bands 2–8 were converted to top of atmosphere (ToA) reflectance using Equations (1) and (2), whereas the DNs of thermal bands 10 and 11 were first converted to ToA radiance (Equation (3)) and later processed to obtain at-satellite brightness temperature in degree Celsius (°C) through Equation (4) (USGS 2013b). The conversion parameter values were obtained from the respective metadata files downloaded along with the satellite data.

(1)

(2)

(3)

(4)

where:

ToA planetary reflectance, without correction for solar angle

Band-specific multiplicative rescaling factor from the metadata (REFLECTANCE_MULT_BAND_x, where x is the band number)

DNs of the band being processed

Band-specific additive rescaling factor from the metadata (REFLECTANCE_ADD_BAND_x, where x is the band number)

ToA planetary reflectance, with correction for solar angle

Local sun elevation angle. The scene centre sun elevation angle in degrees is provided in the metadata (SUN_ELEVATION)

ToA spectral radiance (Watts/(m² × srad × μm))

Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number)

Band-specific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number)

At-satellite brightness temperature (°C)

Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the band number, 10 or 11)

Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the band number, 10 or 11)

Bands 2–7 (excluding panchromatic band of 15 m spatial resolution) and bands 10 and 11 were stacked into two separate Erdas Imagine.img files. Unsigned integer was used as the data type for bands 2–7 (stacked) and panchromatic band 8.

3.2. Spatial enhancement and examination of data

A HPF resolution merging algorithm as proposed by Gangkofner, Pradhan, and Holcomb (2008) was employed to enhance the spatial resolution of spectrally rectified Landsat-8 OLI bands 2–7 (30 m spatial resolution), using panchromatic band 8 (15 m spatial resolution). First, the pixel sizes of both high-spatial resolution band and multispectral bands were interpreted to calculate the value of R, which is the ratio of cell size of multispectral data to that of the high-resolution data. This ratio was 2–1 in case of OLI, thus, a value of 2 was used for R. After that, a high-pass convolution filter kernel (HPK) was applied to filter the high-resolution input data which resulted in an HPF image. The HPK size, which is a function of the R value, was set to 5 × 5 for R = 2. The centre value of the kernel is also determined by the value of R, in this case, it was calculated as 24. The kernel used for producing the HPF image is shown in Figure 2.

Figure 2. High-pass convolution filter kernel.

Multispectral bands were then resampled to the pixel size of the high-resolution image (15 m) using a bilinear 4-nearest neighbour algorithm. Weighting was applied to the HPF image before adding it to the multispectral bands. The value of weight was determined by applying R and the standard deviations (SD) of HPF image and multispectral bands to Equation (5).

(5)

where W is the weight applied to the HPF image, SD(MS) is the standard deviation of the multispectral bands being enhanced, SD(HPF) is the standard deviation of the HPF image and M is the modulating factor determining the crispness of the output image, the value of which depends on R. For R = 2, the recommended value of M was 0.25. Each input multispectral band was then added to the weighted HPF image using Equation (6).

(6)

where x is the multispectral band being processed and W is the multiplicative weight computed by Equation (5). The last step of the HPF resolution merging algorithm was to apply a linear stretch process that rescaled the output DNs to match the mean and standard deviation of the respective input multispectral bands. The output multispectral bands had an improved spatial resolution of 15 m. Figure 3(A) and (B) illustrate the difference between non-enhanced and spatially enhanced images, respectively, at a sample location in the study area. The increase in sharpness and level of detail in the enhanced image are evident. In order to be consistent with the spatial resolution of multispectral bands, thermal bands 10 and 11 were resampled to a resolution of 15 m.

Figure 3. Landsat-8 OLI image in true-colour composite (A) before resolution merging (30 m spatial resolution) and (B) after resolution merging (15 m spatial resolution).

Spectral and thermal variations in built-up, vegetation and water areas were examined, at 10 sample locations each, in spatially enhanced multispectral bands 2–7 and thermal bands 10 and 11. The sample DNs of each of the above classes were separately averaged to plot the graphs shown in Figure 4(A) and (B).

Figure 4. Reflectance of built-up, vegetation and water areas in (A) optical bands 2–7 (DN value) and (B) thermal bands 10–11 (temperature in degree Celsius) of Landsat-8 OLI image.

3.3. Built-up area extraction approach

3.3.1. NDBI and modified NDBI approaches

In an effort towards automated mapping of urban areas using satellite imagery, Zha, Gao, and Ni (2003) proposed the NDBI. This index was primarily developed for Landsat TM satellite data through analysis of the spectral response of built-up areas in different bands of the imagery. The methodology comprised three arithmetic computations (Zha, Gao, and Ni 2003). First, a continuous NDVI imagery was obtained using Equation (7), which was then recoded into a binary image, NDVI_B, according to the conditional argument given in Equation (8).

(7)

(8)

Second, based on the high reflectance of built-up areas in the 1.55–1.75 µm wavelength range (TM band 5) and their low reflectance in the 0.76–0.90 µm wavelength range (TM band 4), NDBI was computed using Equation (9). A continuous NDBI imagery was obtained as a result which was then recoded into a binary image, NDBI_B, using the conditional argument shown in Equation (10).

(9)

(10)

Finally, the built-up areas were extracted using Equation (11). The resultant binary image, BU_B, showed positive DNs for built-up and bare areas, whereas the rest of land-cover types had a value of zero.

(11)

This methodology was modified and improved by He et al. (2010). The binary NDBI_B and NDVI_B in Equation (11) were replaced by the continuous NDBI (Equation (9)) and NDVI (Equation (7)), respectively, as shown in Equation (12).

(12)

In contrast to the binary output of the original technique proposed by Zha, Gao, and Ni (2003), a continuous image, BU, was produced through this modified approach in which a higher value of a pixel indicated a higher possibility that it indicated a built-up area. The output image was later segmented into a binary image (1 showing built-up and 0 showing non-built-up areas) using an optimal threshold value. A systematic approach based on the idea of the DFPS technique of Chen et al. (2003) was devised and used to determine the optimal threshold value (He et al. 2010).

3.3.2. Proposed BAEM

Zha, Gao, and Ni (2003) and He et al. (2010) found in their studies that TM band 4 (0.76–0.90 µm) and TM band 5 (1.55–1.75 µm) were most suitable for discriminating built-up areas from other land-cover types (Equation (9)). As the spectral ranges of Landsat-8 OLI bands differ from those of TM, it was necessary to study the spectral responses of built-up areas in different bands of OLI data before developing a method for extraction. Due to this, the spectral responses of built-up, vegetation and water areas in different bands of spatially enhanced OLI imagery (15 m resolution) of Lahore were examined. Band 5 (0.85–0.88 µm) and band 6 (1.57–1.65 µm) of OLI closely corresponded to band 4 and band 5 of TM, respectively. On examination, it was found that DNs of OLI band 6 highly correlated to those of band 7 of OLI (2.11–2.29 µm) (Figure 5(A)). The relationships between DNs of both bands were further studied for the pixels of built-up, vegetation and water area samples separately. The Pearson’s R values of 0.997, 0.982 and 0.998 for built-up, vegetation, and water area samples, respectively, indicated a very strong positive relationship between OLI bands 6 and 7 at a significance level of 0.01. For an optimal utilization of data from these two highly correlated bands, PCA was performed. PCA is a mathematical method that, while removing redundancy, preserves all the necessary information in correlated data-sets (Jolliffe 2002; Richards 2013). Each principal component (PC) is orthogonal to the other, and the first PC has the largest possible variance. The first PC image of bands 6 and 7, with the data type of unsigned integer, was used in the proposed method.

Figure 5. Scatter plots showing the correlation between (A) band 6 (x-axis) and band 7 (y-axis), and (B) band 10 (x-axis) and band 11 (y-axis) of the Landsat-8 OLI image.

The degree of imperviousness, which is directly related to urban density, has an effect on land surface temperatures (Oke 1973; Weng 2001; Yuan and Bauer 2007). This characteristic was found examining the thermal data. It was observed that the built-up areas had higher temperature values in thermal bands 10 and 11 of OLI compared to vegetation and water (Figure 4(B)). It was also found that band 10 highly correlated with band 11 (Figure 5(B)). A strong positive correlation between these bands was observed for the pixels of built-up, vegetation and water area samples with Pearson’s R values of 0.981, 0.996 and 0.994, respectively, at a significance level of 0.01. Therefore, the first PC image of bands 10 and 11 was used in the proposed method. The data type of the PCA output was set to unsigned integer to match that of the optical bands. To augment the high response of built-up areas in the PCA of bands 6 and 7, the PCA of bands 10 and 11 was added, and the NDBI formula of Zha, Gao, and Ni (2003) was modified as shown in Equation (13).

(13)

TM band 5 in the original NDBI formula was replaced by the sum of PCA of OLI bands 6 and 7 and PCA of OLI bands 10 and 11, while TM band 4 was swapped with OLI band 5 (Equations (9) and (13)). As in the method proposed by Zha, Gao, and Ni (2003) for built-up area extraction, NDVI_OLI was subtracted from NDBI_OLI. NDVI_OLI was computed by replacing TM band 3 and band 4 in Equation (7) by OLI band 4 and band 5, respectively.

The use of an NDWI has been suggested by Ogashawara and Bastos (2012) for reducing the spectral confusion between built-up and bare areas, arguing that bare areas (pervious surfaces) are likely to have a higher water content than built-up areas (impervious surfaces). The modified normalized difference water index (MNDWI) is the revised form of the NDWI and is reported to yield better results (Xu 2005; Chao and Sheng 2001). Hence, to remove water signatures from the output image and increase the spectral contrast among built-up area, bare land and other land-cover classes, an MNDWI is introduced in the proposed method. This MNDWI processes green and mid-infrared spectral bands to extract water areas (Xu 2005). For OLI, MNDWI_OLI was computed using Equation (14).

(14)

The MNDWI_OLI was then subtracted from NDBI_OLI to obtain an output image with water-free areas. Correlations between NDBI_OLI and NDVI_OLI, NDBI_OLI and MNDWI_OLI, and NDVI_OLI and MNDWI_OLI were also tested to examine the spectral distribution of built-up, vegetation and water areas in these data-sets. Finally, continuous NDBI_OLI, NDVI_OLI and MNDWI_OLI images were used to extract a continuous built-up area image, BAEM_OLI (Equation (15)). Higher DN values in BAEM_OLI indicated a higher possibility of those pixels to represent built-up areas, whereas lower values depicted land-cover classes other than built-up.

(15)

Histogram segmentation was performed using a threshold value to segregate built-up from non-built-up areas in the BAEM_OLI image. DFPS was used to determine the optimal threshold value for this segmentation (Chen et al. 2003). Before application, the data type of the BAEM_OLI image was converted to unsigned integer. The DFPS approach was developed by Chen et al. (2003) basically to determine the optimal threshold value between change and non-change pixels while examining land use/cover transitions. However, this method has also been used to find the threshold value for segmenting a thematic imagery to extract built-up areas (He et al. 2010). The process involving three steps was systematically applied to the BAEM_OLI imagery.

First, sample built-up areas were chosen through visual interpretation of the input imagery, subject to the condition that they should only include built-up pixels and had to be ‘built-up islands’ surrounded by non-built-up pixels. Two visual areas were considered here, an inner window consisting of built-up area pixels only, and an outer frame surrounding the inner window and containing the pixels indicating non-built-up areas. In the second step, the histogram of the BAEM_OLI imagery was examined to find out range and pace. The first search pace (increment) P₁ was calculated using Equation (16), whereby the search range was determined by the difference between the minimum (a) and maximum (b) values of BAEM_OLI.

(16)

where m, which represented the number of possible threshold values in the search process, was decided manually. The potential thresholds were assigned within the range of (a, b) as b-P₁, b-2P₁, …. Threshold values were then processed one by one to determine the success rate of built-up area extraction. For a potential threshold value of k, the success rate L_k was computed using Equation (17).

(17)

where A is the total number of pixels in the sample (within the non-built-up frame area), is the number of pixels detected as built-up inside the inner window and is the number of pixels wrongly detected as built-up in A. For all m thresholds in one particular search process, the maximum and minimum values of L_k were calculated and assigned to L_max and L_min, respectively. These two values were examined whether or not they satisfied the exit condition given in the final step of the DFPS process (Equation (18)). If the exit condition was not satisfied, a new search began with a new search range (k_max – P₁, K_max + P₁) and a new search pace that was calculated using Equation (16), where k_max was the threshold value corresponding to L_max in that particular search iteration.

The final step, the termination of the iterative process of steps two and three, was performed when L_max and L_min satisfied the condition given in Equation (18).

(18)

where δ is the acceptable error constant. Its value in this case was 1 since the image used was of an integer, and not floating point, data type. The optimal threshold value to separate built-up from non-built-up areas in the BAEM_OLI image was the one corresponding to L_max in the final iteration of the DFPS process. Using this threshold value, the BAEM_OLI image was finally categorized into two classes, built-up (pixel value = 1) and non-built-up (pixel value = 0).

3.4. Accuracy assessment

The performance of the proposed BAEM was tested through accuracy assessment at 200 locations (pixels) using a stratified random sampling technique (Congalton 1991). Two strata, built-up and non-built-up, were formed, and random samples were taken from each of them based on their proportion in the output image. In addition, the NDBI approach proposed by He et al. (2010) (Section 3.3.1) was also used to generate a map of built-up areas, which was tested using the same accuracy assessment method in order to compare the performance of the two approaches. Producer and user accuracies of both approaches were compared to assess their efficiency in terms of segregating built-up areas from other land covers/uses. Different areas were visually examined in detail to understand the differences in the outputs of both approaches.

4. Results and discussion

4.1. Examination of Landsat-8 OLI data

Spectral responses of built-up, vegetation and water areas were examined in OLI multispectral bands 2–7 and thermal bands 10 and 11. All three land-cover types had a unique signature in optical bands which indicated that they can be segregated from each other through proper application of their respective indices (Figure 4(A)). Compared to the other two land-cover types, the spectral reflectance curve of vegetation had a prominent slope between bands 5 and 4. This indicated that the NDVI could be more efficient in terms of accuracy compared to NDBI or MNDWI. Another relevant finding was the difference between temperatures of the three land-cover types (Figure 4(B)). Both thermal bands 10 and 11 exhibited the highest temperature values for built-up areas (around 50°C), followed by vegetation (around 46°C) and water (around 40°C). The primary reason for this characteristic of built-up areas is the usage of construction materials that absorb heat during daytime and release it at night (Gartland 2010; Goward 1981). Compared to other land-cover types, the cooling process in construction materials is very slow which subsequently increases the overall temperature in built-up areas. Higher temperature readings of built-up areas in the thermal bands implied that temperature can be used as an aiding factor to separate built-up areas from vegetation and water.

4.2. Relationships among NDBI, NDVI and MNDWI

The NDBI_OLI imagery was derived through Equation (13). As illustrated in Figure 6(A), pixels with white and bright grey tones show built-up areas, whereas light grey, dark grey and black tones show areas with other land-cover types. Densely vegetated areas appear dark in this image, whereas water and a few areas with less dense vegetation are shown in bright tones. However, for extracting the built-up areas, both vegetation and water were later filtered using NDVI and MNDWI, respectively. The NDVI_OLI image is presented in Figure 6(B) in which the pixels with white and bright grey tones show vegetation, moderate to dark grey tones show other land covers, and black is associated with water. A very clear separation between vegetation and other land-cover types is evident in this output. A considerable difference in reflectance of vegetation in OLI bands 4 and 5 was the basic reason for this characteristic of the NDVI_OLI image (Figure 4(A)). The MNDWI_OLI imagery generated using Equation (14) is illustrated in Figure 6(C). Water is shown as white pixels in this image, whereas other land-cover types are depicted in the tones of grey and black. Areas with dense vegetation appear dark whereas the built-up areas are shown in variable tones of grey.

Figure 6. The resultant images of (A) NDBI, (B) NDVI and (C) MNDWI derived from the Landsat-8 OLI image.

The NDBI_OLI, NDVI_OLI and MNDWI_OLI images, apparently looking different from one another, were examined in detail for any interrelationship. It was found that NDBI_OLI had a weak positive correlation with NDVI_OLI, whereas a negative correlation was observed between NDBI_OLI and MNDWI_OLI (Figure 7(A) and (B)). NDVI_OLI also negatively correlated with MNDWI_OLI (Figure 7(C)) implying that these data-sets showed different features. The scatter plots shown in Figure 7 also indicate a clear separation in spectral responses of built-up, vegetation and water areas in NDBI, NDVI and MNDWI imageries. The Pearson’s R values revealing the degree of relationship between NDBI_OLI, NDVI_OLI and MNDWI_OLI images are shown in Table 1. These values and the scatter plots confirmed the divergence of NDBI_OLI, NDVI_OLI and MNDWI_OLI images from one another signifying that these data-sets are ready for further processing for built-up area extraction.

Figure 7. Scatter plots showing the correlation between (A) NDBI (x-axis) and NDVI (y-axis), (B) NDBI (x-axis) and MNDWI (y-axis), and (C) NDVI (x-axis) and MNDWI (y-axis).

Table 1. The Pearson’s R values showing the degree of relationship among NDBI, NDVI and MNDWI.

	NDBI	NDVI	MNDWI
NDBI	1	–	–
NDVI	0.130*	1	–
MNDWI	−0.298*	−0.770*	1

Note: *correlation is significant at the level of 0.01.

4.3. Built-up area map and its accuracy

A continuous raster of the built-up area, BAEM_OLI, was derived from Equation (15) through processing the NDBI_OLI, NDVI_OLI and MNDWI_OLI images. The output is illustrated in Figure 8(A) where the pixels in white, bright grey and light grey tones show built-up areas, whereas the dark grey and black tones represent other land covers. The variability in the bright tones indicates the differences in built-up density, where the areas with high built-up density appear the brightest and vice versa. In order to determine the optimum DN value to segregate built-up from non-built-up areas in the BAEM_OLI image, the DFPS method was applied. The image data type was first converted to unsigned integer, and then its histogram was examined to set the initial search range. The search process is shown in Table 2 showing that the search range changed five times with the pace of 50, 20, 8, 2 and 1. In total, 23 threshold values were tested to determine the optimum, which was identified as 83 with a success rate of 93.72%. Pixels with values above 83 were extracted and classified as built-up areas in the BAEM_OLI image (Figure 8(B)).

Figure 8. Map showing (A) a continuous image of built-up areas derived by BAEM, and (B) segmented output of BAEM (built-up area extracted through segmentation of (A) using the threshold value from DFPS).

Table 2. Results of the double-window flexible pace search for the determination of the optimal threshold value.

Range: 255–5 Pace 50		Range: 105–5 Pace 20		Range: 85–45 Pace 8		Range: 85–69 Pace 2		Range: 85–81 Pace 1
Threshold	Success rate	Threshold	Success rate	Threshold	Success rate	Threshold	Success rate	Threshold	Success rate
255	0	105	0	85	76.92	85	76.92	85	76.92
205	0	85	76.92	77	38.46	83	93.72	84	92.31
155	0	65	26.92	69	7.36	81	84.61	83	93.72
105	0	45	3.27	61	7.12	79	65.38	82	88.46
55	3.27	25	3.27	53	3.27	77	38.46	81	84.61
5	3.27	5	3.27	45	3.27	75	19.23
						73	9.61
						71	8.23
						69	7.36

The accuracy of the proposed BAEM and the approach suggested by He et al. (2010) were assessed using a stratified random sampling technique. Results of BAEM were 80.50% accurate, whereas those derived by the method of He et al. (2010) returned an accuracy of 71.50% (Table 3). The commission error was higher than the omission error in both approaches signifying that both outputs contained some areas that were wrongly extracted as built-up; however, the low omission error indicated that both approaches were able to extract built-up areas quite well (Table 4). Nevertheless, its low values of omission and commission errors and the high k-statistic value suggested that the BAEM was more accurate than the established method.

Table 3. Error matrix for NDBI method and BAEM.

Classified data	Reference data
(NDBI method)	Non-built-up area	Built-up area	Classified total	User’s accuracy (%)
Non-built-up area	44	2	46	95.65
Built-up area	55	99	154	64.29
Reference total	99	101	200
Producer’s accuracy (%)	44.44	98.02
Overall accuracy NDBI method = 143/200 (71.50%)
Classified data (BAEM)
Non-built-up area	50	1	51	98.04
Built-up area	38	111	149	74.50
Reference total	88	112	200
Producer’s accuracy (%)	56.82	99.11
Overall accuracy BAEM = 161/200 (80.50%)

Table 4. Summary of the accuracy assessment of the images derived by NDBI method and BAEM.

	Overall accuracy (%)	Omission error (%)	Commission error (%)	k coefficient
NDBI method	71.50	1.98	35.71	0.43
BAEM	80.50	0.89	25.5	0.59

For comparison, the input OLI image and the outputs of the two built-up area extraction approaches showing the same sample location are displayed in Figure 9. It was found that the method of He et al. (2010) over-classified built-up areas as compared to the BAEM (Figure 9). This finding conforms to the higher commission error in the previous approach as compared to the BAEM (Table 4). Besides, the use of temperature data in the proposed method helped in segregating less dense built-up areas from other land-cover types, thus improving the overall accuracy. Vegetation and water were very well separated from built-up areas indicating that both NDVI and MNDWI were correctly processed in the BAEM. Thus, the reason for the high commission error may be attributed to spectral mixing of built-up areas with land-cover types other than vegetation or water.

Figure 9. Comparison of the results from the NDBI approach and BAEM at a sample location (A) OLI image in false-colour composite, (B) output from NDBI approach and (C) output from BAEM. (Black rectangles in each panel show the places where the NDBI approach wrongly extracted other land covers as built-up areas.)

Another interesting finding was the difference between the outputs of the method suggested by He et al. (2010) and BAEM in terms of segregating water from bare land. It was found that the proposed method was better at extracting water as non-built-up area in comparison to the conventional method (Figure 10), which can be attributed to the inclusion of MNDWI in the BAEM. Moreover, both approaches could not completely separate bare land from built-up areas because of the resemblance between spectral characteristics of these two land-cover types (Herold et al. 2004; Lu and Weng 2004; Lwin and Murayama 2013). However, due to the use of thermal-infrared band, MNDWI and DFPS, the BAEM was somewhat better at this segregation which was evident from lower commission error values as compared to the previous approach.

Figure 10. (A) OLI image in false-colour composite, (B) output from NDBI approach and (C) output from BAEM. (Black rectangles in each panel show the places that produced obvious commission errors in both methods.)

5. Conclusions

Both NDBI and modified NDBI methods were able to extract built-up areas with reasonable accuracy when applied to Landsat TM data (Zha, Gao, and Ni 2003; He et al. 2010). However, this study found that the accuracy of the NDBI approach of He et al. (2010) decreased when applied to Landsat-8 OLI data of urban areas of Lahore, Pakistan. This aspect was explained by the differences between spectral wavelength ranges of the different bands of the two data-sets. Therefore, a new method, BAEM, is proposed in this study for extracting built-up areas using OLI data employing NDBI, NDVI and MNDWI. In addition to the modification made to the established NDBI approach by applying PCA image of OLI bands 6 and 7 instead of just band 6, PCA image of thermal bands 10 and 11 was also integrated in the BAEM. The reason for using PCA images was the significant correlations between bands 6 and 7, and bands 10 and 11 of OLI. The use of PCA images consequently helped in capturing the variance and comprehending a broader range of useful data. Since higher temperature values were observed in built-up areas compared to other land-cover types, inclusion of thermal data in BAEM resulted in better separation of built-up from non-built-up areas. Weak correlations among NDBI, NDVI and MNDWI indicated a reduced spectral confusion and produced a better contrast among different land-cover types in the output. The DFPS method was quite efficient at determining the optimal threshold value for segmenting the BAEM output to segregate built-up areas from other land-cover types. Although the use of DFPS reduced the element of automation, the flexibility of this method can help in achieving improved results when applying the BAEM in areas with topographic conditions different from this study area. Another advantageous option could be the double application of DFPS to identify two threshold values, one for separating high-density built-up areas from low-density ones, and the other for separating low-density built-up areas from non-built-up ones. The results of the proposed method revealed a higher accuracy than the one suggested by He et al. (2010); however, it was found that both methods could not completely separate built-up areas from bare land. Since the urban heat island effect is more significant at night when the temperature difference between built-up areas and other land covers is greater (Oke 1982), integrating night-time thermal data in the proposed method could yield better results in separating built-up from bare areas. Night-time remotely sensed light data, which has been helpful in mapping the spatial extent of urban areas (Voogt and Oke 2003), can also be incorporated to improve the overall accuracy of BAEM in future studies. Further application and analysis of the proposed method on data from other satellites is suggested.

Acknowledgements

The authors gratefully acknowledge the support from the Asian Institute of Technology, Thailand, and the Japanese Government for carrying out this research. We are thankful to the Earth Resources Observation and Science Center (EROS), United States Geological Survey (USGS), for providing Landsat-8 OLI data free of charge for this study. The authors would also like to thank the reviewers for their insightful comments and valuable suggestions.