Global urban structural growth shows a profound shift from spreading out to building up

Experimental design

We modified the methodological approach from our two previously published papers^7,8 on assessing urban growth in both upward and outward domains across global cities. The modifications are threefold. Most importantly, we used a time series developed from a longer record of backscatter data from three scatterometers, with a goal of understanding the dynamics of 3D urban growth over three decades. Second, we used a more rigorous definition of the urban domain and a larger number of cities. Finally, we used the WSF-evo dataset for measuring urban extent and outward growth (compared with nightlight and global human settlement layer BF data in earlier papers). Our spatial analysis covers global urban land on a 0.05° lat/lon grid. Our temporal analysis has annual time steps from 1993 to 2021, the period of available microwave backscatter data. We analyze and compare decadal trends, as this is the duration of individual scatterometer datasets and is long enough for observed trends to be important metrics of change.

Urban domain

In this study, we used global Morphological Urban Areas (MUAs)⁴¹ to define the spatial extent of 1,567 urban areas in 151 countries (Extended Data Table 1 and see Supplementary Information Table 1 for a full list of MUAs). In our earlier papers^7,8, city domains were defined by an 11 × 11 grid centered around the city core. That approach had issues of double counting of urban grid cells in adjacent cities and did not account for the actual shape of metropolitan regions (for example, Extended Data Fig. 8).

Taubenböck et al.⁴¹ identify some MUAs as an aggregation of several cities (for example, Los Angeles–Long Beach–Santa Ana–Riverside–San Bernardino–Mission Viejo in the USA); to simplify notation we designate these by the first city name followed by a plus sign (for example, Los Angeles+; for a full list, see Supplementary Table 2). We overlaid all MUA polygons with a 0.05° lat/lon grid and identified all grid cells that were contained within or intersected an MUA polygon with that MUA (see the example in Extended Data Fig. 8). We also extracted MUA total population and polygon area data from this dataset; population was generated from United Nations data for 2014⁴¹. We computed MUA mean population density as the ratio of total (polygon) population to polygon area.

We excluded five MUAs for which we did not have either backscatter data (Honolulu, USA) or WSF-3D data (Al Ain, UAE; Hafar-al-Batin and Medina, Saudi Arabia; Zahedan, Iran). We then masked out all individual MUA grid cells with open water area ≥50%, using the Global 1-km Consensus Land Cover dataset⁴². After that, individual MUA size ranged from 389 grid cells (Guangzhou+, China) to one grid cell (Arkhangelsk, Russia). The total area of all 34,880 urban grid cells in our analysis is ~890,000 km², about 75% larger than the total area of all MUA polygons (~510,000 km²) because all 0.05° lat/lon grid cells both within and intersecting the MUA polygon boundaries were included (see the example in Extended Data Fig. 8).

To summarize our global analysis, we aggregated many results to 12 global regions, and to illustrate the details at the MUA scale, we provide results for a single large city in each region (Extended Data Table 1). We also partitioned the Africa region into five subregions and the Other Asia region into two subregions (Extended Data Table 1).

WSF data

We aggregated the WSF-evo BF data for 1985–2015 (ref. ²⁰) to 0.05° lat/lon (BF as fraction in range 0–1), and then extracted annual values for all MUA urban grid cells. For each grid cell, we computed annual rates of increase in BF at the grid cell level by ‘decade’ as ΔBF_1990s = (BF₂₀₀₀ − BF₁₉₉₃)/7, ΔBF_2000s = (BF₂₀₀₉ − BF₁₉₉₉)/10 and ΔBF_2010s = (BF₂₀₁₅ − BF₂₀₁₀)/5.

