A global map of urban extent from nightlights

The major data used in this study include the DMSP/OLS NTL, high spatial resolution land cover, water mask, and gas flare mask. The DMSP/OLS NTL measures lights on the Earth surface from cities and settlements with persistent lighting, and other sources such as gas flares, fires, and illuminated marine vessels (Zhang et al 2013). The data were recorded as digital number (DN) from 0 to 63 with a 1 km spatial resolution, spanning −180° to +180° in longitude and −65° to +75° in latitude. Annual cloud-free composites were built using the highest-quality data based on a number of constraints (Elvidge et al 2009). The NTL data in the year 2000 were chosen to build the global product of urban areas for validation and comparison with other similar products. The water mask was derived from MODIS 250 m land–water mask (MOD44W), and gas flare data were obtained from the NOAA national geographic data center (Elvidge et al 2009).

High-resolution land cover data for the US in 2006 and China in 2005 were obtained from US Geological Survey National Land Cover Dataset (NLCD) and Resources and Environment Data Center of the Chinese Academy of Science, respectively, both with an original spatial resolution of 30 m (Homer et al 2004, Liu et al 2010). Land cover types mainly include open water, urban, evergreen forest, deciduous forest, shrub, grassland, cropland and wetland. The land-cover data for China were built through visual interpretation of Landsat TM images and aggregated to a 1 km scale map of each land cover type (Liu et al 2010). The US land-cover data layer was also aggregated from 30 m to 1 km spatial resolution, resulting in an urban percentage map.

A unified definition of 'urban' has not been reached, and thus the estimates of global urban land vary widely with different definitions in previous studies (Liu et al 2014). In this study, similar to Zhou et al (2014), we define urban land as those 1 km pixels with urban percentages larger than 20%, consistent with the developed areas in the US NLCD dataset (Fry et al 2011).

2.2.1. Cluster-based method

The threshold techniques were developed to delineate urban extent based on the DMSP/OLS NTL observation to account for the bias of light features in NTL being substantially larger than the lighting sources on the ground. The threshold is defined as the DMSP/OLS DN value above which the pixel is classified as urban area. Moreover, urban pixels tend to group together as clusters. As such, our cluster-based method generally follows this logic by segmenting the NTL into clusters and delineates urban extent in each cluster based on the NTL DN values and a cluster-specific threshold.

Our method for mapping global urban area is an extension of the cluster-based method originally developed in Zhou et al (2014). This revised method includes five major steps: data preprocessing, urban cluster segmentation, parameterization for the logistic model, threshold estimation, and urban extent delineation (figure 1). First, as a pre-processing step, water and gas flare pixels are excluded from the NTL data. Second, using a marker-controlled watershed segmentation algorithm (Zhou et al 2014), the processed NTL data were segmented into potential urban clusters, each consisting of groups of similar and spatially continuous pixels. Not all pixels in each potential cluster are urban. More details about segmentation can be found in Zhou et al (2014).

Zoom In Zoom Out Reset image size

Third, an optimal threshold for each cluster is then estimated to delineate the actual urban area from the NTL data in each potential urban cluster. To estimate the optimal threshold, we developed a logistic model with two parameters a/b and β', as formulated in equations (1) and (2). This model ingests cluster-level NTL information and gives an estimate of optimal threshold for delineating urban extent in each cluster. Building on the work by Zhou et al (2014), we developed a parameterization method to estimate a/b and β' and calculated these two parameters for each country (see next section 2.2.2). Here our parameterization method was developed from the NTL data in years 2006 and 2005 for the US and China, respectively, in order to be temporally consistent with the available high spatial resolution land cover dataset. The optimal threshold is calculated as follows

$$\begin{eqnarray}&&\begin{array}{lllllllllllllll} {\rm NT}{{{\rm L}}_{{\rm thld}}}=\frac{1}{1+{{{\rm e}}^{-\beta ^{\prime} (x^{\prime} -x{{^{\prime} }_{{\rm mean}}})}}} \\ \quad \quad \quad \quad \times \left( {\rm NT}{{{\rm L}}_{{\rm max} }}-{\rm NT}{{{\rm L}}_{{\rm min} }} \right)+{\rm NT}{{{\rm L}}_{{\rm min} }}, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&x^{\prime} ={\rm ln} {\mkern 1mu} \left( {{S}^{\frac{a}{b}}}\cdot {\rm NT}{{{\rm L}}_{{\rm mean}}} \right),\end{eqnarray} \tag{ 2 }$$

where x' is the index for estimating the optimal threshold based on combined effects of both mean lighting magnitude and cluster size, S is the cluster size, NTL_mean is the mean NTL DN in each cluster, and a/b is the parameter for the relative contribution of each of two effects. NTL_thld is the optimal threshold to delineate the urban area in the potential cluster, and x'_mean is the mean value of x'. NTL_min and NTL_max are minimum and maximum NTL DN in the study area, and β' is the parameter in the logistic model.

