A Modified Normalized Difference Impervious Surface Index (MNDISI) for Automatic Urban Mapping from Landsat Imagery

Impervious surface area (ISA) is a key factor for monitoring urban environment and land development. Automatic mapping of impervious surfaces has attracted growing attention in recent years. Spectral built-up indices are considered promising to map ISA distributions due to their easy, parameter-free implementations. This study explores the potentials of impervious surface indices for ISA mapping from Landsat imagery using a case study area in Boston, USA. A modified normalized difference impervious surface index (MNDISI) is proposed, and a Gaussian-based automatic threshold selection method is used to identify the optimal MNDISI threshold for delineating impervious surfaces from background features. To evaluate its effectiveness, comparison analysis is conducted between MNDISI and the original NDISI using Landsat images from three sensors (TM/ETM+/OLI-TIRS) acquired in four seasons. Our results suggest that built-up indices are sensitive to image seasonality, and summer is the best time phase for ISA mapping. With reduced uncertainties from automatic threshold selection, the MNDISI extracts impervious surfaces from all Landsat images in summer with an overall accuracy higher than 87% and an overall Kappa coefficient higher than 0.74. The proposed method is superior to previous index-based ISA mapping from the enhanced thermal integration and automatic threshold selection. The ISA maps from the TM, ETM+ and OLI-TIRS images are not significantly different. With enlarged data pool when all Landsat sensors are considered and automation of threshold selection proposed in this study, the MNDISI could be an effective built-up index for rapid and automatic ISA mapping at regional and global scales.


Introduction
Impervious surfaces are usually defined as anthropogenic features that do not allow water to penetrate through ground, e.g., building roofs, asphalt/cement roads, parking lots, sidewalks and transportation infrastructures [1].Impervious surface area (ISA) is not only an indicator of urbanization process, but also a key sensitive factor for monitoring environmental problems in built-up areas such as high surface runoff [2], urban heat island effect [3], transport of water pollutants [4], degradation in water quality [5,6] and air pollution [7].ISA only accounts for a small proportion of terrestrial lands [8].However, it represents Earth's fastest growing land uses that have dramatically altered climate, environmental condition and hydrology at local, regional, and global scales [9].Facing rapid urbanization all over the world, the above-mentioned concerns have triggered a surge of interest in impervious surface studies.
Due to the relatively low cost, large-area coverage and short revisit cycles, satellite remote sensing has been widely applied to map impervious surfaces.Among these studies, Landsat imagery is the most commonly used because the satellite series provide nearly 45-year data records with wide-swath coverage, free availability and relatively high spatial resolution (e.g., [10][11][12][13]).Various approaches have been employed to extract ISA from Landsat imagery, which can be grouped into five categories: (a) pixel/object-based classification [13][14][15]; (b) spectral mixture analysis (SMA) [16,17]; (c) regression model [18,19]; (d) decision tree [20,21]; and (e) spectral index-based segmentation [22,23].Although classification methods have been widely applied, there are challenges and limitations for applying at regional and global scales, e.g., subjective scene-to-scene data analysis, time consuming, and complicated computing [24].For pixel-based classification, it is difficult to resolve the mixed-pixel problems in Landsat imagery [13,14].For object-based classification, the difficulty in optimizing segmentation parameters hinders its application to a large area [13].The SMA-based methods have proven effective for handling the mixed-pixel problems.However, impervious surfaces are commonly overestimated in areas with low-density urban features and underestimated in high-density urban areas [14].The SMA-based methods are also not suitable for large-area mapping due to difficulties in endmember selection, inter-class variability quantification, and complicated implementation process [24,25].The regression models are effective in mapping national and even global ISA distributions [25].However, their challenges are to select suitable dependent variables from coarse resolution images and independent variables from high-quality ISA reference data.Decision tree methods are primarily rule-based such as the classification and regression tree (CART) algorithm [26].A decision tree can effectively process large, high dimensional and nonlinear data, which is suitable for large-area ISA mapping.However, it is more sensitive to noises than other algorithms and, therefore, relies heavily on the quality of sample data [14].
For large-area ISA mapping, spectral indices are considered promising methods due to their easy implementation, parameter-free and convenience in practical applications [23,24].A set of built-up indices have been developed, including the urban index (UI) [27], impervious surface area index [28], normalized difference built-up index (NDBI) [29], index-based built-up index (IBI) [30], exponential enhancing impervious surface method (P index) [31], normalized difference impervious surface index (NDISI) [22], enhanced built-up and bareness index (EBBI) [32], biophysical composition index (BCI) [24], nightlight-data modified NDISI [33], built-up and bare-land index (BBI) [34], normalized difference impervious index (NDII) [35], and combinational build-up index (CBI) [23].Challenges and limitations, however, still exist in the performance of these indices.For example, most of these indices, such as NDBI, BCI, and CBI, require masking out water bodies or bare soil before extracting ISA.While water could be easily delineated in satellite imagery, bare soil is often mixed with build-ups because of their spectral similarity [26].In addition, it is difficult to select an optimal threshold to segment ISA from background surfaces in these indices [25].Among the published indices, the NDISI takes advantage of thermal infrared (TIR) bands to account for distinctively different emissivity of impervious surfaces, which has distinct advantages for urban mapping.The TIR of all Landsat sensors, however, always has lower resolution than reflective bands of the same image acquisition, which casts a negative effect on identifying ISA from NDISI.
This study aims to address the above-mentioned problems in index-based ISA mapping.We propose a modified normalized difference impervious surface index (MNDISI) by integrating the sharpened TIR band into the NDISI, and develop an automatic threshold selection method based on Generalized Gaussian Model (GGM) to optimally extract ISA from Landsat imagery.Conducting a case study in Boston, Massachusetts, USA, we thoroughly test the applicability and seasonal sensitivity of this index by comparing it with the original NDISI using image sources from three Landsat sensors (TM, ETM+, and OLI-TIRS) acquired in four seasons.The paper is organized as follows.Section 3 describes the study area and datasets.The methodology of the MNDISI development and automatic threshold selection is in Section 3. Results are detailed in Section 4 and discussions are presented in Section 5. Finally, conclusions are summarized in Section 6.

