1 Introduction↩︎

Understanding how people move within and between urban and rural areas is crucial for addressing a wide range of societal challenges. Human mobility data inform research and policy decisions across multiple domains, including transportation planning and city optimization [1]–[3], epidemic modeling [4]–[9], inequalities and segregation [10]–[14], and disaster response [15]–[18]. In each of these settings, mobility is not only a measure of movement volume, but also a relational structure linking places, populations, and opportunities.

Despite the ubiquity of digital traces generated by mobile phones, GPS devices, and connected applications, access to reliable and representative mobility data remains limited [8], [19], [20]. High-resolution datasets are often proprietary, expensive, or subject to restrictive licensing agreements that limit their use in scientific research and public sector applications [20]. When data are publicly available, they are often fragmented, poorly documented, or difficult to integrate into scalable and reproducible workflows [20]. These barriers are especially important for nationwide mobility datasets, where repeated downloading, spatial harmonization, and processing of large daily files can become a substantial obstacle even when the data themselves are open.

To address part of this gap, the Spanish Ministry of Transport and Sustainable Mobility began releasing detailed daily mobility data derived from mobile phone signals [21]–[25]. This initiative began in response to the COVID-19 pandemic, with an initial collection phase spanning from February 2020 to May 2021. A second release phase, starting in January 2022 and continuing to the present, reflects an ongoing effort to improve transparency and enable data-driven decision-making. The released datasets include aggregated origin-destination matrices, trip-count indicators, and overnight stay locations, available at three spatial aggregation levels: districts, municipalities, and greater urban areas. The data are updated daily and published with a short delay [21]–[25].

However, open publication alone does not guarantee scientific usability. Researchers and practitioners still need to identify the correct files, retrieve them repeatedly across dates and spatial scales, harmonize spatial identifiers, handle geometries, convert raw files into analysis-ready formats, and document the workflow in a reproducible way. These tasks are technical but consequential: without standardized access and preprocessing, analyses become harder to replicate, compare, and extend. This is particularly limiting for network-based mobility research, where origin-destination matrices must often be transformed into weighted graphs, stratified by population attributes, and analysed across multiple time periods.

The paper makes two contributions. First, it describes the design, data model, and functionality of pySpainMobility, showing how the package organizes the main spatial units and mobility data types released by the Ministry. Second, it demonstrates the analytical value of the package through a use case on income-stratified mobility inequality in Spain’s inter-province mobility network. Using district-level origin-destination data for four representative weeks in 2023, we construct income-specific mobility layers and analyse how income groups differ in destination concentration, destination hierarchy, spatial reach, and destination-level income selectivity. The use case shows that the package enables more than data retrieval: it supports reproducible, national-scale, network-based analysis of social and spatial inequality using official open mobility data.

The rest of the paper is structured as follows. Section 2 reviews the broader landscape of mobility data sources. Section 3 introduces the study areas used in the Spanish datasets and explains how they are managed in the library. Section 4 details the structure and technical usage of the main data types accessible through the package. Section 5 presents the income-stratified mobility use case, including network construction, measures, results, and discussion. Section 6 concludes with a summary of contributions, limitations, and directions for future work. Code, documentation, and notebooks to reproduce all analyses are available at github.com/pySpainMobility and pySpainMobility.github.io/pySpainMobility/.

2 Mobility Data Sources↩︎

Before presenting the design of pySpainMobility, it is useful to situate the Spanish datasets within the broader landscape of mobility data sources. Human movement can be observed through several technologies, each with distinct trade-offs in spatial accuracy, temporal resolution, population coverage, and accessibility, and these trade-offs shape both the analyses that are feasible and the reliability of their results.

The most widely used source is mobile network data, and in particular Call Detail Records (CDRs), which log the time and serving cell tower of calls and text messages. Because location is inferred from tower positions and recorded only when an event occurs, CDRs are spatially coarse and temporally sparse, which limits detailed trajectory reconstruction [26]–[28]. eXtended Detail Records (XDRs), generated by mobile-internet usage, mitigate this sparsity by sampling locations more frequently and more continuously [29], [30].

GPS data collected through application SDKs or onboard navigation systems are far more accurate in space and time, and several providers released such data during the COVID-19 pandemic [31]–[33]. Their main drawbacks are limited transparency about how and from whom data are collected, and known representativeness issues, since the social profiles of users (e.g., age or income) are rarely disclosed. Aggregated indicators published by technology firms such as Google and Facebook are broadly accessible and were widely used to track behavioral change during lockdowns, but they offer lower spatial and temporal granularity than raw traces [8].

