Nothing Special   »   [go: up one dir, main page]

Next Article in Journal
Variability and climate change trend in vegetation phenology of recent decades in the Greater Khingan Mountain area, Northeastern China
Previous Article in Journal
Developing Theoretical Marine Habitat Suitability Models from Remotely-Sensed Data and Traditional Ecological Knowledge
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Building a Better Urban Picture: Combining Day and Night Remote Sensing Imagery

1
Center for Geospatial Information, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, 1068 Xueyuan Avenue, Shenzhen University Town, Shenzhen, Guangdong 518055, China
2
School of Forestry and Environmental Studies, Yale University, 195 Prospect Street, New Haven, CT 06511, USA
3
Department of Geography, Central Michigan University, Mount Pleasant, MI 48859, USA
4
Google Earth Engine, Google Inc., 1600 Amphitheatre Pkwy, Mountain View, CA 94043, USA
*
Author to whom correspondence should be addressed.
Remote Sens. 2015, 7(9), 11887-11913; https://doi.org/10.3390/rs70911887
Submission received: 30 June 2015 / Revised: 31 August 2015 / Accepted: 1 September 2015 / Published: 16 September 2015
Graphical abstract
">
Figure 1
<p>Clouds and cloud shadows and their effects on NDVI. Left: over vegetated land surface, middle: over water surface, right: over bare land surface. These sites are randomly selected. (<b>A</b>) True color Landsat-7/ETM+ image. (<b>B)</b> Landsat-7/ETM+ NDVI image.</p> ">
Figure 2
<p>Five select testing sites (black rectangles) and two evaluation sites (pink rectangles) in the conterminous USA (overlaid on the NLCD 2006 land cover map).</p> ">
Figure 3
<p>Distributions of maximum, minimum NDVI values within each of the water, vegetation, and bare land categories in the five testing sites. Site #1 and #2 are dominated by vegetation (crops or forests), Site #3 is dominated by water, and Site #4 and #5 are dominated by bare land.</p> ">
Figure 4
<p>The procedure to generate NDVI composites from multi-year Landsat-7/ETM+ imagery with the proposed Mixed NDVI method on Google Earth Engine.</p> ">
Figure 5
<p>Nominal 2006 NDVI composites created with the Maximum NDVI criteria and the Mixed NDVI criteria over a water site (<b>upper row</b>) and a bare land site (<b>bottom row</b>).</p> ">
Figure 6
<p>Landsat-7/ ETM+ NDVI composite of the conterminous US in 2006.</p> ">
Figure 7
<p>Distributions of dates when good quality data were collected that went into the 2000 and 2006 Chicago NDVI composites. (<b>A</b>) Date distribution without year differentiation; (<b>B</b>) cumulative curves of A; (<b>C</b>) date distribution with year differentiation in the 2000 composite; and (<b>D</b>) date distribution with year differentiation in the 2006 composite.</p> ">
Figure 8
<p>Acquisition years and NDVI values in the 2006 composite in a subregion in the Chicago area. (<b>A</b>) Acquisition years: black: 2006, grey: 2007, and white: 2008; (<b>B</b>) NDVI values increase from black to white; (<b>C</b>) Distribution of NDVI; (<b>D</b>–<b>F</b>) Google Earth snapshots.</p> ">
Figure 9
<p>The 2006 NDVI (<b>A</b>), NTL (<b>B</b>), NDUI (<b>C</b>), and NLCD (<b>D</b>) maps in the Salt Lake City area, Utah.</p> ">
Figure 10
<p>A close-up look at the urban areas of the Salt Lake City region. (<b>A</b>) Downtown of Salt Lake City; (<b>B</b>) commercial strip; (<b>C</b>) industrial zone; (<b>D</b>) Salt Lake City International Airport; (<b>E</b>) industrial Zone; (<b>F</b>) sand pit at the Point of the Mountain; (<b>G</b>) Bingham Canyon Open-Pit Copper Mine.</p> ">
Figure 11
<p>The distribution of NDVI values in the 2006 NDVI composite (<b>A</b>) and the distribution of the 2006 NDUI within each NLCD 2006 Developed class and the Cultivated Crops class, excluding water, in the Chicago region (<b>B</b>).</p> ">
Figure 12
<p>The latitudinal profiles of NDUI and NLCD ISA in Salt Lake City (<b>A</b>), and the linear regression model between them (<b>B</b>). Panel (<b>C</b>) shows a Google Earth snapshot, which reveals that the piece of land within the dotted rectangle was bare in 2006. The upper graph of panel A is a false color composite (Red: NDUI, Green: NDUI, Blue: NLCD ISA).</p> ">
Figure 13
<p>Distributions of NDUI (<b>A</b>), BCI (<b>B</b>), NDVI (<b>C</b>), and NDBI (<b>D</b>) within both the NLCD 2006 barren land class and the Developed classes in the two evaluation sites.</p> ">
Figure 14
<p>Some not-well-lit Developed areas cannot be well captured in the 2006 NDUI in the Salt Lake City region.</p> ">
Versions Notes

Abstract

:
Urban areas play a very important role in global climate change. There is increasing need to understand global urban areas with sufficient spatial details for global climate change mitigation. Remote sensing imagery, such as medium resolution Landsat daytime multispectral imagery and coarse resolution Defense Meteorological Satellite Program/Operational Linescan System (DMSP/OLS) nighttime light imagery, has provided a powerful tool for characterizing and mapping cities, with advantages and disadvantages. Here we propose a framework to merge cloud and cloud shadow-free Landsat Normalized Difference Vegetation Index (NDVI) composite and DMSP/OLS Night Time Light (NTL) to characterize global urban areas at a 30 m resolution, through a Normalized Difference Urban Index (NDUI) to make full use of them while minimizing their limitations. We modify the maximum NDVI value multi-date image compositing method to generate the cloud and cloud shadow-free Landsat NDVI composite, which is critical for generating a global NDUI. Evaluation results show the NDUI can effectively increase the separability between urban areas and bare lands as well as farmland, capturing large scale urban extents and, at the same time, providing sufficient spatial details inside urban areas. With advanced cloud computing facilities and the open Landsat data archives available, NDUI has the potential for global studies at the 30 m scale.

Graphical Abstract">

Graphical Abstract

1. Introduction