We also aggregated the WSF-3D data¹² on building volume (as sum), building fraction (as mean) and building height (as mean) to a 0.05° lat/lon grid, using the same aggregation method as with WSF-evo data above, and then extracted values for all MUA urban grid cells. These data were not used for any temporal or change analyses presented here, but for a comparison of spatial patterns in 3D building metrics to backscatter across large cities. For comparison to ASCAT backscatter PR, these gridded WSF-3D data were smoothed with a 3 × 3 mean smoothing grid to more closely match the effective resolution of the ASCAT data⁴³. We chose the year of ASCAT data to compare with WSF-3D based on the best temporal match. WSF-3D building height represents 2012 status (based on TanDEM-X digital elevation model) and WSF-3D building fraction is defined by data from 2019 (it is independent from the WSF-evo product). WSF-3D building volume is the product of building height and building fraction and grid area, and so represents some hybrid of the urban state in 2012 and 2019; we chose 2015 ASCAT PR data for comparison to this 3D metric. We also compared annual time series of WSF-evo BF and ASCAT backscatter PR, 2007–2015, for Damascus, Aleppo and Homs in Syria, which were heavily bombed in the early 2010s²⁵.

Urban microwave backscatter data

We extracted summer season backscatter, σ° in dB, for all the MUA urban grid cells from a global 0.05° seasonal backscatter dataset covering 1993–2020 (ref. ⁴⁴); see ref. ⁴⁰ for dataset details. The summer season was designated as July–August–September for all grid cells in the Northern Hemisphere and January–February–March for all grid cells in the Southern Hemisphere. This backscatter annual time step data consists of data from three scatterometers: ERS (1993–2000), QuikSCAT or QSCAT (1999–2009) and ASCAT (2007–2020). We added an additional year of ASCAT backscatter σ° data (2021) to the time series, using the same data source and processing as described in ref. ⁴⁰. We converted backscatter in dB to power return ratio, PR (=10^(σ° / 10)) for all analyses.

