Global Distribution and Density of Constructed Impervious Surfaces

We present the first global inventory of the spatial distribution and density of constructed impervious surface area (ISA). Examples of ISA include roads, parking lots, buildings, driveways, sidewalks and other manmade surfaces. While high spatial resolution is required to observe these features, the new product reports the estimated density of ISA on a one-km2 grid based on two coarse resolution indicators of ISA – the brightness of satellite observed nighttime lights and population count. The model was calibrated using 30-meter resolution ISA of the USA from the U.S. Geological Survey. Nominally the product is for the years 2000-01 since both the nighttime lights and reference data are from those two years. We found that 1.05% of the United States land area is impervious surface (83,337 km2) and 0.43 % of the world's land surface (579,703 km2) is constructed impervious surface. China has more ISA than any other country (87,182 km2), but has only 67 m2 of ISA per person, compared to 297 m2 per person in the USA. The distribution of ISA in the world's primary drainage basins indicates that watersheds damaged by ISA are primarily concentrated in the USA, Europe, Japan, China and India. The authors believe the next step for improving the product is to include reference ISA data from many more areas around the world.


Introduction
Human beings around the world build, use and maintain constructed impervious surfaces for shelter, transportation and commerce. It is a universal phenomenon -akin to clothing -and represents one of the primary anthropogenic modifications of the environment. Expansion in population numbers and economies combined with the popular use of automobiles has lead to the sprawl of development and a wide proliferation of constructed impervious surfaces. The quantity of development that has occurred in recent decades is staggering. In the USA there are a million new homes and 16,000 kilometers of paved road built each year. The worldwide pattern of sprawl development will continue in the coming decades in response to both population growth and growth in living standards.
Constructed impervious surfaces are both hydrological and ecological disturbances. However, constructed surfaces are different from most other types of disturbances in that recovery is arrested through the use of materials that are resistant to decay and are actively maintained. The same characteristics that make impervious surfaces ideal for use in construction produce a series of effects on the environment [1]. Impervious surfaces alter sensible and latent heat fluxes, causing urban heat islands [2]. In heavily vegetated areas, the proliferation of ISA reduces the sequestration of carbon from the atmosphere [3]. ISA alters the character of watersheds by increasing the frequency and magnitude of surface runoff pulses [4]. The increased overland flow alters the shape of stream channels, raising water temperatures, and sweeping urban pollutants into aquatic environments [5,6]. Hydrologic consequences of ISA include increased flooding, reductions in ground water recharge, and reductions in surface water quality. A widely accepted scale for the impacts of ISA on holds that watershed areas are stressed if they contain 1-10% ISA, impacted if they contain 10%-25% ISA and are degraded if the contain more than 25% ISA [1,7].
Spatial grids depicting the density of impervious surface area (ISA) were developed to address applications requiring specific measures or estimates of constructed areas. There are three basic remote sensing approaches to estimating ISA [8,9]. The first approach is to map constructed areas using high spatial resolution imagery [10,11,12]. Typically ISA products derived from high spatial resolution imagery (in the +/-one meter range) cover small areas and to date there has not been a standardization of methods which would facilitate the merger of products generated by diverse organizations.
The second approach is to use moderate spatial resolution multispectral data (such as Landsat) to estimate the density of ISA. The U.S. Geological Survey (USGS) recently released a 30-meter ISA density grid for the USA derived from Landsat as part of the National Land Cover Database (NLCD-2001) [11]. The USGS used a combination of spectral and spatial methods to estimate the density of ISA. They used a sub-sample of high spatial resolution imagery to establish the methodology and to conduct an accuracy assessment. Noted sources of error for the product include the underestimation of ISA in areas having extensive tree cover and overestimation of ISA in areas with extensive bare ground. The advantage of extracting ISA from Landsat style data is that is possible to cover an entire region or country using a standardized methodology.
The third approach is to use indicators to estimate the density of ISA. One such indicator is population numbers [13] based on data collected by government agencies. The concept behind this is that where there are people -surely there is ISA. One can imagine using gridded population data to make gridded ISA data. However, not all population grids are well suited for this application. In particular, if the population tally is based on residency (where people sleep), the values are quite low in heavily built up commercial centers and transportation areas. Ambient Population products are likely better for this application since they estimate the average population over a typical diurnal cycle. Another indirect method is the estimation of ISA based on coverage coefficients developed for standard land cover classes, such as low density residential, high density residential, commercial / industrial [14].
In 2004 NOAA produced a one km 2 grid of ISA densities for the conterminous USA using indirect data sources to estimate ISA [15]. The input data included radiance calibrated nighttime lights, street and road density (from the U.S. Census Bureau) and three Landsat derived urban land cover classes from the early 1990's [16]. The ISA estimates were calibrated using gridded point counts made of ISA from high spatial resolution aerial photography selected along transects crossing thirteen major urban centers. Using these data sources it was estimated that as of the year 2000-01 the USA held 112,000 square kilometers of constructed impervious surfaces, an area nearly the size of the State of Ohio.
As an extension of the methods used by NOAA in 2004 [15] we have produced the first global grid of ISA density. We have used the product to estimate the world's total ISA, to rank the leading countries in total ISA and to calculate the quantity of ISA per person for individual countries. In addition, we have aggregated the ISA density for the major watershed units of the world to identify those watersheds impacted by the proliferation of ISA.