Other sources are complementary. Geotagged social-media posts and check-ins are sporadic and user-dependent but useful for tourism analysis and event detection [34]; credit-card transactions capture spatially located spending and movement between residence, work, and consumption [35]; satellite imagery supports urban-expansion analysis when combined with mobility indicators [36]; and wearables add physiological and activity data for travel-behavior and health-related studies [37]. Traditional census data remain an essential baseline for population and commuting statistics but, given their low update frequency, are poorly suited to real-time or event-driven analysis [8].

Despite this proliferation of sources, a central limitation persists: the scarcity of open and standardized datasets that are comprehensive in scope and reproducible in methodology [20]. Many datasets are restricted to particular platforms, modes, or regions and rely on opaque preprocessing pipelines whose undocumented choices (stop-detection thresholds, clustering parameters, sampling criteria) materially affect the derived indicators and undermine comparability across studies. The relative absence of shared benchmarks further slows methodological progress compared with fields such as computer vision or natural language processing [20].

The datasets released by the Spanish Ministry of Transport and Sustainable Mobility address part of this gap. They provide daily, nationwide indicators derived from mobile-network data for more than 13 million users, at three levels of spatial aggregation [21], [22], [24], [25]. Data for 2020–2021 (version 1) are based on CDRs, whereas from 2022 onward (version 2) a passive-measurement approach records devices attached to the network even when inactive, not only during calls or data use, which substantially improves data quality [38]. The provider’s market share exceeds 30%, and correction procedures are applied to under-represented groups such as children and the elderly, so that only limited segments of the population are subject to self-selection bias [21]. Open publication, however, does not by itself guarantee scientific usability: users must still locate and repeatedly download the correct files across dates and spatial scales, harmonize spatial identifiers, handle geometries, and convert raw files into analysis-ready formats. pySpainMobility responds to these challenges, a mission shared by parallel community efforts in R [39], through a Python-first workflow designed for scalable, reproducible, and analysis-ready access to Spanish mobility data by:

• generating analysis-ready mobility layers by transforming raw demographic records into standardized, stratified origin–destination matrices, including dimensions such as income, age, and gender;

• supporting national-scale longitudinal workflows through an Apache Arrow backend and compressed Parquet caching, making large daily files and multi-season analyses tractable within standard Python environments;

• unifying spatial and network pipelines by linking mobility flows with official geometries across multiple territorial scales, from districts and municipalities to greater urban areas.

By integrating data retrieval, demographic stratification, efficient storage, and spatial interoperability, pySpainMobility provides a reproducible bridge between official mobility records and network-based urban research.

3 Study Areas↩︎

The Zones module inpySpainMobilityenables programmatic access to the spatial aggregations used in the official Spanish mobility datasets. These study areas divide the national territory into discrete, non-overlapping zones that serve as the fundamental spatial units for aggregating mobility indicators. The three available levels of granularity are districts (distritos), municipalities (municipios), and greater urban areas (grandes áreas urbanas).

These spatial definitions are essential for interpreting the data, and their configuration affects the analysis outcomes. For example, districts are the finest spatial resolution and tend to be larger in rural areas due to lower population density, while municipalities and greater urban areas are more heterogeneous and can better represent urban-rural transitions.

The Zones object supports flexible instantiation. Users must specify the desired spatial granularity using the zones parameter, which supports various aliases (e.g., ‘municipalities’, ‘muni’, ‘gau’). Additionally, the version parameter allows alignment with the temporal availability of the datasets: Version 1 covers data between February 2020 and May 2021 (excluding greater urban areas), while Version 2 includes data from January 2022 to June 2025. The study areas are retrieved from official sources and can be stored locally in GeoJSON format for further use.

To retrieve and load a study area as a GeoDataFrame, the user can run the code shown in Listing 1:

from pySpainMobility import Zones

zones = Zones(zones='municipalities', version=2, output_directory='data')

gdf = zones.get_zone_geodataframe()

gdp.plot()

This command downloads the relevant geometries and returns them as a GeoPandas object, allowing immediate spatial manipulation and visualization. Figure 2 illustrates the differences between the three levels of aggregation and provides a zoomed-in view of the Madrid area, where spatial contrasts are most evident.

The average surface area of the zones varies across study area types. The average area size is 133.56 square kilometers for districts, 193.45 square kilometers for municipalities, and 242.79 square kilometers for Greater Urban Areas

These values underline the substantial variation in spatial granularity across the country. In sparsely populated rural zones, even districts may span large areas, while in cities, the differences in resolution become more pronounced.