Fourth, we estimated the optimal threshold value for each cluster identified from the year 2000 global NTL data at the global level using the logistic model and the derived parameters in step 3. Finally, we mapped the urban extent according to the estimated optimal threshold in each cluster. We chose the NTL data in the year 2000 to build the global urban area product because of the availability of other datasets such as MODIS urban product close to this time period for validation and inter-comparison.

2.2.2. Parameterization of the logistic model

In this study, we developed a parameterization method to estimate the key parameters (a/b and β') in the logistic model using the NTL data only. As reported by Zhou et al (2014), the parameters of a/b and β' in the logistic model show slight difference in the US and China at the country level, although the optimal thresholds are not highly sensitive to these parameters. We calculated these two parameters at the regional level by dividing the US into nine Census regions, and dividing China into three economic regions (figure 2). The optimal thresholds derived from high spatial resolution reference data for the US in 2006 and for China in 2005 were used to estimate these two parameters in the logistic model in each region according to the method by Zhou et al (2014). We performed this analysis at the regional level for several reasons. First, we can collect a number of samples to evaluate the variations of these two parameters. Second, because of differences in factors such as socioeconomic development and geography, these regions can cover a range of urbanization. Third, it can reduce the possibility of separation of potential urban clusters across different regions. The possibility of shared clusters across boundaries is high at finer spatial units, such as province or state.

Zoom In Zoom Out Reset image size

In order to extend the logistic model in other regions without high spatial resolution data, we analyzed the two parameters (a/b and β') in the logistic model for all regions in the US and China. We found that the parameter a/b has a narrow distribution around 0.23 with a standard deviation of 0.01 (figure 3(a)), and that the parameter β' is highly correlated to regional nightlight mean magnitude for each region (figure 3(b)).

Zoom In Zoom Out Reset image size

Based on this analysis, we took the mean value of 0.23 for a/b, for the logistic model in each country. The parameter β' for each country was derived from the regression equation between this parameter and regional nightlight mean magnitude (figure 3(b)). It is worth noting that the parameterization scheme could also be applied at a sub-national level to further improve urban mapping in large countries, particularly for regional scale urban studies. Although this paper is focused on global urban mapping, we evaluated the sensitivity of the optimal thresholds to these two parameters to evaluate the reliability of this parameterization scheme (see the results and discussion section below).

The newly developed 1 km urban product was aggregated to 1/12° as a percentage map (figure 4) to be comparable with other previously published products. Based on this new map, overall urbanization at the global level is about 0.54%, which is close to 0.49% derived from MODIS observations (Schneider et al 2009). Urban areas show large spatial heterogeneity, globally. The global distribution of urban area as depicted in figure 4 appears realistic, and it is not surprising that Europe shows the largest urbanized fraction (2.3%) as compared with other regions, with North America next at 1.2%. The urbanized land in Asia, Latin America, and Former Soviet Union (FSU), and Oceania are 0.7%, 0.4%, 0.2%, and 0.1%, respectively. The urbanized land in Africa is less than 0.1%, the lowest among all regions.

Zoom In Zoom Out Reset image size

At the country level, urban land area varies from lower than 0.01% to higher than 10%. Urban land area in most countries (>70%) is less than 1%, and more than 90% of countries have urban land area of less than 5% (figure 5). In some regions such as FSU, urban land area in all of the countries within the region covers less than 1% of total land area. Urban land area in large countries such as Russia, with a large amount of areas in high latitude zones, is close to 0.1%. The countries with relatively high urban land area are mostly in Europe. Among the large countries, the US has urban land area of more than 1%.

Zoom In Zoom Out Reset image size

Major urban clusters are further uncovered by the longitudinal and latitudinal distributions of both urban area and percentage (figure 6). The latitudinal zones around 30° N have the largest urban area due to contributions from China and the US, while the latitudinal zones with the highest urbanized percentage move North, around 45° N, largely due to the influence of European countries. The longitudinal distribution of urban area shows several large centers. One of the most evident ones is around 80° W, encompassing the 'Boston–New York–Washington' corridor of the US Northeast region, and this longitudinal distribution of urbanized percentage is consistent with that of total urban land area. This type of analysis is useful for developing adaptation and risk management measures for urban infrastructure, transportation, energy, and water systems when considered together with other factors such as climate variability and change, and high impact weather events.

Zoom In Zoom Out Reset image size

Two sensitivity analyses were performed to examine the impact of the parameters of a/b and β' on the selection of the optimal thresholds (figures 7 and 8). In the first analysis, we evaluated the R² and root mean squared error (RMSE) between the optimal thresholds derived from high resolution land cover data and our logistic model by varying the parameters (a/b and β') in a range of possible values (figure 7). The results indicate that the choice of these parameters have small impact on the R² when a/b is lower than 0.3, and that R² values begin to decrease when a/b is larger than 0.3, and the RMSE reaches its minimum value when a/b falls in the range of 0.2–0.25, for most regions. More important, the RMSE is as low as 4 when a/b is within a range of 0.15–0.3, for all regions. The impact of parameter β' on the R² is also minor, and it reaches its maximum R² value around 0.9 in most regions when β' falls in the range of 0.75–1.0. The RMSE is as low as 3 when β'is within the same range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first sensitivity analysis indicates that the sensitivity of the optimal thresholds to the two parameters (a/b and β') is low when these parameters are sufficiently close to their optimal values. Then, in the second analysis, we examined the difference between optimal thresholds from high-resolution land cover data and those using the logistic model with parameters of a/b = 0.23 and β' derived from the regression equation for each region (figure 8). The results indicate that the parameterization method in the logistic model based on a regression equation for the parameter β' together with a mean value (0.23) for the parameter a/b performs well in determining the optimal thresholds for all study regions.