Nighttime Lights
A radiance calibrated nighttime lights product was assembled using data acquired by the U.S. Air Force Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS). The Designed as a cloud imager, the DMSP OLS collects global data in two broad spectral bands -visible and thermal. The DMSP satellites fly in polar orbits and each collects fourteen orbits per day. With a 3000 km swath width, each OLS is capable of collecting a complete set of images of the earth twice a day. At night, a photomultiplier tube (PMT) intensifies the visible band signal to enable the detection of moonlit clouds. The boost in gain enables the detection of lights present at the earth's surface. Most of the lights are from human settlements (cities and towns) and fires, which are ephemeral. Lights from gas flares are also detected and can be readily identified offshore or in isolated areas not impacted by urban lighting.
NGDC serves as the long-term archive for DMSP data and has data holding extending from 1992 to the present. In the standard collection mode the OLS gain is turned up quite high for the detection of moonlit clouds. Under these operating conditions the nighttime visible data are saturated (DN=63) in urban centers. In addition, during the operational collections the system gain can vary, modified by on-board algorithms. It is possible to turn the gain down to avoid saturation in urban centers, but then the dimmer lights are not detected due to the limited dynamic range of the sensor. To overcome these limitations, NGDC requested the data collection in three overlapping fixed gain settings (low, medium and high) on alternating 24-hour periods. The low gain collections provide unsaturated data of bright urban cores. The high gain setting data provide detection of dim lights -but has saturation in the bright urban cores. The medium gain setting covers the gap between the detections achieved with the high and low gain settings. The digital numbers from the low and medium gain settings were converted to the radiance units of the high gain data based on the sensor's preflight calibration. Each orbit was processed with automatic algorithms that identify image features (such as lights and clouds) and the quality of the nighttime data. The basic algorithms have been described in references 17,18,19. The following criteria identified the best nighttime lights data from the set of fixed gain collections: 1. Center half of orbital swath (best geolocation and sharpest features). The unsaturated nighttime data from each individual orbit meeting the above criteria were averaged in a 30 arc second grid (Platte Carree projection) using the radiance increment for a single DN from the highest gain setting of the collections (1.35 x 10 -10 watts/cm 2 /sr). Masking was used to remove the lights from natural gas flares present on land areas in places such as Nigeria and Russia [20].

LandScan 2004
The U.S. Department of Energy, Oak Ridge National Laboratory has produced an evolving series of spatially disaggregated global population count data sets, known as LandScan. The basic concept of the LandScan data sets is to perform a spatial allocation of census reported population numbers based on models developed with spatially disaggregated data. The first LandScan product [21] used DMSP nighttime lights for the mapping of human settlements. However, the nighttime lights were subsequently dropped [22] due to the overt affect of economic development on the extent and brightness of lighting. We used the LandScan 2004 product, which included input from three satellite data sources: NASA MODIS land cover [23], the topographic data from the Shuttle Radar Topography Mission (SRTM) [24], and the high resolution Controlled Image Base (CIB) from the U.S. National Geospatial Intelligence Agency (NGA). The term population count is used instead of population density -which is based on residence. On a population density grid commercial centers and airports have very low numbers, despite the fact that there are substantial numbers of people present during certain hours. Population count products, also referred to as ambient population, attempt to represent the spatial distribution of population based on person hours. The Landscan 2004 product has substantial ambient population numbers in commercial centers and even along some major transportation routes, such as interstate highways. The Landscan 2004 product assigns population values of "zero" to airport lands.