In addition to the geometries, the module provides access to a mapping table that allows interoperability between spatial aggregations and alignment with official statistics. This is particularly important when merging external data sources or when transitioning between datasets at different levels of aggregation.

To access the mapping table, the user can call the get_zone_relations() method, as shown in Listing 2:

from pySpainMobility import Zones

zones = Zones(zones='municipalities', version=2, output_directory='data')

df = zones.get_zone_relations()

The returned DataFrame links identifiers across spatial schemes, allowing users to convert district codes to municipality or urban area codes, align mobility zones with official census references, and aggregate or disaggregate data across spatial hierarchies. This functionality ensures compatibility with external datasets and supports complex geographic analyses requiring hierarchical spatial logic.

4 Mobility Data↩︎

The Mobility module inpySpainMobilityprovides a standardized interface for downloading and preprocessing the main datasets made available by the Spanish Ministry of Transport and Sustainable Mobility. These datasets are published under the national open data strategy and contain anonymized information derived from more than 13 million mobile phone users. The module allows access to daily mobility indicators in a programmatic and reproducible way, offering support for both historical (Version 1: February 2020 – May 2021 based on CDRs) and recent (Version 2: from January 2022 onward based on passive measurements) data collections.

Upon instantiation, the user specifies the version of the dataset, the spatial aggregation level (districts, municipalities, or greater urban areas), and the temporal range of interest. The zones parameter accepts various aliases (e.g., ‘distritos’, ‘municipios’, ‘gau’) for user flexibility. The start date is mandatory, while the end date, if not specified, will be equal to the start date. All data are downloaded directly from the Ministry’s open data portal and stored in Parquet format within the designated output directory. To keep national-scale processing tractable, the module parses the large daily files through an Apache Arrow backend by default, with an automatic fallback to pandas and optional Dask support for out-of-core workloads, and caches downloaded data as compressed Parquet to avoid repeated transfers.

The main functionalities of the module are described below.

4.1 Origin–Destination Matrices↩︎

Origin–destination (OD) matrices represent the volume and characteristics of daily flows between pairs of geographic zones. Each record captures aggregated movements from an origin to a destination and includes auxiliary attributes such as traveled distance, inferred activity at each location (e.g., home, work/study, frequent visit), and demographic indicators. These include age (categorized into 0–24, 25–44, 45–64, 65+, or not disaggregated), gender (male, female), and income level (\(<\euro 10{,}000\), \(\euro 10{,}000\)–\(15{,}000\), and \(>\euro 15{,}000\)).

Only trips longer than 500 meters are included, and all movements refer to Spanish residents. The data are particularly valuable for characterizing commuting behavior, inter-regional connectivity, and socio-demographic patterns of mobility.

To retrieve OD data, the user can rely on the get_od_data() method. An example is shown in Listing 3.

from pySpainMobility import Mobility

mobility_data = Mobility(version=2, zones='municipalities', start_date='2022-03-20', end_date='2022-03-24', output_directory='./data/')

mobility_data.get_od_data(keep_activity=True)

The resulting file is saved as a compressed Parquet dataset, enabling scalable analysis with Python data tools.

4.2 Number of Trips per Person↩︎

In addition to origin–destination data, the Mobility module allows the extraction of records quantifying the number of trips made by individuals residing in a specific zone. Each entry includes the overnight stay zone, the number of individuals belonging to each demographic group, and the count of trips taken (grouped as 0, 1, 2, or more than 2 trips per day).

This dataset supports behavioral mobility profiling and provides insights into activity intensity across populations and regions. It is particularly suited for temporal comparisons, peak detection, and studies on the relationship between mobility and socioeconomic conditions.

The method get_number_of_trips_data() can be used as shown in Listing 4.

from pySpainMobility import Mobility

mobility_data = Mobility(version=2, zones='municipalities', start_date='2022-03-20', end_date='2022-03-24', output_directory='./data/')

mobility_data.get_number_of_trips_data()

The output is similarly saved in a Parquet file and aligned with the other mobility indicators.

4.3 Overnight Stays↩︎

The third available dataset focuses on overnight stays and captures, for each date, the location where individuals spend the night, based on their area of residence. Each record identifies the residence zone, the overnight location, and the number of individuals involved. The data reflect temporary relocation behaviors, such as tourism flows, weekend displacements, or prolonged visits.

This dataset provides a complementary angle to the trip-based indicators. Rather than tracking movement events, it captures the spatial distribution of presence during the nighttime, making it ideal for estimating short-term population concentrations, especially during holidays or seasonal periods.

