Remote Sensing of Environment 205 (2018) 253–275

Contents lists available at ScienceDirect

Remote Sensing of Environment
journal homepage: www.elsevier.com/locate/rse

Using Landsat and nighttime lights for supervised pixel-based image
classiﬁcation of urban land coverT⁎

Ran Goldblatta, , Michelle F. Stuhlmacherb, Beth Tellmanb, Nicholas Clintonc, Gordon Hansona,
Matei Georgescub, Chuyuan Wangb, Fidel Serrano-Candelad, Amit K. Khandelwale,
Wan-Hwa Chengb, Robert C. Balling JrbaSchool of Global Policy and Strategy, UC San Diego, 9500 Gilman Drive La Jolla, CA 92093-0519, USA
School of Geographical Sciences and Urban Planning, Arizona State University, 976 S Forest Mall, Tempe, AZ 85281, USAc
Google, Inc., 1600 Amphitheatre Pkwy, Mtn. View, CA 94043, USAdLaboratorio Nacional de Ciencias de la Sostenibilidad, Universidad Nacional Autónoma de México, Apartado Postal 70-275 Ciudad Universitaria, UNAM 04510 México,
D.F. Circuito exterior s/n anexo al Jardín Botánico Exterior, MexicoeColumbia Business School, Columbia University, New York, NY 10027, USAb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Urbanization
Built-up land coverNighttime light
Image classiﬁcation
Google Earth Engine

Reliable representations of global urban extent remain limited, hindering scientiﬁc progress across a range of
disciplines that study functionality of sustainable cities. We present an eﬃcient and low-cost machine-learning
approach for pixel-based image classiﬁcation of built-up areas at a large geographic scale using Landsat data.
Our methodology combines nighttime-lights data and Landsat 8 and overcomes the lack of extensive groundreference data. We demonstrate the eﬀectiveness of our methodology, which is implemented in Google EarthEngine, through the development of accurate 30 m resolution maps that characterize built-up land cover in threegeographically diverse countries: India, Mexico, and the US. Our approach highlights the usefulness of data
fusion techniques for studying the built environment and is a ﬁrst step towards the creation of an accurate
global-scale map of urban land cover over time.

1. Introduction

urbanization in a consistent manner at high resolution and global scale(Ban et al., 2015; Giri et al., 2013).
The revolution in geospatial data has transformed how we study
cities. Previous approaches leveraged household surveys but these are
expensive to collect, produced infrequently, and subject to measurement problems. Since the 1970s, however, terrestrial Earth-observation
data have been continuously collected in various spectral, spatial andtemporal resolutions. As improved satellite imagery becomes available,
new remote-sensing methods and machine-learning approaches havebeen developed to convert terrestrial Earth-observation data intomeaningful information about the nature and pace of change of urban
landscapes and human settlements (Ban et al., 2015; Chen et al., 2015;
CIESIN, 2005; Gaughan et al., 2013; Pesaresi et al., 2016; Potere et al.,2009; Seto et al., 2011; Taubenböck et al., 2012).
The availability of satellite data has triggered the development of
new methods to map global land cover using remotely-sensed data such
as Landsat (Chen et al., 2015; Gaughan et al., 2013; Goldblatt et al.,2016; Patel et al., 2015), MODIS (Moderate Resolution Imaging Spectroradiometer) (Schneider et al., 2009, 2010; Wan et al., 2015), DMSP-

Urbanization has been a fundamental trend of the past two centuries
and a key force shaping the development of the modern world. Between
1950 and 2014, the share of the global population living in urban areas
increased from 30% to 54%, and in the next few decades is projected to
expand by an additional 2.5 billion urban dwellers, primarily in Asia
and Africa (Seto et al., 2011; UN, 2014). Urban population growth isaccompanied by a dramatic increase in the land area incorporated incities (Seto et al., 2011). While urbanization in rapidly growing nations
is helping lift hundreds of millions of people out of poverty, it is also
creating immense societal challenges by increasing greenhouse-gas
emissions, destabilizing fragile ecosystems, and creating new demands
on public services and infrastructure that impose signiﬁcant burdens on
the environment (Ban et al., 2015). Timely and reliable information on
the extent of urban areas is fundamental for the support of sustainableurban development and management (Ban et al., 2015; Jacob and Ban,2015). Despite the importance of understanding the drivers of urban
growth, we are still unable to quantify the magnitude and pace of⁎Corresponding author.
E-mail addresses: rgoldblatt@ucsd.edu, ranrang@gmail.com (R. Goldblatt).

https://doi.org/10.1016/j.rse.2017.11.026
Received 20 February 2017; Received in revised form 12 November 2017; Accepted 29 November 20170034-4257/ © 2017 Elsevier Inc. All rights reserved.

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Moreover, mapping land cover at a global scale and with high precision
requires eﬀective, eﬃcient and operational approaches to deal with avery large volume of data. For example, it is estimated that over 10,000
Landsat satellite images are required to cover the entire Earth at 30 m
resolution (Chen et al., 2015). Until recently, the majority of studies
that analyze urbanization have been limited in their geographic scalebecause of the lack of extensive high-resolution satellite data, scarcity
of ground-reference data, and computational constraints. Emerging
cloud-based computational platforms now allow for scaling analysis
across space and time. Google Earth Engine (GEE) is one platform thatleverages cloud-computing services to achieve planetary-scale utility.
GEE has been previously used to map population (Patel et al., 2015;
Trianni et al., 2015), urban areas (Goldblatt et al., 2016), and surface
water (Pekel et al., 2016). This paper contributes to this literature bydeveloping a machine-learning methodology for supervised high-resolution image classiﬁcation of built-up areas using GEE's cloud-based
computational platform.

OLS (Elvidge et al., 2014; Liu et al., 2012; Xiao et al., 2014; Zhang and
Seto, 2013) and other spaceborne High-Resolution (HR), Very-HighResolution (VHR) and Synthetic-aperture radar (SAR) radar sensors(Ban et al., 2015; Gamba et al., 2011; Jacob and Ban, 2015). Recent
studies have developed automated and semi-automated classiﬁcation
procedures to map global land cover at a 30 m resolution with highaccuracy (Ban et al., 2015; Chen et al., 2015). Because Landsat satellites
have been collecting data from Earth since 1972, Landsat data are often
used for analysis of urban change (Patel et al., 2015), and are ideal forland cover mapping (Woodcock et al., 2008). Nighttime light data arealso associated with developed land (Elvidge et al., 2014; Levin andDuke, 2012; Sutton, 2003) and can be used to infer the extent of urban
areas (Bagan and Yamagata, 2015; Small and Elvidge, 2013; Zhang and
Seto, 2013), as well as economic activity at the local, regional andnational levels (Elvidge et al., 2014; Henderson et al., 2003; Keola et al.,2015). Sensors on board the Operational Line-scan System of the Defense Meteorological Satellite Program (DMSP-OLS) have captured artiﬁcial lighting since the early 1990's. A pixel's nighttime light valuethat exceeds a speciﬁed threshold, which may vary across regions orcountries, signiﬁes urban development (Henderson et al., 2003; Liuet al., 2016; Small and Elvidge, 2013; Su et al., 2015; Wei et al., 2014;
Zhou et al., 2014, 2015). However, inference using nighttime-light data
are often inaccurate, particularly in low-density urban areas (Zhang and
Seto, 2013). DMSP-OLS can also exaggerate the extent of urban areas(Henderson et al., 2003; Small et al., 2005), while overlooking small ordeveloping settlements. In addition, the extent and intensity of lit areas
cannot directly delimit urban regions due to the “blooming” eﬀect(Imhoﬀ et al., 1997) and “saturation” of pixels (Hsu et al., 2015).
Blooming refers to the identiﬁcation of lit areas as consistently larger
than the settlements with which they are associated (Small et al., 2005);
saturation occurs when pixels in bright areas, such as in city centers,
reach the highest possible digital number (DN) value (i.e., 63) and nofurther details can be recognized (Hsu et al., 2015).
Until recently, most remote sensing studies focused on local settings(Herold, 2009). Mapping land cover at a national or regional scale ischallenging because of the lack of high-resolution global imagery, theheterogeneous and complex spectral characteristics of land, and thesmall and fragmented spatial conﬁguration of many cities (Chen et al.,2015; Herold, 2009). In the case of mapping urbanization, existing
maps of urban land show considerable disagreement on the locationand extent of urbanization (Potere et al., 2009; Seto et al., 2011) andare limited across space and time. These inconsistencies arise in part
because the delineation of urban land depends on the input data
(Schneider et al., 2010), which may capture diﬀerent dimensions ofurbanization, such as built-up land cover or land use and population
density (Bagan and Yamagata, 2014; Stevens et al., 2015; Tatem et al.,2007).

1.2. Research objective and contribution
The use of nighttime remotely-sensed data to map urbanization is
not new to the literature. Remotely-sensed data on artiﬁcial lighting haslong been considered an economical way to map urbanization anddevelopment across the globe (Elvidge et al., 2009). By utilizing the
distribution of vegetation land cover, the combination of nighttime anddaytime data increases the heterogeneity of urban and suburban land
cover (e.g., distinguishing between built-up land cover and vegetationin urban areas) and improves the characterization of inter-urban
variability in nighttime luminosity (Zhang et al., 2013). This, in turn,improves the ability to detect urban features (Lu et al., 2008; Ma et al.,2014) including sub-pixel fractional urban land cover (Huang et al.,2016). Several spectral indices that combine nighttime light and vegetation spectral characteristics have been developed, including the
Vegetation Adjusted NTL Urban Index (VANUI) (Jing et al., 2015;
Zhang et al., 2013), the Normalized Diﬀerence Urban Index (NDUI)(Zhang et al., 2015) and the Normalized Diﬀerence Spectral Vector(NDSV) (Trianni et al., 2015). These indices increase the separability
between urban and non-urban land cover.
We develop a methodology that combines nighttime and daytime
remotely-sensed data. We collect training examples automatically using
DMSP-OLS data, and use them for classiﬁcation of built-up areas withdaytime Landsat 30 m spatial resolution imagery. Previous studies thatcombine nighttime and daytime data have either been limited in their
spatial application (i.e., the ability to generalize the method and to
apply it over regions with heterogeneous land cover) or spatial resolution (i.e., many of the existing approaches, for example those that
rely on MODIS or DMSP-OLS, are limited in their spatial resolution). In
this study, we adopt a hexagonal tessellation mapping approach tohandle large variation across regions (where we refer to each hexagon
in the hexagonal grid as a hex-cell). We collect training examples fromeach hex-cell and classify the hex-cell as an independent unit of analysis.
Our methodology can be applied across heterogeneous land coverand across time and, crucially, does not rely on expensive hand-labeled
examples. It requires minimal manual adjustments for training andclassiﬁcation, and does not require adjustments to local parameters.
This feature makes the approach scalable across space and time.
Importantly, the methodology is time-invariant and can be applied
whenever Landsat and DMSP-OLS data are coincident. The robustnessof the methodology lies in our approach to sample training examples(i.e., according to the relative intensity of the emitted light at night andthe distribution of vegetation land cover) and on the per-hex-cellclassiﬁcation, which allows us to account for regional variations in theland cover. Finally, we assess the accuracy of the methodology using an
extensive dataset of 84,564 hand-labeled polygons characterizing builtup (BU) and not built-up (NBU) pixels for each of the three study areas,

1.1. Detecting urbanization using machine learning
Urban areas can be detected in satellite imagery using various machine-learning approaches (e.g., supervised, unsupervised and semisupervised). These approaches typically rely on reference data thatmark urban features, either for training or validation. Reference data
are fundamental not only for mapping urbanization across space, butalso for classiﬁcation over time (Boucher and Seto, 2009). Some of thereference datasets used for classiﬁcation include Landsat-based urban
maps (Potere et al., 2009), census-based population databases (Stevens
et al., 2015), hand-labeled examples (Goldblatt et al., 2016), and datacollected via crowd-source platforms, such as OpenStreetMap (OSM)(Belgiu and Drǎguţ, 2014; Estima and Painho, 2015). However, because
they are expensive to collect, reference datasets for large geographic
scales are scarce (Miyazaki et al., 2011). Due to the scarcity of groundreference data, it is often necessary to exploit existing global coarse
datasets and classiﬁcation products to create accurate higher-resolution
maps of urban areas (Kasimu et al., 2009; Trianni et al., 2015).
254

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 1. Annual changes in share of urban populationin India, US and Mexico compared to world average(UN, 2014).

100.0

Urban population (% of total)

90.0
80.0
70.0
60.0

50.0
40.0
30.0
20.0

10.0

India

US

2050

2046

2038

2042

2034

2030

2026

2022

2014

Mexico

2018

2006

2010

2002

1998

1994

1986

1990

1982

1978

1974

1966

1970

1962

1954

1958

1950

0

World

2. Conceptual framework: Infused DMSP-OLS/Landsat
methodology

which are described below.
1.3. Study areas

Our methodology combines Earth-observation datasets from two
domains: DMSP-OLS and Landsat. Using a subset of the data, we
identify the optimal parameters to the classiﬁers and then classify thebuilt-up land cover in the entire area of each country. The classiﬁcation
procedure involves ﬁve steps (Fig. 2): (1) Divide each country into a
uniform hexagonal grid; (2) Pre-process Landsat 8 images; (3) Extractlabeled examples from DMSP-OLS; (4) Perform supervised image classiﬁcation; and (5) Validate and test. We next describe these steps indetail.

