Performances of WorldView 3, Sentinel 2, and Landsat 8 data in mapping impervious surface

Remote Sensing Applications: Society and Environment 15 (2019) 100246

Contents lists available at ScienceDirect

Remote Sensing Applications: Society and Environment journal homepage: www.elsevier.com/locate/rsase

Performances of WorldView 3, Sentinel 2, and Landsat 8 data in mapping impervious surface

T

George Xiana,∗, Hua Shib, Jon Dewitza, Zhuoting Wuc a

U.S. Geological Survey (USGS), Earth Resources Observation and Science Center (EROS), Sioux Falls, SD, 57198, USA ASRC Federal ASDP, LLC, Contractor to the USGS EROS, Sioux Falls, SD, 57198, USA c USGS National Land Imaging Program Flagstaff, AZ, 86001, USA b

A R T I C LE I N FO

A B S T R A C T

Keywords: Impervious surface Landsat Sentinel WorldView Urban

Many efforts have been made to map developed impervious surface from remotely sensed information in the last two decades. The U.S. Geological Survey (USGS) developed the National Land Cover Database (NLCD) to provide consistent land cover and change products for the Nation since 2001. Percent impervious surface area (ISA), one of the products in NLCD as a continuous field and estimated with Landsat imagery, represents the fraction of human-made impervious area in a 30 m resolution grid. ISA is used to map urban land cover types and extents for the United States. However, it is still a challenge to quantify highly heterogeneous features in many urban areas using remotely sensed data with spatial and spectral resolutions similar to Landsat and to determine the impacts of remotely sensed data characteristics on ISA mapping.

In this study, we explored using remotely sensed data that have different spatial and spectral resolutions to evaluate their performances in quantifying developed impervious surfaces as a continuous field. Two urban areas were selected as study areas: San Francisco, California, and Dallas, Texas. Satellite images from WorldView-3, Sentinel-2, and Landsat 8 were employed to quantify developed ISAs in the three areas. We compared ISAs estimated with images from these sensors and evaluated the benefits and limitations of spectral and spatial resolutions for mapping ISA. Our study found that WorldView-3 images, especially in the 20 m estimate, outperformed the other two sensors with relatively smaller errors in imperviousness estimation in urban areas that have relatively dry and sparse vegetation cover. Due to the relatively coarse spatial resolution, both area and intensity of imperviousness estimated from Landsat images were relatively lower than from the other two sensors. Comparing ISAs estimated from WorldView and Sentinel images, the underestimate of ISA in high intensity urban centers by Landsat imagery cannot be ignored. 1. Introduction The world population is increasing, and since 2008 more than half of the world's population lives in urban centers (IPCC 2013). Associated with population growth is the expansion of urban land use and land cover around the world. Urban development usually creates various areas of non-evaporating and non-transpiring impervious surfaces ∗

composed of concrete, asphalt, stone, and metal. Developed impervious surface area (ISA) alters vegetation and other land conditions in terms of nature, abundance, pattern, and biodiversity (Cadenasso et al., 2007; Piano et al., 2017; Shi et al., 2005; Xian and Homer 2010). The large human population and intensive commercial, industrial, and transport activities in many urban areas have brought profound changes in environment and landscape structure (Grimm et al., 2008; Li and Zhou 2017; Piano et al., 2017). These altered environmental conditions exist in most urban areas and could have a wide range of environmental challenges resulting in different ecological and environmental consequences. For example, urban runoff, which is highly related to the intensity of ISA, directly affects the hydrological cycle and degrades water quality by loading pollutants from urban centers to streams (Berezowski et al., 2012; Mouri et al., 2012; Wang et al., 2011; Xian et al., 2007). Developed ISA substantially impacts natural landscape characteristics and has been recognized as a prominent stressor for ecological conditions in watersheds. Many previous studies also used ISA to connect to the urban heat island effect and showed influences of ISA on urban land surface temperature change (Dousset and Gourmelon 2003; Fu and Weng 2018; Gallo and Xian 2016; Giovannini et al., 2011; Imhoff et al., 2010; Li et al. 2018a, 2018b; Weng et al., 2004; Xian 2008; Yuan and Bauer 2007). Accurate information about the spatial extent of ISA and its temporal change is critical for urban planning and management to assess variations of urban land cover and associated ecological and environmental effects.

Corresponding author. U.S. Geological Survey (USGS), Earth Resources Observation and Science Center (EROS), Sioux Falls, SD, 57198, USA. E-mail address: [email protected] (G. Xian).

https://doi.org/10.1016/j.rsase.2019.100246 Received 6 February 2019; Received in revised form 5 June 2019; Accepted 11 June 2019 Available online 13 June 2019 2352-9385/ Published by Elsevier B.V.

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

the continuous field data estimate approach. It is poorly understood how different spatial resolutions and spectral information affect ISA mapping performance in terms of spatial distribution and cover intensity. Recent advances in remote sensing technology have made available multiple spectral and spatial resolution data from several satellites including WorldView-3 (WV3), Sentinel-2, and Landsat 8. These datasets have different spectral and spatial resolutions and provide the potential to derive information about the nature and properties of different materials on the urban ground at sub-pixel level (Benedetti et al., 2018; Gong et al., 2019; Shi et al., 2017; Sun et al., 2019; Tavares et al., 2019). Comparing urban impervious surface as a continuous field estimate derived from different sensors will improve our understanding of capabilities and advantages of these data for urban land mapping. The objectives of this study are to 1) systematically analyze performances of remotely sensed images from sensors that have different spectral and spatial resolutions for ISA mapping; 2) evaluate the role of different spectral bands in modeling ISA; and 3) evaluate accuracies of total ISA estimated from continuous field imperviousness. We used the same regression tree modeling algorithms as used for the U.S. Geological Survey (USGS) National Land Cover Database (NLCD) impervious surface mapping (Xian and Homer 2010) to map ISA distributions. Images from WV3, Sentinel-2, and Landsat 8 were used separately to quantify ISA at different spatial scales. We assessed ISAs mapped from different images using different statistic metrics and analyzed which spectral bands in different images are most frequently used for mapping developed ISA based on model output performance parameters.