To download the overnight stays dataset, the user can rely on the following command.

from pySpainMobility import Mobility

mobility_data = Mobility(version=2, zones='municipalities', start_date='2022-03-20', end_date='2022-03-24', output_directory='./data/')

mobility_data.get_overnight_stays_data()

As with other methods, the output is saved to the specified directory and formatted for immediate use in data analysis environments.

5 Use Case: Income-stratified Mobility Networks and Seasonal Inequality↩︎

Urban and regional inequalities are often studied through the spatial distribution of residence, income, employment, or access to services [13], [40]–[43]. Mobility adds a relational dimension to this problem. Places are not only unequal because of their internal socioeconomic composition, but also because they are connected to different populations through flows [44]–[46]. From a network perspective, inequality can therefore be expressed in the structure of mobility itself [10], [31], [40], [47], [48].

Despite the relevance of this perspective, empirical analyses of income-stratified mobility networks remain difficult to reproduce at national scale. They require access to large origin-destination matrices, consistent spatial units, demographic attributes, and a processing pipeline capable of transforming raw daily files into comparable network layers. This is precisely the type of workflow that pySpainMobility is designed to support. By standardizing access to the Spanish open mobility data, the package makes it possible to construct income-specific mobility networks with a transparent and replicable procedure, reducing the technical overhead that often separates open data availability from actual scientific use.

The Spanish mobility data provide a particularly useful setting because origin-destination flows are disaggregated by income category. We exploit this feature to represent the national mobility system as a set of income-specific layers defined over the same geography. In this multilayer representation, each income group generates a directed weighted mobility network. Comparing these layers allows us to ask whether income groups participate in a common spatial system or whether they move through distinct destination geographies.

We use this setting to study income-stratified inter-province mobility across four representative weeks in 2023. Seasonality is not treated as a standalone descriptive object, but as a temporal perturbation of the mobility system. If income-based differences were only a matter of volume, seasonal changes would mainly rescale the layers without altering their relative structure. Conversely, persistent differences in concentration, destination ranking, and destination composition would indicate that income groups occupy distinct positions in the mobility network.

The use case is organized around three linked analytical questions. First, we ask whether income groups differ in the concentration of their inter-province destinations. Second, we examine whether income-specific layers share the same destination hierarchy or instead reveal a separation between low-, middle-, and high-income mobility geographies. Third, we measure whether destination provinces receive income compositions that deviate from the national seasonal baseline, indicating destination-level income selectivity. Together, these analyses show how pySpainMobility can be used not only to retrieve official mobility data, but also to operationalize inequality as a network property: a difference in connectivity, destination structure, and social composition of flows.

5.1 Network construction↩︎

We construct the use case from district-level origin-destination flows for four representative weeks in 2023: 16–22 January, 10–16 April, 17–23 July, and 9–15 October. These weeks are used as seasonal snapshots of winter, spring, summer, and autumn, respectively. For each week, pySpainMobility is used to retrieve the corresponding daily mobility files, harmonize the spatial identifiers, and store the processed data in a format suitable for graph construction.

For each season, we define a directed weighted mobility network \(G_t=(V,E_t)\), where \(V\) is the set of districts and an edge \((i,j)\in E_t\) represents trips from origin district \(i\) to destination district \(j\) during season \(t\). Edge weights correspond to the total number of trips observed during the representative week. Self-loops are excluded, so the analysis focuses on mobility between distinct districts. Indeed, the main analysis focuses on inter-province mobility, defined administratively as flows whose origin and destination districts belong to different provinces.

Income groups are represented as layers of the same mobility system. We construct one aggregate layer and three income-specific layers, corresponding to the income categories available in the source data: low income (\(<10{,}000\) euros), middle income (\(10{,}000\)–\(15{,}000\) euros), and high income (\(>15{,}000\) euros). Each income layer is therefore a directed weighted graph defined over the same spatial system, but with edge weights restricted to trips generated by the corresponding income group.

Although the graphs are constructed at the district level, the main empirical results are reported at province level. We aggregate district-to-district flows into province-to-province flows by summing all edge weights whose origin districts belong to the same origin province and whose destination districts belong to the same destination province. This aggregation improves interpretability for national-scale comparison while preserving the income-layer structure of the mobility network. The resulting province-level network is used to analyse inbound destination shares, income-specific destination hierarchies, and destination-level income selectivity.

The full processing pipeline is retained as part of the reproducible workflow. Supplementary tables report the number of nodes, active edges, and total edge weights by season and income layer, together with additional processing checks (Appendix 7.1).