Due to broad impacts of the concentrated built environment and human activities, cities are now considered agents of global change [1]. Although occupying only about 2% of the global land surface [2], cities worldwide are now hosting more than 50% of the world’s population [3], producing more than 90% of the world gross domestic production (GDP) [4], consuming more than 70% of the available energy [5], and generating more than 71% of the anthropogenic greenhouse gas emissions [6]. As a major factor shaping the Earth system, the importance of cities is increasing rapidly [7]. By 2100, 70%–90% of the world’s population, which is projected to add another three billion, will live in urban regions [8]. In order to harness the growth and development benefits of urbanization while proactively managing its negative effects [9], there is growing interest in monitoring urban areas at regional and global scales to help thinking globally, as well as obtaining fine details of urban spatial structures at the local scale to facilitate acting locally. For example, practical carbon emission mitigation strategies require knowledge of cities down to the street and individual building level [10].
Remote sensing imagery with various spatial resolutions has long been used to monitor and map cities locally and globally. Coarse resolution remote sensing imagery can help monitor urban areas at large scales but cannot provide sufficient spatial details at the street level. Furthermore, existing satellite-based efforts at mapping global urban extents with coarse resolution imagery fail to agree on the size and pattern of urban land use, with estimates ranging from 0.2% to 2.4% of terrestrial land surface circa 2000 [11]. The disagreement comes from three major components: no common definition of urban areas, coarse image resolutions, and the performance of algorithms.
Finer resolution images are very useful in mapping urban details but their application to large scales is often limited. For example, with its medium resolution (30 m) and long historical archive, Landsat makes it possible to map urban extents with considerable accuracy. In fact, it has been the main resource in mapping urban areas, especially since its availability to the public as of 2008 at no charge [12]. In a meta-analysis of global urban land expansion, Seto et al. [13] reported a total of 326 studies that have used remotely sensed images, mainly from Landsat, to map urban land conversion scattered around the world up to the year 2000. Although Landsat has contributed significantly to the understanding of urban land cover and land use change, these efforts have largely focused on individual city or city-region studies, with few comparative or global-scale studies at the Landsat resolution [13]. This is mainly due to the Landsat data availability at the time, the limitation of computing facility, and algorithms to analyze the data. Furthermore, all these studies are very labor-intensive due to the optical spectral complexity of urban materials and their confusion with bare lands, including fallow land. In an urban area, many land cover types/surface materials are spectrally similar, which makes it extremely difficult to analyze an urban scene using a single sensor with a limited spectral range [14,15]. For example, some land cover types are spectrally indistinguishable from each other within a Landsat Thematic Mapper (TM) scene, such as asphalt versus water or shadow, and bare soil versus newly completed concrete sections [16].
Various efforts have reported their success in enhancing urban features in remote sensing images, including the Biophysical Composition Index (BCI) [17] and the multi-source Normalized Difference Impervious Surface Index (MNDISI) [18], to make them more distinguishable from background features. However, due to the confusion between impervious surface and water, water has to be masked out first through an unsupervised classification process before calculating BCI, and constructing BCI involves cross-pixel comparison, which thus leads to heavy computation and non-automation [17]. MNDISI combines Landsat surface temperature and fine resolution International Space Station night light images to highlight urban areas [18]. However, temperature can vary significantly across space; for example, in summertime, areas in lowlands can be warmer than areas in highlands. Furthermore, finer resolution night light images are not available for the entire globe and for multiple times.
Unlike Landsat data, DMSP/OLS Nighttime Light (NTL) imagery can provide thematic information to effectively highlight urban areas out of the background comprising other land cover and land use types, thus minimizing the confusion between urban areas and bare lands, and it can reduce the burden of distinguishing land use types within urban areas. Consequently, DMSP/OLS NTL imagery has been a very important resource to characterize and map urban areas [19,20,21,22]. However, it suffers from saturation and overglow effects, as well as its coarse resolution [23]. Mapping urban areas with only DMSP/OLS NTL thus remains a big challenge.
Combining DMSP/OLS NTL with daytime optical spectral information such as Moderate Resolution Imaging Spectroradiometer (MODIS) can further highlight and increase variation in urban areas [24,25,26] through overcoming limitations in each individual source. An observed fact is that the profiles of vegetation density and luminosity intensity along transects passing through urban cores often exhibit an inverse relationship: vegetation density decreases while luminosity increases from the rural area to the urban core. Such a relationship has been used to generate an index to characterize urban areas and bring back the spatial details into saturated urban cores in the NTL images through a combination of MODIS NDVI and DMSP/OLS luminosity data [25]. However, the coarse resolutions remain the key limitation.
The purpose of this study is to highlight urban areas by integrating the medium resolution Landsat and the coarse resolution DMSP/OLS NTL. We aim to investigate the potential of combining the coarse resolution DMSP nighttime light data with medium resolution Landsat data to enhance urban features at regional and global scales. Considering that both Landsat and DMSP/OLS have long historical archives, we expect our investigation on combining them together will also enable tracking urban expansion over time. Specifically, we focus on two major issues. First, we generate cloud and cloud shadow-free Landsat NDVI composites in growing seasons to maximize the separability between fallow lands and urban areas. Second, we combine Landsat NDVI with DMSP/OLS NTL to highlight urban features at the Landsat scale through constructing and implementing a normalized difference urban index (NDUI) on a cloud computing platform, the Google Earth Engine, for regional and global scale urban area analysis. The Google Earth Engine gathers together the Earth’s raw satellite imagery data—petabytes of historical, present, and future data, including those from Landsat—and makes them easily available for analysis via expert-provided algorithms. With its progressive remote sensing image processing capabilities, the advanced computing system Google Earth Engine now allows efficiently processing and characterizing global-scale Landsat time-series data sets for quantifying land change [27].
We expect combining coarse resolution DMSP/OLS nighttime light imagery with the medium resolution Landsat NDVI composites will maximize the separation between urban areas and barren lands, which are often very easily confused. Furthermore, we aim to make the entire image processing and information extraction highly automated to allow for regional and global scale applications. Evaluation results show that the proposed NDUI can greatly enhance urban areas, which can then be easily distinguished from bare lands including fallows and desserts. With the capability to delineate urban boundaries and, at the same time, to present sufficient spatial details within urban areas, the NDUI has the potential for urbanization studies at regional and global scales.

2. Methodology

2.1. Data and Preprocessing

Google Earth Engine hosts the entire Landsat-7 Enhanced Thematic Mapper Plus (ETM+) data archive from the U.S. Geological Survey Earth Resources Observation and Science (EROS) Center Landsat archive. All the Landsat-7/ETM+ L1T data, regardless of the amount of cloud cover, are used in this study. We only use the L1T data to reduce the impact of misregistration errors, which is critical for temporal compositing. The L1T geo-location error is less than 30 m even in areas with substantial terrain relief [28]. Each Landsat-7/ETM+ L1T scene is approximately 185 km by 185 km and since May 2003 has 22% missing pixels occurring in a repeating along-scan stripe pattern because of the ETM+ scan line corrector (SLC) failure [29]. In this study, we use all the six 30 m reflective Landsat ETM+ wavelength bands: Blue (band 1: 0.45–0.52 µm), Green (band 2: 0.53–0.61 µm), Red (band 3: 0.63–0.69 µm), Near-infrared (band 4: 0.78–0.90 µm), and two Middle-infrared (bands 5 and 7: 1.55–1.75 µm and 2.09–2.35 µm). We use all the available Landsat-7/ETM+ imagery from 2000 to 2002 and from 2006 to 2008 to generate NDVI composites before and after the SLC failure to examine its impacts on the compositing processes.
We convert raw Digital Number (DN) to Top of Atmosphere (TOA) reflectance for Landsat-7/ETM+ using coefficients derived from [30]. Bad pixels with at least one spectral band missing often exist in ETM+ images, especially at the edges of a scene. They must be first identified and masked out before generating the NDVI composite. We generate a bad ETM+ pixel (including gaps due to SLC-off) mask by examining the product of Bands 1–5 and Band 7: if the product equals 0, the pixel is a bad one and must be discarded.
We obtained the version 4 DMSP/OLS NTL data from the website of the National Geophysical Data Center at the National Oceanic and Atmospheric Administration [31] and the data were also hosted on Google Earth Engine. The version 4 DMSP/OLS NTL data consists of the 3000 km swath width NTL scan lines divided into an array of grids with 0.55 km spatial resolution, which are aggregated and composited to 30 arc second (about 1 km) grids. Each grid in the annual composite data contains a DN that indicates the average nighttime light intensity observed within each year. Background noise in the composite were identified and replaced with values of zero during the compositing process, and final data values range from 0 to 63. The version 4 DMSP/OLS NTL contains several products and we use the cleaned stable lights products, which contain lights from cities, towns, and other sites with persistent lighting, including gas flares. Ephemeral events, such as wildfires and lightning, have been discarded.
While the original DMSP/OLS NTL values are recorded in the range [0,63], ETM+ NDVI values are generated in the range [−1.0, 1.0]. To make them comparable, we normalize the DMSP/OLS NTL into the range of [0.0, 1.0] by dividing the raw NTL data by its maximum DN value (63.0). Since the stable DMSP/OLS NTL contain gas flares that are often in remote areas, we mask out gas flares with the gas flare products developed by Elvidge et al. [31,32].
Both ETM+ and NTL images are then automatically reprojected into a common projection (geographic) on the fly in Google Earth Engine. The 1 km NTL is then resampled to a 30 m resolution to match the spatial resolution of ETM+.
We obtain the National Land Cover Database 2006 (NLCD 2006) land cover products and the Impervious Surface Area (ISA) products from the United States Geological Survey (USGS) for developing compositing algorithms and for evaluation purposes. The NLCD products are the only publicly available national land cover products at the Landsat scale and have been strictly evaluated [33,34,35] at the time of the current analysis. To be consistent, we adopt the NLCD definition of urban: areas characterized by a high percentage (30% or greater) of constructed materials (e.g. asphalt, concrete, buildings, etc.) [36].

2.2. ETM+ NDVI Composing

Since its launch in 1999, Landsat-7/ETM+ has accumulated a very dense stack of multiple spectral band imagery at a 30 m spatial resolution. Its opening to the public, free of charge, since 2008 [12] has introduced a new era of global applications using remote sensing data at the 30 m resolution [27,37,38,39]. However, its utilization has been greatly impeded due to cloud and cloud shadow contamination and inadequate algorithms to analyze them. Consequently, the vast majority of the applications are confined to individual small regions [13]. In terms of cloud cover, a Landsat scene (185 km by 185 km) could be in one of the three statuses: clear (no clouds in the scene), partially contaminated (with some clouds but part of the land can still be seen from space), or completely contaminated (with 100% cloud cover). Manually selecting clear Landsat images is the main strategy in most of the studies, which can be very labor-intensive and time-consuming for a large region. However, clear images are often rare, especially in tropical areas with consistent heavy cloud covers [40]. Consequently, good image pixels containing useful information in partially contaminated scenes are often not well exploited.
In order to use clear pixels in partially contaminated images, researchers have developed two major strategies to automatically deal with clouds in remote sensing images: 1) exhausting mapping efforts (direct) and 2) multi-temporal image compositing by selecting only good quality data based on some criteria (indirect). Automated cloud classification methods based on a single Landsat image [41,42,43,44,45,46,47,48] achieved high accuracies in detecting clouds and their shadows. Recent cloud classification efforts based on multi-temporal images [49,50,51,52,53,54,55,56] have been proposed to better detect clouds and cloud shadows. The method proposed in [57] only deals with clouds and ignored cloud shadows, while the method proposed by Zhu et al. [48] is designed to detect clouds and associated shadows simultaneously in Landsat images. Although these existing cloud and cloud shadow classification methods can reach high accuracy, residual clouds, especially thin clouds, still pose challenges. Furthermore, supervised classification methods are often labor-intensive and time-consuming, and thus are not suitable for global applications.
Alternatively, satellite image compositing was developed originally to automatically pick up image pixels free of cloud or cloud shadow contamination in Advanced Very High Resolution Radiometer (AVHRR) time series pixel by pixel. Such a technique aims to produce a composite with good image pixels from different times even if there is no single clear image available during a specific period. It can thus make use of data from all images that are not completely contaminated. The Maximum NDVI Value method and the Maximum Apparent Temperature method [58] are the most commonly used compositing methods. Both methods are based on the observations from coarse resolution AVHRR images that clouds and cloud shadows decrease NDVI values and apparent temperatures over vegetated surfaces. The Maximum NDVI Value method works well over vegetated surfaces but fails over desert and water because clouds can have greater NDVI values than bare land and water [58]. The Maximum Apparent Temperature method assumes that the land surface remains stable during a short period of time, such as 10 days. However, over a longer period, the surface might change dramatically. A pixel can be warmer on a summer day even in a cloud shadow than on a clear winter day. Landsat visits the same place every 16 days, which means that in order to have enough observations from Landsat for a successful composite, a longer time period is required. In our case, we plan to do nominal annual composing using all images over three years. Thus, the Maximum Apparent Temperature method cannot be directly used to generate global cloud and cloud shadow-free Landsat composites.
Despite its limitations, the maximum NDVI compositing method is simple, easy, and fast, and thus has the potential to be applied to generate global NDVI composites and can guarantee good quality data (if it does exist) from growing seasons will be included. This method was modified for generating MODIS vegetation indices by incorporating the smallest view angle constraint to minimize sensor view variations associated with the anisotropic surface reflectance properties, i.e., Bidirectional Reflectance Distribution Function (BRDF) and helps reduce spatial and temporal discontinuities in the composited product [59].
However, the MODIS version does not treat vegetation, bare land, and water separately, either. Here we modify the Maximum NDVI compositing method to generate a global ETM+ NDVI composite on the Google Earth Engine platform through treating vegetation, bare land, and water differently. We begin with a careful examination of the impacts of clouds and cloud shadows over different types of land covers, i.e. vegetated surface, bare land, and water surface, on NDVI.