Study Area
The study area covers one Landsat tile (12/31) in Boston, the capital city of Massachusetts, USA (Figure 1).As a historic municipal city, Boston hosts more than half million people living in an area of 233 km 2 [36].Boston is located on the east coast of the United States with its daily temperature in a range from −16.5 • C to 35.7 • C, having an average of 8.9 • C in spring, 21.7 • C in summer, 12.5 • C in autumn, and −0.13 • C in winter [37].The large dynamics of temperature and greenness between summer and winter is optimal for studying the seasonality impacts on index-based impervious surface extraction.
Remote Sens. 2017, 9, 942 3 of 18 follows.Section 3 describes the study area and datasets.The methodology of the MNDISI development and automatic threshold selection is in Section 3. Results are detailed in Section 4 and discussions are presented in Section 5. Finally, conclusions are summarized in Section 6.

Study Area
The study area covers one Landsat tile (12/31) in Boston, the capital city of Massachusetts, USA (Figure 1).As a historic municipal city, Boston hosts more than half million people living in an area of 233 km 2 [36].Boston is located on the east coast of the United States with its daily temperature in a range from −16.5 °C to 35.7 °C, having an average of 8.9 °C in spring, 21.7 °C in summer, 12.5 °C in autumn, and −0.13 °C in winter [37].The large dynamics of temperature and greenness between summer and winter is optimal for studying the seasonality impacts on index-based impervious surface extraction.

Datasets
Three image sets, acquired from Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI)-Thermal Infrared Sensor (TIRS), were collected in the study area.Four images for each sensor were downloaded to evaluate seasonal impacts on ISA mapping.Their acquisition dates, spectral bands and spatial resolutions are listed in Table 1.Atmospheric correction was performed to each image using the Atmospheric/Topographic Correction for Mountainous Terrain (ATCOR) software developed by German Aerospace Center, Wessling, Germany [38].The Landsat images archived in the U.S. Geological Survey (USGS) data clearinghouse have been geo-rectified [39].All images in Table 1 are geometrically matched to each other.
Training and validation samples were collected from the National Land Cover Database (NLCD) class maps on 2001 and 2011 downloaded from the USGS Data Clearinghouse.Only two classes were considered in this study, impervious surface (IS) and pervious surface (PS).In each NLCD class map, classes 22 (low intensity urban), 23 (medium intensity urban) and 24 (high intensity urban) were considered IS in this study.These NLCD classes represent land covers with

Datasets
Three image sets, acquired from Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI)-Thermal Infrared Sensor (TIRS), were collected in the study area.Four images for each sensor were downloaded to evaluate seasonal impacts on ISA mapping.Their acquisition dates, spectral bands and spatial resolutions are listed in Table 1.Atmospheric correction was performed to each image using the Atmospheric/Topographic Correction for Mountainous Terrain (ATCOR) software developed by German Aerospace Center, Wessling, Germany [38].The Landsat images archived in the U.S. Geological Survey (USGS) data clearinghouse have been geo-rectified [39].All images in Table 1 are geometrically matched to each other.
Training and validation samples were collected from the National Land Cover Database (NLCD) class maps on 2001 and 2011 downloaded from the USGS Data Clearinghouse.Only two classes were considered in this study, impervious surface (IS) and pervious surface (PS).In each NLCD class map, classes 22 (low intensity urban), 23 (medium intensity urban) and 24 (high intensity urban) were considered IS in this study.These NLCD classes represent land covers with 20-100% impervious surfaces [40].All other classes were considered PS.In this way, a reference IS/PS layer was exported from the NLCD class map.The sampling points for each class were randomly selected using a stratified random sampling method.Sample points in 2001 were used for analysis of the 2001 Landsat ETM+ images, and those in 2011 were used for the 2010-2011 TM images and the 2015-2016 OLI-TIRS images.For each class, both training and validation sample sets contained 20% of the class pixels that were randomly selected in the reference layer.

Methodology
This study modified a current built-up index-the normalized difference impervious surface index (NDISI) developed by Xu [22].By sharpening the TIR band assisted with reflective bands, the modified NDISI (MNDISI) was proposed.An automatic generalized Gaussian-based threshold selection method was then developed to identify the optimal threshold for segmenting IS from other land covers.To evaluate the performance of the MNDISI, its outputs were statistically compared with those from the original NDISI.The seasonal sensitivity of both indices in ISA mapping was further examined.