Developed ISA is highly associated with the urban land use condition, such as the size and density of built-up areas. By mapping ISA as a continuous field, the sub-pixel percent ISA can be used to quantitatively determine the spatial extent and development density of highly heterogeneous urban areas (Esch et al., 2009; Hu and Weng 2009; Jensen and Cowen 1999; Mushore et al., 2017; Sexton et al., 2013; Small 2003; Wu 2004; Wu and Murray 2003; Wulder et al., 2018; Xian and Homer 2010). However, mapping impervious surface features is complicated by a mix of widely varied spectral responses representing many varied features. Remotely sensed data acquired from medium resolution satellites (e.g., Landsat) have provided rich observations of the terrestrial surface across large geographic extents for more than four decades. These datasets have been widely used to monitor urban land cover and land use change at both regional and national scales (Lu et al., 2011; Powell et al., 2007; Sexton et al., 2013; Song et al., 2016; Xian et al., 2011; Xie and Weng 2016; Zhou and Wang 2008). However, impervious surface mapping using Landsat imagery is still a challenge, especially in areas that have low ISA intensity and where urban lands are mixed with nonurban lands. These ISA and non-ISA mixed areas usually contain other land cover types, such as bare soil and vegetation, that have similar spectral signatures. Mixed pixels usually dominate most developed areas in many medium resolution images including Landsat. In addition to the spatial resolution challenge for mapping urban land (Jensen and Cowen 1999), the spectral resolution of a sensor could also impact the ability to derive detailed information on the nature and properties of different surface materials on the ground (Herold et al., 2004; Hu and Weng 2009; Rashed et al., 2005; Weng 2012). These urban land cover types and other areas such as concrete roads and some natural rocks may have similar reflectance features and thus are difficult to separate. These mixed classes need to be accountable in ISA estimation with some degree of spectral separability. The use of 30 m resolution Landsat imagery has produced several major land cover mapping products across the globe (Grekousis et al., 2015; Li et al., 2018c). However, impervious surface mapped from Landsat data has been considered too coarse for mapping urban biophysical descriptors (Gómez et al., 2016; Slonecker et al., 2010). In recent years, an increasing number of satellites are collecting data at spatial scales relevant to anthropogenic interactions with the landscape. These data from different sensors with different spatial and spectral resolutions usually are used independently when they are employed for land cover mapping. The differences in mapping land cover, especially in continuous field mapping, using these different sensors have not been systematically evaluated. In recent decades, QuickBird, IKONOS, and WorldView (WV), which have relatively higher spatial or both spatial and spectral resolutions, have been used in urban land cover mapping (Hamedianfar and Shafri 2016; Hamedianfar et al., 2014; Lu and Weng 2009; Salehi 2013; Small 2003; Yang and Li 2015) and impervious surface mapping (Iabchoon et al., 2017). However, most studies using high-resolution imagery for mapping urban land cover focused on using the information for pixel-based thematic land cover classification (Hamedianfar et al., 2014; Iabchoon et al., 2017). Furthermore, these high-resolution images have not been extensively used for mapping impervious surfaces in large areas due to lack of data availability, data cost, as well as time and labor required to process the large volume of data. Little research has been conducted to examine combining the multiple resolutions of sensor data in ISA mapping, especially through

