## Materials and Methods
### Methods
#### 2.2.1 Preprocessing  
Based on metadata for the USGS-HTMC (available from [[https://thor-f5.er.usgs.gov/ngtoc/metadata/misc/](https://thor-f5.er.usgs.gov/ngtoc/metadata/misc/)]([https://thor-f5.er.usgs.gov/ngtoc/metadata/misc/](https://thor-f5.er.usgs.gov/ngtoc/metadata/misc/)), accessed on 10 September 2021), the geographic footprints of each map sheet contained in the archive can be obtained, allowing for reconstruction of the grid (the so-called graticule) in which the USGS-HTMC map sheets are organized. For each quadrangle (i.e., grid cell of the graticule), we identified the earliest available map sheet and its scale within the boundaries of the Boston and Atlanta metropolitan statistical areas in 2010 [87] and automatically downloaded these maps from the AWS S3 archive [81]. By doing so, we obtained 33 maps for the Boston metro area, and 180 maps for the Atlanta metro area. Based on the corner coordinates available for each map sheet, we removed the map collars and generated a seamless mosaic of the maps per study area (see Figure 1c,d). The preprocessing steps for the US study areas are shown in Figure 12a. For the study areas in the UK, we obtained the historical maps from [82,83] and mosaicked them and, in case of the London study area, georeferenced them. Individual map sheets for the Sao Paulo and Lahore study areas were manually georeferenced. All maps and map composites were then spatially aggregated by computing the RGB averages, separately per channel, within blocks of 100 m \(\times\) 100 m (RGB\({}_{100}\)). This 100 m \(\times\) 100 m grid represents the analytical unit for the subsequent analyses. Such a spatial aggregation allows for the fast processing of large amounts of maps and facilitates the integration with other gridded data. Despite a native spatial resolution of the GHSL of 30 m \(\times\) 30 m, we use a grid of 100 m \(\times\) 100 m cells as our analytical unit, since such spatially aggregated RGB signals from the historical maps will allow modelling of \"urban areas\" in a generalized and more robust manner, regardless of how built-up structures are depicted in historical maps (e.g., individual buildings, street blocks, etc. ), as a 100 m \(\times\) 100 m grid cell typically extends across multiple street blocks. The resulting urban areas are assumed to encompass built-up structures, impervious surfaces, and other land use types considered to be \"urban\" in a broader sense, such as smaller intra-urban green spaces or similar. The effect of spatially aggregating RGB information found in the historical maps to create the RGB\({}_{100}\) layers can be seen in Figure 6b,c. Using a finer analytical unit may result in oversampled results, in particular when the scanning resolution of historical maps is low. Moreover, historical maps may suffer from positional inaccuracies introduced by a range of reasons, such as inaccurate topographic measurements due to the instruments used at the time of map creation, cartographic displacements and generalization, humidity or heat-induced (nonlinear) deformations of the paper map, and inaccurate georeferencing [32]. Most of these positional uncertainties are unknown and difficult to quantify. Thus, we use an analytical unit of 100 m \(\times\) 100 m to obtain historical urban extents at a fine spatial granularity (as compared to coarser spatial-historical datasets such as HYDE or HISDAC-US), accounting for and potentially mitigating positional inaccuracies in the signals from the historical maps to a certain degree. The GHSL built-up land surface as well as the HYDE urban area raster data were both clipped to the historical map extent of each study area, and resampled to create a 100 m \(\times\) 100 m grid that is consistent with the aggregated map data (Figure 6a, see Figure 12b for the underlying data processing).