The Modified Normalized Difference Impervious Surface Index (MNDISI)
The NDISI is an effective method to enhance impervious surfaces and to suppress other land covers such as soil, sand, and water body.It is calculated as [22]: where ρ Green , ρ NIR , and ρ SWIR1 are surface reflectance of green, near-infrared, and the first shortwave-infrared band, respectively.T b is the brightness temperature of the thermal band as calculated in Chander et al. [41].The MNDWI is the modified normalized difference water index [42] to enhance water pixels with its calculation described in Equation (2).Both surface reflectance and T b in Equation ( 1) are stretched to [0, 255].For Landsat OLI-TIRS, they are stretched to [0, 65,535] to maintain the 16-bit radiometric resolution.For Landsat imagery, the pixel size of TIR band is much larger than 30 m of the reflective bands (as shown in Table 1).Here we apply the emissivity modulation theory to downscale the TIR image to 30 m pixel size.Land surface emissivity (ε) varies with land covers on ground surfaces.In urban environments, vegetated surfaces have stronger thermal holding capacity and therefore, have higher cooling effects than non-vegetated areas.Sobrino et al. [43] measured the relationships between ε and green covers on a variety of ground surfaces.Linking to Landsat-extracted NDVI, at each 30 m pixel, it is calculated as [43]: where ρ Red is surface reflectance of the red band; P V is the proportion of vegetation in a pixel (also referred to as fractional vegetation cover) with its calculation described in Equation ( 4); and NDV I min and NDV I max represent the minimum and maximum NDVI values in an image, respectively.A pixel is considered bare soil when its NDVI value is less than NDV I min , and full vegetation with its NDVI higher than NDV I max .Pixels with NDVI values in between are a mixture of bare soil and vegetation.For Landsat images acquired in peak growing season (summer), Sobrino et al. [43] suggested NDV I min = 0.2, and NDV I max = 0.5.For those acquired in other seasons, NDV I min is set 0.1-0.2, and NDV I max is 0.4-0.5 for images in winter, spring and autumn, respectively.Brightness temperature (T b ) in Equation ( 1) is extracted from a regular Landsat TIR band at coarser resolution.Contrarily, the ε in Equation ( 3) is extracted from the 30 m reflectance bands, and varies with green covers on ground.In this way, the emissivity-corrected land surface temperature (T s ) is sharpened to 30 m pixel size [44]: where T s is the 30 m land surface temperatures in unit of K; λ is the central wavelength of the thermal band; and ρ = 1.438 × 10 −2 (mK).
Using the sharpened T s in Equation ( 1), this study develops the modified NDISI (MNDISI) that is written as: Similar to Equation ( 1), all input parameters in Equation ( 6) are stretched to [0, 255] for TM and ETM+ and [0, 65,535] for OLI-TIRS.

Automatic Threshold Selection of the MNDISI
Impervious surfaces have higher MNDISI values than pervious surfaces.Figure 2 demonstrates an example MNDISI histogram derived from a TM image acquired on 30 August 2010, in which the IS pixels are located in the right side and PS pixels are in the left.The challenge of ISA mapping is thus to identify an optimal threshold in the MNDISI image.Manual selection of the threshold has been common in past index-based studies, which is time consuming, subjective and may highly affect the accuracies of ISA extraction.Here, we propose an automatic threshold selection method based on the generalized Gaussian model [45,46].Given a threshold value T (as marked in Figure 2), the generalized Gaussian optimization calculates the cost function (J) in the Kittler-Illingworth algorithm [47]: where ℎ( ) expresses the MNDISI histogram, with X varying from to , the minimum and maximum values of MNDISI incrementing at step l (l = 0.01 in this study).
, , and are positive constants related to the shape parameters ( and ) of the generalized Gaussian distribution for IS and PS pixels, respectively.The shape parameter determines the decay rate of the density function [48].and denote the mean MNDISI of IS and PS pixels that are separated at T. and are a priori probabilities associated with IS and PS, respectively.(Ω, ) is the entropy associated with the binary set of classes Ω = , at T. For detailed information about the generalized Gaussian model, please see Bazi et al. [45].
The optimal threshold * for segmenting IS and PS in the MNDISI histogram minimizes the cost function in Equation (7).Once * is determined, it assigns all MNDISI pixels to IS (MNDISI > * ) and PS (MNDISI < * ), respectively.The ISA map is thus extracted.