2. Study areas We selected two study areas in the United States for this study: San Francisco, California, and Dallas, Texas. The San Francisco area has dense, tall buildings in the business district, parks, dense road systems, sandy beaches, parking lots, and residential housing. The Dallas area has an airport, multiple level road systems, residential areas, and business complexes mixed with wetlands. These areas were selected based on the following factors: 1) the diversity of imperviousness surface features; 2) the presence of mixed pixels that contain imperviousness and other land cover classes; and 3) the availability of data from multiple sensors on the same or nearly the same date. 3. Data and methods 3.1. Remote sensing images Data used in the study include images from the National Agriculture Imagery Program (NAIP), WV3, Landsat 8, and Sentinel-2 (Table 1). The NAIP image has four spectral bands with a spatial resolution of 0.6 m–1 m and was degraded to 2 m. The WV3 image, degraded from its original 1.38 m–2 m resolution to match NAIP 2 m data, has eight bands from coastal to near-infrared (NIR). The Sentinel-2 image has three spatial resolutions and thirteen spectral bands. Ten spectral bands from Sentinel-2 were used, excluding coastal aerosol, water vapor, and cirrus bands. The spectral bands were re-organized as four combinations. The first one includes only four spectral bands (blue, green, red, and NIR) at 10 m resolution. The second combination includes six bands (Red edge

Table 1 The remotely sensed data used in the study. NIR: near-infrared; SWIR: shortwave infrared.

Bands Area San Francisco, CA Dallas, TX

NAIP (training)

WV3

Sentinel-2

Landsat 8

Blue, Green, Red, NIR1 (1m) date 06/25/2016 07/24/2014

Coastal, Blue, Green, Yellow, Red, Red edge, NIR-1, NIR-2 (10m, 20m, 30m) date 03/08/2017 06/18/2015

Blue, Green, Red, Red-Edge 1-4, NIR, SWIR1, SWIR-2 (10m, 20m, 30m) date 03/01/2017 07/24/2016

Coastal, Blue, Green, Red, NIR, SWIR-1, SWIR-2, (30m) date 05/13/2017 06/20/2015

2

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

3.3. Mapping impervious surface

1 to 4 and two shortwave infrared (SWIR)) at 20 m resolution. The third and fourth bands combine the first two datasets and have all ten spectral bands degraded to 20 m and 30 m resolutions, separately. Seven spectral bands at 30 m resolution were used from Landsat 8 Collection 1, excluding the pan and thermal infrared (TIR) bands. All NAIP, WV3, Sentinel-2, and Landsat images, the last two having less than 5% cloud cover, were acquired from USGS Earth Resources Observation and Science (EROS) Center (https://www.usgs.gov//centers/ eros). The NAIP images were chosen to develop both training and validation datasets. The other images were used to map imperviousness distributions.

The regression tree modeling algorithms that are similar to the algorithms used for the NLCD continuous field mappings including impervious surface and shrub/grass (Xian and Homer 2010; Xian et al., 2015) were used to estimate percent imperviousness. Regression tree (RT) models are rule-based predictive models in which a multivariate linear model is associated with a rule or a condition and is created according to the relevant training dataset. A training dataset that contains a dependent variable and several independent variables is required to constrain the parameters in each linear regression equation for projection. The regression tree models connect the responding variable (dependent variable) to the controlling variables (independent variables) at a specific condition based on their correlations. Generally, a selected image from a sensor and its derived indices were used as the controlling variable and the training data as the responding variable. Therefore, the frequently used controlling variables in either conditions or linear equations could be found in RT models. To enhance modeling performance, we calculated four indices including Normalized Difference Vegetation Index (NDVI), Soil-adjusted Vegetation Index (SAVI), Normalized Difference Water Index (NDWI), and Normalized Difference Built-Up Index (NDBI) using Sentinel-2 and Landsat 8 images (Table 3). For the WV3 image, four indices including NDWI, NDVI, Normalized Difference Soil Index (NDSI), and Non-Homogeneous Feature Index (NHFD) were calculated. These index images were also used as additional independent variables to create RT models for ISA mapping for comparison.

3.2. Development of training and validation datasets Methods designed for this study consist of several procedures including training and validation data developed using high-resolution NAIP images, impervious surface modeling using WV3, Sentinel-2, and Landsat 8 images, and product validation and comparison using data developed from NAIP and separated from the training dataset. Fig. 1 is a flowchart illustrating the conceptual procedures used to quantify ISA using three different image datasets. The NAIP images were first classified using an unsupervised classification software provided by ERDAS 20161. The image was first classified as 12 different categories including tree, grass, water, soil, agricultural land, and developed impervious surface. After that, a visual check was performed to re-categorize these classes into developed impervious surface, non-developed, and uncertain classes such as shadows from tall buildings in the center of San Francisco. In the last step, a manual editing process characterized the images as impervious and non-impervious by removing most of the uncertain classes to non-impervious. Some of the uncertain classes are related to shadows from buildings and these classes were labeled as impervious class. By comparing with the NAIP image, the final cleanup ensured that the two classes were all correct at 2 m resolution and all impacts of seasonality were not included. The binary image was further upscaled and converted to 10 m, 20 m, and 30 m resolution to represent percentages of imperviousness at these scales by using the area-weighted spatial model that was designed using GDAL and Python to perform the scaling change function (Shi et al., 2018). The method averages all of the original small pixels that make up the new big pixels, and the quality of resampled images were good because more original pixels were taken into account in the computation. Parts of these fractional ISA data were randomly chosen as training datasets, and other small portions of the data were randomly selected as validation samples for the three spatial resolutions. To ensure the representation of the entire range of imperviousness, both training and validation datasets were sampled across each study area with equal sampling numbers of 0–100 percent. The locations of sampling points for both training and validation were kept the same for mapping ISA using different images at different spatial resolutions. Fig. 2 demonstrates the spatial distributions of validation samples at different spatial resolutions. The numbers of sampling points were reduced from ten of thousands to about several thousand following spatial resolution from high to coarse. All pixels that were identified as cloud, shadow, water, or no-data, as well as pixels located at the edge of non-impervious surface, were manually excluded from training and validation populations. The WV3, Sentinel-2, and Landsat 8 images were also resampled and geometrically re-registered to these training images separately so that their grids would overlap exactly at the 10 m, 20 m, and 30 m pixels of the training and validation datasets. Table 2 illustrates the numbers of training and validation samples used for images from different sensors in different study areas.

3.4. Accuracy assessment The accuracy assessments for ISA estimates were conducted using validation datasets developed from high-resolution NAIP images. Three statistic metrics were calculated, including root mean square error (RMSE), mean error (ME), and relative area error (RAE), by the following formulas: n

∑i = 1 (Xobs, i − Xmod, i )2

RMSE =

(1)

n m

ME =

∑i = 1 (Xobs, i − Xmod, i) n

RAE =

(2)

m n

∑i = 1 Aobs, i − ∑i = 1 Amod, i n ∑i = 1

Aobs, i

(3)

where Xobs,i and Xmod,i are a validation measure and a model estimate of percent imperviousness for a sample point i, Aobs,i and Amod,i are total areas of ISA for a sample point i, n is the number of validation samples, and m is the size of the validation samples. RMSE and RAE are used to measure the overall error estimate and overall error for areal estimate, respectively. ME measures mean error. We also evaluated correlations between validation data and ISA estimates from RT models. 4. Results 4.1. Characteristics of impervious surface distribution The study areas (Fig. 3-I) and spatial distributions of ISA estimated from RT models and relevant images from three sensors are displayed in Fig. 3 (II-III). All modeling results shown here and in the accuracy assessment were obtained without using derived indices as listed in Table 3 except for Landsat-derived ISAs. In the San Francisco area, ISAs derived from WV3 contain high intensity ISAs in areas that have high density buildings in all three resolutions (Fig. 3-II a-d). ISAs derived from Sentinel-2 contain less of these high intensity ISAs but are relatively consistent in all three resolutions (Fig. 3-II e-h). The Landsat-

1 Use of any trade, product, or company names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

3

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Fig. 1. Processing workflow of mapping and assessing ISAs using images from different sensors.

confidence level. Several statistical metrics including RMSE, ME, RAE, medians of validation data, and medians of model estimate of percent ISAs were evaluated and are listed in Table 4. In San Francisco, RMSEs at 20 m and 30 m in WV3-derived ISAs are smaller than RMSEs in Sentinel-2 and Landsat 8. RMSE from Landsat-derived ISA is the largest. MEs in San Francisco are positive in both 10 m and 30 m estimates for WV3derived ISAs but negative for 20 m. MEs for both Sentinel- and Landsatderived ISAs are negative. Similar to MEs, the RAEs are less than 3% for ISA estimated from WV3 and Sentinel-2. The Landsat-derived RAE is about -27.01%, the largest in all study areas. In Dallas, the RMSE values are about 6% for ISAs estimated from WV3 and about 9% – 12% for ISAs estimated from Sentinel-2. The ME values are all positive (0.3 – 0.19%) for ISA derived from WV3 and all negative (about -0.44 – 0.49%) for ISAs derived from Sentinel-2. The RAEs are positive (< 0.4%) for WV3-estimated ISA and negative (-0.56 – -0.84%) for Sentinel-estimated ISA. Both RMSE and RAE values for Landsat-estimated ISA are much less than those in San Francisco. In general, ISAs produced from WV3 images have relatively higher similarity with the true ISA than those from Sentinel-2 images in San Francisco and Dallas. However, ISAs derived from Sentinel-2 images have higher close-to-true median values in the 10 m and 20 m estimates than WV3-derived ISAs. Landsat-derived ISAs have the much lower median values than the validation median ISA in all three areas. In all these comparisons, ISAs were obtained without using indices listed in Table 3 for both WV3 and Sentinel-2 data. Comparing to ISAs with use of indices, both statistical metrics and r2 values were not better than those without using indices. However, Landsat-derived ISA shows slight

derived ISA, however, shows less high intensity ISA in the business centers where tall buildings are seen (Fig. 3-II i-j). In the Dallas area, ISAs estimated from WV3 and Sentinel-2 have similar patterns: high ISAs in the airport and business centers and medium intensity ISAs in residential areas (Fig. 3-III a-d, e-h). However, ISA intensities estimated from Landsat are relatively lower even in areas with tall buildings and in the airport (Fig. 3-III i-j). 4.2. Accuracies of ISAs The scatterplots between validation data and model estimates of ISA and associated correlation of determination (r2) are shown in Fig. 4. In San Francisco, r2 value of WV3-derived ISA is the largest at 10 m (Fig. 4a), the second at 30 m (Fig. 4 c), and smallest at 20 m (Fig. 4b). Similarly, r2 values for Sentinel-derived ISA vary from large to small at 10 m, 20 m, and 30 m (Fig. 4 d-f). Both WV3- and Sentinel-estimated ISA show apparent overestimations the low imperviousness density areas. In Dallas, r2 values for WV3-derived ISA (Fig. 4 g-i) show the largest at 20 m and smallest at 30 m. The values for Sentinel-derived ISA (Fig. 4 j-l) vary from large to small as the spatial resolution changes from 10 m to 30 m. These modeled ISAs do not have apparent underestimate issues in the low-density imperviousness areas. WV3 performs slightly better than Sentinel-2 in ISA estimates with a relatively large r2 in the same spatial resolution. For the Landsat-derived ISA, the r2 value in San Francisco is smaller than the r2 derived from WV3 but larger than the r2 derived from Sentinel-2 (Fig. 4 m). However, the r2 value from Landsat-derived ISA in Dallas is smaller than that from in both WV3and Sentinel-derived ISA (Fig. 4 n). All r2 values are significant at 95% 4

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Fig. 2. Examples of spatial distributions of validation samples across datasets at 10 m, 20 m, and 30 m resolutions.

performances of different sensors in quantifying ISA in different ranges could be assessed. In Fig. 5, all MEs are averaged in every 5% of imperviousness. In San Francisco, overestimates of ISA are observed from 5% to about 60% for both WV3 and Sentinel-2 in all three estimates (Fig. 5 a-c). In the above 60% range, ISAs derived from both Sentinel-2 and WV3 are slightly undervalued. However, ISA derived from Landsat underestimates the intensity of true ISA in almost all ranges (Fig. 5 c). In Dallas, intensities of ISA derived from WV3 and Sentinel-2 are slightly underestimated in ranges between 40 and 50% and overestimated in other ranges in the 10 m estimates (Fig. 5 d). The ISAs derived from WV3 are slightly overvalued in ranges between 60% and

Table 2 Numbers of training and validation datasets. Area Sensor WorldView-3

Sentinel-2 Landsat 8

Resolution 10m 20m 30m 10m 20m 30m 30m

San Francisco, CA

Dallas, TX

Training 10,000 8,000 7,000 10,000 8,000 7,000 7,000

Training 8,000 6,000 5,000 8,000 6,000 5,000 5,000

Validation 300 250 230 300 250 230 130

Validation 900 800 600 800 700 600 600

Table 3 Indices used in regression tree models for impervious surface mapping. Letters represent different bands: C for coastal, NIR for near-infrared, R for red, G for green, Y for yellow, RE for red edge. L is a canopy background adjustment factor and 0.5 is used to minimize soil brightness variations and eliminate the need for additional calibration for different soils in SAVI. SN2 and LS8 represent Sentinel-2 and Landsat 8. Index Ratio

Purpose

Formula

Normalized Difference Water Index (NDWI)

Identify areas of standing water in size greater than one pixel

NDWI =

Normalized Difference Vegetation Index (NDVI)

Identify areas of vegetation and determine the health of each vegetation class.

NDVI =

(R − NIR2) (R + NIR2)

WV3

Normalized Difference Soil Index (NDSI)

Identify areas where soil is the dominant background or foreground material

NDSI =

(G − Y ) (G + Y )

WV3

Non-Homogeneous Feature Difference (NHFD) Soil-adjusted Vegetation Index (SAVI)

Classifying areas which contrast against the background, which can be identified as manmade Classifying vegetation areas

Normalized Difference Water Index (NDWI)

Extracting built-up features indices range from -1 to 1.

NDWI =

Normalized Difference Built-Up Index (NDBI)

Extracting built-up features indices range from -1 to 1.

NDBI =

improvements in statistical metrics and the r2 value in Dallas. The MEs for ISAs estimated from the three sensors in the two areas in every 5% of imperviousness were also evaluated so that the

Sensor (C − NIR2) (C + NIR2)

NHFD = SAVI =

(RE − C ) (RE + C )

(1 + L)(NIR − R) (NIR + R + L)

WV3

WV3 SN2, LS8

(G − NIR) (G + NIR)

SN2, LS8

(SWIR − NIR) (SWIR + NIR)

SN2, LS8

90% but slightly undervalued in the above 95% range. In the 20 m estimate, intensities of ISAs derived from both WV3 and Sentinel-2 are overestimated in ranges below 50% and underestimated in ranges

5

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Fig. 3. The study area (I), images and mapped impervious surfaces in San Francisco (II) and Dallas (III). In I, WV3 images (a) are displayed by bands R, NIR-1, G, and derived ISAs at 10 m, 20 m, and 30 m (b-d). Sentinel-2 images (e), displayed by bands SWIR-2, NIR, G, and derived ISA in 10 m, 20 m, and 30 m (f-h). The last row includes Landsat 8 images (i), displayed by bands SWIR-2, NIR, R, and derived ISA at 30 m (j). In II, images and derived ISA are in the same order: a-d are WV3 images and derived ISA; e-h are Sentinel-2 images and derived ISA; i-j are Landsat images and derived ISA. The non-ISA areas are in transparence with background images of WV3, Sentinel-2, and Landsat images at 10 m, 20 m, and 30 m.

6

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Fig. 4. Correlations between modeled ISA and validation data. X-axis and y-axis in all panels represent validation and modeled ISAs. Panels a-c and d-f are WV3- and Sentinel-derived ISAs in the San Francisco area. Panels g-i and j-l are WV3- and Sentinel-derived ISAs in the Dallas area. Panels m and n are Landsat-derived ISAs in the San Francisco and Dallas areas. 7

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Table 4 Validations in the two study areas. Median values of validation (M_v), median values of model estimate (M_m), RMSE, ME, and RAE are in percentages. WV3

San Francisco

Dallas

RMSE ME RAE M_v M_m RMSE ME RAE M_v M_m

10 m 15.70 1.87 2.81 71 78 6.29 0.30 0.42 88 91

Sentinel-2 20 m 14.48 -0.18 -0.28 69 70 5.94 0.18 0.31 62 63

30 m 13.44 0.28 0.46 64 66 6.42 0.19 0.34 57 57

10 m 15.66 -0.99 -1.55 67 69 9.05 -0.44 -0.56 96 94

Landsat 8 20 m 15.18 -1.42 -2.44 67 68 10.59 -0.49 -0.78 67 65

30 m 15.67 -1.46 -2.31 64 65 12.29 -0.49 -0.84 60 60

30 m 20.10 -13.34 -27.01 45 35 13.34 -3.52 -6.17 59 58

Fig. 5. Mean error (ME) of percent impervious surface estimated from images of three sensors at 10 m, 20 m, and 30 m in San Francisco (a-c) and Dallas (d-f). MEs are averaged in every 5% imperviousness.

ranges above 20%. In general, ISA estimated from 20 m WV3 imagery has the smallest ME values. Also, ME values from WV3 images are smaller than those from the other two sensors in areas that imperviousness is above 70%. By evaluating RT model outputs, we also found the most frequently used bands in condition tests and projection models. Table 5 shows the first three most frequently used bands in different models and in different areas. In the RT models created from WV3, the most frequently used bands in condition tests and in models are the coastal and red bands, respectively. In the RT models created from Sentinel-2, the most frequently used bands in condition tests and in models are the NIR band in the 10 m modeling and the blue and green bands in both 20 m and 30 m modeling. In the RT models created from Landsat, the most frequently used band in condition tests becomes NDVI. The most

above 60% for Sentinel-2 (Fig. 5 e). The WV3-derived ISA is close to the true distribution in ranges between 60% and 90% and Sentinel-2 derived ISA is slightly underestimated in the range. In the 30 m estimates, the WV3-derived ISA is close to the true distribution in all ranges. ISAs derived from other two sensors are overvalued in the below 20% range and undervalued in ranges between 40% and 60% for Landsat-derived ISA (Fig. 5 f). The ISAs derived from both Sentinel-2 and Landsat are underestimated for the ranges above 80%. The comparisons of ME values suggest that intensities of ISAs derived from both WV3 and Sentinel-2 are overestimated in the less than 30% range and underestimated in the over 70% range. Such overestimation in the low intensity range and underestimation in high intensity range could balance the overall accuracy in the areal estimate of ISA. However, the intensities of ISA derived from Landsat are undervalued in almost all 8

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Table 5 The most frequently used bands in regression models and band definitions: C for coastal, B for blue, G for green, R for red, Y for yellow, RE for red edge, NIR-1 for near-infrared-1, NIR-2 for near-infrared-2, SWIR-1 for shortwave infrared-1, SWIR-2 for shortwave infrared-2, RE-1 for red-edge-1, RE-2 for red-edge-2, RE-3 for rededge-3, RE-4 for red-edge-4. Indices are the same as in Table 3. Rank-1 is the first most used band in the regression model; Rank-2 is the second most used band in regression tree; and Rank-3 is the third most used band in regression tree model. Study Area

San Francisco, CA

Sensor

Resolution

Rank-1

Rank-2

Rank-3

Rank-1

Rank-2

Rank-3

WV3

10m

C NIR-1 C NIR-1 C Y B NIR B G, R B G SAVI B

NIR-2 Y NIR-2 R NIR-2 C NIR B NIR NIR NIR B SAVI C

NIR-1 B NIR-1 NIR-2 NIR-2 B, NIR-1 NIR R SWIR-2 B SWIR-2 RE-4, R NDBI SWIR-1

C R C R C B B NIR R SWIR-2 B SWIR-1 SAVI SWIR-2

NIR-1 Y, RE, B NIR-2 RE NIR-2 C, RE R NIR SWIR-2 R, SWIR-1 R R, SWIR-2 R R

NIR-2 C R C, G R G NIR R B NIR SWIR-2 NIR SWIR-2 SWIR-1

20m 30m 10m Sentinel-2

20m 30m 30m

Landsat 8

Conds Model Conds Model Conds Model Conds Model Conds Model Conds Model Conds Model

Dallas, TX

measured using RAE. By comparing ISA estimates from images of three sensors, we found that WV3 underestimated imperviousness areas in 20 m had the smallest RAEs (Table 4). ISA derived from Sentinel-2 data had the lowest underestimates at 10 m. However, the RAEs of ISA estimated from Landsat reached -27.013% and -6.17%. The shadows of tall buildings in the central San Francisco area could be the major cause for such a large RAE. Again, the ISAs derived from WV3 data had relatively lower RAE than the other two sensors at the same spatial resolution. We did not use SWIR bands in WV3 images due to the lack of data from these bands in our study areas. The comparisons using ISAs derived from WV3 images may not completely represent the capability of WV3 in estimating imperviousness, especially in areas where ISA is mixed with denser vegetation canopy. For mapping impervious surfaces in urban areas, especially areas with land cover type heterogeneity, Landsat's 30 m spatial resolution inherently has larger groups of spectrally mixed pixels. In modeling practices, the bias toward underestimation can be rectified and accounted for by adjusting and adding training data that provide more distinct spectral signatures for regression modeling to reliably separate mixing features. However, this adjustment could also produce higher commission error outside of target areas. The coarser resolution of Landsat compared to Sentinel-2 and WV3 tended to underestimate the percent imperviousness due to the more mixed and muddled spectral signatures from a large amount of pixels. A higher resolution, such as 20 m, for future Landsat missions can improve accuracy for impervious surface mapping. All solar reflective spectral regions, including visible, NIR, and SWIR bands, are used in impervious surface mapping; therefore, maintaining the full reflective spectral capability is needed for future Landsat. Furthermore, additional spectral bands such as Red edge bands can be helpful in more accurate impervious area mapping and reduction in estimation bias. Additionally, a yellow band can be helpful in certain environments such as dryland ecosystems. Overall, maintaining the continuity of current Landsat spectral bands is critical, and additional improvements in spatial and spectral resolution are needed for impervious surface mapping.

frequently used bands in models are coastal, blue, and SWIR-2 in San Francisco and Dallas, respectively. The second most frequently used spectral bands varied by location, spatial resolution, and sensor. 5. Discussion The goal of this study was to evaluate the performances of imperviousness estimates using WV3, Sentinel-2, and Landsat 8 top of atmosphere reflectance data through regression tree models. While this type of assessment has previously been performed, especially by comparing high resolution images with Landsat, to facilitate the use of satellite images to map ISA distributions, this assessment focused on the data from three sensors to provide recommendations for future Landsat development. The modeling method used here is the same as used for mapping the USGS National Land Cover impervious surface product. One of advantages of continuous field mapping is that it can reveal the density distribution of the target variable in a range from 1 to 100%. The performances of different satellite data in estimating fractional imperviousness can therefore be evaluated in different density ranges as well as the areal estimate. The differences in spatial resolution between sensors were significant and resulted in different performances in mapping ISA. For ISAs derived from WV3 data, 20 m estimates performed better than other ISAs at 10 m and 30 m with lower MEs and RAEs. The break points of ISA overestimates were moved from 80% at 10 m to 65% at 20 m and 30 m in the San Francisco area where some non-impervious surfaces such as sands and bare soil have large albedos for high reflectance (Fig. 4b). Similarly, the ISA at 20 m in Dallas also had a better performance by reducing underestimates in low and overestimates in high imperviousness areas (Fig. 4e). The WV3 image at 20 m is slightly better in RMSE, ME, and RAE but not for r2 for characterizing imperviousness using RT models. However, the Sentinel-derived ISA at 10 m outperformed ISAs in other spatial resolutions. Adding additional spectral bands to the 20 m and 30 m estimates did not improve the accuracy of modeled ISAs. This also explained why ISAs estimated from RT models with additional index inputs did not achieve better results. The Landsatderived ISAs had the lowest accuracies in both study areas, although additional index inputs slightly improved the accuracies of ISAs. The ISAs estimated from Landsat data were underestimated in almost all imperviousness ranges except in very low-density ranges. Developed ISAs estimated in San Francisco and Dallas from WV3, which has higher spatial and spectral resolutions, outperformed ISAs mapped from both Sentinel-2 and Landsat in all spatial resolutions. The performance of total imperviousness area estimate was

6. Conclusions Our RT model results found that WV3 and Sentinel-2 have similar capabilities to detect and quantify developed impervious surface at 10 m, 20 m, and 30 m resolutions, although the overall performance of WV3 is slightly better at 20 m. The spatial distributions of ISA produced from images of these two sensors represent urban land cover intensities 9

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

well in the two study areas. However, Sentinel-2, which has multiple Red edge bands, did not show substantial improvement in mapping ISA but reduced biases associated with non-ISA in the 20 m and 30 m estimates. The ISA derived from images of Landsat 8 can also represent general spatial features of urban land cover. As Landsat 8 has fewer spectral bands and coarser spatial resolution than the other two sensors, ISAs derived from Landsat data are underestimated, especially in high intensity urban centers. Such underperformance results in about 6%–27% overall undervalues in imperviousness area estimates. These ISA underestimates are substantial in city centers where tall buildings cast shadows. Our results also suggest that RAEs of ISA estimated at 20 m using eight bands from WV3 are less than 3%. The ISA estimates using 10 m WV3 or Sentinel-2 images were not better than those using 20 m images in both RMSE and areal estimates.

Gong, P., Liu, H., Zhang, M., Li, C., Wang, J., Huang, H., Clinton, N., Ji, L., Li, W., Bai, Y., Chen, B., Xu, B., Zhu, Z., Yuan, C., Ping Suen, H., Guo, J., Xu, N., Li, W., Zhao, Y., Yang, J., Yu, C., Wang, X., Fu, H., Yu, L., Dronova, I., Hui, F., Cheng, X., Shi, X., Xiao, F., Liu, Q., Song, L., 2019. Stable classification with limited sample: transferring a 30m resolution sample set collected in 2015 to mapping 10-m resolution global land cover in 2017. Sci. Bull. 64, 370–373. Grekousis, G., Mountrakis, G., Kavouras, M., 2015. An overview of 21 global and 43 regional land-cover mapping products. Int. J. Remote Sens. 36, 5309–5335. Grimm, N.B., Faeth, S.H., Golubiewski, N.E., Redman, C.L., Wu, J., Bai, X., Briggs, J.M., 2008. Global change and the ecology of cities. Science 319, 756–760. Hamedianfar, A., Shafri, H.Z.M., 2016. Integrated approach using data mining-based decision tree and object-based image analysis for high-resolution urban mapping of WorldView-2 satellite sensor data. J. Appl. Remote Sens. 10, 025001. Hamedianfar, A., Shafri, H.Z.M., Mansor, S., Ahmad, N., 2014. Improving detailed rulebased feature extraction of urban areas from WorldView-2 image and lidar data. Int. J. Remote Sens. 35, 1876–1899. Herold, M., Roberts, D.A., Gardner, M.E., Dennison, P.E., 2004. Spectrometry for urban area remote sensing—development and analysis of a spectral library from 350 to 2400 nm. Rem. Sens. Environ. 91, 304–319. Hu, X., Weng, Q., 2009. Estimating impervious surfaces from medium spatial resolution imagery using the self-organizing map and multi-layer perceptron neural networks. Rem. Sens. Environ. 113, 2089–2102. Iabchoon, S., Wongsai, S., Chankon, K., 2017. Mapping urban impervious surface using object-based image analysis with WorldView-3 satellite imagery. J. Appl. Remote Sens. 11. Imhoff, M.L., Zhang, P., Wolfe, R.E., Bounoua, L., 2010. Remote sensing of the urban heat island effect across biomes in the continental USA. Rem. Sens. Environ. 114, 504–513. IPCC, 2013. The Physical Science Basis. The Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK and New York, NY, USA, pp. 36. Jensen, J.R., Cowen, D.C., 1999. Remote sensing of urban/suburban infrastructure and socio-economic attributes. Photogramm. Eng. Rem. Sens. 65, 611–622. Li, H., Zhou, Y., Li, X., Meng, L., Wang, X., Wu, S., Sodoudi, S., 2018a. A new method to quantify surface urban heat island intensity. Sci. Total Environ. 624, 262–272. Li, X., Zhou, Y., 2017. Urban mapping using DMSP/OLS stable night-time light: a review. Int. J. Remote Sens. 38, 1–17. Li, X., Zhou, Y., Asrar, G.R., Zhu, Z., 2018b. Creating a seamless 1 km resolution daily land surface temperature dataset for urban and surrounding areas in the conterminous United States. Rem. Sens. Environ. 206, 84–97. Li, X., Zhou, Y., Zhu, Z., Liang, L., Yu, B., Cao, W., 2018c. Mapping annual urban dynamics (1985–2015) using time series of Landsat data. Rem. Sens. Environ. 216, 674–683. Lu, D., Moran, E., Hetrick, S., 2011. Detection of impervious surface change with multitemporal Landsat images in an urban-rural frontier. ISPRS J. Photogrammetry Remote Sens. 66, 298–306. Lu, D., Weng, Q., 2009. Extraction of urban impervious surfaces from an IKONOS image. Int. J. Remote Sens. 30, 1297–1311. Mouri, G., Shinoda, S., Oki, T., 2012. Assessing environmental improvement options from a water quality perspective for an urban–rural catchment. Environ. Model. Softw 32, 16–26. Mushore, T.D., Odindi, J., Dube, T., Matongera, T.N., Mutanga, O., 2017. Remote sensing applications in monitoring urban growth impacts on in-and-out door thermal conditions: a review. Remote Sensing Applications: Society and Environment 8, 83–93. Piano, E., De Wolf, K., Bona, F., Bonte, D., Bowler, D.E., Isaia, M., Lens, L., Merckx, T., Mertens, D., van Kerckvoorde, M., De Meester, L., Hendrickx, F., 2017. Urbanization drives community shifts towards thermophilic and dispersive species at local and landscape scales. Glob. Chang. Biol. 2554–2564. Powell, R., Roberts, D., Dennison, P., Hess, L., 2007. Sub-pixel mapping of urban land cover using multiple endmember spectral mixture analysis: Manaus, Brazil. Rem. Sens. Environ. 106, 253–267. Rashed, T., Weeks, J., Stow, D., Fugate, D., 2005. Measuring temporal compositions of urban morphology through spectral mixture analysis: toward a soft approach to change analysis in crowded cities. Int. J. Remote Sens. 26, 699–718. Salehi, B., 2013. A combined object- and pixel-based image analysis framework for land cover classification of VHR imagery. Photogramm. Eng. Rem. Sens. 79, 999–1014. Sexton, J.O., Song, X.-P., Huang, C., Channan, S., Baker, M.E., Townshend, J.R., 2013. Urban growth of the Washington, D.C.–Baltimore, MD metropolitan region from 1984 to 2010 by annual, Landsat-based estimates of impervious cover. Rem. Sens. Environ. 129, 42–53. Shi, H., Auch, R.F., Vogelmann, J.E., Feng, M., Rigge, M., Senay, G., Verdin, J.P., 2018. Case study comparing multiple irrigated land datasets in Arizona and Colorado, USA. J. Am. Water Resour. Assoc. 54, 505–526. Shi, H., Singh, A., Kant, S., Zhu, Z., Waller, E., 2005. Integrating habitat status, human population pressure, and protection status into biodiversity conservation priority setting. Conserv. Biol. 19, 1273–1285. Shi, L., Ling, F., Ge, Y., Foody, G.M., Li, X., Wang, L., Zhange, Y., Du, Y., 2017. Impervious surface change mapping with an uncertainty-based spatial-temporal consistency model: a case study in Wuhan city using Landsat time-series datasets from 1987 to 2016. Rem. Sens. 1148. Slonecker, T., Milheim, L., Claggett, P., 2010. Landscape indicators and land cover change in the mid-atlantic region of the United States, 1973-2001. GIScience Remote Sens. 47, 163–186. Small, C., 2003. High spatial resolution spectral mixture analysis of urban reflectance. Rem. Sens. Environ. 88, 170–186. Song, X.-P., Sexton, J.O., Huang, C., Channan, S., Townshend, J.R., 2016. Characterizing

Ethical statement I testify on behalf of all co-authors that our article submitted to Remote Sensing Applications: Society and Environment. 1)this material has not been published in whole or in part elsewhere; 2)the manuscript is not currently being considered for publication in another journal; 3)all authors have been personally and actively involved in substantive work leading to the manuscript, and will hold themselves jointly and individually responsible for its content. Declarations of interest None. Acknowledgements The authors would like to thank the USGS National Land Imaging Program for their funding support of this research. We also thank Kevin Gallo and Roger Auch for reviewing the manuscript. Hua Shi's work was performed under USGS contract 140G0119C0001. Appendix ASupplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.rsase.2019.100246. References Benedetti, A., Picchiani, M., Del Frate, F., 2018. Sentinel-1 and sentinel-2 data fusion for urban change detection. Institute of Electrical and Electronics Engineers Inc, pp. 1962–1965. Berezowski, T., Chormański, J., Batelaan, O., Canters, F., Van de Voorde, T., 2012. Impact of remotely sensed land-cover proportions on urban runoff prediction. Int. J. Appl. Earth Obs. Geoinf. 16, 54–65. Cadenasso, M.L., Pickett, S.T., Schwarz, K., 2007. Spatial heterogeneity in urban ecosystems: reconceptualizing land cover and a framework for classification. Front. Ecol. Environ. 5, 80–88. Dousset, B., Gourmelon, F., 2003. Satellite multi-sensor data analysis of urban surface temperatures and landcover. ISPRS J. Photogrammetry Remote Sens. 58, 43–54. Esch, T., Himmler, V., Schorcht, G., Thiel, M., Wehrmann, T., Bachofer, F., Conrad, C., Schmidt, M., Dech, S., 2009. Large-area assessment of impervious surface based on integrated analysis of single-date Landsat-7 images and geospatial vector data. Rem. Sens. Environ. 113, 1678–1690. Fu, P., Weng, Q., 2018. Variability in annual temperature cycle in the urban areas of the United States as revealed by MODIS imagery. ISPRS J. Photogrammetry Remote Sens. 146, 65–73. Gallo, K., Xian, G., 2016. Changes in satellite-derived impervious surface area at US historical climatology network stations. ISPRS J. Photogrammetry Remote Sens. 120, 77–83. Giovannini, L., Zardi, D., de Franceschi, M., 2011. Analysis of the urban thermal fingerprint of the city of trento in the alps. Journal of Applied Meteorology and Climatology 50, 1145–1162. Gómez, C., White, J.C., Wulder, M.A., 2016. Optical remotely sensed time series data for land cover classification: a review. ISPRS J. Photogrammetry Remote Sens. 116, 55–72.