5.2 Measuring Income-stratified Mobility Inequality↩︎

We quantify income-stratified mobility inequality through three complementary network measures. Each measure captures a different aspect of how income layers are organized in the inter-province mobility network: the concentration of destinations, the similarity between income-specific destination hierarchies, and the income selectivity of destination provinces.

First, we measure the concentration of destinations within each income layer. Let \(w_{jgt}^{in}\) denote the total inbound inter-province flow received by destination province \(j\) from income group \(g\) in season \(t\). The destination share of province \(j\) for income group \(g\) is

\[s_{jgt} = \frac{w_{jgt}^{in}}{\sum_{j} w_{jgt}^{in}}.\]

We then compute concentration ratios over the ranked distribution of destination shares [49]. For a given income group and season, the concentration ratio \(CR_k(g,t)\) is defined as

\[CR_k(g,t) = \sum_{j \in \mathcal{T}_k(g,t)} s_{jgt},\]

where \(\mathcal{T}_k(g,t)\) is the set of the \(k\) destination provinces with the largest inbound shares for income group \(g\) in season \(t\). High values of \(CR_k\) indicate that an income layer concentrates a large fraction of its inter-province mobility in a small number of destinations.

Second, we compare the destination hierarchies of different income layers. For each season, provinces are ranked according to their inbound share within each income group. We then compute Spearman rank correlations [50] between pairs of income layers. If two income groups have a high correlation, they prioritize a similar set of destination provinces. Low correlations indicate that the groups are embedded in different destination hierarchies, even though they move through the same national spatial system.

Third, we measure the spatial reach of each income layer, following the notion of a characteristic travel distance in human mobility [26]. For every inter-province origin-destination pair with available geographic matches, we compute the distance between the representative points of the origin and destination districts. Let \(d_{ij}\) be the distance between origin district \(i\) and destination district \(j\), and let \(w_{ijgt}\) be the corresponding number of trips for income group \(g\) in season \(t\). The weighted mean distance of income layer \(g\) in season \(t\) is

\[\bar{d}_{gt} = \frac{\sum_{(i,j)} w_{ijgt} d_{ij}}{\sum_{(i,j)} w_{ijgt}},\]

where the sum is taken over inter-province flows with valid geographic coordinates. This measure captures the average distance traveled by flows, weighting each origin-destination pair by its observed trip volume.

Finally, we introduce a destination income selectivity measure. This measure asks whether a province receives an income composition of inbound travelers that differs from the national seasonal composition of inter-province mobility. For each destination province \(j\), we define the income composition of inbound flows as

\[p_{jgt} = \frac{w_{jgt}^{in}}{\sum_g w_{jgt}^{in}},\]

and the national seasonal income composition as

\[p_{gt}^{nat} = \frac{\sum_j w_{jgt}^{in}}{\sum_j \sum_g w_{jgt}^{in}}.\]

We compare the destination-specific distribution \(p_{jt}\) with the national baseline \(p_t^{nat}\) using Jensen–Shannon divergence [51]:

\[JSD_{jt} = JSD(p_{jt}, p_t^{nat}).\]

A value close to zero indicates that the destination receives an income mix similar to the national seasonal baseline. Higher values indicate stronger destination income selectivity, meaning that the province attracts an inbound mobility composition disproportionately associated with one or more income groups. To characterize the direction of selectivity, we also identify the income group most overrepresented relative to the national baseline.

Finally, we compute a normalized income-mixing entropy [52]–[54] for each destination,

\[H_{jt} = -\frac{1}{\log 3}\sum_g p_{jgt}\log p_{jgt},\]

where \(H_{jt}=1\) corresponds to a perfectly balanced distribution across the three income groups, while lower values indicate stronger dominance by one group. In the main text, entropy is used as a complementary measure to the Jensen–Shannon index. Full tables for all concentration ratios, destination selectivity values, and robustness checks are reported in Appendix 7.2 and Appendix 7.3.

5.3 Results↩︎

The income-stratified networks reveal that seasonal mobility inequality is not only a matter of different mobility volumes. Income groups differ in how concentrated their destinations are, in how similarly they rank destination provinces, and in the income composition of the provinces they reach. We present these results in three steps, moving from the structure of income-specific layers to the social selectivity of destinations.

5.3.1 Unequal Destination Concentration and Hierarchy↩︎