Performance Evaluation
With an independent reference datasets from the NLCD products, accuracy assessment was conducted using the commonly applied error matrix approach [49,50].The accuracy measures, i.e., overall accuracy (OA) and overall kappa coefficient (OK), were calculated.
To further evaluate the MNDISI performance, comparison analysis is conducted between the IS maps extracted from the MNDISI and the original NDISI.For each index, the spectral discrimination index (SDI) is calculated as [51,52]: where and express the mean values of an index for IS and PS, respectively, and and are the standard deviations for the two classes.
The SDI reveals the degree of separation between IS and PS in aspects of their relative locations and spreads in the histograms [24].A higher SDI indicates good separability.With the training data, the SDI values for both indices are compared to examine the effects of seasonality (spring, summer, autumn, and winter) and Landsat sensor types on ISA extraction.Here, we propose an automatic threshold selection method based on the generalized Gaussian model [45,46].Given a threshold value T (as marked in Figure 2), the generalized Gaussian optimization calculates the cost function (J) in the Kittler-Illingworth algorithm [47]: where h(X) expresses the MNDISI histogram, with X varying from L min to L max , the minimum and maximum values of MNDISI incrementing at step l (l = 0.01 in this study).a IS , b IS , a PS and b PS are positive constants related to the shape parameters (β IS and β PS ) of the generalized Gaussian distribution for IS and PS pixels, respectively.The shape parameter determines the decay rate of the density function [48].m IS and m PS denote the mean MNDISI of IS and PS pixels that are separated at T. P IS and P PS are a priori probabilities associated with IS and PS, respectively.H(Ω, T) is the entropy associated with the binary set of classes Ω = {ω PS , ω IS } at T. For detailed information about the generalized Gaussian model, please see Bazi et al. [45].
The optimal threshold T * for segmenting IS and PS in the MNDISI histogram minimizes the cost function in Equation (7).Once T * is determined, it assigns all MNDISI pixels to IS (MNDISI > T * ) and PS (MNDISI < T * ), respectively.The ISA map is thus extracted.

Performance Evaluation
With an independent reference datasets from the NLCD products, accuracy assessment was conducted using the commonly applied error matrix approach [49,50].The accuracy measures, i.e., overall accuracy (OA) and overall kappa coefficient (OK), were calculated.
To further evaluate the MNDISI performance, comparison analysis is conducted between the IS maps extracted from the MNDISI and the original NDISI.For each index, the spectral discrimination index (SDI) is calculated as [51,52]: where µ IS and µ PS express the mean values of an index for IS and PS, respectively, and σ IS and σ PS are the standard deviations for the two classes.
The SDI reveals the degree of separation between IS and PS in aspects of their relative locations and spreads in the histograms [24].A higher SDI indicates good separability.With the training data, the SDI values for both indices are compared to examine the effects of seasonality (spring, summer, autumn, and winter) and Landsat sensor types on ISA extraction.

The Sharpened TIR Imagery
The sharpened 30-m thermal bands are displayed for summer images of ETM+ (Figure 3), TM (Figure 4) and OLI-TIRS (Figure 5).Images in other seasons (as listed in Table 1) were also processed and the ISA distributions were extracted.To better demonstration, in this section we only visually show the results from summer images.The seasonal variations of ISA mapping from all images in this study are described in Section 4.3 below.
The thermal images, including the original and sharpened TIR images, are stretched to [0, 255] for TM and ETM+ and [0, 65,535] for OLI-TIRS.The sharpened TIR images obviously show more details of urban information.Spatial variations in built-up areas are better revealed.The sharpened TIR image in Figure 3 is from the 60 m ETM+ band; that in Figure 4 is from the 120 m TM band; and that in Figure 4 is from the 100 m TIRS band.In original thermal bands, road network information is often lost, hot spots of dense urban cores are broader, and vegetated residential areas have higher values due to mixed-pixel problems.All these problems lead to IS overestimation.The sharpened thermal band, on the other hand, better identifies urban areas and road networks for improved ISA mapping.

The Sharpened TIR Imagery
The sharpened 30-m thermal bands are displayed for summer images of ETM+ (Figure 3), TM (Figure 4) and OLI-TIRS (Figure 5).Images in other seasons (as listed in Table 1) were also processed and the ISA distributions were extracted.To better demonstration, in this section we only visually show the results from summer images.The seasonal variations of ISA mapping from all images in this study are described in Section 4.3 below.
The thermal images, including the original and sharpened TIR images, are stretched to [0, 255] for TM and ETM+ and [0, 65,535] for OLI-TIRS.The sharpened TIR images obviously show more details of urban information.Spatial variations in built-up areas are better revealed.The sharpened TIR image in Figure 3 is from the 60 m ETM+ band; that in Figure 4 is from the 120 m TM band; and that in Figure 4 is from the 100 m TIRS band.In original thermal bands, road network information is often lost, hot spots of dense urban cores are broader, and vegetated residential areas have higher values due to mixed-pixel problems.All these problems lead to IS overestimation.The sharpened thermal band, on the other hand, better identifies urban areas and road networks for improved ISA mapping.

Automatic Threshold Selection for ISA Mapping
With training samples for each image, one optimal threshold was identified to split all MNDISI pixels into IS and PS.Using three summer images as examples, Figure 6 shows the process of automatic threshold selection and its optimization effects for ETM+ (Figure 6a), TM (Figure 6c) and OLI (Figure 6e) sensors.The three peaks in all histograms represent high density IS pixels in the far right, low-density IS in the middle, and PS (all other land covers) in the left.The optimal thresholds, determined from the proposed automatic threshold selection, are −0.08,0.11 and −0.31 (as marked in the red lines), respectively.
Performance of these automatically selected thresholds was evaluated by comparing with randomly selected thresholds.At any given threshold, The IS and PS were segmented from the MNDISI image and the overall accuracy (OA) of the segmentation was calculated using the validation samples.For each image, an OA curve is traced with the threshold iteratively changing at a step of 0.01 (Figure 6b,d,f).For all three images, the automatically selected thresholds (marked in red line) are close to the one with the highest OA value (marked in black line) in all OA curves (Figure 6b,d,f).Because of the higher radiometric resolution of the OLI image, the OA difference between the automatic threshold and the peak threshold seems larger than those of TM and ETM+. Figure 6 indicates that our proposed automatic threshold selection is an effective method to identify the optimal threshold for ISA mapping.