2.2.1. Impact of Clouds and Cloud Shadows on NDVI over Different Land Surfaces

The effects of clouds and their shadows on modifying land surface NDVI values measured from space highly depend on the underlying land cover types: water, vegetation, and bare land (Figure 1). Here we define the bare land class to include all land surfaces that are neither covered by water nor vegetation, such as deserts and dense built-up areas. Over vegetated land surfaces, NDVI values decrease due to both clouds and cloud shadows (Figure 1). Over water, surface clouds cause NDVI values to increase and cloud shadows have little impact on NDVI values or slightly increase NDVI values (Figure 2). Over bare land surface, cloud shadows lead to lower NDVI values, but the effects of clouds on NDVI values vary with the thickness of clouds: thin clouds cause NDVI values to drop while thick clouds cause NDVI values to increase (Figure 3).
Figure 1. Clouds and cloud shadows and their effects on NDVI. Left: over vegetated land surface, middle: over water surface, right: over bare land surface. These sites are randomly selected. (A) True color Landsat-7/ETM+ image. (B) Landsat-7/ETM+ NDVI image.
Figure 1. Clouds and cloud shadows and their effects on NDVI. Left: over vegetated land surface, middle: over water surface, right: over bare land surface. These sites are randomly selected. (A) True color Landsat-7/ETM+ image. (B) Landsat-7/ETM+ NDVI image.
Remotesensing 07 11887 g001
To summarize, over vegetated surface, cloud and cloud shadow-free observations in growing seasons should show the maximum NDVI values in a time series, while over bare land surface they should have the median NDVI values, and over water surface they should exhibit the minimum NDVI values (Table 1). Based on these empirical observations, we can develop mixed NDVI criteria instead of a single maximum NDVI criterion to retrieve cloud and cloud shadow-free observations over different land surfaces. We will therefore retrieve clear-sky observations with a maximum NDVI value over vegetation, with a median NDVI value over bare land, and with a minimum NDVI value over water (Table 1).
Table 1. Summary of the impacts on NDVI from clouds and cloud shadows over different land cover types.
Table 1. Summary of the impacts on NDVI from clouds and cloud shadows over different land cover types.
Land Cover TypeNDVICloud EffectShadow EffectClear-Sky NDVI
VegetationHighMaximum
Bare landLowMedian
WaterVery lowMinimum

2.2.2. Vegetation, Bare Land, and Water Stratification

Applying different composing criteria, as mentioned above, requires first delineating land surfaces into the three categories: water, bare land, and vegetation. Such a delineation must be simple, fast, and easy so that it can be applied to a large area and, eventually, the entire globe. For this purpose, we chose five 2° × 2° regions covering water, bare land, and vegetated areas in the United States as testing sites (Figure 2). We chose these sites in the United States because the NLCD 2006 land cover products can provide detailed ground information at the Landsat scale.
Figure 2. Five select testing sites (black rectangles) and two evaluation sites (pink rectangles) in the conterminous USA (overlaid on the NLCD 2006 land cover map).
Figure 2. Five select testing sites (black rectangles) and two evaluation sites (pink rectangles) in the conterminous USA (overlaid on the NLCD 2006 land cover map).
Remotesensing 07 11887 g002
We examine the distributions of maximum and minimum NDVI values in the three-year period within each of the three categories in these regions. To do that, we first regroup the NLCD 2006 classes into the three categories: water (Water), bare land (Barren), and vegetation (Forest, Shrubland, Herbaceous, Planted/Cultivated, and Wetlands). At this stage, we decide to provisionally ignore the Developed NLCD class, since they are often mixed classes (built-up areas mixed with vegetation or water) with complex spectral signals. From distributions of maximum and minimum NDVI values within each category, we find that vegetated areas often have a maximum NDVI value greater than 0.4, and the maximum NDVI values of bare land and water are very unlikely greater than 0.4 (dashed lines in Figure 3). Thus, we can safely draw a line at 0.4 to separate vegetation from the other two categories. Similarly, water often has a minimum NDVI value less than −0.2 and the other two categories are very unlikely to have minimum NDVI values lower than −0.2 (solid lines in Figure 3). Using the minimum NDVI value of −0.2, we can safely separate water from the vegetation and bare land. Consequently, we design a set of rules based on these observations to delineate vegetation (max. NDVI > 0.4), water (min. NDVI < −0.2) and bare land (anything else) for applying different criteria on the ETM+ time series to retrieve cloud and cloud shadow-free observations. What needs to be pointed out here is that such stratification is not a rigorous land cover classification, but only a rough guess of ground conditions, with the goal to facilitate picking clear acquisitions in a time series. It is not surprising that under such stratification, some sparsely vegetated areas can have low NDVI values and they must be grouped into the bare land category. For example, Site 4 and Site 5 are in dry climate zones and are dominated by the shrub/scrub class, which leads the vegetation category to have a maximum NDVI value less than 0.4 (Figure 3). Furthermore, the High Intensity Developed NLCD classes might have low NDVI values and they can be grouped into the bare land category, while other Developed NLCD classes with lower impervious intensities might have higher NDVI values and can be grouped into the vegetated land category. Such stratification is indeed provisional, since after the compositing procedure each pixel will be assigned with a cloud and cloud shadow-free observation that can be later used to do a more detailed classification if needed.
Figure 3. Distributions of maximum, minimum NDVI values within each of the water, vegetation, and bare land categories in the five testing sites. Site #1 and #2 are dominated by vegetation (crops or forests), Site #3 is dominated by water, and Site #4 and #5 are dominated by bare land.
Figure 3. Distributions of maximum, minimum NDVI values within each of the water, vegetation, and bare land categories in the five testing sites. Site #1 and #2 are dominated by vegetation (crops or forests), Site #3 is dominated by water, and Site #4 and #5 are dominated by bare land.
Remotesensing 07 11887 g003

2.2.3. Building a Nominal Annual NDVI Composite from Multi-Year ETM+ Imagery

In mapping urban areas, fallow land often poses the greatest challenge, due to its spectral similarity to urban features, such as roads, roofs, etc. In order to distinguish urban areas from fallow lands that are often seasonal, growing season images are often required. However, completely cloud-free Landsat images are usually hard to obtain, especially during growing seasons when precipitation is abundant, which means more cloud cover is present. To maximize the degree of separability between fallow lands and urban areas, we built a nominal annual composite with Landsat images taken within a contiguous three-year period. This is based on the fact that on one hand, fallow lands are normally temporal. An observed fact is that a piece of farmland will be reclaimed after being left as fallow for one or, at most, two years. It is very rare that a piece of farmland remains fallow contiguously for more than three years. In rare cases, a piece of farmland is abandoned contiguously for more than three years and wild grasses often will cover its surface, thus making that piece of farmland easily distinguishable from urban areas. On the other hand, urbanization is often path-dependent, i.e. a piece of land once built up is very unlikely to return back to other types of land within a relatively short period, e.g. three years. Furthermore, taking images from three years can maximize the chance of getting cloud and cloud shadow-free observations, especially in regions that have heavy cloud cover and after the ETM+ SLC failure. We built a nominal annual composite for one specific year with a time series comprising all imagery from that year and the following two years. In that way, consistent urban areas can be successfully separated from temporal fallow lands. We summarize the entire procedure to generate ETM+ NDVI composites on the Google Earth Engine platform in a multi-step flow chart (Figure 4).
Figure 4. The procedure to generate NDVI composites from multi-year Landsat-7/ETM+ imagery with the proposed Mixed NDVI method on Google Earth Engine.
Figure 4. The procedure to generate NDVI composites from multi-year Landsat-7/ETM+ imagery with the proposed Mixed NDVI method on Google Earth Engine.
Remotesensing 07 11887 g004

