from pyspainmobility.utils import utils
import pandas as pd
import geopandas as gpd
import os
from os.path import expanduser
from typing import Optional
[docs]
class Zones:
[docs]
def __init__(self, zones: str = 'municipalities', version: int = 2, output_directory: str = None):
"""
Class to handle the zoning related to the Spanish big mobility data. The class is used to download the data and
process it. Selectable granularities are districts (distritos), municipalities (municipios) and large urban areas (grandes áreas urbanas). As a reminder,
mobility data for the COVID-19 period (version 1) are not available for the large urban areas.
Parameters
----------
zones : str
The zones to download the data for. Default is municipalities. Zones must be one of the following: districts, dist, distr, distritos, municipalities, muni, municipal, municipios, lua, large_urban_areas, gau, gaus, grandes_areas_urbanas
``provinces`` dissolves mapped district geometries, excluding zones
without a unique province. These are derived MITMA geometries;
exclusions and geometry repairs are recorded in result ``attrs``.
version : int
The version of the data to download. Default is 2. Version must be 1 or 2. Version 1 contains the data from 2020 to 2021. Version 2 contains the data from 2022 onwards.
output_directory : str
The directory to save the raw data and the processed parquet. Default is None. If not specified, the data will be saved in a folder named 'data' in user's home directory.
Examples
--------
>>> from pyspainmobility import Zones
>>> # instantiate the object
>>> zones = Zones(zones='municipalities', version=2, output_directory='data')
>>> # get the geodataframe with the zones
>>> gdf = zones.get_zone_geodataframe()
>>> print(gdf.head())
name population
ID
01001 Alegría-Dulantzi 2925.0
01002 Amurrio 10307.0
01004_AM Artziniega agregacion de municipios 3005.0
01009_AM Asparrena agregacion de municipios 4599.0
"""
utils.version_assert(version)
utils.zone_assert(zones, version)
self.version = version
self.zones = utils.zone_normalization(zones)
self.source_zones = "distritos" if self.zones == "provinces" else self.zones
self.complete_df = None
self._zoning_links = None
self._downloads_ready = False
# Get the data directory
data_directory = utils.get_data_directory()
self.data_directory = data_directory
# Proper output directory handling
if output_directory is not None:
if os.path.isabs(output_directory):
# Preserve absolute paths even if parent directories do not exist yet.
self.output_path = output_directory
else:
# Treat relative paths as relative to home directory.
home = expanduser("~")
clean_path = output_directory.lstrip("/\\")
self.output_path = os.path.join(home, clean_path)
else:
self.output_path = data_directory
# Ensure directory exists
try:
os.makedirs(self.output_path, exist_ok=True)
except PermissionError as e:
raise PermissionError(f"Cannot create directory {self.output_path}. Please check permissions or use a different path. Error: {e}")
except Exception as e:
raise Exception(f"Error creating directory {self.output_path}: {e}")
def _get_zoning_links(self) -> list:
"""
Resolve available zoning links lazily.
"""
if self._zoning_links is None:
self._zoning_links = (
utils.available_zoning_data(self.version, getattr(self, "source_zones", self.zones))["link"]
.dropna()
.unique()
.tolist()
)
return self._zoning_links
def _ensure_zoning_files_downloaded(self, required=None) -> None:
"""
Download required files only when the user first requests data.
"""
if self._downloads_ready:
return
if required and all(
os.path.isfile(self._resolve_data_file(name))
and os.path.getsize(self._resolve_data_file(name)) > 0
for name in required
):
return
links = self._get_zoning_links()
for link in links:
file_name = link.split("/")[-1]
if required is not None and file_name not in required:
continue
local_path = os.path.join(self.output_path, file_name)
if (
not os.path.exists(local_path)
or os.path.getsize(local_path) == 0
or file_name.endswith((".gz", ".zip"))
):
print("Checking necessary files....")
utils.download_file_if_not_existing(link, local_path)
if self.version == 1 and file_name.endswith(".zip"):
utils.unzip_file(local_path, self.output_path)
if required is None:
self._downloads_ready = True
def _load_zone_geodataframe(self) -> None:
"""
Build/load the zone geodataframe lazily on first access.
"""
if self.complete_df is not None:
return
if self.zones == "provinces":
source = Zones(self.source_zones, self.version, self.output_path)
frame = source.get_zone_geodataframe()
mapping = source.get_province_mapping(source_ids=frame.index, unmapped="exclude")
correspondence = dict(zip(mapping["source_id"], mapping["target_id"]))
excluded = mapping.attrs["unmapped_source_ids"]
retained = frame.loc[frame.index.isin(correspondence), [frame.geometry.name]].copy()
retained["id"] = retained.index.map(correspondence)
missing_geometry = retained.geometry.isna() | retained.geometry.is_empty
if missing_geometry.any():
raise ValueError("Mapped source zones lack geometry: %s" % retained.index[missing_geometry].tolist()[:5])
invalid_geometry = ~retained.geometry.is_valid
if invalid_geometry.any():
retained.loc[invalid_geometry, frame.geometry.name] = retained.loc[invalid_geometry].geometry.make_valid()
if "population" in frame:
retained["population"] = pd.to_numeric(frame.loc[retained.index, "population"], errors="coerce")
retained = retained.reset_index(drop=True)
self.complete_df = retained.dissolve(
by="id", aggfunc=lambda values: values.sum(skipna=False)
) if not retained.empty else retained.set_index("id")
self.complete_df.attrs.update(
derived=True, source_zones=self.source_zones, output_zones="provinces",
unmapped_source_ids=excluded,
method="dissolve_mapped_mitma_zones",
repaired_source_geometries=int(invalid_geometry.sum()),
)
if excluded:
import warnings
warnings.warn("Province geometries exclude %d unmapped source zones; see result.attrs." % len(excluded), RuntimeWarning, stacklevel=2)
return
output_file_path = os.path.join(self.output_path, f"{self.zones}_{self.version}.geojson")
if os.path.exists(output_file_path):
print(f"File {output_file_path} already exists. Loading it...")
self.complete_df = self._canonicalize_geodataframe(
gpd.read_file(output_file_path)
)
return
self._ensure_zoning_files_downloaded()
print("Zones already downloaded. Reading the files....")
if self.version == 2:
def _read_pipe_csv(path, cols):
"""
Read a ‘|’-separated MITMA CSV that may or may not contain a header
and may start with a UTF-8 BOM. Returns a tidy DataFrame.
"""
df = pd.read_csv(
path,
sep="|",
dtype=str,
header=None,
names=cols,
encoding="utf-8-sig",
)
df[cols[0]] = df[cols[0]].str.strip()
if df.iloc[0, 0].upper() == cols[0].upper():
df = df.iloc[1:]
return df
nombre = _read_pipe_csv(
self._resolve_data_file(f"nombres_{self.zones}.csv"),
["ID", "name"],
)
pop = (
_read_pipe_csv(
self._resolve_data_file(f"poblacion_{self.zones}.csv"),
["ID", "population"],
)
.replace("NA", None)
)
zonification = gpd.read_file(
self._resolve_data_file(f"zonificacion_{self.zones}.shp")
)
if zonification.crs is None or zonification.crs.to_epsg() != 4326:
zonification = zonification.to_crs(epsg=4326)
for col in zonification.columns:
if col.lower() in {"id", "id_1", "codigo", "codigoine", "cod_mun"}:
zonification["ID"] = zonification[col].astype(str).str.strip()
break
else:
raise KeyError("No ID-like column found in the shapefile")
complete_df = (
nombre.set_index("ID")
.join(pop.set_index("ID"), how="outer")
.join(zonification.set_index("ID"), how="outer")
)
complete_df = gpd.GeoDataFrame(complete_df, crs="EPSG:4326")
complete_df.reset_index(inplace=True)
complete_df.rename(columns={"ID": "id"}, inplace=True)
complete_df = self._canonicalize_geodataframe(complete_df)
complete_df.to_file(output_file_path, driver="GeoJSON")
self.complete_df = complete_df
return
zonification = gpd.read_file(
os.path.join(self.output_path, f"zonificacion-{self.zones}/{self.zones}_mitma.shp")
)
if zonification.crs is None or zonification.crs.to_epsg() != 4326:
zonification = zonification.to_crs(epsg=4326)
complete_df = zonification
complete_df.rename(columns={"ID": "id"}, inplace=True)
complete_df = self._canonicalize_geodataframe(complete_df)
complete_df.to_file(output_file_path, driver="GeoJSON")
self.complete_df = complete_df
@staticmethod
def _canonicalize_geodataframe(frame: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
"""Return every zoning geometry with the same string ``id`` index.
GeoJSON persists a named index as a regular property. Without this
normalization, a fresh download and a cached GeoJSON exposed different
join contracts to users.
"""
if not isinstance(frame, gpd.GeoDataFrame):
raise TypeError("Zone geometry must be a GeoDataFrame.")
result = frame.copy()
id_column = next(
(column for column in result.columns if str(column).lower() == "id"),
None,
)
if id_column is not None:
identifiers = result[id_column]
elif result.index.name is not None and str(result.index.name).lower() == "id":
identifiers = pd.Series(result.index, index=result.index)
else:
raise KeyError("Zone geometry must provide an 'id' column or index.")
if identifiers.isna().any():
raise ValueError("Zone geometry contains null or empty IDs.")
identifiers = identifiers.astype(str).str.strip()
missing_markers = identifiers.str.lower().isin(
{"na", "nan", "none", "null"}
)
if identifiers.eq("").any() or missing_markers.any():
raise ValueError("Zone geometry contains null or empty IDs.")
if identifiers.duplicated().any():
examples = identifiers[identifiers.duplicated()].head(5).tolist()
raise ValueError("Zone geometry contains duplicate IDs: %s" % examples)
result.index = pd.Index(identifiers, name="id")
if id_column is not None:
result = result.drop(columns=[id_column])
return result
def _resolve_data_file(self, filename: str) -> str:
"""
Resolve a data file path, preferring the instance output path and
falling back to the global default data directory for backward
compatibility.
"""
preferred = os.path.join(self.output_path, filename)
if os.path.exists(preferred):
return preferred
fallback = os.path.join(utils.get_data_directory(), filename)
if os.path.exists(fallback):
return fallback
return preferred
def _read_relation_table(self, filename: str) -> pd.DataFrame:
"""
Read relation CSV files with robust delimiter detection.
"""
path = self._resolve_data_file(filename)
for sep in ("|", ",", ";", "\t"):
try:
df = pd.read_csv(path, sep=sep, dtype=str, encoding="utf-8-sig")
if len(df.columns) > 1:
return df
except Exception:
continue
return pd.read_csv(path, dtype=str, encoding="utf-8-sig")
[docs]
def get_zone_geodataframe(self):
"""
Function that returns the geodataframe with the zones. The geodataframe contains the following columns:
- id: the id of the zone
- name: the name of the zone
- population: the population of the zone (if available)
Parameters
----------
Examples
--------
>>> from pyspainmobility import Zones
>>> # instantiate the object
>>> zones = Zones(zones='municipalities', version=2, output_directory='data')
>>> # get the geodataframe with the zones
>>> gdf = zones.get_zone_geodataframe()
>>> print(gdf.head())
name population
ID
01001 Alegría-Dulantzi 2925.0
01002 Amurrio 10307.0
01004_AM Artziniega agregacion de municipios 3005.0
01009_AM Asparrena agregacion de municipios 4599.0
"""
self._load_zone_geodataframe()
return self.complete_df
[docs]
def get_zone_relations(self):
"""
Return official mapping tables between INE administrative units and
MITMA zoning identifiers.
For version 2, the returned table includes one row per relation entry
with harmonized column names.
For version 1, the method returns one row per MITMA zone id, where
each relation column contains the set of linked INE identifiers.
Parameters
----------
None
Returns
-------
pandas.DataFrame
Relation table between census/municipality identifiers and MITMA
zoning identifiers.
Examples
--------
>>> from pyspainmobility import Zones
>>> zones = Zones(zones='municipalities', version=2, output_directory='data')
>>> rel = zones.get_zone_relations()
>>> rel.columns.tolist()
['census_sections', 'census_districts', 'municipalities',
'municipalities_mitma', 'districts_mitma', 'luas_mitma']
"""
self._ensure_zoning_files_downloaded(
["relacion_ine_zonificacionMitma.csv"] if self.version == 2 else
["relaciones_municipio_mitma.csv", "relaciones_distrito_mitma.csv"]
)
if self.version == 2:
relacion = self._read_relation_table('relacion_ine_zonificacionMitma.csv')
remapping = {
'seccion_ine': 'census_sections',
'distrito_ine': 'census_districts',
'municipio_ine': 'municipalities',
'municipio_mitma': 'municipalities_mitma',
'distrito_mitma': 'districts_mitma',
'gau_mitma': 'luas_mitma'
}
relacion.rename(columns=remapping, inplace=True)
relacion = relacion.replace('NA', None)
return relacion
else:
used_zone = getattr(self, "source_zones", self.zones)[:-1]
relacion = self._read_relation_table(f'relaciones_{used_zone}_mitma.csv')
relacion.rename(columns={f'{used_zone}_mitma': 'id'}, inplace=True)
if used_zone == 'municipio':
temp = self._read_relation_table('relaciones_distrito_mitma.csv')
relacion = relacion.set_index('id').join(temp.set_index('municipio_mitma')).reset_index()
if used_zone == 'distrito':
temp = self._read_relation_table('relaciones_municipio_mitma.csv')
relacion = relacion.set_index('municipio_mitma').join(temp.set_index('municipio_mitma')).reset_index()
to_rename = {
'distrito': 'census_districts',
'distrito_mitma': 'districts_mitma',
'municipio': 'municipalities',
'municipio_mitma': 'municipalities_mitma',
}
relacion.rename(columns=to_rename, inplace=True)
temp_df = pd.DataFrame(relacion['id'].unique(), columns=['id']).set_index('id')
for i in list(relacion.columns):
if i != 'id':
temp_df = temp_df.join(relacion.groupby('id')[i].apply(set))
return temp_df
@staticmethod
def _normalize_relation_identifier(value, *, column: str) -> Optional[str]:
"""Normalize one scalar identifier from a zoning relation table."""
if isinstance(value, (list, tuple, set, frozenset)):
raise ValueError(
f"Relation column '{column}' contains a multi-valued entry. "
"Use the version-2 relation table or resolve it explicitly before "
"building a network mapping."
)
if pd.isna(value):
return None
identifier = str(value).strip()
if identifier.lower() in {"na", "nan", "none", "null"}:
return None
return identifier or None
[docs]
def get_network_mapping(
self,
source_column: str,
target_column: str,
source_ids=None,
) -> pd.DataFrame:
"""Return a verified, deterministic mapping for network aggregation.
The official relationship file can contain several rows for one source
zone. Repeated *identical* pairs are harmless and are removed; two
distinct targets for the same source are ambiguous and raise a
:class:`ValueError`. This prevents an aggregation from silently
choosing an arbitrary territorial correspondence.
Parameters
----------
source_column, target_column : str
Column names returned by :meth:`get_zone_relations`.
source_ids : iterable, optional
IDs occurring in the OD/network being aggregated. When supplied,
the output is restricted to these IDs and missing correspondences
are reported before any flow is transformed.
Returns
-------
pandas.DataFrame
Two string columns, ``source_id`` and ``target_id``, sorted by
``source_id``. It is accepted directly by
:func:`pyspainmobility.network.aggregate_network` and
:func:`pyspainmobility.network.aggregate_od_network`.
"""
mapping = self._relation_pairs(source_column, target_column, source_ids)
mapping["target_id"] = mapping["target_id"].map(
lambda value: self._normalize_relation_identifier(
value, column=target_column
)
)
return self._validate_network_mapping(mapping, target_column)
def _relation_pairs(self, source_column, target_column, source_ids, *, allow_missing=False):
"""Select relation pairs and restrict validation to requested sources."""
relations = self.get_zone_relations()
if relations.index.name == "id" and "id" not in relations.columns:
relations = relations.reset_index()
required = {source_column, target_column}
missing_columns = sorted(required.difference(relations.columns))
if missing_columns:
raise KeyError(
"Unknown relation column(s): %s. Available columns: %s"
% (missing_columns, sorted(relations.columns))
)
mapping = relations.loc[:, [source_column, target_column]].copy()
mapping.columns = ["source_id", "target_id"]
mapping["source_id"] = mapping["source_id"].map(
lambda value: self._normalize_relation_identifier(
value, column=source_column
)
)
# A row without a source cannot describe a source-to-target mapping.
mapping = mapping.loc[mapping["source_id"].notna()].copy()
if source_ids is not None:
if isinstance(source_ids, (str, bytes)):
raise TypeError("source_ids must be an iterable of identifiers, not a string.")
requested = pd.Index(
[
self._normalize_relation_identifier(value, column="source_ids")
for value in source_ids
]
)
if requested.isna().any():
raise ValueError("source_ids contains null or empty identifiers.")
requested = requested.drop_duplicates()
available = pd.Index(mapping["source_id"])
absent = requested.difference(available)
if not absent.empty and not allow_missing:
raise ValueError(
"No relation mapping is available for source IDs: %s"
% absent[:5].tolist()
)
mapping = mapping.loc[mapping["source_id"].isin(requested)]
mapping.attrs["unmapped_source_ids"] = absent.tolist()
return mapping
@staticmethod
def _validate_network_mapping(mapping, target_column):
"""Require exactly one target per source after target normalization."""
missing_target = mapping.loc[mapping["target_id"].isna(), "source_id"]
if not missing_target.empty:
examples = missing_target.drop_duplicates().head(5).tolist()
raise ValueError(
"Relation mapping has source IDs without a target in '%s': %s"
% (target_column, examples)
)
mapping = mapping.drop_duplicates()
target_counts = mapping.groupby("source_id", sort=False)["target_id"].nunique()
ambiguous = target_counts.loc[target_counts > 1].index.tolist()
if ambiguous:
raise ValueError(
"Relation mapping is ambiguous: source IDs map to multiple "
"targets in '%s': %s" % (target_column, ambiguous[:5])
)
return mapping.sort_values("source_id", kind="stable").reset_index(drop=True)
[docs]
def get_province_mapping(
self,
source_column: Optional[str] = None,
municipality_column: str = "municipalities",
source_ids=None,
*,
unmapped: str = "raise",
) -> pd.DataFrame:
"""Map MITMA zones to Spanish provinces through INE municipality IDs.
Province geometries are not a native :class:`Zones` level. Spanish
INE municipality identifiers have five digits and their first two
digits identify the province. This helper exposes that hierarchy as a
validated ``source_id``/``target_id`` mapping, where ``target_id`` is
the two-digit INE province code.
By default, the source column matches this object's selected zoning
level (``id`` for version 1). A zone may contain several municipalities
when all belong to the same province. Missing or invalid municipality
codes and zones spanning multiple provinces raise instead of splitting
their flow. Version-1 sets of municipality IDs are supported.
Code validation checks the five-digit format and province prefix 01–52,
not historical municipality registry membership.
``unmapped='exclude'`` omits entire sources with absent, invalid, or
ambiguous relations and lists them in ``attrs['unmapped_source_ids']``.
The default ``unmapped='raise'`` retains strict mapping validation.
For ``Zones(zones='provinces')`` this maps the underlying district IDs.
"""
if unmapped not in {"raise", "exclude"}:
raise ValueError("unmapped must be 'raise' or 'exclude'.")
if source_column is None:
source_column = "id" if self.version == 1 else {
"distritos": "districts_mitma",
"municipios": "municipalities_mitma",
"gaus": "luas_mitma",
}[getattr(self, "source_zones", self.zones)]
mapping = self._relation_pairs(
source_column, municipality_column, source_ids,
allow_missing=unmapped == "exclude",
).explode("target_id", ignore_index=True)
excluded = set(mapping.attrs.get("unmapped_source_ids", []))
municipalities = mapping["target_id"].map(
lambda value: self._normalize_relation_identifier(
value, column=municipality_column
)
).astype("string")
provinces = municipalities.str[:2]
valid_municipalities = (
municipalities.str.fullmatch(r"[0-9]{5}").fillna(False)
& provinces.isin([f"{code:02d}" for code in range(1, 53)])
)
if not valid_municipalities.all():
examples = mapping.loc[
~valid_municipalities, ["source_id", "target_id"]
].head(5).to_dict("records")
if unmapped == "raise":
raise ValueError(
"Cannot derive province codes: '%s' must contain five-digit INE "
"municipality identifiers with province codes 01–52; invalid relations: %s"
% (municipality_column, examples)
)
excluded.update(mapping.loc[~valid_municipalities, "source_id"])
mapping["target_id"] = provinces
if unmapped == "exclude":
counts = mapping.groupby("source_id")["target_id"].nunique()
excluded.update(counts.index[counts > 1])
mapping = mapping.loc[~mapping["source_id"].isin(excluded)]
mapping = self._validate_network_mapping(mapping, "province")
mapping.attrs["unmapped_source_ids"] = sorted(excluded)
return mapping