Automatic Threshold Selection for ISA Mapping
With training samples for each image, one optimal threshold was identified to split all MNDISI pixels into IS and PS.Using three summer images as examples, Figure 6 shows the process of automatic threshold selection and its optimization effects for ETM+ (Figure 6a), TM (Figure 6c) and OLI (Figure 6e) sensors.The three peaks in all histograms represent high density IS pixels in the far right, low-density IS in the middle, and PS (all other land covers) in the left.The optimal thresholds, determined from the proposed automatic threshold selection, are −0.08,0.11 and −0.31 (as marked in the red lines), respectively.
Performance of these automatically selected thresholds was evaluated by comparing with randomly selected thresholds.At any given threshold, The IS and PS were segmented from the MNDISI image and the overall accuracy (OA) of the segmentation was calculated using the validation samples.For each image, an OA curve is traced with the threshold iteratively changing at a step of 0.01 (Figure 6b,d,f).For all three images, the automatically selected thresholds (marked in red line) are close to the one with the highest OA value (marked in black line) in all OA curves (Figure 6b,d,f).Because of the higher radiometric resolution of the OLI image, the OA difference between the automatic threshold and the peak threshold seems larger than those of TM and ETM+. Figure 6 indicates that our proposed automatic threshold selection is an effective method to identify the optimal threshold for ISA mapping.

Performance Evaluation
Comparison analysis was performed between the proposed MNDISI and the original NDISI using Landsat images from three sensors in four seasons.Automatic threshold selection was employed for both indices.The comparison included seasonality analysis, visual comparison, statistical measures, and accuracy assessment between the two indices.

Seasonality Analysis for ISA Mapping
The separability measures (SDI) of both indices from the three sensors and four seasons are plotted in Figure 7. Generally, IS and PS have higher separability for the MNDISI than the NDISI among all three sensor types: TM (Figure 7a), ETM+ (Figure 7b), and OLI-TIRS (Figure 7c).Seasonal trends of the SDI for the two indices are similar, all having best performance (SDI > 1.0) in summer when vegetation is in peak growth.In spring and autumn, green vegetation is less vigorous and is often exposed with bare soil.Spectral similarity between IS and PS therefore results in lower SDI values in these seasons.Limited greenness and occasional snow disturbance in winter attribute to the lowest SDI values in this season.For all three images in winter, the MNDISI has improved separability compared to NDISI, increasing from around 0.1-0.2 to 0.4-0.5.The relatively low SDI values, however, indicate that winter images cannot be effectively employed to extract impervious surfaces.
In short, Figure 7 suggests that summer imagery was optimal for ISA mapping from both indices.The comparison analyses in the following sections were only performed using the summer images.

Performance Evaluation
Comparison analysis was performed between the proposed MNDISI and the original NDISI using Landsat images from three sensors in four seasons.Automatic threshold selection was employed for both indices.The comparison included seasonality analysis, visual comparison, statistical measures, and accuracy assessment between the two indices.

Seasonality Analysis for ISA Mapping
The separability measures (SDI) of both indices from the three sensors and four seasons are plotted in Figure 7. Generally, IS and PS have higher separability for the MNDISI than the NDISI among all three sensor types: TM (Figure 7a), ETM+ (Figure 7b), and OLI-TIRS (Figure 7c).Seasonal trends of the SDI for the two indices are similar, all having best performance (SDI > 1.0) in summer when vegetation is in peak growth.In spring and autumn, green vegetation is less vigorous and is often exposed with bare soil.Spectral similarity between IS and PS therefore results in lower SDI values in these seasons.Limited greenness and occasional snow disturbance in winter attribute to the lowest SDI values in this season.For all three images in winter, the MNDISI has improved separability compared to NDISI, increasing from around 0.1-0.2 to 0.4-0.5.The relatively low SDI values, however, indicate that winter images cannot be effectively employed to extract impervious surfaces.
In short, Figure 7 suggests that summer imagery was optimal for ISA mapping from both indices.The comparison analyses in the following sections were only performed using the summer images.

Visual Comparison
Figure 8 demonstrates the visual differences of the extracted ISA distributions from the two indices of the TM summer image.Overall, the ISA distributions from both indices are similarly high in urban lands.In areas out of the urban core, the NDISI-extracted IS areas (Figure 8a) have fewer details than those from the MNDISI (Figure 8b), resulting in the loss of small-size impervious surface in these areas.In other words, the MNDISI reduces the omission errors in areas with lower urban densities.