2.3. The Normalized Difference Urban Index (NDUI)

The profiles of vegetation density and luminosity intensity along transects passing through urban cores often exhibit an inverse relationship: vegetation density decreases while luminosity increases from rural areas to the urban core. Such a relationship has been used to generate indices to characterize urban areas by combining MODIS NDVI, a vegetation index, and DMSP/OLS luminosity data [24,25]. We propose here to replace MODIS NDVI with Landsat-7/ETM+ NDVI to capture urban spatial structures at a much finer scale. We construct a Normalized Difference Urban Index (NDUI) by combining ETM+ NDVI with DMSP/OLS luminosity data:
NDUI = (NTL − NDVI) ∕ (NTL + NDVI),   (NDVI ≥ 0)
where NTL is normalized DMSP/OLS nighttime light data, and NDVI is cloud-free ETM+ NDVI. Since water normally has NDVI values less than 0 and urban areas cannot exist under water, we confine NDVI to not be less than 0 in Equation 1 to make the NDUI mathematically more robust and easy to interpret. According to this definition, well-lit urban cores with very sparse vegetation cover will have high NDUI values close to 1, and unlit remote rural areas with abundant vegetation will have low values close to −1.
In the above definition, a cloud and cloud shadow-free NDVI composite from the growing season is critical in order to maximize the separability between urban and vegetated areas, especially agricultural lands experiencing seasonal fallow periods.
To implement the NDUI on the Google Earth Engine platform, both DMSP/OLS NTL and Landsat-7/ETM+ NDVI composites are reprojected into a common geographic coordinate system (geographic projection). We then resample the 1-km DMSP/OLS normalized NTL data down to 30 m by assigning each 30 m sub grid within the same 1-km DMSP/OLS NTL grid the same normalized NTL value. We finally calculate the NDUI using Equation 1 with the 30-m DMSP/OLS normalized NTL and 30-m ETM+ NDVI. The NDUI is a single band image containing information with urban areas enhanced, very similar to NDVI containing information with vegetation enhanced.

2.4. Evaluating NDUI with NLCD 2006

We evaluate the NDUI within each NLCD Developed class (the urban classes) to examine its capability in characterizing urban areas at the Landsat scale at two evaluating sites in the United States (Figure 2A,B). Site A covers most of the Salt Lake City metropolitan area and part of the Salt Lake in Utah, USA. Urban areas in this region are surrounded by densely forested mountains on the east side and by the Salt Lake and desert on the west side. Site B covers the Chicago metropolitan area and its surrounding regions. Heavily urbanized parts of the Chicago metropolitan area are surrounded by Lake Michigan on the east side and by intense agriculture lands on the west and south sides. We chose these two regions as evaluation sites mainly due to their spatial configurations of urban areas, water, and desert. We argue that their spatial configurations can have a wide representativeness. They can test the proposed NDUI’s capability in distinguishing urban built-up areas from not only vegetated areas, but also bare lands and water. In other words, within these two test sites, we will test whether the NDUI can still enhance urban areas when there is no dense vegetation cover right outside of an urban area, since there are a number of cities located in dry lands worldwide. We will focus on evaluating the NDUI’s separability between urban areas and agricultural lands mainly in site B, while evaluating the separability between urban areas and bare lands (desert) mainly in site A.
We investigate the separability between urban areas and bare lands as well as farmland by visually examining their histograms. We further quantify the separability between urban areas and bare lands using the SDI (Spectral Discrimination Index) [60,61]. The SDI measures the degree of separation between the histograms of two classes, with a focus on the relative locations and spreads of two histograms. We calculate SDIs with each of NDUI, BCI, NDVI, and NDBI individually. For a particular index, SDI values greater than 1 indicate good separability between the two classes, while values less than 1 denote poor separability due to large overlaps.

3. Results

3.1. Landsat-7/ETM+ NDVI Composites

To test the algorithm performances, we first chose a small water area within site B and a bare land area within site A and generate nominal 2006 NDVI composites with the proposed Mixed NDVI criteria and with the Maximum NDVI criteria for comparison (Figure 5). It is clear that in both sites our proposed Mixed NDVI criteria successfully generate cloud- and cloud shadow-free NDVI composites, while the Maximum NDVI criteria fails to exclude all clouds and cloud shadows. Furthermore, gaps due to SLC-off in the ETM+ images over both sites are also filled in the composites generated with the proposed Mixed NDVI criteria, while they are still visible in the composites generated with the Maximum NDVI criteria.
Figure 5. Nominal 2006 NDVI composites created with the Maximum NDVI criteria and the Mixed NDVI criteria over a water site (upper row) and a bare land site (bottom row).
Figure 5. Nominal 2006 NDVI composites created with the Maximum NDVI criteria and the Mixed NDVI criteria over a water site (upper row) and a bare land site (bottom row).
Remotesensing 07 11887 g005
For a large region, we show a 2006 ETM+ nominal annual NDVI composite built with three-year ETM+ imagery (2006−2008) in the conterminous US as an example generated with the Mixed NDVI composing framework (Figure 6). In this NDVI composite, the first thing we want to mention is that clouds and cloud shadows have been successfully removed. Second, vegetated areas including agriculture lands, e.g. the Great Plains, have high NDVI values. Third, urban areas, such as the Chicago region and the New York City region, surrounded by vegetated lands and water, are now nicely distinguishable from their vicinity areas. Fourth, dry regions exhibit very small NDVI values, and thus look very similar to urban areas, i.e. it is hard to distinguish urban areas from other bare lands with NDVI only. Fifth, contrasts between neighboring Landsat paths, especially in the eastern part of the country, are discernible.
Figure 6. Landsat-7/ ETM+ NDVI composite of the conterminous US in 2006.
Figure 6. Landsat-7/ ETM+ NDVI composite of the conterminous US in 2006.
Remotesensing 07 11887 g006

3.2. Date Distribution in the NDVI Composites

As mentioned earlier, an important purpose of building the Landsat-7/ETM+ NDVI composites is to collect good quality acquisitions during the growing season to maximize the separation between vegetation and urban built-up areas. An examination of the distribution of dates when clear acquisitions that were entered into the final composite were taken provides a way to check the quality of that composite. We calculated the distribution of dates for all vegetated pixels in the 2006 Chicago composite (Figure 7A). For comparison, we also calculated the distribution of dates for all vegetated pixels in the 2000 Chicago composite (Figure 7A) to examine the effects of gaps due to SLC failure on the composites. It is clear the majority of dates fall in between the end of April (Day of Year (DOY) 120) and early October (DOY 280), roughly the growing season in the Chicago area. A comparison between 2000 and 2006 reveals that the date distribution in 2000 is slightly tighter than in 2006, suggesting that gaps due to SLC failure do reduce the quality of the final ETM+ composite to some degree. A two-sample Kolmogorov–Smirnov test also confirms that the differences between the two date distributions are statistically significant, with a D statistic equal to 0.2855 (the p-value is 2.2e−16 and the total numbers of valid pixels over the evaluation sites are 21,032,235 and 20,732,443). Another fact that can be observed is that the final composites contain data from all three years in each compositing period (Figure 7C,D). Multiple peaks can be observed in the date distributions.
Figure 7. Distributions of dates when good quality data were collected that went into the 2000 and 2006 Chicago NDVI composites. (A) Date distribution without year differentiation; (B) cumulative curves of A; (C) date distribution with year differentiation in the 2000 composite; and (D) date distribution with year differentiation in the 2006 composite.
Figure 7. Distributions of dates when good quality data were collected that went into the 2000 and 2006 Chicago NDVI composites. (A) Date distribution without year differentiation; (B) cumulative curves of A; (C) date distribution with year differentiation in the 2000 composite; and (D) date distribution with year differentiation in the 2006 composite.
Remotesensing 07 11887 g007
In the final 2006 NDVI composite, cloud- and cloud shadow-free data may be acquired in all three years. This is partly because of the gaps due to SLC failure. For example, there are some grey straps can be observed in the date map (Figure 8A) with shapes very similar to those of the gaps caused by the SLC-off. Furthermore, the white patchiness in a subset of the acquisition date map in the Chicago region (Figure 8A) strongly suggests that it might also be a result of the agricultural practices occurring in the field (Figure 8D–E). It is obvious that there is heterogeneity in agricultural practices, such as different cropping, in this region. However, despite that heterogeneity, the Mixed NDVI compositing procedure successfully picked the growing season NDVI values. For example, the majority of the NDVI values for vegetated pixels in each year are larger than 0.8 (Figure 8C).
Figure 8. Acquisition years and NDVI values in the 2006 composite in a subregion in the Chicago area. (A) Acquisition years: black: 2006, grey: 2007, and white: 2008; (B) NDVI values increase from black to white; (C) Distribution of NDVI; (DF) Google Earth snapshots.
Figure 8. Acquisition years and NDVI values in the 2006 composite in a subregion in the Chicago area. (A) Acquisition years: black: 2006, grey: 2007, and white: 2008; (B) NDVI values increase from black to white; (C) Distribution of NDVI; (DF) Google Earth snapshots.
Remotesensing 07 11887 g008