We first examine whether income groups are connected to the same destination geography. For each season and income layer, we rank destination provinces by inbound inter-province flow. We then summarize this distribution in two ways. First, \(CR_{10}\) measures the share of each income group’s mobility captured by its ten largest destination provinces. Higher values indicate a more concentrated mobility layer. Second, Spearman rank correlations compare the destination rankings of different income layers. High correlations indicate that two income groups assign similar relative importance to destination provinces, while low correlations indicate distinct destination hierarchies.

Table 1 shows a clear income gradient in destination concentration. The low-income layer is the most concentrated across its ten largest destinations in every season: its ten largest destination provinces capture between \(73.6\%\) and \(78.3\%\) of inter-province inbound mobility. Middle-income mobility is substantially less concentrated, while high-income mobility occupies an intermediate position in most seasons. This means that low-income inter-province mobility is organized around a narrower set of destinations, rather than being simply a lower-volume version of the other layers.

The rank correlations show that income layers also differ in the destinations they prioritize. The low- and high-income layers are weakly aligned in all seasons, with correlations between \(0.151\) and \(0.228\). Low- and middle-income rankings are more similar, but only moderately so, ranging from \(0.404\) to \(0.566\). By contrast, middle- and high-income hierarchies are strongly aligned, with correlations above \(0.78\) in every season. The destination structure of the low-income layer is therefore distinct from the rest of the mobility system, whereas middle- and high-income mobility follow much more similar spatial hierarchies.

Together, these results indicate that income stratification is embedded in the structure of the inter-province mobility network. Low-income mobility is both more concentrated and organized around a different destination hierarchy, implying a narrower effective geography of inter-province connections. Full concentration ratios, including \(CR_1\), \(CR_3\), \(CR_5\), and \(CR_{10}\), and complete pairwise similarity values are reported in Appendix 7.2.

5.3.2 Unequal Spatial Reach↩︎

The previous results show that income layers differ in the destinations they prioritize. We next examine whether they also differ in spatial reach. For each income layer and season, we compute the weighted mean distance of inter-province trips. Distances are calculated between representative points of origin and destination districts for origin-destination pairs with available geographic matches, and are weighted by trip volume. The measure therefore captures the average distance traveled by observed flows, rather than the average distance between active district pairs.

Table 2 reports weighted mean distances across the four seasonal snapshots. The high-income layer has the greatest spatial reach in every season. In winter, high-income inter-province trips have a weighted mean distance of \(154.6\) km, compared with \(113.5\) km for the middle-income layer and \(98.6\) km for the low-income layer. The same ordering holds in spring, summer, and autumn, with high-income flows consistently traveling farther than middle- and low-income flows.

This result adds a second dimension to the inequality pattern. Low-income mobility is not only more concentrated across destinations; it also has a shorter spatial reach. Conversely, the high-income layer combines a broader destination geography with longer trip distances. Income stratification therefore affects both where groups travel and how far their inter-province mobility extends. The same ordering is visible in the weighted medians reported in Appendix 7.4, which reduces the risk that the result is driven only by long-distance tails. Full distance summaries, including median and interquartile distances, together with geographic match-rate checks, are reported in Appendix 7.4.

5.3.3 Destination Income Selectivity↩︎

The previous results characterize inequality from the perspective of income layers: some groups concentrate their mobility in fewer destinations and follow different destination hierarchies with diverse spatial reaches. We now reverse the perspective and examine destinations themselves. For each season, we first compute the national income composition of inter-province mobility, that is, the share of all inter-province trips generated by the low-, middle-, and high-income groups. We then compute the same income composition separately for the inbound flows received by each destination province. Destination income selectivity is measured as the Jensen–Shannon divergence [51] between these two distributions. A province is therefore more selective when the income mix of its arrivals differs more strongly from the national seasonal baseline.

Figure 3 maps this measure across seasons. Each province is colored according to the income group most overrepresented among its inbound inter-province arrivals relative to the national seasonal baseline, while shade intensity reflects the magnitude of the Jensen–Shannon divergence. The geography of selectivity is not uniform. Some provinces repeatedly appear as low-income selective destinations, while others attract a disproportionately high-income composition of arrivals. For example, Almeria and Murcia display strong low-income overrepresentation in several seasons, whereas Basque provinces such as Gipuzkoa, Araba/Alava, and Bizkaia stand out as high-income selective in winter.

At the national level, destination income selectivity is stronger in winter and spring than in summer and autumn. The weighted mean Jensen–Shannon divergence is \(0.052\) in winter and \(0.054\) in spring, compared with \(0.035\) in summer and \(0.039\) in autumn. The complementary income-mixing entropy follows the opposite pattern, reaching its highest value in summer. This suggests that the summer expansion of inter-province mobility is associated with a more socially mixed geography of arrivals, whereas winter and spring mobility are more socially selective.