Visual Comparison
Figure 8 demonstrates the visual differences of the extracted ISA distributions from the two indices of the TM summer image.Overall, the ISA distributions from both indices are similarly high in urban lands.In areas out of the urban core, the NDISI-extracted IS areas (Figure 8a) have fewer details than those from the MNDISI (Figure 8b), resulting in the loss of small-size impervious surface in these areas.In other words, the MNDISI reduces the omission errors in areas with lower urban densities.To better demonstrate the effectiveness of our proposed MNDISI, subsets with representative urban features are selected to compare the effectiveness of the two indices with summer images (Figure 9): central business district (CBD), residential area, road infrastructure, and suburban-rural area.One Landsat image type is used in a subset.In the CBD area (Figure 9I-a), one can observe that vegetation and impervious surfaces are often mixed.In the NLCD-extracted IS reference layer (Figure 9I-b), almost all pixels are classified as IS except a few vegetative clusters.Given the 60 m resolution of the ETM+ thermal band, both the NDISI (Figure 9I-c) and the MNDISI (Figure 9I-d) generate better ISA results than the NLCD reference.In residential areas (Figure 9II-a), neighborhood streets and building clusters are clear.The NLCD reference layer still overestimates the IS class (Figure 9 II-b).The NDISI can only obtain agglomerated urban clusters (Figure 9II-c).Due to the low resolution (100 m) of the OLI-TIRS thermal band, street network information is lost.With sharpened thermal image, the MNDISI fairly extracts both the agglomerated IS clusters and detailed road network in the neighborhood (Figure 9II-d).Road infrastructures such as highways are a significant component of urban impervious surfaces (Figure 9III-a).The linear feature is easily classified in the NLCD product (Figure 9III-b).The NDISI is unable to extract road information because of the low resolution (120 m) of the TM thermal band.With sharpened thermal image, the MNDISI fairly extracts the highway that is lost in the NDISI-extracted map.In suburban-rural areas, bare surfaces are observed in areas such as newly cleared construction fields and fallow agricultural lands (Figure 9IV-a).Bare fields are often confused with impervious surfaces in Landsat imagery due to their similar spectral characteristics, resulting in noisy, salt-and-pepper appearances in the NLCD reference map (Figure 9IV-b).The IS is better delineated from bare soil in the NDISI (Figure 9IV-c) and MNDISI (Figure 9IV-d  To better demonstrate the effectiveness of our proposed MNDISI, subsets with representative urban features are selected to compare the effectiveness of the two indices with summer images (Figure 9): central business district (CBD), residential area, road infrastructure, and suburban-rural area.One Landsat image type is used in a subset.In the CBD area (Figure 9I-a), one can observe that vegetation and impervious surfaces are often mixed.In the NLCD-extracted IS reference layer (Figure 9I-b), almost all pixels are classified as IS except a few vegetative clusters.Given the 60 m resolution of the ETM+ thermal band, both the NDISI (Figure 9I-c) and the MNDISI (Figure 9I-d

Statistical Comparison
Figure 10 compares the histograms of IS and PS for the NDISI and MNDISI images from the three sensors.Visually, the two classes can be delineated, although there is a large portion of overlay.More specifically, the histograms of IS derived from the MNDISI have two crests corresponding to low density IS and high density IS, respectively.Therefore, the MNDISI can completely identify high density IS pixels, and improve the separability between low density IS and PS pixels.In addition, for both classes, the MNDISI has wider data ranges than the NDISI, which also indicates the MNDISI's improved performance for separating impervious surface from other land covers.All SDI values are around 1.0 and above, indicating good separability between IS and PS classes from both built-up indices.For the same Landsat image, the SDI value of MNDISI-derived maps is higher than that of NDISI, although the SDI values from the OLI-TIRS are similar.All aspects above reveal that our proposed MNDISI has higher capability of extracting impervious surfaces from Landsat imagery.

Statistical Comparison
Figure 10 compares the histograms of IS and PS for the NDISI and MNDISI images from the three sensors.Visually, the two classes can be delineated, although there is a large portion of overlay.More specifically, the histograms of IS derived from the MNDISI have two crests corresponding to low density IS and high density IS, respectively.Therefore, the MNDISI can completely identify high density IS pixels, and improve the separability between low density IS and PS pixels.In addition, for both classes, the MNDISI has wider data ranges than the NDISI, which also indicates the MNDISI's improved performance for separating impervious surface from other land covers.All SDI values are around 1.0 and above, indicating good separability between IS and PS classes from both built-up indices.For the same Landsat image, the SDI value of MNDISI-derived maps is higher than that of NDISI, although the SDI values from the OLI-TIRS are similar.All aspects above reveal that our proposed MNDISI has higher capability of extracting impervious surfaces from Landsat imagery.

Accuracy Assessment
With validation samples, the two accuracy measures (OA and OK) are compared between the ISA outputs of the two indices with summer images from the three Landsat sensors (Table 2).For each index, the OA and OK derived from the three images are similar.The MNDISI improves the ISA extraction with higher accuracies.For the three sensors, the ETM+ image has the lowest increase of accuracies, followed by OLI-TIRS then TM.Coincidently, this sequence agrees with the spatial resolutions of their thermal bands, i.e., 60, 100 and 120 m, respectively.It should be noted, however, the improvements among these three sensors are not significantly large, which also confirms that Landsat imagery from the three sensors has similar performance in ISA mapping.

Discussion
This study aims to extract impervious surfaces from the built-up index with improved thermal integration and automatic threshold selection.Several innovative experiments are carried out.Compared with the original NDISI, our proposed MNDISI takes advantage of sharpened Landsat thermal channels based on emissivity modulation theories.Using reflective bands to refine

Accuracy Assessment
With validation samples, the two accuracy measures (OA and OK) are compared between the ISA outputs of the two indices with summer images from the three Landsat sensors (Table 2).For each index, the OA and OK derived from the three images are similar.The MNDISI improves the ISA extraction with higher accuracies.For the three sensors, the ETM+ image has the lowest increase of accuracies, followed by OLI-TIRS then TM.Coincidently, this sequence agrees with the spatial resolutions of their thermal bands, i.e., 60, 100 and 120 m, respectively.It should be noted, however, the improvements among these three sensors are not significantly large, which also confirms that Landsat imagery from the three sensors has similar performance in ISA mapping.

Discussion
This study aims to extract impervious surfaces from the built-up index with improved thermal integration and automatic threshold selection.Several innovative experiments are carried out.
Compared with the original NDISI, our proposed MNDISI takes advantage of sharpened Landsat thermal channels based on emissivity modulation theories.Using reflective bands to refine emissivity that varies with green covers, it enhances the identification of heterogeneous urban land covers, especially in low impervious surface areas.As shown in the histograms of impervious surfaces and pervious surfaces in Figure 10, the distributions of low-density IS pixels are overlaid with PS pixels.The MNDISI maximally stretches the high-density and low-density IS classes and is thus superior for ISA mapping.Furthermore, subjective threshold selection has been considered a key problem in index-based ISA mapping [25].In this study, a generalized Gaussian-based automatic threshold selection method is proposed to automatically identify the optimal threshold to segment impervious surfaces from other land covers in the MNDISI image.Automation of this process reduces the uncertainties from the subjective manual selection of thresholds that vary from scene to scene.The reduced time and labor cost also makes it more feasible for regional and even global ISA mapping with Landsat imagery.
This study contributes to the literature of ISA mapping with its comprehensive comparison analysis among different image sources (three Landsat sensors) and seasonality (four seasons), which has not been explicitly documented in past studies.The findings of this study suggest that satellite images acquired in summer are most suitable for ISA mapping, and those from all three Landsat sensors (TM/ETM+/OLI-TIRS) are good image sources.Although Landsat data collection is limited by its 16-day revisit cycle, and image quality is frequently contaminated by cloud cover, it is more reliable to collect optimal images in growing seasons from all Landsat sensors.In this sense, ISA mapping at regional even global scales also becomes practically feasible.
While built-up indices have been considered promising for large-area ISA mapping from optical satellite imagery at local, national, regional and global scales, our results indicate that there still exists limitations and challenges for index-based extraction of impervious surfaces.Firstly, the mixed-pixel problem remains a major challenge for ISA mapping using medium resolution Landsat imagery.Secondly, spectral similarity between different urban land covers always raises uncertainties in heterogeneous urban environments.For example, past studies have shown that bare soil can be easily confused with impervious surfaces [14].In the histograms (Figure 10) of this study, the overlaid areas between IS and PS clearly demonstrate the confusion between the brighter PS pixels (such as bare soil) and darker IS pixels (such as low-density urban).Although the MNDISI improved the separation by enhancing the IS values, the overlay remains unsolved.Further investigation is needed to tackle the problem of spectral similarity of urban land covers, especially impervious surfaces and bare soil.More recently, several indices specifically dealing with bare soil surfaces have been published [53][54][55].These bare soil indices could be coupled with built-up indices for improved ISA mapping.Finally, thermal image is utilized in our proposed MNDISI.However, the TIR band is not always available in the same image acquired from one sensor (e.g., Landsat Multispectral Scanner (MSS)).The feasibility of fusing thermal and reflective bands from different sensors will be tested in the future.
Other remote sensing categories have also been examined in mapping urban environments.For example, Synthetic Aperture Radar (SAR) imagery is found promising due to its all-weather conditions and strong double bounce scattering of urban features [56][57][58].In the future work, with the emergence of advanced satellite systems, index-based urban mapping with automatic threshold section could be a promising solution for rapidly and automatically mapping of ISA distribution at a regional and global scale.

Conclusions
This study proposes an improved impervious surface index, the MNDISI, by integrating the sharpened thermal band to reflective bands of Landsat imagery.An automatic threshold selection method is developed to optimally segment impervious surfaces from other land covers in the MNDISI image.Several findings are drawn: (1) Our results show that built-up indices are sensitive to seasonal changes in urban environments; and imagery acquired in summer is the best for ISA mapping.The three Landsat sensors, TM, ETM+ and OLITIRS, do not have significant differences in extracting impervious surfaces.(2) By downscaling thermal information using green cover-related emissivity, the proposed MNDISI enhances urban features especially, for low density impervious surfaces.The OA and OK values derived from Landsat imagery of all three sensors are higher than 87% and 74%, respectively.(3) With the rich set of Landsat imagery from its multiple sensors in the past 40 years, the proposed MNDISI together with automatic threshold selection approach could become an efficient tool for rapid mapping of impervious surfaces at regional and global scales, and therefore provides vital information to study broad impacts of human activities across the globe.

Figure 1 .
Figure 1.Study area and an example Landsat 8 OLI image acquired on 13 July 2016.

Figure 1 .
Figure 1.Study area and an example Landsat 8 OLI image acquired on 13 July 2016.

Figure 2 .
Figure 2.An example of automatically selected threshold (T) in the MNDISI histogram for segmenting impervious surface (IS) and pervious surface (PS) pixels.

Figure 2 .
Figure 2.An example of automatically selected threshold (T) in the MNDISI histogram for segmenting impervious surface (IS) and pervious surface (PS) pixels.

Figure 3 .
Figure 3. Sharpening of the ETM+ image acquired on 29 August 2001: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.An example subset (marked in the box) is shown in the second line of each column.

Figure 3 .
Figure 3. Sharpening of the ETM+ image acquired on 29 August 2001: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.An example subset (marked in the box) is shown in the second line of each column.

Figure 4 .
Figure 4. Sharpening of the TM image acquired on 30 August 2010: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.An example subset is shown in the second line of each column.

Figure 5 .
Figure 5. Sharpening of the OLI-TIRS image acquired on 13 July 2016: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.Two example subsets are shown in the second and third lines of each column.

Figure 4 .Figure 4 .
Figure 4. Sharpening of the TM image acquired on 30 August 2010: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.An example subset is shown in the second line of each column.

Figure 5 .
Figure 5. Sharpening of the OLI-TIRS image acquired on 13 July 2016: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.Two example subsets are shown in the second and third lines of each column.

Figure 5 .
Figure 5. Sharpening of the OLI-TIRS image acquired on 13 July 2016: the first column is the standard false color display; the second column represents the original thermal band; and the last column denotes the sharpened thermal band.Two example subsets are shown in the second and third lines of each column.

Figure 6 .
Figure 6.The MNDISI histograms and the overall accuracy (OA) curves for summer images of: ETM+ (a,b); TM (c,d); and OLI (e,f).Their optimal thresholds are marked in each histogram, and the OA value using this threshold (red line) is positioned in the corresponding OA curve.The black line represents the threshold reaching the highest accuracy of ISA mapping.

Figure 6 .
Figure 6.The MNDISI histograms and the overall accuracy (OA) curves for summer images of: ETM+ (a,b); TM (c,d); and OLI (e,f).Their optimal thresholds are marked in each histogram, and the OA value using this threshold (red line) is positioned in the corresponding OA curve.The black line represents the threshold reaching the highest accuracy of ISA mapping.
) indices.Similarly, the MNDISI results in more detailed IS distributions than the NDISI.
) generate better ISA results than the NLCD reference.In residential areas (Figure 9II-a), neighborhood streets and building clusters are clear.The NLCD reference layer still overestimates the IS class (Figure 9II-b).The NDISI can only obtain agglomerated urban clusters (Figure 9II-c).Due to the low resolution (100 m) of the OLI-TIRS thermal band, street network information is lost.With sharpened thermal image, the MNDISI fairly extracts both the agglomerated IS clusters and detailed road network in the neighborhood (Figure 9II-d).Road infrastructures such as highways are a significant component of urban impervious surfaces (Figure 9III-a).The linear feature is easily classified in the NLCD product (Figure 9III-b).The NDISI is unable to extract road information because of the low resolution (120 m) of the TM thermal band.With sharpened thermal image, the MNDISI fairly extracts the highway that is lost in the NDISI-extracted map.In suburban-rural areas, bare surfaces are observed in areas such as newly cleared construction fields and fallow agricultural lands (Figure 9IV-a).Bare fields are often confused with impervious surfaces in Landsat imagery due to their similar spectral characteristics, resulting in noisy, salt-and-pepper appearances in the NLCD reference map (Figure 9IV-b).The IS is better delineated from bare soil in the NDISI (Figure 9IV-c) and MNDISI (Figure 9IV-d) indices.Similarly, the MNDISI results in more detailed IS distributions than the NDISI.