3.3. Characterizing Urban Areas by NDUI

We chose the Salt Lake City region as our testing site to examine the capability of the NDUI in characterizing and mapping urban areas (Figure 9). An example NDUI map from the Chicago region is shown in Figure 9C. In this map, spatial details can be easily identified. For example, major roads radiating out from downtown are very distinct from the surroundings. It is clear that these spatial details, which are lacking in the NTL image (Figure 9B), come from the Landsat-7/ETM+ NDVI composite (Figure 9A). Visually, the NDUI map is very comparable to the NLCD land cover map (Figure 9D), suggesting that the NDUI successfully enhances urban areas.
The NDUI captures not only regional urban areas but also local urban areas with fine details. For example, road arterials can be easily identified within Salt Lake City (Figure 10). The downtown area, the commercial strip, the two industrial zones, and the Salt Lake City International Airport all pop up and can be easily distinguished from their surroundings (Figure 10).
Figure 9. The 2006 NDVI (A), NTL (B), NDUI (C), and NLCD (D) maps in the Salt Lake City area, Utah.
Figure 9. The 2006 NDVI (A), NTL (B), NDUI (C), and NLCD (D) maps in the Salt Lake City area, Utah.
Remotesensing 07 11887 g009
Figure 10. A close-up look at the urban areas of the Salt Lake City region. (A) Downtown of Salt Lake City; (B) commercial strip; (C) industrial zone; (D) Salt Lake City International Airport; (E) industrial Zone; (F) sand pit at the Point of the Mountain; (G) Bingham Canyon Open-Pit Copper Mine.
Figure 10. A close-up look at the urban areas of the Salt Lake City region. (A) Downtown of Salt Lake City; (B) commercial strip; (C) industrial zone; (D) Salt Lake City International Airport; (E) industrial Zone; (F) sand pit at the Point of the Mountain; (G) Bingham Canyon Open-Pit Copper Mine.
Remotesensing 07 11887 g010

3.4. Separating Urban Area from Agricultural Land

The Chicago region in focus (Figure 2B) can be roughly grouped into three categories: cultivated crops that occupy more than 30% of the region, water body (part of Lake Michigan), and urban areas (the Chicago Metro area). This can be well revealed by the Landsat NDVI composite (Figure 11A). The Mixed NDVI compositing procedure ensures NDVI values from growing seasons and removes fallow land as much as possible. Consequently, there exists a distinct break between bare land (mainly urban areas in this region) and vegetated land (mainly agricultural land) (Figure 11A). Analysis of the NDUI further confirms the separability between the urban area and agricultural land (Figure 11B). All NLCD Developed classes concentrate in the region with positive NDUI values, while the Cultivated Crops class concentrates in the region with NDUI values around −0.8. It is noticed that parts of the Developed, Open Space class and the Developed, Low Intensity class appear in the region with negative NDUI values. This might be because, in the rural areas, there are unlit rural roads, highways, and farmhouses, which have been classified as Developed classes in the NLCD 2006 land use land cover product.
Figure 11. The distribution of NDVI values in the 2006 NDVI composite (A) and the distribution of the 2006 NDUI within each NLCD 2006 Developed class and the Cultivated Crops class, excluding water, in the Chicago region (B).
Figure 11. The distribution of NDVI values in the 2006 NDVI composite (A) and the distribution of the 2006 NDUI within each NLCD 2006 Developed class and the Cultivated Crops class, excluding water, in the Chicago region (B).
Remotesensing 07 11887 g011
A close look at the NDUI values reveals some interesting patterns. The bimodal distribution of the NDUI of the Developed, Medium Intensity class is mainly due to the bimodal distribution of NDVI within that class. For example, street trees are common in Chicago, covering parts of streets and leaving others bare across the entire city. This leads to very different NDVI values even along the same individual street, especially within residential neighborhoods. Many residential streets are much narrower than 30 meters, which leads to mixed Landsat pixels covering street trees, street pavement, and possible rooftops. Consequently, some developed pixels can have relatively high NDVI values.

3.5. Quantitative Assessment of NDUI

To quantitatively examine the NDUI’s capability in enhancing urban features, we compare it with the NLCD2006 ISA product, which is a measure of built-up area density. We extract data values of the NDUI and NLCD ISA from a latitudinal transect in Salt Lake City and plot them together (Figure 12A). We also build a linear model to see if a statistically strong linear relationship between the two exists (Figure 12B). Visually, the latitudinal profiles of the NDUI and NLCD ISA generally match well, while the quantitative assessment only shows a fairly strong positive relationship between them with a R2 value of 0.227. A close look at the profiles reveals that there are some discrepancies between them. The largest discrepancy appears inside the black dotted rectangle (Figure 12A), where most of that piece of land was bare (not built) in 2006, as confirmed by high resolution Google Earth images (Figure 12C). The bare proportion has low NDVI values but was influenced by the night lights from its neighboring built-up land, which led to a high NDUI value while the ISA value was 0. The R2 value increased to 0.372 (y = 0.6998x + 0.1032, R2 = 0.3722) after we removed data points in that region. Pikes and valleys in the two profiles generally match; however, their spatial locations can shift by one or two pixels in many cases. That might be attributed to the low geo-referencing accuracy of some Landsat ETM+ scenes [39], as the NDUI was calculated with the NDVI composite generated from multi-year scenes, and thus there is a chance the geometrical errors in some scenes got transferred into the final product. Another factor attributing to the relatively low R2 value is the spatial distribution of vegetation within urban areas, as discussed earlier. Consequently, each Developed class exhibits a bimodal distribution, which means there is a chance that some pixels with the same ISA value might have different NDUI values.
Figure 12. The latitudinal profiles of NDUI and NLCD ISA in Salt Lake City (A), and the linear regression model between them (B). Panel (C) shows a Google Earth snapshot, which reveals that the piece of land within the dotted rectangle was bare in 2006. The upper graph of panel A is a false color composite (Red: NDUI, Green: NDUI, Blue: NLCD ISA).
Figure 12. The latitudinal profiles of NDUI and NLCD ISA in Salt Lake City (A), and the linear regression model between them (B). Panel (C) shows a Google Earth snapshot, which reveals that the piece of land within the dotted rectangle was bare in 2006. The upper graph of panel A is a false color composite (Red: NDUI, Green: NDUI, Blue: NLCD ISA).
Remotesensing 07 11887 g012

3.6. Comparison with Other Indices in Terms of Separating Urban Areas from Bare Lands

We calculate the distributions of the NDUI, BCI, NDVI, and the Normalized Difference Built Up Index (NDBI) [62] within each of the barren land classes and the Developed classes in the two evaluation sites in the Salt Lake City region (Figure 9D 1−2). The Developed classes are well separated from the barren land class by the NDUI, with the majority of barren land class having NDUI values around −1, while the majority of Developed classes have NDUI values larger than 0 (Figure 13A). Furthermore, SDI values between the barren land class and all the Developed classes are all greater than 1 (Table 2), confirming the good separability between the two groups. The two groups are also well separated by BCI, with the majority of the Barren Land class having BCI values less than −0.25, while the majority of Developed classes have BCI values greater than −0.25. However, a detailed examination of SDIs shows that there exists some degree of confusion between the barren land class and the Developed, Open Space class with a SDI value of 0.78, as well as the Developed, Low Intensity class with a SDI value of 0.94 (Table 2). In contrast, the Developed classes cannot be separated from the barren land class by either NDVI or NDBI (Figure 13C,D). The barren land class totally overlaps the Developed, High Intensity class, with SDI values lower than 1 (Table 2).
The right, smaller peaks of the barren land class curve in both NDUI and BCI are due to the Bingham Canyon Open-Pit Copper Mine southwest of Salt Lake City (Figure 10G) and a large sand pit at the Point of the Mountain between Salt Lake City and Provo (Figure 10F). The mine has been in production since 1906, and has resulted in the creation of a pit over 0.6 miles (0.97 km) deep, 2.5 miles (4 km) wide, and covering 1900 acres (770 ha). Both sites are still used for industrial purposes and are well lit at night though they are classified as barren land in the NLCD 2006 land use land cover product.
Figure 13. Distributions of NDUI (A), BCI (B), NDVI (C), and NDBI (D) within both the NLCD 2006 barren land class and the Developed classes in the two evaluation sites.
Figure 13. Distributions of NDUI (A), BCI (B), NDVI (C), and NDBI (D) within both the NLCD 2006 barren land class and the Developed classes in the two evaluation sites.
Remotesensing 07 11887 g013
Table 2. Separability between the barren land class and the Developed classes as measured by Spectral Discrimination Index (SDI) generated from different indices.
Table 2. Separability between the barren land class and the Developed classes as measured by Spectral Discrimination Index (SDI) generated from different indices.
IndexOpen SpaceLowMediumHigh
NDUI1.271.541.812.17
CI0.780.941.091.30
NDBI0.480.560.490.12
NDVI1.411.441.020.36
We also noticed that some Developed pixels have negative NDUI values (Figure 13A). A further check of these pixels reveals that these are not well lit, with either low or high NDVI values (Figure 14), indicating a failure of the NDUI to capture built-up areas under such conditions.
Figure 14. Some not-well-lit Developed areas cannot be well captured in the 2006 NDUI in the Salt Lake City region.
Figure 14. Some not-well-lit Developed areas cannot be well captured in the 2006 NDUI in the Salt Lake City region.
Remotesensing 07 11887 g014