10

Remote Sensing Applications: Society and Environment 15 (2019) 100246

G. Xian, et al.

Xian, G., Crane, M., Su, J., 2007. An analysis of urban development and its environmental impact on the Tampa Bay watershed. J. Environ. Manag. 85, 965–976. Xian, G., Homer, C., 2010. Updating the 2001 national land cover Database impervious surface products to 2006 using Landsat imagery change detection methods. Rem. Sens. Environ. 114, 1676–1686. Xian, G., Homer, C., Dewitz, J., Fry, J., Hossain, N., Wickham, J., 2011. Change of impervious surface area between 2001 and 2006 in the conterminous United States. Photogramm. Eng. Rem. Sens. 77, 756–762. Xian, G., Homer, C., Rigge, M., Shi, H., Meyer, D., 2015. Characterization of shrubland ecosystem components as continuous fields in the northwest United States. Rem. Sens. Environ. 168, 286–300. Xie, Y., Weng, Q., 2016. Updating urban extents with nighttime light imagery by using an object-based thresholding method. Rem. Sens. Environ. 187, 1–13. Yang, J., Li, P., 2015. Impervious surface extraction in urban areas from high spatial resolution imagery using linear spectral unmixing. Remote Sensing Applications: Society and Environment 1, 61–71. Yuan, F., Bauer, M.E., 2007. Comparison of impervious surface area and normalized difference vegetation index as indicators of surface urban heat island effects in Landsat imagery. Rem. Sens. Environ. 106, 375–386. Zhou, Y., Wang, Y.Q., 2008. Extraction of impervious surface areas from high spatial resolution imagery by multiple agent segmentation and classification. Photogramm. Eng. Rem. Sens. 74, 857–868.

