## Materials and Methods
### Methods
#### 2.2.2 Urban Area Extraction  
As no training data on urban vs non-urban areas are available that could be used for a supervised extraction approach, we developed a simple unsupervised method to extract the urban areas from the spatially aggregated RGB\({}_{100}\) surfaces. In topographic map processing, color space clustering based on techniques such as k-means clustering [88] is commonly used for color reduction and color segmentation of scanned maps [30]. Thus, in order to reduce the color complexity, and to group the map color in meaningful ways, we performed k-means clustering on these surfaces, for a range of k \(\epsilon\)[2,10]. For map composites (mosaics of individually scanned map sheets), such as in the Boston and Atlanta study areas (Figure 1c,d), we conducted a separate clustering analysis for each map sheet (Figure 6d) in order to account for potential differences in contrast or color tone (see Figure 6b,c). Moreover, we used the Elbow method [89] to identify the optimum number of clusters per study area. Herein, we used a rule-based decision mechanism to determine which of the obtained color clusters represents the urban areas contained in the historical maps. This mechanism works as follows. We calculated the area proportion of each detected cluster within the built-up areas reported in the GHSL in 1975, aggregated to 100 m \(\times\) 100 m grid cells. Assuming that the urban areas in the historical map (dated earlier than 1975) are contained within the 1975 built-up areas from the GHSL (BUA\({}_{1975}\)), we identified the cluster of the highest area proportion within the BUA\({}_{1975}\) as the cluster likely to represent the urban areas in the historical maps. Moreover, we tested whether the average R, G, and B values of the  
Figure 6: Illustrating the historical urban area extraction method using historical maps and the GHSL. (**a**) GHSL multi-temporal built-up areas (resampled to 100 m \(\times\) 100 m), (**b**) original historical map sheets from approximately 1900, (**c**) generated 100 m \(\times\) 100 m RGB aggregates (averages per channel), (**d**) color clustering results for k = 4, (**e**) target clusters likely representing urban areas, identified by a rule-based decision mechanism taking into account the GHSL areal proportions per cluster and the cluster brightness, (**f**) target clusters within GHSL 1975 built-up areas only, (**g**) post-processed target cluster areas, and (**h**) back-casted urban areas (i.e., extracted historical urban areas integrated with the GHSL multi-temporal built-up areas). RGB\({}_{100}\) cells within that cluster were less than a given threshold value (e.g., <200). Since urban areas in historical maps are typically depicted in dark or saturated colors (black, grey, red), the use of such a simple brightness-based criterion helps to robustly identify the correct target cluster (i.e., the cluster identified as urban area). The target cluster may still contain a considerable number of false positives (e.g., dark text elements or major roads) (Figure 6e). To reduce these artefacts, we implemented a two-step spatial refinement procedure: First, we exclude grid cells of the target cluster located outside of the BUA\({}_{1975}\) extents (Figure 6f), given that this approach only detects urban growth, not urban shrinkage, which is consistent with the implemented GHSL modeling strategy and built-up concept [69]. By using the BUA\({}_{1975}\) as spatial constraint, as opposed to the later epochs from the GHSL, we reduce the temporal gap between the historical map date and the constraining areas, ensuring that a maximum of false positives is eliminated. Second, we implemented a morphological post-processing strategy, removing further segments of the target cluster that are below a specific area threshold t (10, 50, 100 pixels). This is based on the assumption that settlements require a minimum size to be mapped at all.