Country level urban land and its variability
Large variabilities in the degree of present-day urban land (for the year 2019 or 2020; see Methods) are seen across countries (Fig. 1a) based on eight global datasets, with Vatican City and Singapore showing the highest values (eight-product mean urban percentages of 79.87% and 53.7%, respectively). Ignoring uninhabited territories, on the low end, there are several overseas island territories, such as Cocos, Midway, and Pitcairn Islands, with negligible urban land detected by these satellite-derived products. Overall, China is the country with the most urban land (264403 sq. km covering 2.82% of its total area based on the eight-product mean), followed by the United States (183735 sq. km; 1.94%), India (85760 sq. km; 2.77%), and Russia (59311 sq. km; 0.35%) (Fig. 1b). The present-day global urban land percentage varies between 0.52% in the World Settlement Footprint (WSF) 2019 dataset22 to a four times higher estimate of 2.07% according to the Esri land Cover product (for the year 2020; 1.93% for 2019)23. The eight-product mean global urban percentage is ~ 0.95%.
Since countries have different baseline levels of urbanization (Fig. 1a), we calculate the coefficient of variation or Normalized Root-Mean-Square Deviation (standard deviation across products divided by their mean expressed as a percentage) to standardize the degree of disagreement between these eight data products (Fig. 1c). Larger disagreements are seen for Greenland, countries in East Africa (Ethiopia, Kenya, Uganda, and Tanzania), Russia, countries in south Asia (Afghanistan, Pakistan, India, and Myanmar), Paraguay, etc. Better agreements between datasets are seen for Brazil, Argentina, Japan, most countries in western Europe, parts of Central and South Africa, Canada, etc.
Regional- to city-scale differences across datasets
We also compare the present-day estimates of urban land for four distinct regions in the world: the Great Lakes and Mid Atlantic regions in North America, and the Indo-Gangetic and Yangtze River Basins in Asia (Fig. 2a). While all four of these regions are heavily urbanized, in the last few decades, the first two have shown stable urbanization levels, and the latter two have shown significant urban growth. For present-day urban percentage, higher variabilities are seen for the Indo-Gangetic and Yangtze River Basins (coefficients of variation of 88.1% and 83.3%, respectively) than for the Great Lakes and Mid Atlantic regions (56.5% and 62.8%, respectively). The differences between these products are evident at all spatial scales, from global to local. For example, large local scale differences are evident for the highly urbanized Shanghai Metropolitan Area in China (Fig. 2b) and over Delhi, India (Fig. S1). Since illustrating these differences between datasets for all regions of the world is not possible here, we developed a web app for this purpose (https://ee-tc25.projects.earthengine.app/view/urbancomparison).
The disagreements among the data products reflect differences in methods, inputs, and native resolutions24,25. However, a key difficulty in making any apples-to-apples comparisons between these datasets is that, while all the products represent some aspect of physical urbanization, they define ‘urban’ differently. These specific definitions are already baked into the training data (for supervised learning methods), accuracy estimates, and pre- and post-processing steps. For example, among the four ~ 10 m resolution present-day estimates of urban land – WSF 201922, Esri Land Cover26, European Space Agency (ESA) WorldCover27, and Dynamic World28 – the Dynamic World dataset calls this class ‘Built Area’ and includes urban vegetation and green space in that definition, while the ESA WorldCover calls the class ‘Built-up’ and explicitly excludes urban vegetation in their class definition (Table S1). Interestingly, while the Esri Land Cover dataset does not mention inclusion of urban trees, it generally shows much higher urban land percentages across scales than ESA WorldCover (Figs. 1b, 2a, 3), even though both are based on Sentinel-2 data. Some of the differences between these three datasets are related to methodology, in that the Dynamic World and Esri Land Cover datasets use convolutional neural networks that consider contextual information in the classification through the use of convolution kernels, while the ESA WorldCover uses a random forest approach with each pixel classified independently17. Another difference is due to the choice of minimum mapping unit, which is 50 × 50 m for Dynamic World28 and therefore necessitates a mosaic of built and natural surfaces in areas labelled as ‘Built Area’. Finally, the WSF 2019 dataset22 is for human settlements and excludes roads. For readability, and in line with how some of these products have been used in the scientific literature as a proxy for physical urbanization29–31, we refer to all of them as ‘urban’ in the present study.
Urban growth over time
The explosion of medium-resolution global urban products, and global land cover datasets in general, has been largely made possible due to the free release of the Landsat archive in 200832. Consequently, there are several long-term estimates of urban land at the Landsat resolution (30 m) starting from the 1980s. In contrast, the first-generation global urban land cover products were generally limited to the Moderate Resolution Imaging Spectroradiometer (MODIS) resolution of 250 to 500 m resolution and starting from the year 200119. Some of these multi-year urban land cover products do not extend till 2019/2020 and thus were not included in the earlier comparison of present-day urbanization. In total, we examine twelve global data products, including the complete time series (when multiple years are available) of the eight products considered earlier. The four new datasets considered are the Global Artificial Impervious Area (GAIA)33, World Settlement Evolution (WSE)22, the Copernicus Global Land Service (CGLS) product34, and the Global Urban Footprint (GUF)35. All long-term urban datasets show large global urban growth over time during their respective time spans (Fig. 3a). For, GISA (Global Impervious Surface Area)36, GAIA, and WSE – the three datasets with longest time series – global urban percentage increased by 297.4%, 123.4%, and 111.2%, respectively (three-product mean of 177.3%), for the 1985-2015 common period. This pattern of rapid urban expansion is consistently echoed across all continents (Fig. 3; Fig. S2a), with Africa, Asia, and South America showing a notable rise (three-product means of 226.2%, 425.3%, and 186.5%, respectively), although from a lower baseline compared to North America and Europe.
The impact of urban definition and methodology also reflects on the variability in the change in urban land percentage over time across datasets. For instance, the percentage of urban land in the WSF 2019 dataset is much lower than the values in WSF 201537 for all continents. This is for two reasons: 1) the WSF 2019 dataset uses Sentinel-2 instead of Landsat 8, the latter being much finer (10 m versus 30 m); the scale effect25 and 2) the WSF 2019 uses ancillary data to mask out roads to focus only on pixels where people live (Table S1). Another evident difference in time series arises when comparing the MODIS data38 with the others. The global percentage of urban land increases by only 5.5% between 2001 and 2015 according to the MODIS Land Cover; yet the GISA/GAIA/WSE mean change for the same period is around 40.2%, almost an order of magnitude higher. The low estimate of urban expansion in MODIS is a function of its definition of urban as a minimum of 30% impervious at the 500 m scale39. Conceptually, this means that a MODIS pixel starts being classified as urban at a lower percentage than other datasets, which generally consider the dominant land cover (which can exceed 30% of the area) as the class of a pixel or use higher impervious percentage thresholds (50% for GAIA33, for instance), and that a pixel remains urban over a much larger range of values (from 30–100%). Other key specifics of these differences between the global data products are provided in Table S1.
Implications for observational and modeling applications
Global estimates of urban land have become critical for both science and applications. However, most use cases of these datasets do not simultaneously consider multiple estimates due to a combination of legacy, convenience, and potential redundancies. Here we examine how the choice of dataset may lead to biases for some common use cases. These use cases are divided into: 1) direct incorporation of these products to generate derived datasets, 2) combining global urban datasets with estimates of hazards to quantify environmental risks, and 3) using these products as surface constraints in process-based models.
First, global urban land cover datasets are used as inputs for other derived products. For instance, the most commonly used MODIS land surface temperature (LST) products for urban climate studies (MOD11 and MYD1140) use a classification-based emissivity method41 with the pixel emissivity taken from a look-up table and the class of the pixel according to the MODIS Land Cover product. This needs to be done because LST and surface emissivity cannot be analytically separated using only thermal observations42. Similarly, the MODIS evapotranspiration products mask out any pixels that are classified as urban in the MODIS Land Cover data since the empirical model used is not calibrated for urban surfaces43. There are also some inter-dependencies between different global urban land cover datasets. The CGLS product uses WSF 2015 as the training data for its ‘Urban / Built up’ class (Table S1)34, while ‘Urban Areas’, as classified within the ESA CCI product, are identified based on the GUF dataset35 as well as the Global Human Settlement Layer (GHSL)44 datasets (Fig. S3). Other composite urban datasets, such as global annual urban dynamics (GAUD) dataset30, have also been generated by combining various existing estimates (GUF, GAIA, GHSL, etc.)
Second, various urban land cover datasets are used as inputs for examining urban climate impacts and city-level environmental risks. The choice of dataset influences the magnitude of these estimates. For the surface urban heat island (SUHI) intensity, the impact of urbanization on local surface warming11, (Fig. 4a; see Methods), larger values of absolute coefficient of variation are seen for urban clusters in the middle East, parts of India, southern and eastern Africa, and the southwest United States. Although most datasets capture the well-established impacts of background climate on SUHI and its seasonality45, the choice of dataset can have larger impacts in arid regions during summer and for polar climate in winter i.e. when the actual SUHI signal is small, with inconsistencies in the sign of the SUHI seen (Fig. S5). Long-term changes in urban land are often combined with ancillary datasets to examine land use transitions31 and exposure to environmental hazards over time29,11,8,9. The rates of change over time would depend on the choice of dataset (Fig. 3), while most studies typically use a single product. For instance, ESA CCI Land Cover, GHSL, MODIS Land Cover, WSE have been individually used in these types of studies29,11,8,9. Andreadis et al. (2022)29 and Rentschler et al. (2023)8 both examined increased urbanization in flood-prone areas using two different urban land cover datasets (GAUD and WSE, respectively); therefore finding different magnitudes of these changes. We replicate a comparison of urban growth in flood plains46 between 1985 and 2015 based on GAIA, WSE, and GISA here (Fig. 4b), with particularly large differences seen for Asia and Oceania. Sometimes the choice of dataset can lead to artifacts due to mismatch between two products. For instance, Mentaschi et al. (2022)9 combined the GHSL 2018 dataset with the MYD11 LST product to estimate intra-urban SUHI extremes. However, as noted earlier, this LST product is constrained by the MODIS Land Cover through the classification-based emissivity method41. As such, we would expect artifacts in LST for a proportion of pixels due to the GHSL data considering a pixel as urban while the emissivity in the LST product being defined for a rural surface. Similar artifacts would be expected for other combinations of MODIS LST with non-MODIS urban land cover estimates47,48 or when MODIS LST is used to validate simulations from models that use different urban emissivity constraints42,45,49.
Third, urban land cover products are incorporated into process-based models, including weather and climate models, as surface input datasets2,14,23,50. Since different land cover types in land models use distinct prescribed radiative, thermodynamic, and aerodynamics properties, the land cover data used strongly modulates crucial variables like the components of the surface energy budget5,45,49 and thus the lower boundary conditions for the atmosphere in coupled model simulations. One of the most common mesoscale models used for urban climate research – the Weather Research and Forecasting (WRF) model12,51,52 – uses the MODIS urban land cover as the default surface dataset. Newer versions of the urban components of this model can also use the local climate zone (LCZ) classification system53, with a recent global 100 m dataset planned to become the default urban representation for future releases of WRF15. Earth system models (ESMs) rarely resolve urban areas, but one of the few such ESMs with an urban representation – the Community Earth System Model (CESM)54–56 – uses a circa 2001 estimate of urbanization14. This urban dataset is also used in other ESMs that have branched off from CESM57,58 and has also been incorporated into regional models59. Large differences between these three products (MODIS Land Cover, Demuzere et al. (2022)15, and Jackson et al. (2010)14) as well as other present-day estimates of urban land are evident (Figs. 5a, S4). Note that except the MODIS Land Cover data, the other two are not pure estimates of physical urbanization. Jackson et al. (2010)14 actually uses population-based thresholds of urban density while several of the LCZ classes represent different mixes between built and natural surfaces. For example, LCZ9, the sparely built class, is characterized by a high abundance of natural land cover and behaves thermally like a natural land-cover, and is often excluded as a built up class60. As such, there are potential mismatches here that should be kept in mind. For example, since CESM does not use the low-density urban class within the urban model61, it is implicitly assumed that anything up to medium-density would be an appropriate representation of physical urbanization at the grid-scale, which leads to massive overestimation for regions of the world with high population density and low physical urbanization such as Asia and Africa (Fig. 5b). Similarly, since the MODIS Land Cover can be as low as 30% impervious, the model may overestimate urban thermodynamic/radiative/aerodynamic impacts on climate2,49 for urban-to-rural transition zones since urban models have not traditionally accounted for urban vegetation50. As the newer urban land cover datasets are incorporated into these models, it is important to be aware of consistency in definitions. Using Dynamic World, which includes some vegetation in the urban class, is not appropriate if the model treats the entire urban grid as an impervious surface, such as in CESM, its offshoots57,58, and most versions of WRF. Similarly, since WSF 2019 removes urban roads, incorporating this dataset into a process-based model will capture only a fraction of the physical impact of urbanization on weather and climate.