There is an offset in backscatter PR between the three sensors⁴⁰ (also see Extended Data Fig. 1), so we did not construct a continuous PR time series. To reduce impacts of short-term PR variability not likely to be due to urban development, we computed grid cell level annual backscatter trends (R command 1m) separately for each sensor and decade: the 1990s: ERS (1993–2000); the 2000s: QuikSCAT (1999–2009 in the Northern Hemisphere, 2000–2009 in the Southern Hemisphere); and the 2010s: ASCAT (2010–2021). The sensor offset leads to a bias in annual trends, which is evident in the comparison of ASCAT and QSCAT trends for their three overlapping years of data, 2007–2009 (Extended Data Fig. 9a). To intercalibrate ASCAT to QSCAT, we scaled the ASCAT trends by the 2007–2009 mean QSCAT PR divided by the 2007–2009 mean ASCAT PR (Extended Data Fig. 9b). This grid cell level intercalibration factor was then applied to all ASCAT trends for the 2010s. Evaluating the ERS–QSCAT intercalibration is more problematic, as there are only 2 years of overlapping data for Northern Hemisphere cities (1999 and 2000) and only 1 year (2000) for Southern Hemisphere cities; 2-year ‘annual trends’ are just differences (with R² identically 1.0) and 1-year trends are undefined. Since most city-level backscatter PR trends were small in the 1990s (for example, Extended Data Fig. 1) a 2-year difference can be more noise than signal. We applied the same 3-year intercalibration method to ERS PR trends for the 1990s, that is, multiplying by the 1999–2001 QSCAT mean PR divided by the 1998–2000 ERS PR (using 2000–2002 mean QSCAT PR for the 12% of our target cities that are in the Southern Hemisphere (Extended Data Fig. 9c,d).

Building floor area data

In statistical datasets, building floor area is often computed as the building footprint times the number of stories in a building, and thus the floor area will be approximately proportional to the total building volume. We retrieved global building floor area data from online reports of the International Energy Agency (IEA). Specifically, we used the IEA global annual total floor area for 2000–2020 ( and this same global floor area for the years 2010, 2015, 2019, 2020 and 2021 but disaggregated into two subsets, ‘advanced economy’ countries and ‘emerging and developing economy’ countries (hereafter ‘other’ countries); country lists for this disaggregation come from IEA global energy analyses⁴⁵ and are listed in Extended Data Table 1).

We retrieved the national area of completed (new) building floor area for China for 1985–2021 (ref. ⁴⁶). This national building footprint data was reported every 5 years from 1985 to 2005, and annually thereafter. We linearly interpolated between the 5-year reporting intervals before 2005 to get an estimate of annual new floor area. We also retrieved annual completed (new) building floor area for two municipalities: Beijing for 1990–2020 (ref. ⁴⁷) and Shanghai for 1978–2019 (ref. ⁴⁸). Since these national and municipal data represent only floor area of new completed buildings, with no information of total existing floor area, we computed accumulated new floor area for each region since the year before the particular data record begins, as a measure of increasing building volume at an annual time step. These datasets do not account for loss of floor area from demolition or restructuring of existing buildings.

Methodology

As backscatter strength is very sensitive to buildings as corner reflectors, construction of new buildings (infilling or lateral extension) or taller buildings can similarly influence backscatter. Therefore, we consider building volume, which would increase with either type of construction, to be the best proxy for backscatter strength, and combine the backscatter data with independent built fraction data to better discriminate different growth types.

We constructed annual time series of MUA-scale means of grid cell backscatter PR (1993–2021) and WSF-evo BF (1990–2015). Across the spatial variation of any large city, QSCAT backscatter has been shown to correlate with building height for nine US cities¹⁸. ASCAT backscatter PR has been shown to correlate with 2015 building volume estimates¹³ for China, Europe and the USA⁴⁰. We extended this analysis to global coverage, comparing 2012 ASCAT backscatter with WSF-3D building height, 2019 backscatter with WSF-3D building fraction and the mean of 2012 and 2019 backscatter with WSF-3D building volume, with WSF-3D data smoothed over 3 × 3 grids. Extended Data Fig. 9 shows these correlations for 12 MUAs, one per region.

To test whether temporal development in urban microwave PR and urban BF track temporal urban building growth, we computed aggregate annual microwave backscatter PR as the sum over all 0.05° urban MUA grid cells (1) globally, (2) in ‘advanced economy’ and ‘other economy’ countries, (3) in China, (4) in Beijing+ and (5) in Shanghai+, and compared these results with annual building floor area data (global and disaggregated into two ‘economies’) and to annual newly completed building floor area data (China, Beijing+ and Shanghai+). Similarly, we computed the mean annual WSF-evo BF for these same domains. We note that the IEA building floor area data is not restricted to urban domains, nor is the China new building floor area data, and that the Beijing and Shanghai statistical floor area data are for these municipality provinces, whose spatial extent does not fully align with the Beijing+ and Shanghai+ MUAs.

We constructed gridded maps of decadal bivariate trends in urban backscatter PR and WSF-evo BF. We then binned grid cells by their values of these decadal bivariate trends of MUA change into four broad categories: low or high rates of change (decadal trends) in backscatter PR coupled with low or high rates of change (decadal trends) in WSF-evo BF; these are (1) slow growth (low ΔPR/Δt and low ΔBF/Δt), (2) area-dominant or outward growth (low ΔPR/Δt and high ΔBF/Δt), (3) rapid 3D or up-and-out growth (high ΔPR/Δt and high ΔBF/Δt) and (4) height-dominant or upward growth (high ΔPR/Δt and low ΔBF/Δt), where t is time in years. We also quantified whether decadal growth rate trends in urban BF and PR accelerated or decelerated from the 1990s to the 2000s to the 2010s, aggregated at the regional scale and for individual MUAs.

Finally, we implemented a k-means analysis to group the global set of MUA grid cells into clusters. We used all MUA grid cells (n = 34,880), and performed separate, independent clustering for each decade, each based on the decade’s ‘initial state’ (for example, year 2000 QSCAT PR and WSF-evo BF) and growth rate (for example, mean annual rates of increase in QSCAT PR and WSF-evo BF in the 2000s). Since these clusters were generated for each decade independently, we use the actual PR values, not intercalibrated. Following ref. ⁸, we generated five clusters or typologies for each decade. We then identified which clusters had similar characteristics between decades and through this identified unique clusters across the three decades. Finally, we tracked grid cell trajectories from their 1990s typology to their 2000s typology to their 2010s typology, to study shifts and changes in urban growth over three decades at global and regional scales.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.