4. Discussion

Our assessment shows that the NDUI successfully minimizes confusion between urban areas and bare lands, including fallows and barren lands, thus enhancing urban areas at the regional scale. This success comes from the fact that the NDUI combines the composited Landsat NDVI and DMSP/OLS NTL. The composited Landsat NDVI ensures that agricultural lands all have maximum NDVI values taken in growing seasons to minimize the impact on urban detection from seasonal fallow lands. Since deserts and other barren lands not related to urban uses are normally not lit, DMSP/OLS NTL can effectively separate them from urban areas. Furthermore, we show that the NDUI can not only characterize urban areas at the regional scale, but also reveals their fine spatial details at a 30 m resolution, indicating the success of combining medium resolution Landsat-7/ETM+ with coarse resolution DMSP/OLS NTL. The fine spatial details lacking in DMSP/OLS NTL are well compensated by the medium resolution Landsat NDVI. Although there are other indices such as BCI that can similarly enhance urban areas, their requirement for first masking out water makes them less competitive than the proposed NDUI, which can be highly automated and scalable.
The Mixed NDVI value compositing procedure is a value-adding process. We have shown that it can effectively remove not only clouds and cloud shadows, but also gaps caused by SLC-off in the ETM+ imagery since May 2003. This strategy ensures the picking of growing season pixels as much as possible, which thus relieves the burden of finding cloud-free images in the growing season which normally has more rainy days. First, the growing season composite forms a three-mode NDVI distribution (Figure 11). The three peaks roughly correspond to the three broad land categories: water, vegetated land, and non-vegetated land. The implication of this is that it simplifies surface characterization that could otherwise be very complicated in an image taken at a time not in the growing season. Consequently, non-vegetated lands such as urban areas can be easily separated from water and vegetated lands. Second, it fills gaps due to SLC failure. The reason that these gaps can be filled by the Mixed NDVI compositing procedure may be attributed to slight positional shifting of the scan lines in the orbit at different revisits. Although there are about 22% pixels missing in each ETM+ scene since May 2003, the locations of the gaps can be slightly different during different visits. This gives the sensor the precious opportunity that every single piece of the global land surface has the chance to be scanned during a long period, such as three consecutive years. Furthermore, the Mixed NDVI value compositing method is simple, fast, and can be easily extended to generate cloud- and cloud shadow-free Landsat NDVI composites at the global scale.
We generated Landsat-7/ETM+ NDVI composites from original Landsat imagery taken during three consecutive years to: 1) minimize the influence from fallows; and 2) improve the quality of the final NDVI composites. The gaps due to SLC-off in Landsat-7/ETM+ imagery are successfully filled when we use all the imagery from three years to increase the chance of getting cloud- and cloud shadow-free and gap-free observations. The comparison between the date distribution of the 2000 and 2006 Landsat NDVI composites shows that the quality of the final NDVI composite will be improved if we have more gap-free observations (Figure 7). The number of total observations in a specific period will increase if we can also include Landsat-5/TM imagery or other Landsat-like imagery, such as that from the Chinese Ziyuan-3 satellite. However, that will require inter-satellite calibration and accurate cross-imagery geo-referencing. As we have shown, slight geo-referencing errors in some scenes have affected the NDUI’s capability to quantitatively characterize urban areas.
Land surface stratification is critical to the success of the proposed composing method. We roughly grouped land surface into three categories. We have to point out that such stratification is not intended to do a very accurate land cover classification. Instead, we only need a rough guess of the ground condition. Furthermore, the vegetated category might not include land surface with sparse vegetation cover. Sparse vegetation is treated as bare land under the current composing framework, due to the low maximum NDVI values they can have. We then have to apply the “median NDVI” strategy to retrieve clear-sky observations over those surfaces. Examples are Site 4 and Site 5 (Figure 2), where sparse vegetation is dominant.
This Mixed NDVI compositing method can also be applied to time series data obtained from other remote sensors, such as MODIS, to generate cloud- and cloud shadow-free composites. MODIS covers the entire globe twice a day. Such a dense time series is ideal for generating cloud- and cloud shadow-free satellite image composites.
Although the NDUI is simple and easy to implement, there are still challenges. First, edges between scenes are discernible in the final Landsat-7/ETM+ NDVI composites, mainly due to BRDF effects. Second, gaps caused by SLC-off can be filled, but with data from different dates, which means there can exist a phenology difference in the final Landsat NDVI composites. Landsat-7 needs at least 16 days to completely scan the entire globe, which means that even under perfect weather conditions (no clouds globally during a consecutive 16-day period), a global Landsat-7 composite will still contain a phenology difference. Biases caused by vegetation phenology differences can be even larger than those caused by BRDF effects and are much harder to correct. Due to this reason, we did not correct the BRDF effects during our Landsat-7/ETM+ NDVI compositing procedure, but we will focus on this issue in our future research. Third, we did not perform atmospheric correction on the Landsat-7/ETM+ imagery. Although the compositing method based on NDVI values can effectively minimize atmospheric effects, there is still a chance that residual aerosols can exist in the final NDVI composites. Future efforts to correct all these biases are highly needed. Finally, substantial quantities of impervious surfaces can exist within rural areas, including rural roads, highways, and farmhouses, which might not all be well lit and thus can be missed by the NDUI (Figure 14). Some urban areas are in darkness due to the lack of electricity and humanitarian disasters [63,64,65] and non-urban bare land can be influenced by night lights, which will also lead to underestimation or overestimation accordingly.
Our evaluation results show the NDUI’s potential for regional and global scale urban characterization. However, we have selected sites only in the US to build and evaluate our models. Although these sites were chosen with the goal to have broad representativeness, land cover and land use worldwide can vary dramatically. Furthermore, energy consumption patterns and urban spatial configurations outside the US can also be very different from those inside the US. To extend NDUI use to other regions, it will be necessary to have more testing and evaluation of sites outside the US.

5. Conclusions

The proposed Normalized Difference Urban Index (NDUI) combines medium resolution Landsat NDVI composites with coarse resolution DMSP/OLS nighttime light (NTL) images to better characterize urban areas, easily distinguishing them from non-urban land features at a much finer resolution. To facilitate the construction of the NDUI, we propose a new compositing method to generate cloud- and cloud shadow-free Landsat NDVI composites, by treating water, vegetation, and non-vegetation land surfaces individually with different compositing criteria. Our assessments of the NDUI confirm its ability to bring NTL down to the 30-m scale to reveal urban spatial structures with finer details. Moreover, its simplicity enables rapid characterization and monitoring of urban areas globally. With more than 20 years of the DMSP/OLS NTL data archive and more than 30 years of the Landsat data archive available, the NDUI provides an efficient practice for characterizing historical urban form and spatial structures. The NDUI is a response to urgent calls for simple, timely, and relevant remote sensing measures that help inform practical decision-making about the built environment.
Nevertheless, there are still challenges associated with the NDUI. The NDUI’s capability to characterize urban areas highly depends on night lighting. Consequently, built-up areas that are not lit during the night cannot be detected by the NDUI, and non-urban bare land that is lit during the night will be hard to be distinguished from urban areas by the NDUI. Furthermore, some Landsat scenes with less geo-referencing accuracy in a time series will downgrade the quality of the final NDVI composites and eventually the quality of the NDUI.
The new generation night light sensor, the Visible/Infrared Imager/Radiometer Suite (VIIRS), aboard the Suomi Preparatory Project Satellite (NPP) has been collecting data in the Day/Night Bands (DNB) spanning the visible and infrared regions since October 28 2011 [66]. Compared with DMSP/OLS, VIIRS has a broader dynamic radiometric range and advanced onboard calibration facilities, which enable it to take continuous and consistent measurements of nighttime lights free of saturation. As a successor to Landsat-7/ETM+, the Landsat-8/OLI (Operational Land Imager) has improved sensor signal-to-noise performance and associated improvements in radiometric resolution, and an improved duty cycle that allows the collection of a significantly greater number of images per day [67]. The availability of high quality data from both VIIRS and OLI will improve the NDUI and extend its application in the future.

Acknowledgments

This study was supported by Shenzhen Key Basic Research Project (JC201005270334A) and the 2015 Scientific and Technological Activities Preferential Funding Program for Returning Scholars, Ministry of Human Resources and Social Security of the People’s Republic of China. The authors would like to thank Conghe Song for his valuable comments on quantitative evaluation of NDUI. We thank the three anonymous reviewers whose suggestions helped to improve the clarity of the manuscript.

Author Contributions