the magnitude, timing and duration of urban growth from time series of Landsatbased estimates of impervious cover. Rem. Sens. Environ. 175, 1–13. Sun, Z., Xu, R., Du, W., Wang, L., Lu, D., 2019. High-resolution urban land mapping in China from sentinel 1A/2 imagery based on Google Earth Engine. Rem. Sens. 11. Tavares, P.A., Beltrão, N.E.S., Guimarães, U.S., Teodoro, A.C., 2019. Integration of sentinel-1 and sentinel-2 for classification and LULC mapping in the urban area of Belém, eastern Brazilian Amazon. Sensors 19. Wang, M., Shi, W., Tang, J., 2011. Water property monitoring and assessment for China's inland Lake Taihu from MODIS-Aqua measurements. Rem. Sens. Environ. 115, 841–854. Weng, Q., 2012. Remote sensing of impervious surfaces in the urban areas: requirements, methods, and trends. Rem. Sens. Environ. 117, 34–49. Weng, Q., Lu, D., Schubring, J., 2004. Estimation of land surface temperature–vegetation abundance relationship for urban heat island studies. Rem. Sens. Environ. 89, 467–483. Wu, C., 2004. Normalized spectral mixture analysis for monitoring urban composition using ETM+ imagery. Rem. Sens. Environ. 93, 480–492. Wu, C., Murray, A.T., 2003. Estimating impervious surface distribution by spectral mixture analysis. Rem. Sens. Environ. 84, 493–505. Wulder, M.A., Coops, N.C., Roy, D.P., White, J.C., Hermosilla, T., 2018. Land cover 2.0. Int. J. Remote Sens. 39, 4254–4284. Xian, G., 2008. Satellite remotely-sensed land surface parameters and their climatic effects for three metropolitan regions. Adv. Space Res. 41, 1861–1869.

11