We compared the global map of urban extent from NTL to five previously published and widely used global urban area products approximately for the same period. These include (1) MODIS 500 m urban map (2001–2002); (2) GlobCover (2005); (3) GLC2000 (2000); (4) global rural–urban mapping project (GRUMP) (1995); and (5) NOAA's impervious surface area map (IMPSA) (2000–2001), at the pixel level for two selected sites (figure 5), and regionally (figure 9). The derived urban area from MODIS data is defined as built environment including all non-vegetative and human-constructed facilities that cover greater than 50% of a given landscape unit (Schneider et al 2009). The urban area estimates in GlobCover and GLC2000 are defined as artificial surfaces and associated areas (Bartholomé and Belward 2005, Bicheron et al 2008). The urban area in GRUMP is defined as urban extents with the total population greater than 5000 persons (CIESIN 2011). The urban areas in IMPSA is based on the density of impervious surface area (Elvidge et al 2007b).

Zoom In Zoom Out Reset image size

For the first site in Japan with highly urbanized land, our newly developed product shows the greatest similarity in spatial pattern and urbanization magnitude with the MODIS product. The extent of urban areas in GlobCover and GLC2000 is much smaller compared to other datasets, while the urban extent in GRUMP is the largest among the six products. The urban extent in IMPSA is similar to our estimate and the MODIS product, but its magnitude is smaller. For the second site in Brazil with medium urbanized land, the urban extent and magnitude in our data, the MODIS product, and IMPSA are similar, while the urban extent is still smaller in GlobCover and GLC2000, and larger in GRUMP compared to the other products.

Similar to the comparison at the two selected sites, total urban land area in our newly developed product is closest to the MODIS product at both global and regional levels (figure 9). Urban areas in GlobCover and GLC2000 are generally lower than all other products, and urban area in GRUMP is the highest at the global and regional levels. Urban area in IMPSA is generally lower than our estimation and those from MODIS product, except for Asia. There are several factors, such as definition of urban extent and data used, contributing to the difference between products. Different from these products, often limited by the temporal coverage of used geospatial data, the proposed method in this study provides the possibility to build temporally updateable urban area maps using the long time series of DMSP/OLS NTL data.

We evaluated our urban area product at the pixel level using the five global urban products discussed above as reference data (figure 10). A kappa coefficient of 0.41–0.60 indicates moderate agreement, and a value of 0.61–0.80 indicates substantial agreement (Viera and Garrett 2005). Based on this analysis, we find our newly developed urban area product to be in substantial agreement with the IMPSA product for most regions, globally. One reason for this high degree of agreement is that both products use the same input data, including NTL and NLCD, in mapping urban areas. As the IMPSA method requires auxiliary data such as gridded population, it needs more effort to develop a global product. The requirement for auxiliary data also limits its capability in mapping temporal urban dynamics. It would be useful in future work to examine if the methodologies developed here could be extended to also produce information on impervious surface area, given its importance for applications such as hydrology modeling.

Zoom In Zoom Out Reset image size

The agreement for North America is also high between our product and all other products except for GlobCover. For all regions, our product is in moderate agreement with MODIS product, and especially in substantial agreement for North America, South America, and Oceania. At the global level, our product is in almost substantial agreement with MODIS product. Our newly derived urban area map achieves an overall high comparability with other five products for all regions. According to the measurement of omission error, our product varies the most compared with the GRUMP product, while the omission error is the lowest compared with IMPSA product. The omission error for urban area in our product is lowest in North America, as compared with all five products. According to the measurement of commission error, our product varies the most compared with the GlobCover product, while the commission error is lowest compared with GRUMP product. The commission error for urban area in our product is lowest in Oceania region, compared with all five products.

In addition to the evaluation against other products at the global level, we chose Europe, a region not used in our model development, to assess the accuracy of the derived urban extent using the 100 m spatial resolution Corine land cover 2000 data set (Büttner et al 2004). While our global product is targeted for large-scale urban area mapping and monitoring, and it is also limited by the spatial resolution of NTL data, we believe that the resulting pixel level accuracy (overall accuracy of 87%, Kappa coeeficient of 0.55, producer accuracy of 96%, and user accuracy of 88%) for Europe is reasonably good.