ISA Estimation Model
A model was developed to estimate the density of ISA based on the radiance calibrated nighttime lights and the Landscan population count. The model was developed using the 30-meter USGS Landsat derived ISA (from the NLCD-2001) as the reference data source. The 30-meter ISA was aggregated to a one-kilometer equal area grid in an Albers projection ( Figure 1). The Landscan and radiance calibrated nighttime lights were reprojected to the same one kilometer grid (Figures 2 and 3). Linear regression defined an equation for estimating the density of ISA. Only grid cells with population count values of three or greater were included in the regression. This excluded airport lands from the regression. Also excluded from the regression were outliers having extremely high population counts (greater than 3000 persons / km 2 ) and extremely bright lights (digital numbers greater than 800). The regression included 470,894 grid cells. The equation developed through the regression is as follows:   The regression had an r 2 squared of 0.59 and is highly significant p < .0001. Figure 4 shows the resulting ISA grid of the USA.

Global ISA Grid
The initial global ISA density grid was produced on a 30 arc second grid since that is the native grid of both the Landscan and nighttime lights. This was then converted to a 1 km equal area grid in a Mollweide projection ( Figure 5). A threshold of 0.4 percent was applied to eliminate the salt and pepper noise present at the very low end of the ISA scale. ISA values over 100 were reset to 100. Extractions of the digital values were run to tally the quantity of ISA for countries, sub-national units (states / provinces) and major watersheds. We used the 30 arc second resolution Hydro-1K drainage basins of the world [25] as a representation of the major watersheds. We analyzed the ISA cover present in the major watersheds of the world to tally the percent area of each watershed that was pristine (zero ISA), stressed (1-10% ISA), impacted (10-25%), and damaged (more than 25% ISA).

Results and Discussion
In order to examine the qualities of the product we generated a scattergram of the global ISA values against the USGS ISA aggregated to 1 km resolution ( Figure 6). The scattergram shows that there is a substantial amount of scatter in the NGDC estimates of ISA relative to the reference USGS estimates. However, this scatter diminishes with aggregation when ISA is considered at the sub-national (e.g. state) level. Figure 7 shows the scattergram of the ISA estimates for the forty-eight (plus D.C.) states of the conterminous USA. The USGS ISA estimate is that 90,037 km 2 of the conterminous U.S. is impervious; whereas the NGDC estimate is slightly lower at 83,338 km 2 . The difference between the two estimates is 6,699 km 2 . Both of these ISA estimates for the conterminous USA are lower than the previous estimate from NOAA [14] of 112,000 km 2 , which used gridded point counts of ISA made on high resolution aerial photography as reference data. It is likely that the there are systematic differences in the reference ISA methods used in the previous NOAA [15] and USGS [11] products.
In comparing the global (NGDC) ISA versus the USGS ISA ( Figure 8) the most profoundly overpredicted states in the NGDC product were small and highly urbanized (The District of Columbia, Delaware, New Jersey, and Maryland); whereas the most significantly under-predicted state was large and rural (Montana). The over-predictions occurred in heavily treed regions with widespread detection of dim lighting outside of urban centers in states such as New Jersey and Maryland. The underestimations are concentrated in rural areas in the western two-thirds of the country. Some of the underestimation in the NGDC product is associated with the detection of roads in the Landsat data. In other cases it appears that the USGS product overestimated ISA in areas having extensive bare soil. For the conterminous USA the aggregate area of ISA in the USGS product that had zero ISA in the NGDC product is 10,842 km 2 -more than enough to account for the 6,699 km 2 difference in the total ISA from the USGS and NGDC products. This implies that on an aggregate level the two products are largely equivalent outside of the most sparsely populated areas where no lighting was detected.
The hundred leading countries in terms of total ISA are shown in Table 1. Also listed is the quantity of ISA per person in square meters. At the bottom of the table we list the total ISA for all countries and the average amount of ISA found per person worldwide. Clearly, the countries that measure high on total ISA are either large in area extent and/or total population, or have high levels of economic development.    The countries with particularly high values of ISA per person according to our estimation are almost universally affluent (United States, Canada, Norway, Sweden, Finland, Spain, France, Bahrain, Brunei, Qatar, and the United Arab Emirates). With the exception of Brunei, these countries cluster in the northern hemisphere ( Figure 9). It is interesting to note that Japan and Mexico had 114 m 2 of ISA per person). Japan's moderate level of ISA (relative to their GDP per capita) can be attributed to the topographic and agricultural constraints on development present in that country. Clearly, ISA is a function of population, level of economic development, and the availability of surfaces suitable for building. High levels of ISA per capita calculated at nationally aggregate levels generally identify wealthy countries (e.g. high GDP per capita). However, exploration of ISA as a percentage of watershed surface area identifies watersheds that are likely to be suffering detrimental hydrologic and ecologic effects. Figure 10 shows the percentage of watershed areas that are in the stressed category -with 1-10% ISA. Pristine watersheds having little or no ISA (yellow) are concentrated at northern latitudes, central Asia, portions of Africa, the Amazon Basin and the southern tips of South America and the Arabian Peninsula. Watersheds with more than half of their areas stressed by ISA are prevalent in Japan, China, India, Europe and the eastern half of the USA. Figure  11 shows the percentage of watershed areas falling in the impacted category (10-25% ISA cover).
There are very few watersheds having a substantial percentage of their area in the impacted range. Watersheds having 3% or more of their area in the impacted range are concentrated in Japan, China, Java, India, Europe and the eastern USA. A similar situation was found with the percentage of watershed areas in the damaged category (more than 25% ISA cover), as shown in Figure 12.   To summarize the results on watersheds, Figure 13 shows the percentage of watersheds having half or more of their area in the stressed category. South America and Africa have a very low percentage (3-4%) of their watersheds with half or more of their area in the stressed category. North America (10%) and Asia (18%) are slightly higher. It is clear that at the continental scale the vast majority of watersheds still have very low levels of ISA. At the scale of individual countries, the story is different. In the conterminous USA 30% of the watersheds have half their area or more in the stressed category. On the same measure, Europe has 40% and Japan, China and India (as a group) have 45%. The results suggest that the adverse effects of ISA are concentrated into a small segment of the world's watersheds but are more pervasive in Japan, China, India, Europe and the USA than in other landmasses. Figure 13. Percentage of watersheds having 50% or more of their surface in the stressed category (1-10% constructed impervious surface area).