Qingling Zhang designed the study, performed most of the analysis, and led the writing. Bin Li, David Thau, and Rebecca Moore equally contributed to some of the analysis and writing.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Mills, G. Cities as agents of global change. Int. J. Climatol. 2007, 27, 1849–1857. [Google Scholar] [CrossRef]
  2. Akbari, H.; Menon, S.; Rosenfeld, A. Global cooling: Increasing world-wide urban albedos to offset CO2. Clim. Chang. 2009, 94, 275–286. [Google Scholar] [CrossRef]
  3. World Urbanization Prospects: The 2011 Revision. Available online: http://esa.un.org/unpd/wpp/ (accessed on 30 June 2015).
  4. Gutman, P. Ecosystem services: Foundations for a new rural–urban compact. Ecol. Econ. 2007, 62, 383–387. [Google Scholar] [CrossRef]
  5. Nakićenović, N. Summary for Policy Makers. In Global Energy Assessment: Toward a Sustainable Future; Cambridge University Press: Hong Kong, China, 2012; pp. 16–18. [Google Scholar]
  6. Hoornweg, D.; Bhada, P.; Freire, M.; Trejos, C.; Sugar, L. Cities and Climate Change: An Urgent Agenda; World Bank: Washington, DC, USA, 2010. [Google Scholar]
  7. Seitzinger, S.P.; Svedin, U.; Crumley, C.; Steffen, W.; Adbullah, S.A.; Alfsen, C.; Broadgate, W.J.; Biermann, F.; Bondre, N.R.; Dearing, J.A. Planetary stewardship in an urbanizing world. Ambio 2012, 41, 787–794. [Google Scholar] [CrossRef] [PubMed]
  8. World Population Prospects: The 2010 Revision. Available online: http://www.un.org/en/development/desa/publications/world-population-prospects-the-2010-revision.html (accessed on 30 June 2015).
  9. Scott, A.J. World development report 2009: Reshaping economic geography. J Econ. Geogr. 2009, 9, 583–586. [Google Scholar] [CrossRef]
  10. Gurney, K.R.; Razlivanov, I.; Song, Y.; Zhou, Y.; Benes, B.; Abdul-Massih, M. Quantification of fossil fuel CO2 emissions on the building/street scale for a large US city. Environ. Sci. Technol. 2012, 46, 12194–12202. [Google Scholar] [CrossRef] [PubMed]
  11. Potere, D.; Schneider, A. A critical look at representations of urban areas in global maps. GeoJournal 2007, 69, 55–80. [Google Scholar] [CrossRef]
  12. Woodcock, C.E.; Allen, R.; Anderson, M.; Belward, A.; Bindschadler, R.; Cohen, W.; Gao, F.; Goward, S.N.; Helder, D.; Helmer, E. Free access to landsat imagery. Science 2008, 320. [Google Scholar] [CrossRef] [PubMed]
  13. Seto, K.C.; Fragkias, M.; Güneralp, B.; Reilly, M.K. A meta-analysis of global urban land expansion. PLoS ONE 2011, 6. [Google Scholar] [CrossRef] [PubMed]
  14. Forster, B. An examination of some problems and solutions in monitoring urban areas from satellite platforms. Int. J. Remote Sens. 1985, 6, 139–151. [Google Scholar] [CrossRef]
  15. Hepner, G.F.; Houshmand, B.; Kulikov, I.; Bryant, N. Investigation of the integration of AVIRIS and IFSAR for urban analysis. Photogramm. Eng. Remote Sens. 1998, 64, 813–820. [Google Scholar]
  16. Evaluation of Thematic Mapper Data for Determining Urban Land Cover. Available online: http://www.alwelaie.com/website/universitytheses_details_print.php?theses_id=18148 (accessed on 30 June 2015).
  17. Deng, C.; Wu, C. BCI: A biophysical composition index for remote sensing of urban environments. Remote Sens. Environ. 2012, 127, 247–259. [Google Scholar] [CrossRef]
  18. Liu, C.; Shao, Z.; Chen, M.; Luo, H. Mndisi: A multi-source composition index for impervious surface area estimation at the individual city scale. Remote Sens. Lett. 2013, 4, 803–812. [Google Scholar] [CrossRef]
  19. Small, C.; Pozzi, F.; Elvidge, C.D. Spatial analysis of global urban extent from DMSP-OLS night lights. Remote Sens.Enviro. 2005, 96, 277–291. [Google Scholar] [CrossRef]
  20. Imhoff, M.L.; Lawrence, W.T.; Stutzer, D.C.; Elvidge, C.D. A technique for using composite dmsp/ols “city lights” satellite data to map urban area. Remote Sens. Environ. 1997, 61, 361–370. [Google Scholar] [CrossRef]
  21. Small, C.; Elvidge, C.D.; Balk, D.; Montgomery, M. Spatial scaling of stable night lights. Remote Sens. Environ. 2011, 115, 269–280. [Google Scholar] [CrossRef]
  22. Zhou, Y.; Smith, S.J.; Elvidge, C.D.; Zhao, K.; Thomson, A.; Imhoff, M. A cluster-based method to map urban area from dmsp/ols nightlights. Remote Sens. Environ. 2014, 147, 173–185. [Google Scholar] [CrossRef]
  23. Elvidge, C.D.; Cinzano, P.; Pettit, D.; Arvesen, J.; Sutton, P.; Small, C.; Nemani, R.; Longcore, T.; Rich, C.; Safran, J. The nightsat mission concept. Int. J. Remote Sens. 2007, 28, 2645–2670. [Google Scholar] [CrossRef]
  24. Lu, D.; Tian, H.; Zhou, G.; Ge, H. Regional mapping of human settlements in southeastern China with multisensor remotely sensed data. Remote Sens. Environ. 2008, 112, 3668–3679. [Google Scholar] [CrossRef]
  25. Zhang, Q.; Schaaf, C.; Seto, K.C. The vegetation adjusted ntl urban index: A new approach to reduce saturation and increase variation in nighttime luminosity. Remote Sens. Environ. 2013, 129, 32–41. [Google Scholar] [CrossRef]
  26. Shao, Z.; Liu, C. The integrated use of DMSP-OLS nighttime light and modis data for monitoring large-scale impervious surface dynamics: A case study in the Yangtze River delta. Remote Sens. 2014, 6, 9359–9378. [Google Scholar] [CrossRef]
  27. Hansen, M.; Potapov, P.; Moore, R.; Hancher, M.; Turubanova, S.; Tyukavina, A.; Thau, D.; Stehman, S.; Goetz, S.; Loveland, T. High-resolution global maps of 21st-century forest cover change. Science 2013, 342, 850–853. [Google Scholar] [CrossRef] [PubMed]
  28. Lee, D.S.; Storey, J.C.; Choate, M.J.; Hayes, R.W. Four years of landsat-7 on-orbit geometric calibration and performance. IEEE Trans. Geosci. Remote Sens. 2004, 42, 2786–2795. [Google Scholar] [CrossRef]
  29. Markham, B.L.; Storey, J.C.; Williams, D.L.; Irons, J.R. Landsat sensor performance: History and current status. IEEE Trans. Geosci. Remote Sens. 2004, 42, 2691–2694. [Google Scholar] [CrossRef]
  30. Chander, G.; Markham, B.L.; Helder, D.L. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sens. Environ. 2009, 113, 893–903. [Google Scholar] [CrossRef]
  31. National Geophysical Data Center at the National Oceanic and Atmospheric Administration. Available online: http://www.ngdc.noaa.gov/dmsp/downloadV4composites.html (accessed on 24 January 2013).
  32. Elvidge, C.D.; Ziskin, D.; Baugh, K.E.; Tuttle, B.T.; Ghosh, T.; Pack, D.W.; Erwin, E.H.; Zhizhin, M. A fifteen year record of global natural gas flaring derived from satellite data. Energies 2009, 2, 595–622. [Google Scholar] [CrossRef]
  33. Homer, C.; Huang, C.; Yang, L.; Wylie, B.; Coan, M. Development of a 2001 national landcover database for the United States. Photogramm. Eng. Remote Sens. 2004, 70, 829–840. [Google Scholar] [CrossRef]
  34. Wickham, J.D.; Stehman, S.V.; Gass, L.; Dewitz, J.; Fry, J.A.; Wade, T.G. Accuracy assessment of NLCD 2006 land cover and impervious surface. Remote Sens. Environ. 2013, 130, 294–304. [Google Scholar] [CrossRef]
  35. Xian, G.; Homer, C.; Fry, J. Updating the 2001 national land cover database land cover classification to 2006 by using landsat imagery change detection methods. Remote Sens. Environ. 2009, 113, 1133–1147. [Google Scholar] [CrossRef]
  36. Vogelmann, J.E.; Sohl, T.L.; Campbell, P.V.; Shaw, D.M. Regional land cover characterization using landsat thematic mapper data and ancillary data sources. Environ. Monit. Assess. 1998, 51, 415–428. [Google Scholar] [CrossRef]
  37. Townshend, J.R.; Masek, J.G.; Huang, C.; Vermote, E.F.; Gao, F.; Channan, S.; Sexton, J.O.; Feng, M.; Narasimhan, R.; Kim, D. Global characterization and monitoring of forest cover using Landsat data: Opportunities and challenges. Int. J. Digit. Earth 2012, 5, 373–397. [Google Scholar] [CrossRef]
  38. Loveland, T.R.; Dwyer, J.L. Landsat: Building a strong future. Remote Sens. Environ. 2012, 122, 22–29. [Google Scholar] [CrossRef]
  39. Gong, P.; Wang, J.; Yu, L.; Zhao, Y.; Zhao, Y.; Liang, L.; Niu, Z.; Huang, X.; Fu, H.; Liu, S. Finer resolution observation and monitoring of global land cover: First mapping results with Landsat TM and ETM+ data. Int. J. Remote Sens. 2013, 34, 2607–2654. [Google Scholar] [CrossRef]
  40. Ju, J.; Roy, D.P. The availability of cloud-free Landsat ETM+ data over the conterminous United States and globally. Remote Sens. Environ. 2008, 112, 1196–1211. [Google Scholar] [CrossRef]
  41. Huang, C.; Goward, S.N.; Masek, J.G.; Thomas, N.; Zhu, Z.; Vogelmann, J.E. An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks. Remote Sens. Environ. 2010, 114, 183–198. [Google Scholar] [CrossRef]
  42. Huang, C.; Thomas, N.; Goward, S.N.; Masek, J.G.; Zhu, Z.; Townshend, J.R.G.; Vogelmann, J.E. Automated masking of cloud and cloud shadow for forest change analysis using Landsat images. Int. J. Remote Sens. 2010, 31, 5449–5464. [Google Scholar] [CrossRef]
  43. Irish, R.R.; Barker, J.L.; Goward, S.N.; Arvidson, T. Characterization of the Landsat-7 ETM+ automated cloud-cover assessment (ACCA) algorithm. Photogramm. Eng. Remote Sens. 2006, 72, 1179–1188. [Google Scholar] [CrossRef]
  44. Masek, J.G.; Vermote, E.F.; Saleous, N.E.; Wolfe, R.; Hall, F.G.; Huemmrich, K.F.; Gao, F.; Kutler, J.; Lim, T.K. A landsat surface reflectance dataset for north America, 1990–2000. IEEE Geosci. Remote Sens. Lett. 2006, 3, 68–72. [Google Scholar] [CrossRef]
  45. Oreopoulos, L.; Wilson, M.J.; Várnai, T. Implementation on Landsat data of a simple cloud-mask algorithm developed for MODIS land bands. IEEE Geosci. Remote Sens. Lett. 2011, 8, 597–601. [Google Scholar] [CrossRef]
  46. Roy, D.P.; Ju, J.; Kline, K.; Scaramuzza, P.L.; Kovalskyy, V.; Hansen, M.; Loveland, T.R.; Vermote, E.; Zhang, C. Web-enabled Landsat Data (WELD): Landsat ETM+ composited mosaics of the conterminous United States. Remote Sens. Environ. 2010, 114, 35–49. [Google Scholar] [CrossRef]
  47. Scaramuzza, P.L.; Bouchard, M.A.; Dwyer, J.L. Development of the landsat data continuity mission cloud-cover assessment algorithms. IEEE Trans. Geosci. Remote Sens. 2012, 50, 1140–1154. [Google Scholar] [CrossRef]
  48. Zhu, Z.; Woodcock, C.E. Object-based cloud and cloud shadow detection in landsat imagery. Remote Sens. Environ. 2012, 118, 83–94. [Google Scholar] [CrossRef]
  49. Derrien, M.; le Gléau, H. Improvement of cloud detection near sunrise and sunset by temporal-differencing and region-growing techniques with real-time SEVIRI. Int. J. Remote Sens. 2010, 31, 1765–1780. [Google Scholar] [CrossRef]
  50. Goodwin, N.R.; Collett, L.J.; Denham, R.J.; Flood, N.; Tindall, D. Cloud and cloud shadow screening across Queensland, Australia: An automated method for Landsat TM/ETM+ time series. Remote Sens. Environ. 2013, 134, 50–65. [Google Scholar] [CrossRef]
  51. Hagolle, O.; Huc, M.; Pascual, D.V.; Dedieu, G. A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENμS, LANDSAT and SENTINEL-2 images. Remote Sens. Environ. 2010, 114, 1747–1755. [Google Scholar] [CrossRef] [Green Version]
  52. Jin, S.; Homer, C.; Yang, L.; Xian, G.; Fry, J.; Danielson, P.; Townsend, P.A. Automated cloud and shadow detection and filling using two-date Landsat imagery in the USA. Int. J. Remote Sens. 2013, 34, 1540–1560. [Google Scholar] [CrossRef]
  53. Liu, R.; Liu, Y. Generation of new cloud masks from MODIS land surface reflectance products. Remote Sens. Environ. 2013, 133, 21–37. [Google Scholar] [CrossRef]
  54. Lyapustin, A.; Wang, Y.; Frey, R. An automatic cloud mask algorithm based on time series of MODIS measurements. J. Geophys. Res. Atmos. 2008, 113. [Google Scholar] [CrossRef]
  55. Tseng, D.C.; Tseng, H.T.; Chien, C.L. Automatic cloud removal from multi-temporal spot images. App. Math. Comp. 2008, 205, 584–600. [Google Scholar] [CrossRef]
  56. Wang, B.; Ono, A.; Muramatsu, K.; Fujiwarattt, N. Automated detection and removal of clouds and their shadows from Landsat TM images. IEICE Trans. Inf. Syst. 1999, 82, 453–460. [Google Scholar]
  57. Roy, D.P.; Ju, J.; Lewis, P.; Schaaf, C.; Gao, F.; Hansen, M.; Lindquist, E. Multi-temporal MODIS–Landsat data fusion for relative radiometric normalization, gap filling, and prediction of Landsat data. Remote Sens. Environ. 2008, 112, 3112–3130. [Google Scholar] [CrossRef]
  58. Holben, B.N. Characteristics of maximum-value composite images from temporal AVHRR data. Int. J. Remote Sens. 1986, 7, 1417–1434. [Google Scholar] [CrossRef]
  59. Huete, A.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213. [Google Scholar] [CrossRef]
  60. Kaufman, Y.J.; Remer, L.A. Detection of forests using mid-IR reflectance: An application for aerosol studies. IEEE Trans. Geosci. Remote Sens. 1994, 32, 672–683. [Google Scholar] [CrossRef]
  61. Pereira, J.M.C. A comparative evaluation of NOAA/AVHRR vegetation indexes for burned surface detection and mapping. IEEE Trans. Geosci. Remote Sens. 1999, 37, 217–226. [Google Scholar] [CrossRef]
  62. Zha, Y.; Gao, J.; Ni, S. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. Int. J. Remote Sens. 2003, 24, 583–594. [Google Scholar] [CrossRef]
  63. Li, X.; Li, D. Can night-time light images play a role in evaluating the Syrian crisis? Int. J. Remote Sens. 2014, 35, 6648–6661. [Google Scholar] [CrossRef]
  64. Li, X.; Zhang, R.; Huang, C.; Li, D. Detecting 2014 northern Iraq insurgency using night-time light imagery. Int. J. Remote Sens. 2015, 36, 3446–3458. [Google Scholar] [CrossRef]
  65. Zhang, Q.; Seto, K. Can night-time light data identify typologies of urbanization? A global assessment of successes and failures. Remote Sens. 2013, 5, 3476–3494. [Google Scholar] [CrossRef]
  66. Miller, S.D.; Lee, T.F.; Turk, F.J.; Kuciauskas, A.P.; Hawkins, J.D. Shedding New Light on Nocturnal Monitoring of The Environment with the VIIRS Day/Night Band. In Proceedings of the Atmospheric and Environmental Remote Sensing Data Processing and Utilization: Numerical Atmospheric Prediction and Environmental Monitoring, San Diego, California, USA, 31 July 2005.
  67. Roy, D.P.; Wulder, M.A.; Loveland, T.R.; Woodcock, C.E.; Allen, R.G.; Anderson, M.C.; Helder, D.; Irons, J.R.; Johnson, D.M.; Kennedy, R. Landsat-8: Science and product vision for terrestrial global change research. Remote Sens. Environ. 2014, 145, 154–172. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Zhang, Q.; Li, B.; Thau, D.; Moore, R. Building a Better Urban Picture: Combining Day and Night Remote Sensing Imagery. Remote Sens. 2015, 7, 11887-11913. https://doi.org/10.3390/rs70911887

AMA Style

Zhang Q, Li B, Thau D, Moore R. Building a Better Urban Picture: Combining Day and Night Remote Sensing Imagery. Remote Sensing. 2015; 7(9):11887-11913. https://doi.org/10.3390/rs70911887

Chicago/Turabian Style

Zhang, Qingling, Bin Li, David Thau, and Rebecca Moore. 2015. "Building a Better Urban Picture: Combining Day and Night Remote Sensing Imagery" Remote Sensing 7, no. 9: 11887-11913. https://doi.org/10.3390/rs70911887

APA Style

Zhang, Q., Li, B., Thau, D., & Moore, R. (2015). Building a Better Urban Picture: Combining Day and Night Remote Sensing Imagery. Remote Sensing, 7(9), 11887-11913. https://doi.org/10.3390/rs70911887

Article Metrics

Back to TopTop