Figure 9 .
Figure 9.The ISA maps from the summer images in the four subsets of the study area: CBD (I); residential area (II); highway (III); and suburban-rural area (IV).The rows give: the standard false color image displays (a); the NLCD-extracted IS reference layers (b); the NDISI-extracted ISA maps (c); and the MNDISI-extracted maps (d).

Figure 9 .
Figure 9.The ISA maps from the summer images in the four subsets of the study area: CBD (I); residential area (II); highway (III); and suburban-rural area (IV).The rows give: the standard false color image displays (a); the NLCD-extracted IS reference layers (b); the NDISI-extracted ISA maps (c); and the MNDISI-extracted maps (d).

Figure 10 .
Figure 10.The IS and PS histograms and SDI values of the NDISI (left column) and MNDISI (right column) from summer images of: ETM+ acquired on 29 August 2001 (the first row); TM acquired on 30 August 2010 (the second row); and OLI/TIRS acquired on 13 July 2016 (the third row).

Figure 10 .
Figure 10.The IS and PS histograms and SDI values of the NDISI (left column) and MNDISI (right column) from summer images of: ETM+ acquired on 29 August 2001 (the first row); TM acquired on 30 August 2010 (the second row); and OLI/TIRS acquired on 13 July 2016 (the third row).

Table 1 .
Satellite images from the three Landsat sensors (TM, ETM+, OLI-TIRS) used in this study.

Table 2 .
Accuracy comparison between the NDISI and MNDISI classifications with the three summer images.

Table 2 .
Accuracy comparison between the NDISI and MNDISI classifications with the three summer images.