Conclusion
We have produced the first global grid of ISA at a resolution of one kilometer. The GIS readable product and related tables are freely available at: http://www.ngdc.noaa.gov/dmsp/download.html. The total ISA of the world is estimated to be 579,703 km 2 . This is nearly the same size as the country of Kenya (584,659 km 2 ) and larger than Spain (505,735 km 2 ) or France (546,962 km 2 ). The country with the most ISA is China (87,182 km 2 ) followed closely by the United States (83,881 km 2 ), and India (81,221 km 2 ). China and India's ISA footprints are population driven whereas the United States ISA footprint is more driven by affluence. Explorations of ISA per capita show generally expected patterns in that countries with high population densities (e.g. 'big denominators') show lower levels of ISA per capita. The global average of ISA per capita was estimated to be 93 m 2 per person. Examinations of ISA at the watershed level support ideas that there are both economic and demographic forces contributing to changes in the hydrologic and ecologic functioning of watersheds around the world.
The estimate of ISA is derived solely from the brightness of satellite observed nighttime lights and population count. Both of the input sources (nighttime lights and population count) are produced as 30 arc second grids (~1 km 2 resolution) and could potentially be updated on an annual basis. These two data sources are complementary in that the nighttime lights are generally brightest in the commercial and industrial areas -which are generally not well defined in the population count data. In areas with no detected lighting -the ISA estimate is based solely on population count. Based on equation 1, in the absence of detected lighting the population count at which 100% ISA is reached is slightly over 11,000 persons per square kilometer. At this density each person is directly associated with 91 m 2 of ISA. It should be noted that none of the 100% ISA values in the product were produced by population alone since all areas having population counts of 5000+ had lighting detected.
There are likely systematic errors in that the product calibration was conducted with USA data only. It is well established [26,27] that most developing countries have less lighting at night than countries such as the USA. This results in a lower ISA estimation, which is likely the case since most developing countries have less ISA per person than the USA. In addition, comparison of this product to USGS estimates of ISA in the conterminous USA with NOAA's previous estimate [15] suggests that the new product may be biased towards underestimation globally. The best way to improve the product would be to include ISA reference data from a broad suite of countries -instead of the USA only.
In recent years, porous asphalt and concrete products have been developed [28]. These materials substantially reduce the overland flow of water, reducing the hydrologic impacts of ISA. To date however, the spatial extent of porous roadways and parking lots is quite small.
As the world economy and population expands, the quantity of constructed surfaces of the earth will expand significantly. For the moment, this product stands as the only global ISA grid. We offer the current product as a pathfinder, recognizing that the mapping of constructed area density at both the global and local scale will continue to improve over time.