2.1 INSPIRE building data processing

2.1.1 Data collection

For most parts of Spain, cadastral building data are available through an ATOM interface. The ATOM feed format is an XML-based language that allows automated, web-based content retrieval by providing a machine-readable web interface (IBM, 2023). The ATOM XML files are organized in a hierarchical manner (see Table 1 for examples) and allow the building data to be accessed in GML format, at the municipality level. We created a Python script to automatically download these GML files (https://github.com/johannesuhl/hisdac-es, last access: 10 October 2023). In some cases (e.g., Basque Country, Navarre region, which have their own cadastral systems) we manually downloaded the data available as a Web Feature Service (WFS) (Table 1). All downloaded data were projected in UTM zone 30N (EPSG:25830).

Table 1Source data overview.

^a Last access: 10 October 2023. ^b Last access: 15 October 2021. In the meantime, the WFS has been replaced by a repository available at https://filescartografia.navarra.es/2_CARTOGRAFIA_TEMATICA/2_7_CATASTRO/2_7_1_SHAPEFILE/ (last access: 10 October 2023).

Download Print Version | Download XLSX

2.1.2 Data preprocessing

After downloading and gathering the building data for over 8100 municipalities, covering all regions of Spain, we first calculated and attached the building footprint area (obtained after reprojecting to EPSG:3035) and converted the polygonal building footprint data to centroids, retaining all relevant attributes, to reduce the computational effort for the subsequent data processing. Despite the common INSPIRE framework, attribute names and building function classes differed slightly between the different data sources. Thus, we harmonized the data by renaming columns, and by applying a common building function classification scheme, including the six building function classes “residential”, “commercial”, “industrial”, “agricultural”, “public services”, and “office”. “Public services” is probably the broadest of these categories, including governmental buildings, but also health-related buildings and cultural institutions (e.g., churches or museums). Specifically, building function ontologies differed slightly for the data from the region of Navarre and the province of Araba and were consistent across the other regions/provinces. For example, commercial buildings in Navarre are labeled “trade” instead of “commercial”. The applied mapping scheme can be accessed on the HISDAC-ES GitHub repository (https://github.com/johannesuhl/hisdac-es/blob/main/landuse_mapping.csv, last access: 10 October 2023).

The gridded surfaces of HISDAC-ES are provided in three different spatial reference systems: (a) ETRS89 UTM zone 30N (EPSG:25830) for the Spanish mainland, Balearic Islands, and the exclaves Ceuta and Melilla located in northern Africa; (b) REGCAN-95 (EPSG:4083) for the Canary Islands; and (c) the reference grid of the European Environmental Agency (EEA), which is based on the ETRS89 Lambert azimuthal equal area projection (LAEA; EPSG:3035) and is commonly used for pan-European statistical mapping (e.g., CORINE Land Cover data), for all Spanish territory. This way, users can refer to the data in UTM/REGCAN for mapping purposes (north-oriented, angle-preserving), while the datasets in area-preserving LAEA projection can be used for statistical modeling and integration with other datasets (e.g., gridded statistical data from Eurostat (https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/grids, last access: 10 October 2023) or gridded Spanish census data (INE grid (https://www.ine.es/censos2011_datos/cen11_datos_resultados_rejillas.htm, last access: 10 October 2023)). Thus, we reprojected the building centroids into these reference systems, yielding two sets of harmonized building centroid data: (a) in UTM/REGCAN and (b) in LAEA projection.

After an initial examination of the attribute coverage and completeness, we decided to focus on six well-covered attributes, measuring different aspects of the built environment, including the construction year, building function, number of dwellings, number of building units, building indoor area, and building footprint area (Fig. 1). For clarity, the number of dwellings describes the number of housing units in residential buildings, whereas the number of building units also includes the number of units within non-residential buildings (e.g., number of commercial businesses within a building complex). Furthermore, the building indoor area represents the attribute “official area” and measures the gross indoor area (across all stories) within a building.

Figure 1Input data: INSPIRE-conforming building footprint data provided by the Spanish authorities, including several attributes. (a) Year built, (b) building use, (c) number of dwellings, (d) number of building units, (e) building indoor area, and (f) building footprint area. Data shown for the city of Valencia.

2.1.3 Data aggregation

Based on the preprocessed, harmonized building data and the six selected thematic attributes, we created a range of different aggregated representations of the data. These aggregations include (a) spatial aggregation into grid cells within regular spatial grids of 100 m × 100 m, (b) aggregation to the municipality level by calculating zonal statistics per municipality polygon, and (c) temporal aggregation by stratification into different temporal classes. The combination of spatial aggregation and temporal stratification applied to the different thematic attributes yields a range of different sets of gridded surfaces. For example, we calculated the sum and the mean of the building units (BUNITS) and dwellings (DWEL) over all buildings within a given grid cell, as well as both the sum and mean building indoor area (BIA) and building footprint area (BUFA), based on the building centroids located within a grid cell. The resulting gridded surfaces represent physical features of the built environment. Similarly, we calculated the minimum, maximum, mean, median, mode, and the variety of construction years (COY) per grid cell, which measures settlement age (heterogeneity) and quantifies construction/remodeling activity within each grid cell. Thus, COY statistics describe the age-related features of the built environment.

We stratified the building records by their construction year into temporal classes (epochs) based on 5-year intervals (e.g., built-up before 1900, before 1905) and calculated the number (or density) of buildings (BUDENS) and the total building footprint area (BUFA) per grid cell in each of these epochs from 1900 to 2020. We further binarized these grid cells to measure developed area (DEVA) (i.e., grid cells containing at least one building) and undeveloped areas in each epoch. These gridded surfaces measure the physical evolution of the built environment in Spain. Similarly, we thematically disaggregated the building count surfaces per epoch based on the building function attribute of the buildings in each grid cell and epoch, yielding time series of function-specific building density layers, for six types of building functions (i.e., residential, commercial, industrial, agricultural, public services, and offices) as a proxy measure for built-up land use evolution from 1900 to 2020. Table 2 provides an overview of the gridded surfaces and spatial variables generated by these data processing steps. These surfaces quantify, for example, the building density (i.e., number of buildings per grid cell, BUDENS), the built-up surface density (i.e., building footprint area per grid cell, BUFA), or the built-up intensity (i.e., the total building indoor floor area per grid cell, BIA).

Table 2Overview of the gridded surfaces of the HISDAC-ES, measuring different components of the built environment, including the underlying variables and statistics.

Download Print Version | Download XLSX

As mentioned before, all gridded surfaces are available in UTM Zone 30N for the Iberian Peninsula, in REGCAN-95 for the Canaries, and in Lambert azimuthal equal area (LAEA) projection for the whole extent. For selected variables, we also provide a time series of zonal statistics, aggregated to the municipality boundaries² (see Sect. 3.5). We created these zonal statistics based on (a) spatially joining the municipality identifier and area to each building centroid of our harmonized building dataset (point-in-polygon query) and (b) deriving statistics (sums, densities) for each municipality. All data processing, as well as the evaluation experiments and data visualizations, was implemented in Python 3.8, using libraries such as NumPy, SciPy, GDAL, GeoPandas, pandas, Matplotlib, PIL, and Esri ArcPy. The core component of our grid cell aggregation procedure is the “binned statistic 2D” function in SciPy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic_2d.html, last access: 10 October 2023). The overall processing workflow from the building data to the spatial layers of the HISDAC-ES is shown in Fig. 2.

Figure 2Data processing workflow to create the spatial data layers of the HISDAC-ES, measuring multiple dimensions of the built environment in Spain from 1900 to 2020, obtained from cadastral building footprint data available via the INSPIRE Spatial Data Infrastructure (SDI).

Download

2.2 Evaluation data and agreement assessments

As historical spatial data are generally scarce, the evaluation of the produced historical data is difficult. In order to evaluate the quality of the produced spatial layers in the HISDAC-ES as thoroughly as possible, we employed a range of independent datasets that exhibit coherence to the spatio-temporal processes measured by HISDAC-ES and carried out different evaluations and cross-comparisons. Specifically, we used three types of spatial data for these experiments: (a) recent, remote-sensing-derived datasets (i.e., the Global Human Settlement Layer (GHSL), CORINE Land Cover); (b) spatial–historical land use models (i.e., History Database of the Global Environment); and (c) historical cartographic data (i.e., historical maps and urban atlases) and orthoimagery. For most of these experiments, we implemented the following strategy: (1) when comparing HISDAC-ES to other gridded datasets, we downsampled the dataset of higher resolution to the dataset of lower resolution. This way, additional uncertainty introduced by resampling is kept to a minimum. (2) We conducted agreement experiments at the grid cell level, i.e., based on cell-by-cell map comparison or correlation analysis. (3) From these cell-by-cell-level comparisons, calculated agreement metrics within local or regional strata, defined by administrative boundaries or other classifications, are granular enough that we could assess regional variations of agreement, and large enough to ensure statistical robustness within each local stratum.

2.2.1 Global Human Settlement Layer (GHSL)

To evaluate the plausibility and reliability of the developed area (DEVA) layers from 1975 to 2015, we used built-up areas from the Global Human Settlement Layer (GHS-BUILT R2018, Florczyk et al., 2019) for comparison. The GHS-BUILT surfaces are derived from multispectral remote sensing data (Landsat sensors, Sentinel-2) and map built-up areas globally from 1975 to 2014, at a spatial resolution of 30 m (Fig. 3a, b). They are accompanied by a rural–urban classification (settlement model GHS-SMOD; Fig. 3c, d). GHS-SMOD is available at a resolution of 1 km and classifies each location on Earth into one of seven classes of urbanness, ranging from sparse rural settlements to high-density urban centers (Florczyk et al., 2019). While more recent versions of GHS-BUILT are available at the time of writing, we decided to use GHS-BUILT R2018A due to its fine spatial resolution and because a lot of work has been done and published to quantify the accuracy of GHS-BUILT R2018A across the rural–urban continuum and over time (e.g., Liu et al., 2020; Uhl and Leyk, 2022b, c), whereas little information is available on the accuracy of newer, multi-temporal GHS-BUILT datasets. For example, it has been reported that GHS-BUILT R2018A yields an average Intersection over Union of around 0.35 in 1975, in rural areas, to around 0.65 in 2018, in urban areas, respectively, for selected study areas in the United States (Uhl and Leyk, 2022b), and correlation coefficients of built-up surface fraction >0.7, compared to reference data for selected cities in China (Liu et al., 2020).

Figure 3Remote-sensing- and model-based evaluation data: built-up areas from GHS-BUILT R2018A in (a) 1975 and (b) 2014; GHS-SMOD rural–urban classes in (c) 1975 and (d) 2014; CORINE Land Cover data in (e) 1990 and (f) 2018; and urban area fraction from the HYDE v3.2 dataset in (g) 1900 and (h) 2015. Panels (a)–(f) show the city of Valencia, and panels (g) and (h) show mainland Spain and the Balearic Islands.

We resampled GHS-BUILT to the HISDAC-ES grid (i.e., upsampling from 30 to 100 m spatial resolution) for the epochs 1975, 1990, 2000, and 2015 and, thus, obtained binary grid cells (i.e., built-up vs. not built-up) for each epoch. To reduce spatial misalignment effects, we first upsampled the 30 m GHS-BUILT data to a 10 m × 10 m grid, nesting within the 100 m × 100 m HISDAC-ES grid, and then downsampled to the target grid, assigning 1 (built-up) if at least one 10 m grid cell within the target cell was labeled as built-up. We then quantified the agreement between the resampled, binary GHS-BUILT and the DEVA layers using Precision, Recall, and F1 score for each epoch and for each municipality. Specifically, we overlaid the binary raster surfaces of DEVA and GHSL and calculated the number of true positive (TP), false positive (FP), and false negative (FN) grid cells within each municipality polygon. These zonal statistics of binary agreement categories were then used to calculate municipality-level Precision, Recall, and F1 score. Moreover, we expected the agreement to vary across the rural–urban gradient (see Leyk et al., 2018; Uhl and Leyk, 2022a, b). Hence, we calculated these agreement metrics for each year, within each of the seven GHS-SMOD rural–urban classes, based on the zonal counts of TP, FP, and FN per SMOD class. Moreover, we generated a continuous rural–urban index for each municipality based on the GHS-SMOD layers, constructed from the weighted histogram of SMOD class instances within each municipality polygon (see Uhl et al., 2023c, for details), and assessed the municipality-level agreement trends across this rural–urban continuum (Sect. 4.1). Here, it is worth noting that we qualitatively compared the built-up areas from GHS-BUILT to the World Settlement Footprint (WSF) Evolution dataset (Marconcini et al., 2020a) and found high levels of agreement between these two datasets (Fig. A1). Thus, herein, we compared HISDAC-ES to GHS-BUILT only.

2.2.4 Historical maps and orthoimagery

While the evaluation approaches described in the previous sections are based on measured or modeled data, they suffer from measurement error, resampling errors, and other incompatibilities that may bias the comparative analyses. Thus, we decided to include alternative, historical data sources in our evaluation analyses, allowing for a more unbiased evaluation of HISDAC-ES layers in early years. These data sources include (a) historical planimetric maps (shown for the city of Alicante in Fig. 4a), (b) aerial photographs from 1956 (Fig. 4b), and (c) an urban atlas (Remírez et al., 1998) depicting different urban development phases (Fig. 4c). While no quantitative information on the accuracy of these data sources is available, the underlying maps are handcrafted and based on manual interpretation of orthophotos, topographic measurements, or local domain knowledge and can be deemed to be relatively accurate. Specifically, we manually digitized the areas developed in different time periods for the city of Alicante and Madrid and obtained a similar vector dataset, depicting different historical urban development phases for the city of Valencia (courtesy of Carmen Zornoza-Gallego, 2022a). We quantitatively assessed the agreement between these historical urban extents and the MINCOY/DEVA layers (Sect. 4.4). Moreover, we visually compared the built-up areas depicted in planimetric topographic maps from around 1900³ to the HISDAC-ES DEVA layer for several urban and rural places (Sect. 4.4). Finally, we manually delineated the approximate urban boundaries for the cities of Santiago de Compostela, Madrid, and Alicante based on aerial imagery from 1956⁴ and compared them qualitatively to the developed areas in HISDAC-ES in the same year (Sect. 4.5).

Figure 4Evaluation data based on historical cartography and orthoimagery: (a) Planimetría historical topographic map from 1910 (Minutas catastrones at scale 1:50 000; Instituto Nacional de Geografía, 2022), (b) historical aerial image from 1956 (Ortofotos AMS(B) 1956–1957), and (c) historical urban atlas illustrating different urban development phases (Remírez et al., 1988), and (d) urban areas for different time periods, manually digitized from (c). All datasets shown for the city of Alicante.

3.1 Physical characteristics

Figure 5 displays the gridded surfaces measuring selected contemporary, physical features (in the year 2020) of the Spanish built environment, exemplarily shown for the city of Valencia. Surfaces showing the sums of these features exhibit interesting spatial patterns of the density of building indoor area (BIA) and footprint area (BUFA), decreasing from the city center towards the outskirts, whereas the density of dwellings (i.e., housing units) is higher in the outskirts than in the center part. The grid cell means of these variables are measures of (vertical/horizontal) building size and exhibit different patterns, illustrating the presence of large multi-apartment complexes in the outskirts and small, historical buildings in the center part of the city.

Figure 5Examples of the created gridded surfaces of the HISDAC-ES, quantifying physical characteristics of the built environment: building indoor area (BIA), building footprint area (BUFA), number of dwellings (DWEL), and number of building units (BUNITS). HISDAC-ES contains the sum (top row) and the mean (bottom row) of these variables per grid cell. All datasets shown for the city of Valencia.

3.2 Age-related characteristics

Age-related statistics of the Spanish building stock are measured by different statistics calculated using the construction year attribute of the buildings within each grid cell. For example, the minimum construction year (MINCOY) impressively shows the settlement age patterns in the metropolitan area of Valencia (Fig. 6), depicting the historical city core, as well as recently developed suburban areas in the urban fringe, and older settlements in the surrounding villages. Similar patterns can be observed in the mean (MEANCOY), median (MEDCOY), and the most frequent (i.e., mode) construction year (MODECOY). The maximum construction year (MAXCOY) per grid cell measures the year of the last modification of the building stock and, alongside the construction year variety (i.e., the number of unique construction years per grid cell, VARCOY), is a measure of construction activity, highlighting areas characterized by heavy urban renewal processes. Here it is worth noting that the construction year on record may also represent the year of the last building reformation, which introduces certain bias in the created surfaces (see Discussion section).

Figure 6Examples of the created gridded surfaces of the HISDAC-ES, describing building-age-related characteristics of the built environment: minimum construction year, maximum construction year, mean construction year, median construction year, mode (i.e., most frequent) construction year, and the variety (i.e., number of unique construction years) per grid cell. All datasets shown for the city of Valencia.

3.3 Evolutionary characteristics

Grid-cell-level statistics (i.e., sums of BUFA, counts/densities of buildings, BUDENS) stratified by the construction year attribute yield a time series of gridded surfaces measuring the long-term evolution of cities, towns, and villages. For example, the BUFA and BUDENS surfaces show how the built-up intensity (as measured by built-up surface density and building density) has increased from 1900 to 2020 (Fig. 7). These gridded surfaces uniquely document the long-term urban growth processes, measured at fine spatial granularity and over long time periods. The derived DEVA surfaces show developed/undeveloped land for each point in time, facilitating quantitative, multi-temporal analysis of urban form, e.g., using landscape metrics (Uhl et al., 2021b). The high temporal resolution (i.e., 5 years) of these multi-temporal layers enables the measurement of the evolution of urban extents and building density at fine spatial and temporal detail, as illustrated in a complete time series for the city of Valencia (Appendix Fig. B1). The additional stratification of the BUDENS layers by building function disaggregates the building stock spatially, by age, and by function. As an example shown in Fig. 7 (right column), industrial land use has heavily increased in suburban areas, especially in the northern suburbs (between 1900 and 1960) and later in the southern suburbs (after 1960).

Figure 7Examples of the created gridded surfaces of the HISDAC-ES, quantifying evolutionary characteristics of the built environment: multi-temporal layers of developed area (DEVA), building footprint area (BUFA), and building density (BUDENS), as well as industrial land use in 1900, 1960, and 2020 (as one exemplary land use class). All datasets shown for the city of Valencia.

While the examples in Figs. 5–7 show the city of Valencia, we would like to emphasize the country-wide coverage of HISDAC-ES. For example, the minimum construction year surface (MINCOY) reveals commonalities and differences between settlement age patterns in different cities in Spain (Fig. 8), including polycentric development (e.g., Barcelona, Seville) and monocentric development (most other cities shown). Moreover, HISDAC-ES not only measures urban development across different cities, e.g., by means of the MINCOY surface (Fig. 8), but also long-term land development processes in rural areas, including towns, villages, and scattered, unincorporated settlements, as exemplified by the DEVA surfaces (Fig. 9). The DEVA layers reveal further detail on the spatial configuration of cities in early years, allowing, e.g., for the computation of historical, urban morphological indicators (Appendix Fig. B2). We also provide several supplementary animated data visualizations illustrating the value of HISDAC-ES for quantifying long-term urbanization processes (see Video supplement).

Figure 8HISDAC-ES settlement age surfaces for 30 cities in Spain: minimum construction year (MINCOY) shown for a selection of 30 cities, arranged in an approximate, quasi-geographic space (upper right shows the northeast; lower left shows the southwest).

Figure 9HISDAC-ES multi-temporal surfaces of developed areas: DEVA layers for 1900 and 2020, shown (a) for Greater Barcelona, (b) for the island of Mallorca, and (c) for the greater León area. Basemap: Esri, USGS, NOAA.

3.4 Built-up land use surfaces

Lastly, we show the building density layers stratified by building function, measuring the spatial distribution of different built-up land use classes (Fig. 10). These gridded surfaces not only highlight the dominance of residential land use but also illustrate peri-urban clusters of industrial land use, as well as spatial patterns of commercial land use, which has a mixed clustered and scattered spatial pattern, or agricultural land use, mostly occurring in peri-urban areas to the northeast of Valencia. These surfaces, along with the corresponding multi-temporal land use surfaces (Fig. 7d) enable the quantitative assessment of land use specific evolution of the built environment and add a unique thematic component to the HISDAC-ES data layers. Besides the building density layers stratified per building function, we also provide the total building indoor area and building footprint area per grid, for residential buildings only (RES_BIA, RES_BUFA), facilitating the integration with historical population data and population disaggregation. We created RES_BIA and RES_BUFA for each decade from 1900 to 2020 (approximately in line with the decennial census), and in the LAEA grid only, as these layers are intended for statistical use.

Figure 10Gridded surfaces describing the spatial distributions of different building function classes or building-related land use classes, shown for the city of Valencia in 2020.

3.5 Municipality-level statistics

We provide zonal statistics of building footprint data for over 8100 municipalities in Spain as tabular data and geospatial vector data. These datasets contain the zonal sums of selected variables (i.e., building counts, as well as BUFA, BIA, DWEL, BUNITS, RES_BUFA, RES_BIA) as well as corresponding densities (per municipality area) and allow for coarse-scale analyses, and for the joint analysis with historical population data, available at the municipality level. The visualizations in Fig. 11 illustrate the usefulness of such aggregated statistics to observe and quantitatively assess broad-scale settlement and building stock age patterns. These patterns can be interpreted in the context of historical settlement development but also provide insight into the contemporary building stock age and its spatial variation. As the absolute counts per municipality may be affected by regional trends of municipality area (Fig. 11, top row), we also provide these statistics as densities normalized by the municipality area, which show a different picture (Fig. 11, bottom row) (see Video supplement).

Figure 11Selected multi-temporal municipality-level statistics: building counts and built-up surface density per municipality in 1900, 1960, and 2020. See Appendix Fig. C1 for corresponding maps of the Canary Islands.

4.1 Multi-temporal built-up area evaluation (1975–2014)

The comparison of the developed area (DEVA) to built-up areas from GHS-BUILT reveals several trends: (1) the agreement between DEVA and GHS-BUILT changes from rural to urban areas. Specifically, the precision of DEVA is high in urban areas and low in rural areas (Fig. 12a). A low precision in an agreement assessment would indicate high commission errors, which in this case implies that DEVA labels much more grid cells as built-up than GHS-BUILT. This is encouraging as previous work has shown that GHS-BUILT tends to underreport built-up areas in rural regions (Leyk et al., 2018; Uhl and Leyk, 2022b), and thus, the DEVA layers appear to account for this shortcoming. Similarly, recall of DEVA is slightly lower in urban areas such as the Madrid region (Fig. 12a). As previous work revealed, GHS-BUILT tends to overestimate built-up areas in urban settings (i.e., roads are often classified as built-up). The DEVA layers are not affected by this type of misclassification, resulting in lower recall values. Hence, both low precision (rural settings) and low recall (urban settings) imply high accuracy in the DEVA layers, as the reference data (i.e., GHS-BUILT) suffer from the described shortcomings. (2) We observe an increase of precision over time (Fig. 12b), which is likely due to increasing completeness of built-up areas in GHS-BUILT, particularly in rural areas. Recall, however, shows a different trend over time (Fig. 12b), maximizing, on average, across all municipalities in 1990 and decreasing towards recent epochs. This is likely a combined effect of (a) increasing incompleteness in GHS-BUILT as we go back in time due to poorer quality of underlying Landsat data and (b) increasing incompleteness in DEVA as we go back in time due to a survivorship bias in the INSPIRE building footprint data. Specifically, new buildings that replace an existing (old) building will be attributed with the construction year of the replacement, and the building that existed prior to the replacement is not contained in our data. Thus, urban renewal causes this bias in our data, and this bias manifests in lower recall values towards early points in time.

Figure 12Evaluation of the HISDAC-ES developed area (DEVA) in comparison to built-up areas from the Global Human Settlement Layer per municipality, over time, and across the rural–urban continuum, by means of map comparison. (a) Maps of municipality-level precision, recall, and F1 score in 1975 and 2014; (b) overall temporal trends of municipality-level precision, recall, and F1 score; (c) temporal trends of precision, recall, and F1 score calculated globally (i.e., at the country-level) within strata of the seven GHS-SMOD rural–urban classes; and (d) distributions of municipality-level agreement metrics within strata of a GHS-SMOD-based rural–urban index calculated per municipality. See Appendix Fig. D1 for corresponding maps of the Canary Islands. All agreement metrics are obtained by map comparison on a cell-by-cell basis at a spatial resolution of 100 m × 100 m, calculated locally within municipality boundaries in (a), (b), and (d) and calculated globally (i.e., overall metrics for the whole country) within areas delineated by GHS-SMOD classes in (c).

Looking at the agreement trends over time and across the GHS-SMOD rural–urban classes (Fig. 12c), we observe a sharp increase in agreement from rural to urban areas, and a slight increase over time, implying that the reliability of the DEVA layers is highest in urban centers. Lastly, the distributions of municipality-level agreement metrics across rural–urban strata confirm this trend (Fig. 12d). The peaks in recall in the low-density and rural cluster strata, across all years, indicate that the effects of incompleteness in the reference data (caused by omitting rural settlements) and in the DEVA (caused by survivorship bias) are of similar magnitude and, thus, cause higher levels of agreement. Here it is noteworthy that the more recent GHS-BUILT v2022 is likely to perform better in rural areas, and thus, precision of the HISDAC-ES is expected to increase in such areas, and the agreement gradients across the rural–urban continuum are expected to be less steep (Uhl and Leyk, 2023). However, as mentioned above, we prefer to use GHS-BUILT R2018A because its accuracy has been well studied and makes our interpretations more robust.

4.2 Land use evaluation 1990–2020

We cross-tabulated the land-use-specific building counts from HISDAC-ES in 1990 and 2020, within the land cover classes from CLC for the years 1990 and 2018 (Fig. 13). The absolute building counts per land cover class in Fig. 13a indicated that most buildings underlying the HISDAC-ES fall into urban fabric, industrial, or commercial areas; agriculturally used areas; and, to a lesser extent, areas characterized by forest (Fig. 13a). When plotting the proportions of buildings per HISDAC-ES land use class (Fig. 13b) or per CLC class (Fig. 13c), we observe more interesting patterns. For example, most buildings of any land use (except agriculture) are in areas of continuous urban fabric. Agriculturally used buildings are mostly located in areas classified as “complex cultivation patterns” in CLC. This indicates that the agricultural land use as reported in the INSPIRE building data is highly accurate. Moreover, in 2020, the proportion of buildings in “discontinuous urban fabric” has increased, as compared to 1990, which may be an effect of suburbanization and increasing low-density built-up areas. Finally, the cross-tabulation relative to the CLC classes shows that residential land use is the most dominant across all land cover classes, with a few interesting exceptions: industrial land use also has high proportions in CLC “industrial or commercial units”, and buildings attributed as “public services” in the cadastral data have a peak in “port areas” and “airport” CLC classes. Agricultural buildings peak in CLC classes “rice fields”, “annual crops”, and “agro-forestry areas”, which confirms the high levels of coherence between the two datasets. The peak of agricultural buildings in the “inland marshes” class may indicate higher levels of confusion between CLC classes “inland marshes” and “rice fields”, which may be difficult to differentiate. All these observations confirm that the land use information from the cadastral building data and the derived HISDAC-ES land use layers seems highly plausible. Note that there is a slight temporal gap between the two datasets, as the most recent CLC data are from 2018. However, we expect this discrepancy to be of minor importance.

Figure 13Comparison of the HISDAC-ES land use data (columns) to land cover classes from CORINE Land Cover (rows) for the years 1990 and 2018. The heatmaps show the number of buildings (BUDENS) assigned to grid cells within each combination of HISDAC-ES land use and CLC classes. These bidimensional frequency maps are shown in three variants showing (a) absolute building counts, (b) proportions of buildings per HISDAC-ES land use class, and (c) proportion of buildings per CLC class. Note that the 2020 HISDAC-ES data were compared to the 2018 CLC data.

Download

4.3 Long-term trajectory evaluation (1900–2015)

While the comparisons of HISDAC-ES to GHS-BUILT and CORINE Land Cover data focus on recent decades, the comparison to HYDE's urban fraction estimates examines the long-term agreement with the BUFA and DEVA evolution layers. We observe that the correlation between urban area fraction and BUFA is high after 1980 (i.e., >0.8 for most regions) and decreases as we go back in time but is still at around 0.6 in 1900, in most regions (Fig. 14a). This drop in correlation could be explained by the previously mentioned survivorship bias in HISDAC-ES but could also be attributed to uncertainties in the model-based urban area fractions in HYDE. The correlation to DEVA (Fig. 14b) shows a similar trend but slightly lower levels of correlations. Interestingly, correlations are highest in the southern region (i.e., Andalucia), possibly due to low levels of building stock renewal and, thus, a weaker effect of the survivorship bias in the HISDAC-ES data. Interestingly, correlations reach an early peak for the Madrid region (1960s for BUFA, 1930s for DEVA) and then drop. Such a decreasing agreement towards recent epochs could be related to heavy (peri-)urban renewal in the Madrid region, which would be less well captured in HISDAC-ES (cf. Sect. 4.5). Besides this comparison of quantitative measures per grid cell, we also compared the agreement between developed and undeveloped grid cells in DEVA and HYDE, as measured by the time series of F1 scores in Fig. 14c. We observe extremely high agreement (>0.95) in recent decades and just slightly lower F1 scores in the beginning of the 20th century. These high levels of agreement of the HISDAC-ES and model-based urban area estimations from HYDE underline the high quality of the HISDAC-ES evolutionary layers. It is noteworthy that the high F1 scores may be an effect of the relatively low spatial resolution of HYDE (5^′ × 5^′ grid cells).

Figure 14Regionally stratified agreement analysis of HISDAC-ES developed area (DEVA) and building footprint area (BUFA) variables against HYDE urban area estimates. (a) Time series of Pearson's correlation coefficient of BUFA and HYDE, (b) time series of Pearson's correlation coefficient of DEVA and HYDE, and (c) time series of F1 scores based on developed vs. non-developed grid cells in the HISDAC-ES and HYDE urban area fraction (5^′ × 5^′ grid cells). Regional stratification was done based on the seven NUTS-1 regions (https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts, last access: 10 October 2023).

4.4 Comparison to historical maps

The previous evaluations are based on either remotely sensed data or model-based reference data. Thus, those datasets are limited in their temporal coverage or suffer from uncertainty themselves. Hence, we used multi-temporal urban areas manually digitized from historical maps, covering the period from approximately 1900 to present, for the cities of Alicante, Madrid, and Valencia (Fig. 15). These extents were manually digitized for Alicante and Madrid from an urban atlas (Remírez et al., 1988; Valencia data are courtesy of Carmen Zornoza-Gallego), resulting in increments of urban area newly added in a given time period (see Fig. 4a, b). We rasterized the resulting vector data in the HISDAC-ES grid, encoding the earlier year of each time period (Fig. 15a). For each city, we reclassified the MINCOY surface to match the settlement age categories from the digitized urban areas (Fig. 15b) and calculated agreement measures for each city, based on both the cumulative urban area per point in time and the newly added built-up area in each period (Table 4). As the MINCOY surface only encodes the year of earliest settlement per grid cell, disregarding the settlement density in that year, we also used the BUDENS layers for the respective years (shown for 1900 and 2015 in Fig. 15c, d), to weight each grid cell by its building density. We did this because we assumed that misclassifications are more likely in sparsely populated areas, likely not contained in the urban extents from the historical urban atlases due to generalization. Thus, agreement metrics based on building counts rather than grid cell counts would be more representative and realistic for a comparison between these two datasets. Based on this evaluation strategy, we observe the following:

1. Agreement levels are generally relatively low (F1 between 0.28 and 0.74). This may be attributed to the survivorship bias in HISDAC-ES but also due to definitional differences; i.e., HISDAC-ES measures built-up area, while the urban extents derived from historical atlases measure urban area and, thus, are already a generalized representation of the developed land at a given point in time. They are likely to omit low-density, scattered, peri-urban settlements but include roads, impervious surfaces, and intra-urban green spaces (e.g., parks and cemeteries), which are not directly measured in HISDAC-ES, which is based on the presence of built structures only.

2. Recall is higher than precision. Low precision indicates high commission error (i.e., overestimation), likely because peri-urban, rural settlements are contained in HISDAC-ES but not in the urban extents due to generalization and the (arbitrary) definition of the urban boundaries in the historical urban atlases, similar to what we observed in the comparison to the GHSL (Sect. 4.1).

3. Agreement for cumulative urban extents is higher than for incremental time slices. This effect is to be expected, as the confusion between historical increments is irrelevant when comparing the total built-up/urban area at each point in time. As the urban areas (and the increments) in early time periods can be small, misclassification is more likely, also due to higher levels of survivorship bias in HISDAC-ES for early time periods.

4. Precision based on the number of buildings is higher than precision based on the number of “occupied” grid cells. This indicated that grid cells label ed as “built-up” in HISDAC-ES but not in the historical urban areas tend to have low building density, confirming our aforementioned assumption that low-density settlements are not mapped in the historical urban extents.

5. For cumulative urban extents, precision and recall increase over time. This is a direct effect of the survivorship bias, manifesting in higher omission errors (and thus lower recall) as we go back in time. Moreover, the proportion of scattered low-density settlements (which are not contained in the historical urban atlases) in relation to dense, urban settlements was higher in early than in recent epochs, resulting in an increase of precision over time.

Moreover, the agreement levels are relatively similar across the three cities under study, implying that these observations are likely to be generalizable at least across the larger cities in Spain. Despite the low absolute numbers of the agreement metrics reported in Table 4, which are likely due to definitional differences, the observed trends are in line with theoretical expectations (e.g., survivorship bias decreases over time) and with evaluation results of the HISDAC-US (Uhl et al., 2021c), indicating similar characteristics of the historical settlement layers derived from cadastral/property data in the United States and in Spain. Confusion matrices underlying the agreement metrics shown in Table 4 can be found in the Appendix Fig. E1.

Figure 15Data used for the quantitative comparison of the HISDAC-ES to historical urban extents, derived from historical maps for the cities of Alicante (top row), Madrid (middle row), and Valencia (bottom row): (a) historical urban areas digitized from historical maps after rasterization, (b) the HISDAC-ES MINCOY layer, and the HISDAC-ES building density layers for (c) 1900 and (d) 2020. Data source of the historical urban extents for Valencia: courtesy of Carmen Zornoza-Gallego. Figure in the lower left corner adapted from Zornoza-Gallego (2022a).

Table 4Map comparison results of urban areas digitized from historical maps and HISDAC-ES settlement age surface (MINCOY), for both increments (i.e., newly developed grid cells between two points in time) and cumulatively (i.e., for the total developed land at a given point in time. Note that there is a slight temporal gap in the earliest Madrid epoch, where MINCOY from 1900 is compared to urban extents from 1892.

Download Print Version | Download XLSX

The qualitative comparison of the DEVA layers and planimetric historical maps (predecessors of the Minutas catastrones from 1870–1950, at scale 1:50 000) confirms the previously made observations. As shown in Fig. 16, the DEVA layer (for 1900) mimics the urban areas as depicted in the maps quite well, with some omission errors mainly in peri-urban areas, e.g., for Madrid, Seville, and Terrassa. This could indicate that building replacements (causing survivorship bias in the HISDAC-ES) tend to occur least in the medieval city centers, which are subject to monument protection (assuming concentric growth patterns). Moreover, these disagreement patterns may also be due to temporal inconsistencies between DEVA (1900) and the planimetric maps, which may have slightly different temporal references. In the small villages around Hornillos del Camino (Burgos) (Fig. 16 right column) we rather observe over- than underestimation. This could indicate that in rural, economically less prosperous areas, where less building remodeling occurs, the original building stock is still dominating, and, thus, the HISDAC-ES is less affected by survivorship bias. This observation may imply higher levels of data quality in small, rural places, a promising perspective for long-term settlement modeling in the often under-studied rural settings. The bottom row (Fig. 16c) shows the contemporary (i.e., 2020) building densities, illustrating a positive association between building density and settlement age.

Figure 16Visual comparison of HISDAC-ES developed areas (i.e., DEVA layer in 1900) with historical maps. (a) “Planimetría” maps (b) overlaid with DEVA in 1900 and (c) contemporary building density layer (BUDENS) in 2020), illustrating the change between 1900 and 2020. Each layer is shown for Madrid, Seville, Terrassa (Catalonia), and Hornillos del Camino (Burgos). Historical map source: Instituto Nacional de Geografía (2022).

4.5 Comparison to historical orthophotos from 1956

The urban extents manually digitized from historical orthophotos acquired in 1956 show the urban boundary in those years in Alicante, Madrid, and Santiago de Compostela (Fig. 17a, b, c). When we overlay these boundaries to the built-up areas from HISDAC-ES in the same year, we observe varying levels of agreement: In Alicante the agreement is relatively high (Fig. 17d), as well as in Santiago (Fig. 17f), whereas we observe higher levels of disagreement in the Madrid data (Fig. 17e), mainly occurring in suburban areas. While some of the disagreement may be attributed to differences in definitions (i.e., the urban boundaries drawn in the orthophotos only include dense, urban settlements), in the Madrid case there are additional issues, related to notable activities of urban renewal in previously informal settlements (e.g., the Entrevías neighborhood; example 2 in Fig. 17e). See Fig. F1 for a more detailed discussion of historical reasons for these discrepancies.

Figure 17Visual comparison of HISDAC-ES to aerial imagery from 1956. Aerial image overlaid with manually digitized urban extents for (a) Alicante, (b) Madrid, and (c) Santiago de Compostela and (d)–(f) overlaid with the HISDAC-ES built-up grid cells after and prior to 1955. Dashed black boxes are selected areas of discrepancies (i.e., turquoise areas within the digitized urban boundaries or red areas without the urban boundary), which are enlarged and discussed in the Appendix Fig. A6.

4.6 Attribute completeness and coverage

Lastly, we report the INSPIRE building attribute completeness at the municipality level, compared to the total number of buildings available in each municipality (Fig. 18a). We observe very high levels of completeness of the construction year attribute (Fig. 18b). The building function attribute has a high coverage, except in the Gipuzkoa and Bizkaia provinces on the northern coast (Fig. 18c). Moreover, there is a building function attribute, “other”, that is only available in the Navarre region, which we excluded from the HISDAC-ES dataset (Fig. 18d). Thus, the uncertainty in the (historical) land use layers in HISDAC-ES in Navarre is slightly higher, as it is unknown what building function the “other” class encompasses. The indoor area, number of dwellings, and number of building units attributes also have lower levels of completeness in some areas of the Basque Country (Fig. 18e, f, h).The completeness of “number of dwellings” is higher in buildings labeled as “residential” (Fig. 18g) than across all buildings (Fig. 18f), as this attribute is semantically linked to (mostly) residentially used buildings. Conversely, information on the number of floors is highly complete in the Navarre region (Fig. 18j) but otherwise not covered in the remaining provinces and, thus, has not been used in this version of HISDAC-ES.

Figure 18Attribute completeness and construction year coverage at the municipality level. See the Appendix Fig. G1 for corresponding maps of the Canary Islands.

We also assessed the temporal coverage of construction year information at the municipality level, in order to better understand potential survivorship bias in the data. As can be seen in Fig. 18k, most municipalities have the earliest construction year on record <1850 or <1900. Only a few regions have minimum construction years between 1900 and 1925 (e.g., regions around San Sebastián and Bilbao), whereas very few, scattered municipalities have earliest construction years between 1925 and 1950. In these municipalities, data users should be careful when conducting long-term analyses, as survivorship bias may be high. Likewise, we mapped out the maximum construction year per municipality (Fig. 18l), indicating the recentness of the cadastral building data underlying HISDAC-ES, as in most municipalities, the most recent construction year on record is between 2015 and 2020. Generally, these high levels of attribute completeness and temporal coverage are encouraging and indicate that the layers derived from these attributes are expected to be highly reliable at least for recent points in time. We made these municipality-level attribute completeness statistics available in the data repository.

Table 5HISDAC-ES dataset overview.

Download Print Version | Download XLSX