This result adds a destination-level dimension to the multilayer interpretation. Income-stratified mobility inequality is not only visible in the structure of the layers generated by each income group. It is also visible in the social composition of flows received by places. Provinces differ in whether they attract a mobility mix close to the national baseline or a more selective composition dominated by specific income groups. In this sense, destination income selectivity provides a network-based measure of mobility sorting: it identifies where the mobility system connects places to socially distinctive populations.

Full province-level selectivity values, including the most overrepresented income group and income-specific deviations from the national baseline, are reported in Appendix 7.3.

6 Discussion and Conclusion↩︎

The income-stratified use case shows that pySpainMobility supports more than data retrieval: it enables a complete, reproducible pipeline from official open files to comparable network layers. We discuss the implications along two axes. Section 6.1 considers the methodological contribution, what the use case demonstrates about the package as an instrument for national-scale, network-based analysis. Section 6.2 turns to the substantive findings on income, mobility, and spatial inequality in Spain.

6.1 Methodological Contribution↩︎

Beyond its substantive findings, the use case illustrates the methodological value of pySpainMobility. Constructing income-specific mobility layers at national scale requires retrieving large daily origin-destination files across four seasonal windows, harmonizing district identifiers against the official geometry, filtering by income category, and assembling directed weighted graphs that are comparable across time and groups. Performed manually, this pipeline is error-prone and difficult to reproduce; the package reduces it to a small, documented set of calls and stores analysis-ready outputs that feed directly into standard network tools.

This reproducibility is what makes the multilayer, multi-season design feasible in the first place. Because access, spatial harmonization, and preprocessing are standardized, the same workflow can be re-run for the other demographic dimensions available in the data (age and gender), for other spatial aggregations (municipalities or greater urban areas), for other time windows, and for alternative network measures of accessibility, exposure, or segregation. The contribution is therefore not a single analysis but a replicable template: it lowers the technical barrier that usually separates the release of official open mobility data from its substantive scientific use, and it makes national-scale, network-based studies accessible to researchers without a bespoke data-engineering pipeline.

6.2 Income, Mobility, and Spatial Inequality↩︎

The use case shows that income-stratified mobility inequality is not reducible to differences in aggregate trip volume. When income groups are represented as layers of the same inter-province mobility network, inequality becomes visible in the organization of connections: the low-income layer is more destination-concentrated and shorter-range, low- and high-income layers follow weakly aligned destination hierarchies, and several provinces receive income compositions that deviate from the national seasonal baseline. Income groups thus participate in partially distinct mobility geographies rather than simply moving more or less within a common one.

These findings connect to the literature in the social sciences that treats spatial inequality as a structural property of where groups are located and how they are distributed across space. Research on residential segregation has established both the multiple dimensions along which groups can be separated [54], [55] and the tight link between rising income inequality and the residential sorting of households by income [56]. Our results echo this literature from a relational angle: the pronounced concentration and shorter reach of the low-income layer, together with the broader and longer-range geography of high-income mobility, are the flow-based counterpart of the residential “segregation of affluence” [56], and our divergence-based selectivity measure adapts the same information-theoretic logic used to quantify residential segregation [54].

At the same time, our analysis contributes to a rapidly growing body of work arguing that segregation is not only a residential condition but also an experienced and relational one, produced through everyday movement [57], [58]. Studies using large-scale mobility data have shown that daily travel connects—or fails to connect—neighborhoods of different socioeconomic composition [59], [60], that experienced income segregation is embedded in urban mobility patterns [13], that higher-status individuals typically reach more distant and more diverse destinations [61], and that large, dense cities can paradoxically sustain higher exposure segregation by offering more differentiated estinations [48]. Our results extend this agenda in two respects. First, they move the analysis from the intra-urban, single-city, and predominantly US-based settings that dominate this literature to a national, inter-province scale built entirely on official open data. Second, they show that the weak alignment between low- and high-income destination hierarchies persists across all four seasons: although summer broadens the social mix of arrivals (higher income-mixing entropy and lower selectivity), it does not dissolve the structural separation between income layers. Seasonality modulates the intensity of mobility sorting without eliminating it.