To illustrate our methodology and its applicability in heterogeneousand diverse geographic conditions, we map built-up land cover in three
countries that are characterized by distinct geographic conditions (e.g.,land cover, topography, climate, soil, landform, and fauna): India,
Mexico, and the United States (US).
1.3.1. India
The share of India's population living in urban areas in 2015 was
33%, which is much lower than the corresponding values of 79% and
82% for Mexico and the US, respectively. India is, however, urbanizing
at a relatively rapid rate (Fig. 1). For example, between 2010 and 2015,
India's average annual rate of change of the urban population was
1.14%, compared to 0.36% and 0.21% in Mexico and in the US, respectively. By 2050, half of India's population is likely to be urban.

2.1. Divide each country into a uniform hexagonal grid (“mapping zones”)
Mapping heterogeneous land cover at a large geographic scale requires partitioning the region of interest into a ﬁnite number of relatively homogenous sub-regions that are characterized by similar landform, soil, vegetation, spectral reﬂectance, and image footprints(Homer et al., 2004). This practice is often referred to as ‘zone mapping’(Homer and Gallant, 2001). The partition can be according to diﬀerent
criteria, such as land cover and land use, socio-political deﬁnition, or
size (Hunsaker et al., 1994; O'Neill et al., 1996; Turner, 1989), or bymeans of an artiﬁcial grid system where each element in the grid istreated as an independent region of interest. Examples of such grids inpractice include the Global Grid (GG) system (Theobald, 2016) and the
ISEA DGGs of Gong et al. (2013) for global land cover mapping. There
are many advantages for using an artiﬁcial grid system over administrative division (e.g., census areas or political boundaries) or any otherthematic division (e.g., ecological or climatic regions). First, administrative divisions are often delineated according to criteria that are unrelated to physical characteristics (e.g., the boundary between Washington, D.C., and Maryland). Second, the geographic scale and thesize of the administrative division varies within and between countries(e.g., the division into census tracts and blocks). Unlike a uniform
global artiﬁcial grid system, an administrative division does not allow
for a multi-geographic-scale classiﬁcation approach using recursivelydivided units of analysis.
Here we partition each of our three countries of interest into an
equal-area hexagonal grid (or a hexagonal tessellation). Hex-cells in ahexagonal tessellation are arranged in a contiguous global lattice,which is the most compact arrangement of many equal circles (Birchet al., 2007). Hexagonal grids are advantageous because they arecharacterized by elements that do not have gaps or overlaps and the1.3.2. Mexico
Mexico has gone through three major phases of urbanization. In theﬁrst phase, 1900–1940, urban growth slowly incorporated 10% to 20%of the population. In the second phase, 1940–1980, rapid urban expansion, particularly in Mexico City, increased the share of urban population to 55%. Since the 1980s, more dispersed moderate urban
growth increased the urban share of the population to over 70%(Consejo Nacional de Población, 2012). The current phase of urbanization is characterized by informal urbanization on the city periphery,
representing 65% of all new housing construction in Mexico City, and iseven higher for small to mid-sized cities (Connolly, 2014). Informal
settlements tend to be marginalized in terms of lower socio-economic
development, access to services like water and electricity, and are more
vulnerable to risks like water scarcity and ﬂooding (Aguilar, 2008;
Aguilar and Guerrero, 2013; Consejo Nacional de Población, 2012;
Eakin et al., 2016).
1.3.3. US
In the context of the US, the initial urban growth occurred from1790 to 1890 and the country has become increasingly urban since (USCensus Bureau, 2012). In 1910, the Census Bureau deﬁned an urban
area as one with a population above 2500, and the 1920 census marked
the ﬁrst time that 50% of the US population lived in an urban area (USCensus Bureau, 2016).
255

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 2. Schematic illustration of the methodology.

2.2. Pre-process Landsat 8 images (the classiﬁer's inputs)

center-to-center distances between adjacent grid cells are approximately equal. Moreover, hex-cells have a topology that is symmetrical,
invariant, and of equal area, and can be recursively partitioned into
smaller divisions of grids (Richards et al., 2000) (Fig. 3, right). Since thesize of each hex-cell in the grid may aﬀect the accuracy of the classiﬁcation, we use a hexagonal tessellation that in successive classiﬁcations treats each pixel as belonging to a hex-cell of a diﬀerent size. We
examine the following sizes of hex-cells (diﬀerent center-to-center
distances): a distance of 1 decimal degree (1-degree), 4 decimal degrees(4-degrees) and 8 decimal degrees (8-degrees) from center to center(see Fig. 3 for illustration).

We use Landsat 8 data as classiﬁer inputs (predictors). We apply astandard Top-of-Atmosphere (TOA) calibration on all USGS Landsat 8
Raw Scenes in one year. Since processed DMSP-OLS is only available up
through 2013, our analysis begins that year. We assign a “cloud score”
to each pixel and select the lowest possible range of cloud scores. The
cloud score is a simple cloud-likelihood score, intended to compare
multiple looks at the same point for ‘relative’ cloud likelihood (acrosstime). The score, ranging between 0 and 100, is determined according
to the brightness of a pixel, its temperature, and a computed per-pixel
Normalized Diﬀerence Snow Index (NDSI). Similar to Trianni et al.
(2015), we retrieve only those scenes with less than 10% cloud256Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 3. The three examined hexagon levels: 1, 4, and 8° from center to center, for (a) Mexico, (b) India, and (c) the US (left); The cells in the grid can be recursively partitioned into smaller
divisions creating a hierarchical grid; (d) each pixel is classiﬁed multiple times, each time as part of a diﬀerent hexagonal scale.

2.3. Extract “built-up” and “not built-up” labeled examples from DMSPOLS

coverage.
We compute per-band median values from the accepted pixels. Perpixel median values are insensitive to extreme values, such as very
bright or dark pixels, in a stack of scenes (Flood, 2013). The decimalvalue of each pixel is scaled to an 8-bit range (0–255) for display
purposes. We also calculate six additional per-pixel spectral indices andadd them as additional inputs to the classiﬁer. These indices are commonly used to identify bodies of water (Normalized Diﬀerence Water
Index, NDWI), built-up areas (Normalized Diﬀerence Built-up Index,
NDBI; Urban Index, UI, NDSV), and vegetation (Normalized Diﬀerence
Vegetation Index, NDVI; Enhanced Vegetation Index, EVI). NDWI
(McFeeters, 1996) is commonly used to identify open water or tomeasure vegetation water content. This index and subsequent versionsof it (e.g., the Modiﬁcation of Normalized Diﬀerence Water Index(NDWI)) (Xu, 2006) have been shown to correlate negatively withurban areas and impervious surface (Chen et al., 2006; Xu, 2008). NDVI(Tucker, 1979) and its variant, EVI (Huete et al., 2002), produce
spectral indices that tend to negatively correlate with urban areas.
Previous studies have used these indices to identify urban areas and tomeasure urban characteristics such as heat islands, changes in net primary productivity, and urban form (Chen et al., 2006; Clinton andGong, 2013; Estoque and Murayama, 2015; Li et al., 2013; Myint et al.,2011; Xu, 2008; Yuan and Bauer, 2007; Zhao et al., 2016). UI(Kawamura et al., 1996) and NDBI (Zha et al., 2003) are both used to
identify built up areas and correlate positively with urban land coverand impervious surface (Chen et al., 2006; Estoque and Murayama,
2015; Li et al., 2013). NDSV (Normalized Diﬀerence Spectral Vector)(Angiuli and Trianni, 2014; Trianni et al., 2015) is used to identify
consistent human settlement and artiﬁcial materials. This index, whichis composed out of all possible normalized indexes combining two
bands with the same spatial resolution, has been shown to produce
robust and globally consistent measures of urban land cover (Angiuli
and Trianni, 2014). We describe these indices in detail in Appendix A.

DMSP-OLS sensors detect visible and near-infrared (VNIR) emissions at night. Because highly lit pixels are associated with man-madestructures that emit light, we assume that pixels with DN values thatexceed a threshold represent areas with built-up land cover and manmade structures. In this study we use DMSP-OLS (the “stable light”
band of the ‘F182013’ satellite) to identify highly lit pixels. The “stable
light” band discards ephemeral events, such as ﬁres, but DMSP-OLS
sensors are sensitive to persistent lighting, including from othersources, such as gas ﬂares or volcanoes. Elvidge et al. (2009)'s Global
Gas Flaring Estimates showed that while in India and Mexico there areno gas ﬂares on land, the US has several gas ﬂares on land, primarily inremote locations in Texas and Wyoming. We use their data to remove
pixels that were classiﬁed as “built-up” from these locations during
post-processing.
We sample “built-up” examples from highly lit DMSP-OLS pixels.
We characterize a pixel as highly lit if its DN value exceeds a given
threshold. We deﬁne this threshold as a percentile value (the 50th,
75th, 90th, 95th or the 99th percentile) of all DMSP-OLS pixels within ahex-cell. Additionally, we sample “not-built-up” examples from low-lit
pixels. We characterize a DMSP-OLS pixel as low-lit if its DN value islower than a given threshold (which we deﬁne as the 10th, 25th, 50th
or the 75th percentile of all DMSP-OLS pixels within a hex-cell). We do
not sample examples from DMSP-OLS pixels that are located in theintermediate buﬀer between highly-lit and low-lit pixels(Fig. 4).
Due to the spatial resolution of DMSP-OLS and the blooming eﬀect,areas identiﬁed as highly lit may potentially include non-built-up land
cover. We therefore exclude these types of land cover from lit pixels
based on Landsat's per-pixel NDVI and NDWI values. We characterize apixel as vegetation if its NDVI value is higher than a given threshold.
Due to the heterogeneity of the land cover, we determine this thresholdseparately for each hex-cell. We calculate the 50th, 75th, 90th, 95th
and the 99th percentile NDVI value of all pixels within the hex-cell andremove from the DMSP-OLS highly-lit pixels all Landsat pixels whose257

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 4. The method to identify highly-lit and low-lit pixels for our classiﬁcation training data. A pixel is characterized as highly-lit and as lowly-lit if its DN value exceeds and is below agiven threshold, respectively. The threshold is deﬁned as a percentile value of all DMSP-OLS pixels in the hex-cell. This example presents three 1-degree hex-cells in India. The thresholds
to characterize highly-lit pixels are (a) 27, (b) 30 and (c) 17; The threshold to characterize low-lit pixels is 9. Note, that we do not sample examples from the intermediate buﬀer between
highly-lit and low-lit pixels (marked in light blue). (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)scales and when the data are noisy and high dimensional (Gislason
et al., 2006; Guan et al., 2013). Random Forests are computationallyless intensive than other tree ensemble methods (Jean et al., 2016;Rodriguez-Galiano et al., 2012) and can eﬀectively incorporate many
covariates with minimum tuning and supervision (Stevens et al., 2015).
We set the number of trees in the Random Forest classiﬁer to 20. Previous studies show mixed results as for the optimal number of trees inthe decision tree. The number ranges between 10 trees (Zhang et al.,2012) and 100 trees (Rodriguez-Galiano et al., 2012). Goldblatt et al.
(2016) show that although the performance of Random Forest improvesas the number of trees increases, this pattern holds only up to 10 trees;
the performance remains stable with up to 50 or even 100 decision
trees.
We map each hex-cell as an independent unit of analysis. We only
sample and classify hex-cells that include at least one DMSP-OLS lit
pixel (i.e., DN value above 0). These hex-cells, however, (i.e., the 1degree and 4-degrees hexagonal grids) may also include isolated ruraland suburban settlements that do not emit much light, which we stillwish to capture. To ensure that they are not omitted, pixels in the excluded hex-cells are mapped according to the classiﬁer trained withlarger hex-cells (i.e., the 8-degree hexagonal grid). We clip the classiﬁcation at the 8-degrees hex-cells to the extent of the smaller hex-cells
that were excluded from the analysis.

NDVI value exceeds this threshold. Additionally, we remove bodies ofwater from DMSP-OLS highly-lit pixels. We characterize a pixel as bodyof water if its NDWI value is positive (see Gao, 1996).
Finally, we randomly sample 100,000 pixels in each hex-cell. We
create a point at the center of each Landsat pixel and associate eachpoint with the spectral values of the Landsat composite and derived
spectral indices. Each example includes a label (i.e., highly-lit or not)and the spectral values from Landsat. These training examples are used
to build a classiﬁcation model for each hex-cell, independently.
2.4. Perform supervised image classiﬁcation
We train and classify each hex-cell using a Random Forest classiﬁer.
Random Forests are tree-based classiﬁers that include k decision trees
and p randomly chosen predictors for each recursion. When classifying
an example, its variables are run through each of the k tree predictors,
and the k predictions are merged by voting on the most popular class.
The learning process of the forest involves some level of randomness.
Each tree is trained over a random sample of examples from the training
set and each node's binary question in a tree is selected from a random
subset of p input variables. We use Random Forest because previous
studies ﬁnd that its performance is superior to other classiﬁers(Goldblatt et al., 2016), especially when applied to large geographic
258Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

et al. (2009) and Schneider et al. (2010), propose a similar deﬁnition
for urban. For India, we use Goldblatt et al.'s (2016) ground-reference
dataset, which includes 21,030 examples (30 m by 30 m in size) labeled
as “built-up” or as “not built-up” that span the entire country (4682
polygons labeled as built-up and 16,348 labeled as not built-up). The
polygons are manually labeled by two graduate students (each polygon
is labeled by one student) who were provided with extensive training
and supervised by the researchers. Each student was provided with anequal proportion of samples. The students labeled each polygon either
as “built-up” or as “not built-up” based on a visual interpretation of themost recent available high-resolution satellite image in Google Earth(typically from 2014 to 2015).
For the US and Mexico, we construct a manually labeled groundreference dataset of 20,000 examples (polygons, 30 m by 30 m in size)for each country. To ensure representative “built-up” and “not built-up”
examples, we construct examples using a random stratiﬁed sampling
procedure based on the intensity of nighttime light (the stable lightband of DMSP-OLS). In each country, we identify the DMSP-OLS pixels
whose value is 63 (the highest possible value). Then, we calculate a 5pixel radius circular kernel (a buﬀer of approximately 5 km). In eachcountry, we select a random sample of points: 15,000 from the DMSPOLS pixels with a DN value 63 and 5000 from the 5 km buﬀer zone. A
30 × 30 m polygon is drawn around each sample point. We overlay
these polygons onto high-resolution imagery in Google Earth. The US
and Mexico examples are labeled by six undergraduate students (three
students per country), each polygon is labeled either by one (in Mexico)
or by three students (in the US). The students were provided with extensive training and were supervised by the researchers. The students
were instructed to ﬂag ambiguous polygons where the deﬁnition of thepolygon as built-up or not built up is not absolutely clear; we exclude
these polygons from the dataset. We randomly assign half of the dataset
in each country as the test set (to determine optimal parameters to theclassiﬁer) and half as the validation set (to validate the ﬁnal classiﬁedmap). Between 22% and 27% of the polygons are labeled as built-up(Table 1).
Previous studies (Liu and Zhou, 2004; Rwanga and Ndambuki,
2017; Schneider et al., 2009; Zhang and Seto, 2011) rely on hand-

The output of the classiﬁcation is a per-pixel posterior probability
value (i.e., the probability that a pixel is built-up). We characterize apixel as “built-up” if the probability it is “built-up” is higher than agiven threshold. Previous studies determine this threshold based on
visual assessment (Schneider et al., 2010) or according to the natural
breaks in the data values (Mertes et al., 2015). Here, we determine thisthreshold for each hex-cell using the Otsu algorithm, a clustering-based
image thresholding method (Otsu, 1979). Otsu is a nonparametric andunsupervised method for automatic threshold selection that was originally developed for picture segmentation. The method uses a discriminant criterion to identify an optimal threshold that maximizes the
between-class variance. Speciﬁcally here, the threshold is chosen to
maximize the variance between “built-up” and “not-built-up” classes:2∑ (pk − p )2

(1)

k=1

where the threshold partitions the pixels into classes k = 1 (built-up)
and k = 2 (not built-up). In the remote-sensing domain, Otsu has beenused for classiﬁcation of mangrove forests (Chen et al., 2013) and urban
areas (Sirmacek and Unsalan, 2010).
Finally, we tile the classiﬁcations of the hex-cells to the extent ofeach country and post process the classiﬁcation maps. We clip them tothe extent of the countries' borders and remove misclassiﬁed built-up
pixels from gas ﬂaring areas using the Elvidge et al. (2009) database.
2.5. Validate and test
In each country, we assess the accuracy of the classiﬁcation using
hand-labeled examples. We manually label these examples (polygons,
30 m by 30 m in size) as “built-up” or as “not built-up”. We deﬁne apolygon (30 m by 30 m in size) as “built-up” if the majority of its area
(50% or greater) is paved or covered by human-made surfaces and used
for residential, industrial, commercial, institutional, transportation, orother non-agricultural purposes. Otherwise, it is deﬁned as “not builtup”. Fig. 5 presents examples of “built-up” and “not built-up” labeledpolygons. Previous studies, including Goldblatt et al. (2016), Potere

a

b

c

d

e

f

g

h

Fig. 5. Ground-reference labeled examples. Each example (a polygon 30 m × 30 m in size) is labeled as “not built-up” (top) and as “built-up” (bottom). To illustrate, the examples inFig. 5a, b and d are labeled as “not built-up” (the majority of the polygon in Fig. 5a is an agriculture ﬁeld, the majority of the polygon in Fig. 5b is a backyard covered with grass, and the
majority of the polygon in Fig. 5d is covered with grass). The examples in Fig. 5e and g are labeled as “built-up” (the majority of the polygon in Fig. 5e is a residential building and a paved
road; the majority of the polygon in Fig. 5g is a paved parking lot).

259

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

and below a constant percentile value (the 99th percentile) of all
DMSP-OLS pixels in a hexagon, respectively. We remove fromhighly-lit pixels bodies of water (i.e. Landsat pixels with a positive
NDWI value) and vegetation. We deﬁne a Landsat pixel as vegetation if its NDVI value exceeds either the 50th, 75th, 90th, 95th or the
99th percentile of the NDVI value of all Landsat pixels in a hex-cell.
We perform classiﬁcation with Landsat spectral bands, together with
ﬁve spectral indices (NDVI, NDWI, NDBI, EVI, UI) as inputs to theclassiﬁer. Based on these results, we determine the optimal percentile value to characterize a Landsat pixel as vegetation.
b. Classify highly-lit pixels: We deﬁne a DMSP-OLS pixel as highly-lit
if its DN value exceeds either the 50th, 75th, 90th, 95th or the 99th
percentile value of all DMSP-OLS pixels within a hex-cell (below thisvalue, the pixel is deﬁned as low-lit). We remove bodies of water(i.e., a positive NDWI value) and vegetation (according to the optimal NDVI percentile value found in step (a)) and perform classiﬁcation with Landsat's band values together with ﬁve spectral indices (NDVI, NDWI, NDBI, EVI, UI) as input (we do not include
NDSV here). These results indicate the optimal percentile thresholdto characterize a DMSP-OLS pixel as highly lit. Pixels below thethreshold are classiﬁed as lowly lit.
c. Determine inputs to the classiﬁer: We deﬁne a DMSP-OLS pixel as
highly-lit according to procedure in (b) and remove bodies of water(i.e., a positive NDWI value) and vegetation (according to the deﬁnition of a Landsat's vegetation pixel from step (a)). We assess theperformance of the classiﬁers with ﬁve combinations of inputs: (1)only Landsat band values (without additional spectral indices);
Landsat bands together with (2) NDVI; (3) NDBI; (4) NDVI, NDWI,NDBI, EVI, UI; and (5) NDSV. The result of this step indicate theoptimal inputs to the classiﬁer.
d. Classify lowly-lit pixels: In the procedure above, we deﬁne a pixel
as lowly lit if it is not highly lit. However, this deﬁnition may resultin mislabeled pixels as “not built-up”, especially in proximity tohighly-lit pixels in the core of settlements. In the ﬁnal step of theparameter evaluation, we sample “not built-up” examples only from
DMSP-OLS pixels whose DN value is below either the 10th, 25th,
50th and 75th percentile value of all DMSP-OLS pixels in a hex-cell.
We perform classiﬁcation using the optimal parameters found insteps (a) through (c). Based on these results, we determine thethreshold to identify a DMSP-OLS pixel as lowly-lit.

Table 1the distribution of the built-up and not built-up points for each country.
BU points

NBU points

Total

Hand-labeled polygons
India
4682 (22.26%)
US4386 (21.6%)
Mexico5477 (27.4%)

16,348 (77.74%)
15,898 (78.4%)
14,523 (72.6%)

21,03020,28420,000Tomnod-labeled polygons
India
2276 (23%)
US
3305 (30.7%)
Mexico829 (31.8%)7602 (77%)
7456 (69.3%)
1782 (68.2%)

9878
10,761
2611
84,564labeled examples for accuracy assessment of classiﬁcation maps.
However, a key challenge in any machine learning task is “class noise”
which may be caused by a subjective interpretation of the target label(Foody, 2002). To assure the validity of our hand-labeled examples, we
add an additional source of reference data to the ﬁnal accuracy of ourclassiﬁcation. We utilize DigitalGlobe's Tomnod platform to label
23,250 examples (referred to as Tomnod examples). The examples arelabeled as “built-up” or as “not built-up” according to Tomnod's
CrowdRank algorithm (Wood and Har-Noy, 2013), a measure designed
to estimate crowd consensus and to assure that a suﬃcient number oflabeling individuals reach a consensus for each labeled examples. The
rank is also used to verify the accuracy and reliability of individual
contributors. Together with the hand-labeled examples, the dataset we
use for test and for validation consists of 84,564 labeled examples intotal (Table 1).
We use half of the hand-labeled dataset (the test set) to assess alternative parameters for the classiﬁers, and the other half (the validation set), in addition to the Tomnod examples, to evaluate and report its
performance and the accuracy of the ﬁnal classiﬁcation map. We estimate the performance of the classiﬁers with several performance estimators. We refer to the class “built-up” (BU) as positive and to the class“not built-up” (NBU) as negative: (1) Producer's accuracy, referred to as
True-Positive Rate (TPR), the percentage of actual BU examples classiﬁed correctly as BU; (2) True-Negative Rate (TNR), the percentage ofactual NBU examples classiﬁed correctly as NBU; (3) Balanced
Accuracy, the average of TPR and TNR; and (4) User's accuracy, thepercentage of correctly classiﬁed BU examples of all examples predicted
as BU.

Producer′s accuracy (TPR) = TP/(TP + FN)(2)

TNR = TN/(TN + FP)(3)

Balanced Accuracy = (TPR + TNR)/2(4)

User′s accuracy = TP/(TP + FP)(5)

3. Results
3.1. Optimal parameters to classiﬁer
As explained above, we ﬁrst select in each country 50 1-degree hexcells in order to identify the optimal parameters to the classiﬁer. We
classify these hex-cells and assess the accuracy of the classiﬁcation withthe test set (the results of the accuracy assessment with each of theexamined parameter are presented in Appendix B).
The results show variations between countries in the optimal
parameters to the classiﬁers. First, we ﬁnd that the countries diﬀer inthe threshold to characterize a DMSP-OLS pixel as highly lit. In the US
and in Mexico the optimal threshold is the 90th percentile value of all
DMSP-OLS pixels in a hex-cell, while in India it is the 95th percentile.
We ﬁnd that the optimal threshold to deﬁne a pixel as lowly lit is the75th percentile value in the three countries. Additionally, the countries
diﬀer in the optimal threshold to deﬁne a Landsat pixel as vegetation. In
India and the US, the best performance is achieved when Landsat pixels
are deﬁned as vegetation if their NDVI value exceeds the 90th percentile of all pixels in a hex-cell, while in Mexico this threshold is lower(the 75th percentile). We remove vegetation Landsat pixels fromhighly-lit DMSP-OLS pixels. Moreover, the results show variations between countries in the optimal inputs to the classiﬁers. Classiﬁcation
with NDVI, NDWI, NDBI, EVI, UI as inputs to the classiﬁer results in thebest performance (determined by the highest balanced accuracy) for thewhere TP is the number of the actual BU examples predicted to be BU;
TN is the NBU examples predicted as NBU; FN is the actual BU examples
predicted as NBU and FP is the actual NBU examples predicted as BU.
2.5.1. Optimal parameters to the classiﬁer
We use half of the hand-labeled dataset (the test set) to classify 50 1degree hex-cells per country and to assess alternative parameters for theclassiﬁers. We select hex-cells for which we have the largest number ofhand labeled examples (at least 35 examples). This accuracy assessment
procedure is designed to identify (1) the optimal percentile values tocharacterize a DMSP-OLS pixel as highly-lit and as lowly-lit; (2) the
optimal percentile value to characterize a Landsat pixel as containing
vegetation (according to NDVI); and (3) the optimal input features tothe classiﬁer. We proceed in the following steps, also show in Fig. 6:
a. Deﬁne a Landsat pixel as vegetation: First, we deﬁne highly-lit
and low-lit DMSP-OLS pixels as pixels whose DN values are above
260

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 6. The procedure to identify the optimal parameters to theclassiﬁer. We select in each country 50 hexagon (1-decimal-degreefrom center to center) and identify (1) the optimal percentile value
to characterize a DMSP-OLS pixel as highly-lit and as lowly-lit; (2)the optimal percentile value to characterize a Landsat pixel as
containing vegetation (according to NDVI); and (3) the optimal
input features to the classiﬁer. Each experiment (B-D) relies on the
optimal parameters found in the previous experiments. (* The ﬁve
indices are NDVI, NDBI, UI, NDWI, EVI).

a

b

c

d

3.2. The Otsu method

Table 2The optimal parameters to the classiﬁers, per country (India, US and Mexico): percentile
NDVI value to characterize a Landsat pixel as vegetation, percentile DN value to characterize a DMSP-OLS pixel as highly-lit or low-lit and inputs to the classiﬁer.
NDVI threshold(percentiles)

Highly-lit pixels(percentiles)

Inputs toclassiﬁer

Low-lit pixels(percentiles)

India
US

90th
90th

95th
90th

75th
75th

Mexico

75th

90th

NDSV
NDVI, NDWI,NDBI, EVI,UINDVI, NDWI,NDBI, EVI,
UI

Random forests predict for each new example the probability it is apositive example (in this case, a built-up pixel). This is a posterior
probability that ranges between 0 and 1. A pixel is deﬁned as built-up if
the probability it is built-up exceeds a given threshold, which we
identify for each hex-cell using the Otsu method. The Otsu method has
shown to provide satisfactory results for thresholding an image when
the histogram of the pixels values is characterized by a non-normal
distribution, especially bimodal (Ng, 2006).
As clearly illustrated in Fig. 7 (the distribution of the posterior
probability values of all pixels within 12 hexagons), the distribution of
the posterior probability values is symmetric (U shaped), with two highpeaks in each histogram (hex-cell) that are separated by deep valleys.
These conditions allow for an eﬀective performance of the Otsu
threshold methods in separating the left side of the histogram (lowposterior probability values) from the right side of the histogram (highposterior provability values).
To summarize, classiﬁcation of 50 hex-cells in each country shows
variations in the optimal parameters to the classiﬁer. The countries
diﬀer in the optimal deﬁnition of DMSP-OLS pixels as highly-lit, as wellas in the deﬁnition of a Landsat pixel as vegetation and in the optimal
inputs to the classiﬁers. Based on these results, we next map the builtup land cover in the entire area of each country.

75th

Table 3
Performance measures for classiﬁcation of 50 1-degree hex-cells (classiﬁcation with the
optimal parameters to the classiﬁer): True positive rate (TPR), True negative rate (TNR),
balanced accuracy and User's accuracy.

India
US
MexicoTPR (producer's
accuracy)

TNR

Balancedaccuracy

User's accuracy

77.5%
82.5%
88.7%

84.1%
75.1%
78.9%

80.8%
78.8%
83.8%

60.1%
67.3%
67.6%

3.3. Classiﬁcation of the entire area of each country

US and Mexico, while classiﬁcation with NDSV as input to the classiﬁer
results in the best performance of the classiﬁers for India. With theseparameters (Table 2), classiﬁcation of the 50 hex-cells yields an accurate classiﬁcation of built-up land cover, indicated by a high balanced
accuracy of 79%, 81% and 84% for classiﬁcation in the US, India and
Mexico, respectively (Table 3). User's accuracy ranges between 60% in
India to 67% and 68% in the US and Mexico, and the TNR and TPR
measures are above 75% in all countries.

In the experiments described above, we treated each hex-cell (1degree distance from center to center) as the unit of analysis for training
and for classiﬁcation. To determine the optimal scale of the hex-cells in
the hexagonal grid division, we next classify the entire area of eachcountry using three sizes of hexagonal grid divisions: 1, 4, and 8 decimal degrees from center to center. We calculate the thresholds tocharacterize highly-lit and lowly-lit DMSP-OLS pixels within each (1, 4or 8°) hex-cell. As shown in Fig. 8 (the spatial distribution of the DN
threshold to determine a highly-lit pixel, per hex-cell scale and per261

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

b

c

d

200000180000160000140000120000100000
80000
60000
40000
20000
0

200000180000160000140000120000100000
80000
60000
40000
20000
0

100000
90000
80000
70000
60000
50000
40000
30000
20000
10000
0
0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

e

f

g

h

100000
90000
80000
70000
60000
50000
40000
30000
20000
10000
0

350000300000250000300000250000200000250000200000200000150000150000150000100000100000100000

50000

0

0

0
0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

50000
0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

50000

i

j

k

l

100000
90000
80000
70000
60000
50000
40000
30000
20000
10000
0

350000300000250000300000250000200000250000200000200000100000100000100000
50000

50000

50000

0

0

0
0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

150000150000150000

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

0
0.125
0.25
0.375
0.5
0.625
0.75
0.875
1

a
200000180000160000140000120000100000
80000
60000
40000
20000
0

Fig. 7. The histograms of the posterior probability values of all pixels in 12 1-degree hex-cell in India (a–d), US (e–h) and Mexico (i–l). The dashed line indicates the posterior probability
threshold according to Otsu method.

When classiﬁcation is performed with 8-degree hex-cells, the balanced
accuracy falls by up to 4.6% and the Producer's accuracy drops moresigniﬁcantly by up to 13.4% (for example, from 63.2% to 49.8% inIndia). The countries also diﬀer in the level of hex-cell that results in thehighest User's accuracy rate. For example, we ﬁnd the highest User's
accuracy in India for classiﬁcation with hex-cell level of 4-degree while
in Mexico for classiﬁcation with hex-cell level of 1-degree (70.4% and76.7%, respectively).

country), there is large variation across and within countries. Thisvariation is also clearly illustrated in Fig. 9, the distribution histogramof the DN threshold to determine highly-lit and low-lit pixels.
Next, we mask out vegetation and bodies of water from highly-lit
pixels, select built-up and not built-up examples, and perform supervised image classiﬁcation with NDVI, NDWI, NDBI, EVI, UI (in the
US and Mexico) and NDSV (in India) as additional inputs to the classiﬁer. We classify each hex-cell as an independent unit of analysis andpredict, for each pixel, the probability (posterior probability) it is builtup. The mean value of the posterior probability ranges between 0.29
and 0.43 and the distribution of the posterior probability thresholdwithin each country is skewed left. Fig. 10 presents the distribution of
the per hex-cell posterior probability threshold in each country, foreach hex-cell scale.
As shown in Table 4, there are also variations between countries in
the optimal hex-cell scale. Classiﬁcation with the smallest level of hexcell (1-degree) yields the best performance in India and Mexico, whileclassiﬁcation with 4-degrees hex-cell scale yields the best performancein the US (though in Mexico, classiﬁcation with 1-degree hex-cell only
marginally exceeds classiﬁcation with 4-degree hex-cells, by 0.4%). Of
the three countries, we ﬁnd the highest accuracy performance forclassiﬁcation in the US, indicated by a balanced accuracy of 82% and
Producer's accuracy of 73% (for classiﬁcation with 4-degree hex-cells).
The balanced accuracy in India and in Mexico is a bit lower (around77%) and the Producer's accuracy ranges between 61% to 63% (in India
and Mexico, respectively, for classiﬁcation with 1-degree hex-cells).

3.4. Assessment of DMSP-OLS characterization of highly- and low-lit pixels
To evaluate the accuracy of our characterization of highly-lit andlow-lit DMSP-OLS pixels, we perform an additional per-hex-cell accuracy assessment. In each country, we adopt the optimal parameters to
the classiﬁer (i.e., hex-cell scale, threshold to characterize a Landsat
pixel as vegetation and threshold to characterize a DMSP-OLS pixel as
highly-lit or as lowly-lit) to identify highly-lit DMSP-OLS pixels (i.e.,
locations of built-up examples) and lowly-lit pixels (i.e., locations of not
built-up examples). We compare this characterization with our validation set. This analysis is intended to assess our deﬁnition of highly-lit
and low-lit regions from which we sample built-up and not built-up
examples (we relate to the characterization of DMSP-OLS pixels as thepredicted class). We ﬁnd that the balanced accuracy ranges between
67.2% in Mexico to 77% and 79.3% in India and the US, respectively.
Furthermore, when vegetation and bodies of water are not excluded
from highly-lit pixels the balanced accuracy drops in all countries, by262

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 8. The threshold to identify a pixel as highly-lit, per level of hex-cells (1, 4, and 8°) and per country (India, Mexico and the US).

up to 7.5%. By their nature, DMSP-OLS data are relatively “noisy” andsuﬀer from many limitations (e.g., coarse spatial resolution, blooming
eﬀect), yet by relying on local thresholds to characterize highly-lit andlowly-lit pixels as well as by removing vegetation land cover and bodies
of water, we improve the identiﬁcation of built-up and not built-up
locations.

3.5. Accuracy assessment with validation set
The results indicate that the smallest level of hex-cells (1-degree)
results in the best performance of the classiﬁers in India and in Mexico,while classiﬁcation with 4-degrees hex-cells results in the best performance in the US. Before evaluating these classiﬁcation maps with thevalidation set and with the Tomnod examples, we post-process these263

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.
70%

70%

70%

60%

60%

60%

50%

50%

50%

40%

40%

40%

30%

30%

30%

20%

20%

20%

10%

10%

10%

0%

0%

a

0%
5 10 15 20 25 30 35 40 45 50 55 60 63

5 10 15 20 25 30 35 40 45 50 55 60 63b

5 10 15 20 25 30 35 40 45 50 55 60 63

c

70%

70%

70%

60%

60%

60%

50%

50%

50%

40%

40%

40%

30%

30%

30%

20%

20%

20%

10%

10%

10%

0%

0%
5 10 15 20 25 30 35 40 45 50 55 60 63

d0%
5 10 15 20 25 30 35 40 45 50 55 60 63e

5 10 15 20 25 30 35 40 45 50 55 60 63

fFig. 9. The distribution of the threshold to identify a pixel as highly-lit ((a)India; (b)US; (c) Mexico), and low-lit ((d) India; (e) US; (f) Mexico) per level of hex-cells (1, 4, and 8°).

4. Discussion

maps. We clip them at the extent of each country's border, and remove
misclassiﬁed built-up pixels from areas with gas ﬂares using the Elvidge
et al. (2009) database. The ﬁnal post-processed classiﬁcation maps span
the full area of the three countries (Fig. 11).
The confusion matrix for the accuracy assessment with the validation set is presented in Table 5. The results show a high balanced accuracy rate of around 75% in India and Mexico and around 81% in the
US (Table 5). The accuracy measures for assessment with our validation
set and Tomnod examples are very close (Table 5). User's accuracy
ranges between 66% and 77% and the overall accuracy measure is
above 80% in all countries, (Table 5).

We present a machine learning approach to map built-up areas at alarge geographic scale. Our methodology utilizes nighttime-light data
(derived from DMSP-OLS) as a source for training examples of built-up
and non-built-up areas. These examples are used for supervised image
classiﬁcation in Landsat 8 imagery. To our knowledge, this is the ﬁrst
study to develop a practical and simple pixel-based image classiﬁcationof built-up areas that utilizes nighttime light data as a source of training
examples and high-resolution daytime imagery as the input source forthe classiﬁer. This methodology can be applied to map built-up areas
across space. Although many classiﬁcation products map urban land,
they are typically limited in their temporal and/or spatial resolution

70%

70%

70%

60%

60%

60%

50%

50%

50%

40%

40%

40%

30%

30%

30%

20%

20%

20%

10%

10%

10%

0%

0%
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

a0%
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

b

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1cFig. 10. The distribution of the posterior probability threshold to determine a pixel as “built-up” in India, the US and Mexico.

264

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

(OSM) can also be used to map urban areas (Belgiu and Drǎguţ, 2014;
Estima and Painho, 2015). OSM is a valuable source for ground-reference data because of its geographic extent and availability, but OSM
is subject to concerns about completeness and reliability (Schlesinger,2015). The use of OSM for supervised image classiﬁcation remains
challenging due to the risk of imbalanced distribution of class labels(including spatial coverage), the presence of errors or missing classassignments (“class-noise”), and inaccurate polygon boundary delineations (Johnson and Iizuka, 2016). Our methodology overcomes thelack of such data by utilizing low-resolution DMSP-OLS data for classiﬁcation of built-up areas in Landsat imagery. By partitioning countries
into smaller regions, we allow the parameters of the classiﬁcation
model to vary across space. The results indicate that the three examinedcountries diﬀer in the optimal hex-cell scale, or the size of the hex-cellthat is used as the unit of analysis (for training and for classiﬁcation).
Classiﬁcation at a local scale (a hex-cell) allows one to account forheterogeneity in environmental conditions. Although many studies
address the eﬀect of the classiﬁers' hyper-parameters on their performance, here we show that classiﬁers also have an optimal spatial scalethat can be and should be discovered through techniques similar to
those we employ here. We leave the choice of the classiﬁer scale foreach pixel for future research.

Table 4
Performance measures as a factor of the hex-cell level for India, the US, and Mexico.
Accuracy assessment with test set (classiﬁcation of all hexagons in each country).
Overallaccuracy

TPR (producer's
accuracy)

TNR

Balanced

User'saccuracyIndia
8° 84.2%
4° 84.7%
1° 83.8%

49.8%
53.9%
63.2%

94.0%
93.5%
89.7%

71.9%
73.7%
76.5%

70.3%
70.4%
63.9%

US
8° 87.2%
4° 86.8%
1° 85.9%

68.8%
72.8%
71.1%

92.7%
91.0%
90.4%

80.8%
81.9%
80.7%

73.6%
70.6%
68.7%

Mexico8° 80.5%
4° 80.9%
1° 82.4%

64.1%
64.1%
61.4%

87.7%
88.3%
91.7%

75.9%
76.2%
76.6%

69.6%
70.6%
76.7%

which limits their use to track urbanization processes over time.
Mapping built-up areas at a global scale is challenging because of
the scarcity of extensive ground-reference data for supervised classiﬁcation and validation. Crowd-sourced datasets, such as OpenStreetMap

Fig. 11. Classiﬁed built-up land cover (in red) across the full countries: (a) Mexico, (b) India and (c) the US. Post-processed classiﬁcation maps, using 1-degree hex-cells. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

265

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Table 5
Performance measures as a factor of the hex-cell level for India, the US, and Mexico after post processing. Accuracy assessment with validation set and with Tomnod examples.
Overall accuracy

TPR (producer's accuracy)

TNR

Balanced

User's accuracy

Post processed maps – validation with validation set
India (1-dgree)
84.0%
US (4-degree)
86.2%
Mexico (1-degree)
82.0%

59.6%
70.7%
58.6%

91.0%
90.8%
92.3%

75.3%
80.7%
75.4%

65.6%
69.0%
77.1%

Post processed maps – validation with Tomnod set
India (1-dgree)
83.6%
US (4-degree)
85.9%
Mexico (1-degree)
80.2%

56.5%
72.1%
55.0%

91.7%
92.1%
91.9%

74.1%
82.1%
73.5%

67.2%
80.3%
76.0%

Confusion matrix for Post processed maps (validation with validation set)⁎
India

US

MexicoPredicted

Actual

⁎

BU
NBU
Sum

BU

NBU

SUM

BU

NBU

SUM

BU

NBU

SUM

1521
799
2320

1031
8110
9141

2552
8909
11,461

1734
778
2512

719
7648
8367

2453
8426
10,879

1961
584
2545

1385
6993
8378

3346
7577
10,923

Confusion matrices for classiﬁcation with examples in the validation set. Comparison between predicted and actual built-up (BU) and not-built-up (NBU) examples.

the ﬁne boundaries of built-up areas with high precision (Fig. 13). We
note that although the ground-reference data are not free of measurement error, this methodology allows us to collect millions of labeled
examples. This large sample is valuable for training a classiﬁer to map
built-up land cover.
We used 60,000 labeled examples (polygons) to assess the accuracyof our classiﬁcation compared to the deﬁnition of urban and built-up
land cover according to MCD12Q1 UMD MODIS classiﬁcation scheme.
The results show that the accuracy of our classiﬁcation exceeds
MCD12Q1 UMD MODIS classiﬁcation scheme, indicated by a higher
balanced accuracy of between 2.5%–5.6% (Table 6). As illustrated inFig. 14, the extent of the built-up land cover that we detect with our
methodology is larger than the extent of the land cover that is classiﬁed
as urban and built-up by the MCD12Q1 classiﬁcation.
Our classiﬁcation also exceeds other national high-resolution landcover and land-use maps. To illustrate, a comparison between our examples and the US National Land Cover Database (NLCD) classiﬁcation
map showed a lower balanced accuracy rate in the NLCD product(72.3%) (because we deﬁne a polygon as “built-up” if 50% or greater of
its area is built, we compare our classiﬁcation to the NLCD classes“Developed, Medium Intensity” and “Developed High Intensity” as“built-up”). Fig. 15 presents a comparison between our classiﬁcation
and NLCD's classiﬁcation of medium and high- intensity in three casestudies in the US.
Finally, classiﬁcation of urban land cover can be performed using
Soft Classiﬁcation (SC) techniques. This technique is the creation of acontinuous surface where each pixel is characterized by the fraction ofurban land, or impervious surface that is incorporated within it (Huang
et al., 2016; MacLachlan et al., 2017; Song et al., 2016; Van de Voorde
et al., 2011). Alternatively, one can use Hard Classiﬁcation (HC) techniques, where the extent of urban land is characterized as discrete perpixel binary class (presence/absence of urban land). Many of the existing global urban maps (e.g., the Global Rural-Urban Mapping Project(GRUMP), MODIS Urban Land cover and GLC00) represent the surface
on Earth according to a discrete thematic classiﬁcation (Potere et al.,2009). In part, this is because many of the common classiﬁers, such as
Random Forests, predict for each pixel the probability (posterior
probability) it belongs to a given class (e.g., a pixel is characterized asurban if the probability it is urban exceeds a given threshold). However,a discrete map can further be translated into percentage of urban land,for example, by aggregating high-resolution binary classiﬁcation to a

We assess the validity of our approach using an external validationthat uses a total of around 84,500 hand-labeled examples. The results
demonstrate the robustness of our approach and its applicability inheterogeneous regions. Yet, the degree to which the “localization” ofour classiﬁers aﬀects their performance varies between regions anddepends on the characteristics and heterogeneity of the mapped landscape. For example, because of the relatively similar spectral characteristics of built-up areas and bare land, it is often hard to diﬀerentiate between these land cover types through remote sensing,
especially in regions characterized by heterogenous land cover (e.g.
complex pattern of built-up areas, bare land and vegetation) (He et al.,2010). Thus, diﬀerences in the complexity of the land cover within andbetween countries, as well as between diﬀerent hexagonal scales, may
aﬀect the accuracy of the classiﬁcation as we observe in this study.
Applying this method at a global scale will require automatic methods
for selecting, in any pixel, the scale of classiﬁcation that maximizes
accuracy. Of course, the accuracy of the classiﬁcation also depends on
the quality of the ground-reference data. Our methodology performs
well especially when built-up and not built-up examples can be suﬃciently diﬀerentiated in highly-lit areas. Because we rely on NDVI andNDWI indices to separate between built-up and not built-up examples inhighly-lit DMSP-OLS pixels, our methodology may work best in regionsthat are characterized by vegetation land cover. Future research could
rely on our method to improve accuracy in settings that are characterized by bare land or bare soil, for example, by utilizing the Enhanced Built-Up and Bareness Index (EBBI) (As-Syakur et al., 2012). In
addition, due to the lack of on-board calibration and unstable radiometric performance of the DMSP-OLS sensors, the absolute radiance oflight cannot directly represent temporal changes in the intensity of thelight, and thus, inter-sensor calibration is required to make our approach operational in time.
Our methodology overcomes the need for expensive hand-labeleddata for supervised classiﬁcation as well as many of the limitations
associated with DMSP-OLS data. First, by identifying local per-hex-cell
thresholds to characterize DMSP-OLS pixels as highly-lit, we accountfor the blooming eﬀect and identify the boundaries of highly-lit settlements. As illustrated in Fig. 12, due to the blooming eﬀect, lit areas
are consistently larger than the built-up land cover they are associatedwith. Second, by identifying vegetation land cover and bodies of wateraccording to NDVI and NDWI indices, we minimize the eﬀect of theover saturation of DMSP-OLS in city centers, which allows us to capture266

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

a

b

c

d

e

f

g

h

i

Fig. 12. A comparison between our classiﬁcation of built-up areas and lit pixels according to DMSP-OLS in (a) Ahmedabad, Gujarat, India; (b) New Delhi, Delhi, India; (c) Kolkata, West
Bengal, India; (d) Phoenix, Arizona, US; (e) Washington DC, US; (f) Denver, Colorado, US; (g) Mexico City, Mexico; (h) Guadalajara, Mexico; (i) Puebla, Mexico. The top ﬁgure in each citypresents our classiﬁcation; the bottom ﬁgure presents the DN values of DMSP-OLS stable lights band: 30–55 (blue), 56–61 (green), 63 (red). (For interpretation of the references to colourin this ﬁgure legend, the reader is referred to the web version of this article.)(according to our classiﬁcation) in each cell. Fig. 16 presents, as an
example, per-pixel fractions of built-up land cover in nine cities in ourinvestigated countries. The ﬁgure demonstrates variations in built-up
density, both within and between cities. Some cities are characterized
by a clear gradient from the center of the city towards the periphery

continuous coarser resolution surface representing per pixel percentage
urban land (Schneider et al., 2009, 2010; Zhou et al., 2015). To illustrate how our classiﬁcation maps can be translated into a continuous
characterization of built-up land cover, we create an artiﬁcial grid of
500 m ∗ 500 m cells and calculate the fraction of built-up land cover267

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

a

b

c

d

Fig. 13. Classiﬁcation of built-up areas (in orange) in (a) Nagpur, India; (b) Hyderabad, India; (c) Houston US; (d) Mexico City, Mexico.

(e.g., Hyderabad, India and Culiacán, Mexico) while other cities (e.g.,
Atlanta, Georgia, US), are less dense and the built-up density is more
evenly distributed throughout the city. Future research could also build
on our methodology to estimate sub-pixel fraction of built up land
cover, for example, by adopting the FCM and PCM algorithms to derive
relative and absolute measures of class membership within a pixel(Foody, 2002) or by applying Spectral Unmixing methods (e.g., Mitraka
et al., 2016) that are especially useful for classiﬁcation in complex
urban/peri-urban settings (Poursanidis et al., 2015).
To summarize, we have developed a framework for pixel-based
image classiﬁcation of built-up areas that uses hexagonal tessellation tocombine nighttime and daytime light data. The methodology yields
high-resolution, high quality depictions of built-up areas across threehighly diverse countries.
The study contributes to the existing literature in several aspects.
First, previous studies that combine nighttime and daytime data havebeen limited either in their spatial generalization or in their spatial

Table 6
Accuracy assessment (balanced accuracy) of our methodology for classiﬁcation, MODISMCD12Q1 (MODIS), DMSP-OLS highly-lit areas (deﬁned as pixels with a DN value above
the 95th or 90th percentile), Accuracy assessment with 60,000 labeled examples.
Our BU classiﬁcation

MODIS UMD

India
TPR (producer's accuracy)
TNR
Balanced accuracy

59.9%
91.2%
75.5%

61.7%
84.3%
73.0%

US
TPR (producer's accuracy)
TNR
Balanced accuracy

71.8%
90.9%
81.4%

64.4%
87.3%
75.8%

MexicoTPR (producer's accuracy)
TNR
Balanced accuracy

62.7%
92.6%
77.6%

53.4%
94.1%
73.7%

a

b

c

d

Fig. 14. A comparison between areas classiﬁed as built-up using our methodology (top) and areas classiﬁed as built up and urban by MCD12Q1 UMD MODIS classiﬁcation scheme(bottom), in (a) Ahmedabad, Gujarat, India; (b) New Delhi, Delhi, India; (c) Phoenix, Arizona, US; (d) Mexico City, Mexico.

268

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

NLCD

Our classification

Reference

1.

1a

1b

3a

2b

3a

3b

2.

3.

Fig. 15. A comparison of our classiﬁcation (column a) and NLCD's classiﬁcation* (column b), in (1) Columbus, Ohio; (2) Las Vegas, Nevada; (3) Atlanta, Georgia.
*Classes 23–24: Developed, Medium and Intensity (Impervious surfaces account for 50% to 100% of the total cover).

In today's era of big data, a globally consistent and data-driven
method of deﬁning and classifying urban areas has extensive applications. Economics, urban planning, climate modeling, water-resource
management, hazard-response eﬀorts, and urban-ecosystem assessments all use geographic data on urban areas. Having information on
urban extent that is spatially and temporally consistent and deﬁned athigh resolution is both essential to a wide range of disciplines and necessary for helping society better understand the drivers of urbanization, therefore promoting sustainable urban development.

resolution. To cope with regional variations, we adopt a hexagonal
tessellation mapping approach to handle heterogeneous land cover withcomplex spectral characteristics. Second, despite signiﬁcant progress inapplying machine-learning techniques to process satellite data, there isa paucity of reference datasets to detect urban areas. The approach we
present uses DMSP-OLS and Landsat data to collect training examples
and Landstat data for the input to the classiﬁer. Third, to determine thebinary-class of a pixel (i.e., “built-up” or “not built-up”) we use the Otsu
algorithm to determine the threshold of the posterior probability foreach hex-cell. Finally, our approach requires minimal manual input fortraining and for classiﬁcation and does not rely on local adjustableparameters to the classiﬁers, making it scalable across space and time,
which is the subject of ongoing work.

Acknowledgements
MG was supported by the National Science Foundation Grants EAR269

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

Fig. 16. percentage of built-up land cover in: India (a) New Delhi, (b) Ahmedabad, (c) Hyderabad; US (d) Atlanta, (e) Phoenix, (f) Kansas City; and Mexico (g) Culiacán, (h) Monterrey, (i)Gómez Palacio. The value of each 500 m*500 m pixel represent the percentage of built-up land cover (according to Landsat) that is incorporated in it.

Centre (IGC project number 89448). RG and GH were supported by theCenter on Global Transformation at UCSD. We thank the Institute ofGeography, National Autonomous University of Mexico, Dr. Armando
Peralta Higuera, who provided support for two students, Juan Alberto
Guerra Moreno and Kimberly Mendez Gomez for validating the algorithm in Mexico.

1204774 and DMS-1419593, the U.S. Department of Agriculture NIFA
grant 2015-67003-23508. MG and MS were also supported by the
National Science Foundation Sustainability Research Network (SRN)
Cooperative Agreement 1444758, the Urban Water Innovation Network(UWIN). MS was also supported by Dr. B. L. Turner II's Gilbert F. White
Environment and Society Fellowship at Arizona State University. AK,
RG and GH were supported by funds from the International Growth
Appendix A. Indices used as inputs to the classiﬁers

Table A1presents the bands and the indices used as inputs to the classiﬁers. We now turn to describe the indices in details.
I. NDVI (Normalized Diﬀerence Vegetation Index)
270

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

NDVI expresses the relation between red visible light (which is typically absorbed by a plant's chlorophyll) and near-infrared wavelength (whichis scattered by the leaf's mesophyll structure) (Pettorelli et al., 2005). It is computed as:(NIR − RED)/(NIR + RED)

(A1)

where NIR is the near infra-red wavelength and RED is the red wavelength. The values of NDVI range between (− 1) and (+ 1). An average NDVI
value was calculated for each pixel.
II. NDBI (Normalized Diﬀerence Built-up Index)
NDBI expresses the relation between the medium infra-red and the near infra-red wavelengths (Zha et al., 2003). It is computed as:(A2)(MIR − NIR)/(MIR + NIR)

where MIR is the medium infra-red and NIR is the near infra-red wavelength. The index assumes a higher reﬂectance of built-up areas in the medium
infra-red wavelength range than in the near infra-red.
III. NDWI (Normalized Diﬀerence Water Index)
NDWI expresses the relation between the green (G) and the NIR (near infra-red) wavelength, with a scaling of − 1 to +1 (McFeeters, 1996). It is
computed as:(A3)(G − NIR)/(G + NIR)
The positive values typically represent open water areas, while negative values represent non-water features, like soil and vegetation.
IV. UI (Urban Index)
UI is the diﬀerence between the short infra-red (SWIR) and the near-infrared wavelengths (Kawamura et al., 1996). It is computed as:(A4)(SWIR−NIR)/(SWIR+NIR)
Like NDBI, UI assumes high brightness in SWIR in urban areas as opposed to the NIR.
V. EVI (Enhanced Vegetation Index)

EVI is an improved vegetation index with higher sensitivity in high biomass regions where NDVI tends to saturate, reduces atmospheric inﬂuences, and removes the canopy background signal (Huete et al., 2002). It is computed as:(2.5∗ ((NIR − R)/(NIR + 6R − 7. 5B + 1))(A5)This is a similar formula to the NDVI, which takes advantage of high reﬂectance of vegetation in the NIR band as opposed to the R band. The blue
band (B) is used (with a coeﬃcient of 7.5) to correct for aerosol inﬂuences in the red band (R). There is an integer of 1 added to the denominator toadjustment for nonlinear NIR and R radiant transfer through canopies. 2.5 is applied as a gain to the index.
VI. NDSV (Normalized Diﬀerence Spectral Vector)
NDSV (Normalized Diﬀerence Spectral Vector) (Trianni et al., 2015) is used to identify consistent human settlement and artiﬁcial materials. The
index is composed out of all possible normalized indexes combining two bands with the same spatial resolution. Each pixel is characterized by aspectral vector of normalized indexes, which maps the original multispectral pixel into a higher dimensional space (Trianni et al., 2015). The index is
calculated with each of Landsat's 30 m resolution bands as follows (Trianni et al., 2015):

NDSVij =Bi − Bj
Bi + Bj

where Bi and Bj are two generic bands.
Table A1
Bands in Landsat and indices used.

Landsat 8
B1
B2
B3
B4
B5
B6
B7
B8
B10

Spectral band

Wavelength (micrometers)

Resolution (meters)

Band 1 – Ultra blue
Band 2 - Blue
Band 3 - GreenBand 4 – Red
Band 5 - Near Infrared (NIR)
Band 6 - SWIR 1
Band 7 - SWIR 2
Band 8 - PanchromaticBand 10 - Thermal Infrared (TIRS) 1

0.43–0.45
0.45–0.51
0.53–0.59
0.64–0.67
0.85–0.88
1.57–1.65
2.11–2.29
0.50–0.68
10.60–11.19

30
30
30
30
30
30
30
15
100 (resampled to 30)271

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

B11
NDVI
NDWI
NDBI
EVI
UI

Band 11 - Thermal Infrared (TIRS) 2(B5 − B4) / (B5 + B4)(B3 − B5) / (B3 + B5)(B6 − B5) / (B6 + B5)
2.5 ∗ ((B5 / B4) / (B5 + 6 ∗ B4 − 7.5 ∗ B2 + 1))(B7 − B5) / (B7 + B5)

11.50–12.51

100 (resampled to 30)
30
30
30
30
30

Appendix B. Performance measures with ﬁfty hex-cells
Table B1
Performance measures with ﬁfty-hex-cells per country.
NDVI percentile

Highly lit pixels

Inputs to classiﬁer

Low-lit pixels

TPR

TNR

Balanced

User's accuracy

Optimal NDVI threshold
India
50
75
90
95
99
US
50
75
90
95
99
Mexico
50
75
90
95
99

99
99
99
99
99
99
99
99
99
99
99
99
99
99
99

All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*

99
99
99
99
99
99
99
99
99
99
99
99
99
99
99

56.16%
58.24%
61.27%
58.75%
59.04%
65.84%
72.42%
76.14%
76.76%
76.22%
69.94%
72.96%
72.40%
72.22%
72.61%

93.23%
92.25%
91.18%
91.27%
90.67%
84.51%
80.71%
79.85%
78.79%
77.06%
94.63%
93.55%
92.90%
92.87%
92.61%

74.69%
75.25%
76.23%
75.01%
74.85%
75.18%
76.57%
77.99%
77.78%
76.64%
82.28%
83.26%
82.65%
82.55%
82.61%

71.89%
69.86%
67.33%
67.49%
66.13%
72.53%
69.99%
70.11%
69.20%
67.35%
86.59%
84.87%
83.50%
83.39%
82.97%

Highly-lit pixels
India
90
90
90
90
90
US
90
90
90
90
90
Mexico
75
75
75
75
75

99
95
90
75
50
99
95
90
75
50
99
95
90
75
50

All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*
All*

99
99
99
99
99
99
99
99
99
99
99
99
99
99
99

58.89%
72.71%
78.98%
87.62%
91.00%
76.14%
77.15%
80.64%
87.45%
94.72%
72.96%
72.71%
86.01%
92.90%
95.25%

91.18%
85.57%
77.75%
60.56%
58.23%
79.85%
78.55%
76.43%
64.26%
37.42%
93.55%
92.66%
82.58%
58.65%
39.42%

75.04%
79.14%
78.36%
74.09%
74.61%
77.99%
77.85%
78.53%
75.86%
66.07%
83.26%
82.69%
84.29%
75.77%
67.34%

67.33%
60.84%
52.26%
40.66%
40.19%
70.11%
69.07%
67.99%
60.31%
49.81%
84.87%
83.09%
71.00%
52.69%
43.81%

Inputs to classiﬁer
India
90
90
90
90
90
US
90
90
90
90
90
Mexico
90
90
90
90
90

95
95
95
95
95
90
90
90
90
90
90
90
90
90
90

No Indices
NDVI
NDBI
All*
NDSVI
No Indices
NDVI
NDBI
All*
NDSVI
No Indices
NDVI
NDBI
All*
NDSVI

99
99
99
99
99
99
99
99
99
99
99
99
99
99
99

71.13%
71.85%
70.34%
72.71%
73.87%
80.25%
78.93%
80.09%
80.64%
80.17%
86.11%
84.81%
85.48%
86.01%
85.41%

86.52%
86.28%
86.41%
85.57%
86.05%
74.51%
75.32%
75.61%
76.43%
76.48%
82.44%
82.72%
82.48%
82.58%
82.69%

78.83%
79.06%
78.37%
79.14%
79.96%
77.38%
77.13%
77.85%
78.53%
78.32%
84.28%
83.77%
83.98%
84.29%
84.05%

61.94%
61.76%
61.49%
60.84%
62.03%
66.16%
66.51%
67.10%
67.99%
67.91%
70.86%
70.88%
70.75%
71.00%
70.98%

Low-lit pixels

272

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

India

US

Mexico

90
90
90
90
90
90
90
90
90
90
90
90

95
95
95
95
90
90
90
90
90
90
90
90

NDSV
NDSV
NDSV
NDSV
All*
All*
All*
All*
All*
All*
All*
All*

10
25
50
75
10
25
50
75
10
25
50
75

84.16%
82.72%
79.91%
77.54%
93.57%
91.87%
88.30%
82.49%
93.28%
92.65%
91.42%
88.75%

69.04%
74.82%
80.70%
84.12%
44.49%
56.18%
67.48%
75.08%
59.50%
64.07%
71.64%
78.87%

76.60%
78.77%
80.31%
80.83%
69.03%
74.02%
77.89%
78.79%
76.39%
78.36%
81.53%
83.81%

45.61%
50.33%
56.09%
60.10%
51.14%
56.56%
62.78%
67.28%
53.32%
56.11%
61.51%
67.56%

Observation, Remote Sensing Applications Series. CRC Press, pp. 97–120.
Estima, J., Painho, M., 2015. Investigating the Potential of OpenStreetMap for Land Use/
Land Cover Production: A Case Study for Continental Portugal. In: Arsanjani, J.J.,
Zipf, A., Mooney, P., Helbich, M. (Eds.), OpenStreetMap in GIScience, Lecture Notes
in Geoinformation and Cartography. Springer International Publishing, pp. 273–293.
Estoque, R.C., Murayama, Y., 2015. Classiﬁcation and change detection of built-up lands
from Landsat-7 ETM + and Landsat-8 OLI/TIRS imageries: a comparative assessmentof various spectral indices. Ecol. Indic. 56, 205–217. http://dx.doi.org/10.1016/j.
ecolind.2015.03.037.
Flood, N., 2013. Seasonal composite Landsat TM/ETM + images using the Medoid (amulti-dimensional median). Remote Sens. 5, 6481–6500. http://dx.doi.org/10.3390/
rs5126481.
Foody, G.M., 2002. Status of land cover classiﬁcation accuracy assessment. Remote Sens.
Environ. 80, 185–201. http://dx.doi.org/10.1016/S0034-4257(01)00295-4.
Gamba, P., Aldrighi, M., Stasolla, M., 2011. Robust extraction of urban area extents in HR
and VHR SAR images. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 4, 27–34.
http://dx.doi.org/10.1109/JSTARS.2010.2052023.
Gao, B.C., 1996. NDWI—A normalized diﬀerence water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 58 (3), 257–266.
Gaughan, A.E., Stevens, F.R., Linard, C., Jia, P., Tatem, A.J., 2013. High resolution population distribution maps for Southeast Asia in 2010 and 2015. PLoS One 8, e55882.
http://dx.doi.org/10.1371/journal.pone.0055882.
Giri, C., Pengra, B., Long, J., Loveland, T.R., 2013. Next generation of global land covercharacterization, mapping, and monitoring. Int. J. Appl. Earth Obs. Geoinf. 25,
30–37. http://dx.doi.org/10.1016/j.jag.2013.03.005.
Gislason, P.O., Benediktsson, J.A., Sveinsson, J.R., 2006. Random forests for land coverclassiﬁcation. Pattern Recogn. Lett. 27, 294–300. http://dx.doi.org/10.1016/j.
patrec.2005.08.011.
Goldblatt, R., You, W., Hanson, G., Khandelwal, K.A., 2016. Detecting the boundaries ofurban areas in India: a dataset for pixel-based image classiﬁcation in Google earth
engine. Remote Sens. 8. http://dx.doi.org/10.3390/rs8080634.
Gong, P., Wang, J., Yu, Le, Zhao, Yongchao, Zhao, Yuanyuan, Liang, L., Niu, Z., Huang,X., Fu, H., Liu, S., Li, C., Li, X., Fu, W., Liu, C., Xu, Y., Wang, X., Cheng, Q., Hu, L.,
Yao, W., Zhang, Han, Zhu, P., Zhao, Z., Zhang, Haiying, Zheng, Y., Ji, L., Zhang, Y.,
Chen, H., Yan, A., Guo, J., Yu, Liang, Wang, L., Liu, X., Shi, T., Zhu, M., Chen, Y.,
Yang, G., Tang, P., Xu, B., Giri, C., Clinton, N., Zhu, Z., Chen, Jin, Chen, Jun, 2013.
Finer resolution observation and monitoring of global land cover: ﬁrst mapping results with Landsat TM and ETM + data. Int. J. Remote Sens. 34, 2607–2654. http://dx.doi.org/10.1080/01431161.2012.748992.
Guan, H., Li, J., Chapman, M., Deng, F., Ji, Z., Yang, X., 2013. Integration of orthoimagery
and lidar data for object-based urban thematic mapping using random forests. Int. J.
Remote Sens. 34, 5166–5186. http://dx.doi.org/10.1080/01431161.2013.788261.
He, C., Shi, P., Xie, D., Zhao, Y., 2010. Improving the normalized diﬀerence built-up index
to map urban built-up areas using a semiautomatic segmentation approach. IEEE
Geosci. Remote Sens. Lett. 1, 213–221.
Henderson, M., Yeh, E.T., Gong, P., Elvidge, C., Baugh, K., 2003. Validation of urban
boundaries derived from global night-time satellite imagery. Int. J. Remote Sens. 24,
595–609. http://dx.doi.org/10.1080/01431160304982.
Herold, M., 2009. Some recommendations for global eﬀorts in urban monitoring andassessments from remote sensing. In: Global Mapping of Human Settlement, Remote
Sensing Applications Series. CRC Press. http://dx.doi.org/10.1201/
9781420083408-c2.
Homer, C., Gallant, A., 2001. Partitioning the Conterminous United States into Mapping
Zones for Landsat TM Land Cover Mapping. USGS.
Homer, C., Huang, C., Yang, L., Wylie, B.K., Coan, M., 2004. Development of a 2001
National Land-Cover Database for the United States. Photogramm. Eng. Remote.
Sens. 7, 829–840.
Hsu, F.-C., Baugh, K.E., Ghosh, T., Zhizhin, M., Elvidge, C.D., 2015. DMSP-OLS radiancecalibrated nighttime lights time series with intercalibration. Remote Sens. 7,
1855–1876. http://dx.doi.org/10.3390/rs70201855.
Huang, X., Schneider, A., Friedl, M.A., 2016. Mapping sub-pixel urban expansion in China
using MODIS and DMSP/OLS nighttime lights. Remote Sens. Environ. 175, 92–108.
http://dx.doi.org/10.1016/j.rse.2015.12.042.
Huete, A., Didan, K., Miura, T., Rodriguez, E.P., Gao, X., Ferreira, L.G., 2002. Overview of
the radiometric and biophysical performance of the MODIS vegetation indices.
Remote Sens. Environ. 83, 195–213. http://dx.doi.org/10.1016/S0034-4257(02)00096-2.

References
Aguilar, A.G., 2008. Peri-urbanization, illegal settlements and environmental impact in
Mexico City. Cities 25, 133–145. http://dx.doi.org/10.1016/j.cities.2008.02.003.
Aguilar, A.G., Guerrero, F.L., 2013. Poverty in peripheral informal settlements in Mexico
City: the case of Magdalena Contreras, Federal District. Tijdschr. Voor Econ. En Soc.
Geogr. 104, 359–378. http://dx.doi.org/10.1111/tesg.12012.
Angiuli, E., Trianni, G., 2014. Urban mapping in Landsat images based on normalized
diﬀerence spectral vector. IEEE Geosci. Remote Sens. Lett. 11, 661–665. http://dx.
doi.org/10.1109/LGRS.2013.2274327.
As-Syakur, A.R., Adnyana, I., Arthana, I.W., Nuarsa, I.W., 2012. Enhanced built-up andbareness index (EBBI) for mapping built-up and bare land in an urban area. Remote
Sens. 4 (10), 2957–2970.
Bagan, H., Yamagata, Y., 2014. Land-cover change analysis in 50 global cities by using a
combination of Landsat data and analysis of grid cells. Environ. Res. Lett. 9, 064015.
http://dx.doi.org/10.1088/1748-9326/9/6/064015.
Bagan, H., Yamagata, Y., 2015. Analysis of urban growth and estimating population
density using satellite images of nighttime lights and land-use and population data.
GISci. Remote Sens. 52, 765–780. http://dx.doi.org/10.1080/15481603.2015.
1072400.
Ban, Y., Jacob, A., Gamba, P., 2015. Spaceborne SAR data for global urban mapping at
30m resolution using a robust urban extractor. ISPRS J. Photogramm. Remote Sens.
103, 28–37. http://dx.doi.org/10.1016/j.isprsjprs.2014.08.004.
Belgiu, M., Drǎguţ, L., 2014. Comparing supervised and unsupervised multiresolution
segmentation approaches for extracting buildings from very high resolution imagery.
ISPRS J. Photogramm. Remote Sens. 96, 67–75. http://dx.doi.org/10.1016/j.
isprsjprs.2014.07.002.
Birch, C.P.D., Oom, S.P., Beecham, J.A., 2007. Rectangular and hexagonal grids used forobservation, experiment and simulation in ecology. Ecol. Model. 206, 347–359.
http://dx.doi.org/10.1016/j.ecolmodel.2007.03.041.
Boucher, A., Seto, K., 2009. Methods and challenges for using high-temporal resolution
data to monitor urban growth. In: Global Mapping of Human Settlement, Remote
Sensing Applications Series. CRC Press. http://dx.doi.org/10.1201/9781420083408c16.
Chen, X.L., Zhao, H.M., Li, P.X., Yin, Z.Y., 2006. Remote sensing image-based analysis of
the relationship between urban heat island and land use/cover changes. Remote Sens.
Environ. 104, 133–146. http://dx.doi.org/10.1016/j.rse.2005.11.016.
Chen, C.-F., Son, N.-T., Chang, N.-B., Chen, C.-R., Chang, L.-Y., Valdez, M., Centeno, G.,
Thompson, C.A., Aceituno, J.L., 2013. Multi-decadal mangrove forest change detection and prediction in Honduras, Central America, with Landsat imagery and a
Markov chain model. Remote Sens. 5, 6408–6426. http://dx.doi.org/10.3390/
rs5126408.
Chen, Jun, Chen, Jin, Liao, A., Cao, X., Chen, L., Chen, X., He, C., Han, G., Peng, S., Lu, M.,
Zhang, W., Tong, X., Mills, J., 2015. Global land cover mapping at 30 m resolution: a
POK-based operational approach. ISPRS J. Photogramm. Remote Sens. 103, 7–27.
http://dx.doi.org/10.1016/j.isprsjprs.2014.09.002.
CIESIN, 2005. Gridded Population of the World, Version 3 (GPWv3) Data Collection.
Clinton, N., Gong, P., 2013. MODIS detected surface urban heat islands and sinks: global
locations and controls. Remote Sens. Environ. 134, 294–304. http://dx.doi.org/10.
1016/j.rse.2013.03.008.
Connolly, P., 2014. Latin Amercian Informal Urbanism: Contexts, Concepts, andContributions with Speciﬁc Reference to Mexico. In: Becerra, A., Hernandez, F.
(Eds.), (Re)viewing Latin American Cities. Cambridge Scholars Publishing,Cambridge.
Consejo Nacional de Población, 2012. Sistema Urbano Nacional. Col. Juárez, Del.
Cuauhtémoc, Mexico.
Eakin, H., Lerner, A.M., Manuel-Navarrete, D., Hernández Aguilar, B., Martínez-Canedo,A., Tellman, B., Charli-Joseph, L., Fernández Álvarez, R., Bojórquez-Tapia, L., 2016.
Adapting to risk and perpetuating poverty: Household's strategies for managing ﬂood
risk and water scarcity in Mexico City. Environ. Sci. Pol. 66, 324–333. http://dx.doi.
org/10.1016/j.envsci.2016.06.006.
Elvidge, C.D., Ziskin, D., Baugh, K.E., Tuttle, B.T., Ghosh, T., Pack, D.W., Erwin, E.H.,
Zhizhin, M., 2009. A ﬁfteen year record of global natural gas ﬂaring derived fromsatellite data. Energies 2, 595–622. http://dx.doi.org/10.3390/en20300595.
Elvidge, C., Hsu, F.-C., Baugh, K., Ghosh, T., 2014. National trends in satellite-observed
lighting: 1992–2012. In: Global Urban Monitoring and Assessment through Earth273

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

A.M., Kemper, T., Soille, P., Syrris, V., 2016. Operating Procedure for the Productionof the Global Human Settlement Layer From Landsat Data of the Epochs 1975, 1990,
2000, and 2014. Publications Oﬃce of the European Union, Ispra (VA), Italy.
Pettorelli, N., Vik, J.O., Mysterud, A., Gaillard, J.M., Tucker, C.J., Stenseth, N.C., 2005.
Using the satellite-derived NDVI to assess ecological responses to environmental
change. Trends Ecol. Evol. 20 (9), 503–510.
Potere, D., Schneider, A., Angel, S., Civco, D.L., 2009. Mapping urban areas on a global
scale: which of the eight maps now available is more accurate? Int. J. Remote Sens.
30, 6531–6558. http://dx.doi.org/10.1080/01431160903121134.
Poursanidis, D., Chrysoulakis, N., Mitraka, Z., 2015. Landsat 8 vs. Landsat 5: a comparison based on urban and peri-urban land cover mapping. Int. J. Appl. Earth Obs.
Geoinf. 35, 259–269. http://dx.doi.org/10.1016/j.jag.2014.09.010.
Richards, T., Gallego, J., Achard, F., 2000. Sampling for forest cover change assessment at
the pan-tropical scale. Int. J. Remote Sens. 21, 1473–1490. http://dx.doi.org/10.
1080/014311600210272.
Rodriguez-Galiano, V.F., Ghimire, B., Rogan, J., Chica-Olmo, M., Rigol-Sanchez, J.P.,
2012. An assessment of the eﬀectiveness of a random forest classiﬁer for land-coverclassiﬁcation. ISPRS J. Photogramm. Remote Sens. 67, 93–104. http://dx.doi.org/10.
1016/j.isprsjprs.2011.11.002.
Rwanga, S.S., Ndambuki, J.M., 2017. Accuracy assessment of land use/land cover classiﬁcation using remote sensing and GIS. Int. J. Geosci. 8, 611.
Schlesinger, J., 2015. Using Crowd-Sourced Data to Quantify the Complex Urban Fabric OpenStreetMap and the Urban–Rural Index. In: Arsanjani, J., Zipf, A., Mooney, P.,
Helbich, M. (Eds.), OpenStreetMap in GIScience. Springer International Publishing,
Switzerland, pp. 295–315.
Schneider, A., Friedl, M.A., Potere, D., 2009. A new map of global urban extent from
MODIS satellite data. Environ. Res. Lett. 4, 044003. http://dx.doi.org/10.1088/
1748-9326/4/4/044003.
Schneider, A., Friedl, M.A., Potere, D., 2010. Mapping global urban areas using MODIS
500-m data: new methods and datasets based on “urban ecoregions”. Remote Sens.
Environ. 114, 1733–1746. http://dx.doi.org/10.1016/j.rse.2010.03.003.
Seto, K.C., Fragkias, M., Güneralp, B., Reilly, M.K., 2011. A meta-analysis of global urban
land expansion. PLoS One 6, e23777. http://dx.doi.org/10.1371/journal.pone.
0023777.
Sirmacek, B., Unsalan, C., 2010. Urban area detection using local feature points andspatial voting. IEEE Geosci. Remote Sens. Lett. 7, 146–150. http://dx.doi.org/10.
1109/LGRS.2009.2028744.
Small, C., Elvidge, C.D., 2013. Night on earth: mapping decadal changes of anthropogenicnight light in Asia. Int. J. Appl. Earth Obs. Geoinf. 22, 40–52. http://dx.doi.org/10.
1016/j.jag.2012.02.009.
Small, C., Pozzi, F., Elvidge, C.D., 2005. Spatial analysis of global urban extent from
DMSP-OLS night lights. Remote Sens. Environ. 96, 277–291. http://dx.doi.org/10.
1016/j.rse.2005.02.002.
Song, X.-P., Sexton, J.O., Huang, C., Channan, S., Townshend, J.R., 2016. Characterizing
the magnitude, timing and duration of urban growth from time series of Landsatbased estimates of impervious cover. Remote Sens. Environ. 175, 1–13. http://dx.doi.
org/10.1016/j.rse.2015.12.027.
Stevens, F.R., Gaughan, A.E., Linard, C., Tatem, A.J., 2015. Disaggregating census data for
population mapping using random forests with remotely-sensed and ancillary data.
PLoS One 10, e0107042. http://dx.doi.org/10.1371/journal.pone.0107042.
Su, Y., Chen, X., Wang, Chongyang, Zhang, H., Liao, J., Ye, Y., Wang, Changjian, 2015. A
new method for extracting built-up urban areas using DMSP-OLS nighttime stable
lights: a case study in the Pearl River Delta, southern China. GISci. Remote Sens. 52,
218–238. http://dx.doi.org/10.1080/15481603.2015.1007778.
Sutton, P.C., 2003. A scale-adjusted measure of “urban sprawl” using nighttime satellite
imagery. Remote Sens. Environ. 86, 353–369. http://dx.doi.org/10.1016/S00344257(03)00078-6.
Tatem, A.J., Noor, A.M., von Hagen, C., Gregorio, A.D., Hay, S.I., 2007. High resolution
population maps for low income nations: combining land cover and census in East
Africa. PLoS One 2, e1298. http://dx.doi.org/10.1371/journal.pone.0001298.
Taubenböck, H., Esch, T., Felbier, A., Wiesner, M., Roth, A., Dech, S., 2012. Monitoring
urbanization in mega cities from space. Remote Sens. Environ. 117, 162–176. http://dx.doi.org/10.1016/j.rse.2011.09.015.
Theobald, D.M., 2016. A general-purpose spatial survey design for collaborative scienceand monitoring of global environmental change: the global grid. Remote Sens. 8, 813.
http://dx.doi.org/10.3390/rs8100813.
Trianni, G., Lisini, G., Angiuli, E., Moreno, E.A., Dondi, P., Gaggia, A., Gamba, P., 2015.
Scaling up to national/regional urban extent mapping using Landsat data. IEEE J. Sel.
Top. Appl. Earth Obs. Remote Sens. 8, 3710–3719. http://dx.doi.org/10.1109/
JSTARS.2015.2398032.
Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring
vegetation. Remote Sens. Environ. 8, 127–150. http://dx.doi.org/10.1016/00344257(79)90013-0.
Turner, M.G., 1989. Landscape ecology: the eﬀect of pattern on process. Annu. Rev. Ecol.
Syst. 20, 171–197.
UN, 2014. World Urbanization Prospects: The 2014 Revision. https://esa.un.org/unpd/
wup/CD-ROM/, Accessed date: 17 February 2001.
US Census Bureau, 2012. Increasing Urbanization: Population Distribution by City Size,1790 to 1890. https://www.census.gov/dataviz/visualizations/005/, Accessed date:17 February 2001.
US Census Bureau, 2016. Urban and Rural Areas. https://www.census.gov/history/
www/programs/geography/urban_and_rural_areas.html, Accessed date: 17 February
2001.
Van de Voorde, T., Jacquet, W., Canters, F., 2011. Mapping form and function in urban
areas: an approach based on urban metrics and continuous impervious surface data.
Landsc. Urban Plan. 102, 143–155. http://dx.doi.org/10.1016/j.landurbplan.2011.

Hunsaker, C.T., O'Neill, R.V., Jackson, B.L., Timmins, S.P., Levine, D.A., Norton, D.J.,
1994. Sampling to characterize landscape pattern. Landsc. Ecol. 9, 207–226. http://dx.doi.org/10.1007/BF00134748.
Imhoﬀ, M., Lawrence, W.T., Stutzer, D.C., Elvidge, C.D., 1997. A technique for using
composite DMSP/OLS “City Lights” satellite data to map urban area. Remote Sens.
Environ. 61, 361–370. http://dx.doi.org/10.1016/S0034-4257(97)00046-1.
Jacob, A., Ban, Y., 2015. Sentinel-1A SAR Data for Global Urban Mapping: Preliminary
Results. 2015 IEEE International Geoscience and Remote Sensing Symposium(IGARSS). Presented at the 2015 IEEE International Geoscience and Remote Sensing
Symposium (IGARSS)pp. 1179–1182. http://dx.doi.org/10.1109/IGARSS.2015.
7325982.
Jean, N., Burke, M., Xie, M., Davis, W.M., Lobell, D.B., Ermon, S., 2016. Combining satellite imagery and machine learning to predict poverty. Science 353, 790–794.
http://dx.doi.org/10.1126/science.aaf7894.
Jing, W., Yang, Y., Yue, X., Zhao, X., 2015. Mapping urban areas with integration of
DMSP/OLS nighttime light and MODIS data using machine learning techniques.
Remote Sens. 7, 12419–12439. http://dx.doi.org/10.3390/rs70912419.
Johnson, B.A., Iizuka, K., 2016. Integrating OpenStreetMap crowdsourced data andLandsat time-series imagery for rapid land use/land cover (LULC) mapping: casestudy of the Laguna de Bay area of the Philippines. Appl. Geogr. 67, 140–149. http://dx.doi.org/10.1016/j.apgeog.2015.12.006.
Kasimu, A., Tateishi, R., Hoan, N.T., 2009. Global urban characterization using population density, DMSP and MODIS data. In: 2009 Joint Urban Remote Sensing Event.
Presented at the 2009 Joint Urban Remote Sensing Event, pp. 1–7. http://dx.doi.org/
10.1109/URS.2009.5137493.
Kawamura, M., Jayamana, S., Tsujiko, Y., 1996. Relation between social and environmental conditions in Colombo Sri Lanka and the urban index estimated by satellite
remote sensing data. Int. Arch. Photogramm. Remote. Sens. Spat. Inf. Sci. 31,
321–326.
Keola, S., Andersson, M., Hall, O., 2015. Monitoring economic development from space:using nighttime light and land cover data to measure economic growth. World Dev.
66, 322–334. http://dx.doi.org/10.1016/j.worlddev.2014.08.017.
Levin, N., Duke, Y., 2012. High spatial resolution night-time light images for demographic
and socio-economic studies. Remote Sens. Environ. 119, 1–10. http://dx.doi.org/10.
1016/j.rse.2011.12.005.
Li, C., Li, J., Wu, J., 2013. Quantifying the speed, growth modes, and landscape patternchanges of urbanization: a hierarchical patch dynamics approach. Landsc. Ecol. 28,
1875–1888. http://dx.doi.org/10.1007/s10980-013-9933-6.
Liu, H., Zhou, Q., 2004. Accuracy analysis of remote sensing change detection by rulebased rationality evaluation with post-classiﬁcation comparison. Int. J. Remote Sens.
25, 1037–1050. http://dx.doi.org/10.1080/0143116031000150004.
Liu, Z., He, C., Zhang, Q., Huang, Q., Yang, Y., 2012. Extracting the dynamics of urban
expansion in China using DMSP-OLS nighttime light data from 1992 to 2008. Landsc.
Urban Plan. 106, 62–72. http://dx.doi.org/10.1016/j.landurbplan.2012.02.013.
Liu, Y., Delahunty, T., Zhao, N., Cao, G., 2016. These lit areas are undeveloped: delimiting
China's urban extents from thresholded nighttime light imagery. Int. J. Appl. Earth
Obs. Geoinf. 50, 39–50. http://dx.doi.org/10.1016/j.jag.2016.02.011.
Lu, D., Tian, H., Zhou, G., Ge, H., 2008. Regional mapping of human settlements insoutheastern China with multisensor remotely sensed data. Remote Sens. Environ.
112, 3668–3679. http://dx.doi.org/10.1016/j.rse.2008.05.009.
Ma, L., Wu, J., Li, W., Peng, J., Liu, H., 2014. Evaluating saturation correction methods
for DMSP/OLS nighttime light data: a case study from China's cities. Remote Sens. 6,9853–9872. http://dx.doi.org/10.3390/rs6109853.
MacLachlan, A., Roberts, G., Biggs, E., Boruﬀ, B., 2017. Subpixel land-cover classiﬁcation
for improved urban area estimates using Landsat. Int. J. Remote Sens. 38,
5763–5792. http://dx.doi.org/10.1080/01431161.2017.1346403.
McFeeters, S.K., 1996. The use of the normalized diﬀerence water index (NDWI) in thedelineation of open water features. Int. J. Remote Sens. 17, 1425–1432.
Mertes, C.M., Schneider, A., Sulla-Menashe, D., Tatem, A.J., Tan, B., 2015. Detectingchange in urban areas at continental scales with MODIS data. Remote Sens. Environ.
158, 331–347. http://dx.doi.org/10.1016/j.rse.2014.09.023.
Mitraka, Z., Frate, F.D., Carbone, F., 2016. Nonlinear spectral unmixing of Landsat
imagery for urban surface cover mapping. IEEE J. Sel. Top. Appl. Earth Obs. Remote
Sens. 9, 3340–3350. http://dx.doi.org/10.1109/JSTARS.2016.2522181.
Miyazaki, H., Iwao, K., Shibasaki, R., 2011. Development of a new ground truth databasefor global urban area mapping from a gazetteer. Remote Sens. 3, 1177–1187. http://dx.doi.org/10.3390/rs3061177.
Myint, S.W., Gober, P., Brazel, A., Grossman-Clarke, S., Weng, Q., 2011. Per-pixel vs.
object-based classiﬁcation of urban land cover extraction using high spatial resolution imagery. Remote Sens. Environ. 115, 1145–1161. http://dx.doi.org/10.1016/j.
rse.2010.12.017.
Ng, H.-F., 2006. Automatic thresholding for defect detection. Pattern Recogn. Lett. 27,
1644–1649. http://dx.doi.org/10.1016/j.patrec.2006.03.009.
O'Neill, R.V., Hunsaker, C.T., Timmins, S.P., Jackson, B.L., Jones, K.B., Riitters, K.H.,
Wickham, J.D., 1996. Scale problems in reporting landscape pattern at the regional
scale. Landsc. Ecol. 11, 169–180. http://dx.doi.org/10.1007/BF02447515.
Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE Trans.
Syst. Man Cybern. 9, 62–66. http://dx.doi.org/10.1109/TSMC.1979.4310076.
Patel, N.N., Angiuli, E., Gamba, P., Gaughan, A., Lisini, G., Stevens, F.R., Tatem, A.J.,
Trianni, G., 2015. Multitemporal settlement and population mapping from Landsat
using Google Earth Engine. Int. J. Appl. Earth Obs. Geoinf. 35 (Part B), 199–208.
http://dx.doi.org/10.1016/j.jag.2014.09.005.
Pekel, J.-F., Cottam, A., Gorelick, N., Belward, A.S., 2016. High-resolution mapping of
global surface water and its long-term changes. Nature 540, 418–422. http://dx.doi.
org/10.1038/nature20584.
Pesaresi, M., Ehrlich, D., Ferri, S., Florczyk, A., Carneiro Freire, S.M., Halkia, S., Julea,274

Remote Sensing of Environment 205 (2018) 253–275

R. Goldblatt et al.

mapping urban areas from TM imagery. Int. J. Remote Sens. 24, 583–594. http://dx.
doi.org/10.1080/01431160304987.
Zhang, Q., Seto, K.C., 2011. Mapping urbanization dynamics at regional and global scales
using multi-temporal DMSP/OLS nighttime light data. Remote Sens. Environ. 115,2320–2329. http://dx.doi.org/10.1016/j.rse.2011.04.032.
Zhang, Q., Seto, K.C., 2013. Can night-time light data identify typologies of urbanization?
A global assessment of successes and failures. Remote Sens. 5, 3476–3494. http://dx.
doi.org/10.3390/rs5073476.
Zhang, H., Zhang, Y., Lin, H., 2012. Urban land cover mapping using random forest
combined with optical and SAR data. In: 2012 IEEE International Geoscience andRemote Sensing Symposium. Presented at the 2012 IEEE International Geoscienceand Remote Sensing Symposium, pp. 6809–6812. http://dx.doi.org/10.1109/
IGARSS.2012.6352600.
Zhang, Q., Schaaf, C., Seto, K.C., 2013. The vegetation adjusted NTL urban index: a new
approach to reduce saturation and increase variation in nighttime luminosity.
Remote Sens. Environ. 129, 32–41. http://dx.doi.org/10.1016/j.rse.2012.10.022.
Zhang, Q., Li, B., Thau, D., Moore, R., 2015. Building a better urban picture: combining
day and night remote sensing imagery. Remote Sens. 7, 11887–11913. http://dx.doi.
org/10.3390/rs70911887.
Zhao, S., Liu, S., Zhou, D., 2016. Prevalent vegetation growth enhancement in urban
environment. Proc. Natl. Acad. Sci. 113, 6313–6318. http://dx.doi.org/10.1073/
pnas.1602312113.
Zhou, Y., Smith, S.J., Elvidge, C.D., Zhao, K., Thomson, A., Imhoﬀ, M., 2014. A clusterbased method to map urban area from DMSP/OLS nightlights. Remote Sens. Environ.
147, 173–185. http://dx.doi.org/10.1016/j.rse.2014.03.004.
Zhou, Y., Smith, S.J., Zhao, K., Imhoﬀ, M., Thomson, A., Bond-Lamberty, B., Asrar, G.R.,
Zhang, X., He, C., Elvidge, C.D., 2015. A global map of urban extent from nightlights.
Environ. Res. Lett. 10, 054011. http://dx.doi.org/10.1088/1748-9326/10/5/
054011.

03.017.
Wan, B., Guo, Q., Fang, F., Su, Y., Wang, R., 2015. Mapping US urban extents from MODIS
data using one-class classiﬁcation method. Remote Sens. 7, 10143–10163. http://dx.
doi.org/10.3390/rs70810143.
Wei, Y., Liu, H., Song, W., Yu, B., Xiu, C., 2014. Normalization of time series DMSP-OLS
nighttime light images for urban growth analysis with pseudo invariant features.
Landsc. Urban Plan. 128, 1–13. http://dx.doi.org/10.1016/j.landurbplan.2014.04.
015.
Wood, S., Har-Noy, S., 2013. A Better World with Enhanced Analytics Capabilities. 16.
GeoInformatics Emmeloord, pp. 6–8.
Woodcock, C.E., Allen, R., Anderson, M., Belward, A., Bindschadler, R., Cohen, W., Gao,F., Goward, S.N., Helder, D., Helmer, E., Nemani, R., Oreopoulos, L., Schott, J.,
Thenkabail, P.S., Vermote, E.F., Vogelmann, J., Wulder, M.A., Wynne, R., 2008. Free
access to Landsat imagery. Science 320, 1011. http://dx.doi.org/10.1126/science.
320.5879.1011a.
Xiao, P., Wang, X., Feng, X., Zhang, X., Yang, Y., 2014. Detecting China's urban expansion
over the past three decades using nighttime light data. IEEE J. Sel. Top. Appl. Earth
Obs. Remote Sens. 7, 4095–4106. http://dx.doi.org/10.1109/JSTARS.2014.
2302855.
Xu, H., 2006. Modiﬁcation of normalised diﬀerence water index (NDWI) to enhance open
water features in remotely sensed imagery. Int. J. Remote Sens. 27, 3025–3033.
http://dx.doi.org/10.1080/01431160600589179.
Xu, H., 2008. A new index for delineating built-up land features in satellite imagery. Int.
J. Remote Sens. 29, 4269–4276. http://dx.doi.org/10.1080/01431160802039957.
Yuan, F., Bauer, M.E., 2007. Comparison of impervious surface area and normalized
diﬀerence vegetation index as indicators of surface urban heat island eﬀects in
Landsat imagery. Remote Sens. Environ. 106, 375–386. http://dx.doi.org/10.1016/j.
rse.2006.09.003.
Zha, Y., Gao, J., Ni, S., 2003. Use of normalized diﬀerence built-up index in automatically275