Limitations should be noted. The income categories are those provided in the official files and are therefore coarse, allowing broad stratification but not a detailed analysis of the income distribution. The analysis uses four representative weeks rather than a continuous time series, which is appropriate for demonstrating the library and identifying clear patterns but leaves the temporal robustness of the results for future work. Finally, the measures describe observed flows, not the mechanisms that generate them: differences between income layers may reflect constraints, preferences, employment structures, family networks, or tourism, which cannot be disentangled from the mobility matrices alone. These limitations also mark natural extensions—finer socioeconomic strata, longer time series, intra-urban aggregations, and the integration of mobility layers with residential segregation measures—for which the standardized pySpainMobility workflow provides a direct foundation. Together, these results demonstrate how official open mobility data can be operationalized as a multilayer network for the national-scale study of spatial inequality, reinforcing the broader contribution of pySpainMobility: lowering the technical barrier to mobility analysis while enabling substantive research on networks and forms of inequality.

List of Abbreviations↩︎

Declarations↩︎

Ethics approval and consent to participate↩︎

Not applicable.

Consent for publication↩︎

All authors have read and approved the final manuscript and consent to its publication.

Availability of data and material↩︎

The dataset analysed in this work is publicly available via https://github.com/pySpainMobility/pySpainMobility.

Competing interests↩︎

The authors declare that they have no competing interests.

Funding↩︎

B.L. has been supported by PRECRISIS (PRotECting public spaces thRough Integrated Smarter Innovative Security), funded by the European Union Internal Security Fund (ISFP-2022-TFI-AG-PROTECT-02-101100539). B.L. and M.L. have been supported by the PNRR ICSC National Research Centre for High Performance Computing, Big Data and Quantum Computing (CN00000013), under the NRRP MUR program funded by the NextGenerationEU. B.L. also acknowledges the support of the PNRR project FAIR - Future AI Research (PE00000013), under the NRRP MUR program funded by the NextGenerationEU. B.L. and M.L. also acknowledge the support of the European Union’s Horizon Europe research and innovation programme under grant agreement No. 101120237 (ELIAS).

Authors’ contributions↩︎

C.B. designed and developed the package and performed the experiments. M.L. designed and developed the package and coordinated the writing of the paper. C.B., B.L., T.G.M.R., and M.L. contributed to the interpretation of the results and to writing the manuscript. All authors read and approved the final manuscript.

Acknowledgements↩︎

Not applicable.

7 Supplementary material↩︎

This appendix reports supplementary material for the use case presented in Section 5. It documents the graph-construction audit, provides full concentration and destination-hierarchy tables, reports supplementary destination income-selectivity results, and gives additional statistics for the spatial-reach analysis.

7.1 Processing audit↩︎

Table 3 summarizes the seasonal graph objects created for the use case before the province-level aggregation. The audit reports the number of active nodes, active directed edges, and total edge weight for each representative week.

A small number of origin-destination identifiers appearing in the official mobility files could not be matched to the released district geometry layer. Since this inconsistency originates from the official source data, we do not alter the origin-destination matrices. We retain the official flows for network construction and use geographic attributes only when a match is available.

7.2 Destination concentration and hierarchy↩︎

Table 4 extends the concentration results reported in Table 1. It reports \(CR_1\), \(CR_3\), \(CR_5\), and \(CR_{10}\), where \(CR_k\) is the percentage of each income group’s inbound inter-province mobility captured by its top \(k\) destination provinces.

Table 5 reports the full pairwise Spearman rank correlations between income-specific destination hierarchies. These values correspond to the similarity block in Table 1.

7.3 Destination income selectivity↩︎

Table 6 reports the national seasonal summary of destination income selectivity and income mixing. The weighted mean Jensen-Shannon divergence summarizes how strongly destination provinces deviate from the national seasonal income composition of inter-province mobility. The weighted mean entropy summarizes the income mixing of destination provinces.

Tables 1–3 report the destination provinces where each income layer is most overrepresented relative to the national seasonal baseline. Provinces are ranked by overrepresentation, measured as the difference between the income group’s share among arrivals to the destination and the same income group’s national share of inter-province mobility. The JSD column reports the overall income selectivity of the destination province.

7.4 Spatial reach by income layer↩︎

Table 10 reports supplementary distance statistics for the inter-province flows analysed in Section 5.3.2. Distances are computed between representative points of origin and destination districts for origin-destination pairs with available geographic matches. All statistics are weighted by trip volume, so they describe the distance distribution of observed mobility flows rather than the unweighted distribution of active origin-destination pairs. The table uses only the inter-province subset of flows, defined as trips whose origin and destination districts belong to different provinces.

The supplementary statistics confirm the pattern reported in the main text. In every season, the high-income layer has the largest weighted mean and median distance, while the low-income layer has the shortest spatial reach. The interquartile range also shows that high-income mobility has a longer upper tail, especially in spring, summer, and autumn.

References↩︎