{ "cells": [ { "cell_type": "markdown", "id": "d24a49d3", "metadata": {}, "source": [ "# 2. PART TWO\n", "Looking at the facts and figures above, integrated urban development must promote sustainable transportation and infrastructure. Therefore, public transit services are one of the prerequisites for transforming the world toward a net-zero future. Let us explore a national-scale dataset related to urban transportation infrastructure services in Germany. \"[ioerDATA](https://data.fdz.ioer.de/)—the Institutional Research Data Repository of the IOER—is a curated dataset on [urban public transit service frequency](https://doi.org/10.71830/ABPCUS).\"\n", "\n", "\n", "## Urban Public Transit Frequency Indicator — Germany\n", "\n", "Three analyses:\n", "\n", "| # | Analysis | Source |\n", "|---|---|---|\n", "| 1 | Country-level point density of transit stations | `xx_am_transit_frequency_station.gpkg` |\n", "| 2 | City-level summary statistics | `de_transit_frequency_indicator_aggregated_statistics_all_cities.csv` (file id 6215) |\n", "| 3 | City-level frequency EDA + spatial statistics | `data/..._grid_single_cities//*.gpkg` |\n", "\n", "Analyses 3 follows the method of\n", "[1_exploratory_data_analysis](https://github.com/ssujit/public_transit_germany/blob/main/jupyter_notebook/1_exploratory_data_analysis.ipynb) and\n", "[3_spatial_statistics_neighborhood](https://github.com/ssujit/public_transit_germany/blob/main/jupyter_notebook/3_spatial_statistics_neighborhood.ipynb).\n", "\n", "### Running this on Colab\n", "\n", "Run the cells in order; section 0 installs what Colab lacks and picks up your API token.\n", "\n", "To open it, either **File -> Upload notebook** and pick this `.ipynb` (works regardless of\n", "repository visibility), or, if you have linked Colab to GitHub and can see this private\n", "repo, **File -> Open notebook -> GitHub** and search `ssujit/urban_mobiltiy4netzero`.\n", "Once the repository is public this badge works directly:\n", "\n", "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/ssujit/urban_mobiltiy4netzero/blob/main/scripts/gtfs.ipynb)\n", "\n", "**Your token:** add it once as a Colab secret named `IOER_API_TOKEN` (key icon in the left\n", "sidebar, with \"Notebook access\" enabled) and it is picked up automatically. Otherwise you\n", "are prompted for it, hidden as you type. Either way the token is never written into the\n", "notebook, so the file stays safe to commit.\n", "\n", "**Access:** the account behind the token needs download rights on this DOI. The dataset\n", "custodian's own account has them; a participant account generally does not until access is\n", "granted (see section 1, which checks this before anything else runs).\n", "\n", "### Two things to know before running\n", "\n", "1. **Almost every file in this dataset is access-restricted.** You need an ioerDATA API\n", " token (ioerDATA → your account → API Token). Set it as `IOER_API_TOKEN` in your\n", " environment, or paste it when prompted below.\n", "2. **Column naming differs from the reference notebooks.** The per-city grid files use\n", " `fi` for the frequency indicator (the reference EDA notebook says `pti`, because it\n", " reads a differently-named copy from GitHub). The INSPIRE grid CSV uses `pti_`.\n", " This notebook uses `fi`." ] }, { "cell_type": "markdown", "id": "fd25ae76", "metadata": {}, "source": [ "## 0 · Environment" ] }, { "cell_type": "code", "execution_count": null, "id": "4cebee42", "metadata": {}, "outputs": [], "source": [ "# Dependencies. Safe to re-run; already-satisfied packages are skipped.\n", "# Colab already has pandas/matplotlib/seaborn and usually geopandas, so only the\n", "# spatial-statistics stack and the fast OGR reader normally need installing.\n", "import importlib.util # NB: plain `import importlib` does not expose .util\n", "import sys\n", "\n", "IN_COLAB = \"google.colab\" in sys.modules\n", "\n", "PIP_NAME = {\"sklearn\": \"scikit-learn\"} # import name -> pip name\n", "\n", "missing = [m for m in (\"geopandas\", \"pyogrio\", \"mapclassify\",\n", " \"esda\", \"libpysal\", \"splot\", \"sklearn\")\n", " if importlib.util.find_spec(m) is None]\n", "\n", "if missing:\n", " pkgs = \" \".join(PIP_NAME.get(m, m) for m in missing)\n", " print(\"installing:\", pkgs)\n", " %pip install -q {pkgs}\n", " print(\"\\nIf an import below fails, restart the runtime \"\n", " \"(Runtime -> Restart session) and re-run from here.\")\n", "else:\n", " print(\"all dependencies already present\")" ] }, { "cell_type": "code", "execution_count": null, "id": "57bf36c1", "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "import warnings\n", "from pathlib import Path\n", "from getpass import getpass\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import requests\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "sns.set_theme(style=\"whitegrid\")\n", "plt.rcParams[\"figure.dpi\"] = 110\n", "\n", "IN_COLAB = \"google.colab\" in sys.modules # restated so this cell stands alone\n", "\n", "print(f\"python {sys.version.split()[0]} | geopandas {gpd.__version__} | \"\n", " f\"pandas {pd.__version__} | matplotlib {plt.matplotlib.__version__}\")\n", "print(f\"running on Colab: {IN_COLAB}\")" ] }, { "cell_type": "markdown", "id": "8af179f2", "metadata": {}, "source": [ "## 1 · Connect to ioerDATA and index the dataset\n", "\n", "Rather than hard-coding file ids (which change between dataset versions), pull the file\n", "listing once and build a lookup table. Every download below resolves its id through this\n", "index." ] }, { "cell_type": "code", "execution_count": null, "id": "f61690b6", "metadata": {}, "outputs": [], "source": [ "BASE_URL = \"https://data.fdz.ioer.de\"\n", "DOI = \"doi:10.71830/ABPCUS\"\n", "\n", "# Token from the environment if present, otherwise prompt (input stays hidden).\n", "\n", "\n", "def get_api_token():\n", " \"\"\"Colab secret -> environment variable -> hidden prompt.\"\"\"\n", " if IN_COLAB:\n", " try:\n", " from google.colab import userdata\n", " token = userdata.get(\"IOER_API_TOKEN\")\n", " if token:\n", " print(\"token loaded from Colab secrets\")\n", " return token.strip()\n", " except Exception:\n", " pass # secret not set, or access not granted to this notebook\n", " token = os.environ.get(\"IOER_API_TOKEN\")\n", " if token:\n", " print(\"token loaded from the environment\")\n", " return token.strip()\n", " return getpass(\"ioerDATA API token: \").strip()\n", "\n", "\n", "api_token = get_api_token()\n", "headers = {\"X-Dataverse-key\": api_token}\n", "\n", "resp = requests.get(\n", " f\"{BASE_URL}/api/datasets/:persistentId/\",\n", " params={\"persistentId\": DOI},\n", " headers=headers,\n", " timeout=60,\n", ")\n", "resp.raise_for_status()\n", "\n", "files = resp.json()[\"data\"][\"latestVersion\"][\"files\"]\n", "\n", "index = pd.DataFrame(\n", " [\n", " {\n", " \"id\": f[\"dataFile\"][\"id\"],\n", " \"filename\": f[\"dataFile\"][\"filename\"],\n", " \"directory\": f.get(\"directoryLabel\", \"\"),\n", " \"mb\": round(f[\"dataFile\"].get(\"filesize\", 0) / 1e6, 2),\n", " \"restricted\": f.get(\"restricted\", False),\n", " }\n", " for f in files\n", " ]\n", ")\n", "\n", "print(f\"{len(index)} files, {index.mb.sum():.0f} MB total, \"\n", " f\"{index.restricted.sum()} restricted\")\n", "\n", "# --- Permission check -------------------------------------------------------\n", "# A valid token is not the same as access. Every analysis file in this dataset is\n", "# restricted, and the token only authenticates you; the data custodian still has to\n", "# grant download rights. Check before doing anything else.\n", "me = requests.get(f\"{BASE_URL}/api/users/:me\", headers=headers, timeout=30).json()\n", "print(f\"authenticated as {me['data']['displayName']} (id {me['data']['id']})\")\n", "\n", "probe = index[index.filename ==\n", " \"de_transit_frequency_indicator_aggregated_statistics_all_cities.csv\"].iloc[0]\n", "perm = requests.get(f\"{BASE_URL}/api/access/datafile/{probe.id}/userPermissions\",\n", " headers=headers, timeout=30).json()[\"data\"]\n", "\n", "if perm[\"canDownloadFile\"]:\n", " print(\"access granted - the analyses below will run\")\n", "else:\n", " print(\n", " \"\\nNO DOWNLOAD ACCESS to the restricted files.\\n\"\n", " \"Neither downloading nor streaming can work until this is granted: both hit the\\n\"\n", " \"same endpoint and both return 403. Request access per file with\\n\\n\"\n", " f' requests.post(f\"{{BASE_URL}}/api/access/datafile/{{probe.id}}/requestAccess\",\\n'\n", " \" headers=headers, timeout=30)\\n\\n\"\n", " \"or use the Request Access button on the file's ioerDATA page, then wait for the\\n\"\n", " \"custodian (s.sikder@ioer.de) to approve.\"\n", " )\n", "\n", "index[index.directory.isin([\"\", \"data\"])]" ] }, { "cell_type": "code", "execution_count": null, "id": "bbcc1091", "metadata": {}, "outputs": [], "source": [ "# Colab runtimes are wiped when the session ends. Set USE_DRIVE = True to keep the\n", "# downloads in your Google Drive so a later session reuses them instead of refetching.\n", "USE_DRIVE = False\n", "\n", "if IN_COLAB and USE_DRIVE:\n", " from google.colab import drive\n", " drive.mount(\"/content/drive\")\n", " DATA_DIR = Path(\"/content/drive/MyDrive/ioer_transit/data/raw\")\n", "else:\n", " DATA_DIR = Path(\"data/raw\")\n", "\n", "\n", "def fetch(filename, directory=None, overwrite=False):\n", " \"\"\"Download one dataset file by name, cached under data/raw//.\n", "\n", " Returns the local Path. Re-running is free: an existing file is reused.\n", " \"\"\"\n", " if \"index\" not in globals():\n", " raise RuntimeError(\"Run section 1 first - it authenticates and builds `index`.\")\n", "\n", " hits = index[index.filename == filename]\n", " if directory is not None:\n", " hits = hits[hits.directory == directory]\n", " if len(hits) == 0:\n", " raise FileNotFoundError(f\"{filename!r} is not in this dataset version\")\n", " if len(hits) > 1:\n", " raise ValueError(\n", " f\"{filename!r} is ambiguous; pass directory= one of \"\n", " f\"{hits.directory.tolist()}\"\n", " )\n", " row = hits.iloc[0]\n", "\n", " out = DATA_DIR / row.directory / row.filename\n", " if out.exists() and not overwrite:\n", " return out\n", " out.parent.mkdir(parents=True, exist_ok=True)\n", "\n", " print(f\"downloading {row.filename} ({row.mb} MB) ...\", end=\" \", flush=True)\n", " r = requests.get(\n", " f\"{BASE_URL}/api/access/datafile/{row.id}\",\n", " headers=headers,\n", " stream=True,\n", " timeout=300,\n", " )\n", " if r.status_code in (401, 403):\n", " raise PermissionError(\n", " f\"No access to {row.filename}. This file is restricted — check that your \"\n", " \"token is valid and that access has been granted for this dataset.\"\n", " )\n", " r.raise_for_status()\n", "\n", " with open(out, \"wb\") as fh:\n", " for chunk in r.iter_content(1024 * 1024):\n", " fh.write(chunk)\n", " print(\"done\")\n", " return out\n", "\n", "\n", "print(f\"cache: {DATA_DIR.resolve()}\")" ] }, { "cell_type": "markdown", "id": "4d7c3365", "metadata": {}, "source": [ "### Streaming instead of downloading\n", "\n", "ioerDATA serves `206 Partial Content` and advertises `Accept-Ranges: bytes`, so GDAL can\n", "read a remote GeoPackage through `/vsicurl` and pull only the pages it needs — no file on\n", "disk. The helpers below do that.\n", "\n", "**When streaming is worth it:** reading the schema, or pulling a bbox/column subset out of\n", "a large layer. Fetching 1,200 features from the 106 MB national station layer touches a few\n", "MB instead of all of it.\n", "\n", "**When it is not:** reading a whole layer. A GeoPackage is a SQLite database, so a full\n", "table scan becomes thousands of small range requests — measurably slower than one bulk\n", "download. For the per-city grids (100–370 KB) always just `fetch()` them." ] }, { "cell_type": "code", "execution_count": null, "id": "59d093c5", "metadata": {}, "outputs": [], "source": [ "import io\n", "\n", "# GDAL passes these to libcurl on every /vsicurl request.\n", "os.environ[\"GDAL_HTTP_HEADERS\"] = f\"X-Dataverse-key: {api_token}\"\n", "os.environ[\"CPL_VSIL_CURL_USE_HEAD\"] = \"NO\" # Dataverse HEAD omits Content-Length\n", "os.environ[\"GDAL_HTTP_MULTIRANGE\"] = \"YES\" # coalesce reads into one request\n", "os.environ[\"GDAL_HTTP_MERGE_CONSECUTIVE_RANGES\"] = \"YES\"\n", "os.environ[\"VSI_CACHE\"] = \"TRUE\"\n", "os.environ[\"VSI_CACHE_SIZE\"] = str(64 * 1024 * 1024)\n", "\n", "\n", "def _row(filename, directory=None):\n", " hits = index[index.filename == filename]\n", " if directory is not None:\n", " hits = hits[hits.directory == directory]\n", " if len(hits) != 1:\n", " raise ValueError(f\"{filename!r} matched {len(hits)} files; pass directory=\")\n", " return hits.iloc[0]\n", "\n", "\n", "def remote_url(filename, directory=None):\n", " \"\"\"GDAL virtual path for a dataset file — readable without downloading it.\"\"\"\n", " return f\"/vsicurl/{BASE_URL}/api/access/datafile/{_row(filename, directory).id}\"\n", "\n", "\n", "def read_remote(filename, directory=None, **kwargs):\n", " \"\"\"Read a layer straight off the server.\n", "\n", " Pass bbox=(minx, miny, maxx, maxy) and/or columns=[...] to pull a subset —\n", " that is the case where streaming beats downloading.\n", " \"\"\"\n", " return gpd.read_file(remote_url(filename, directory), engine=\"pyogrio\", **kwargs)\n", "\n", "\n", "def read_remote_csv(filename, directory=None, **kwargs):\n", " \"\"\"Stream a CSV into a DataFrame without writing it to disk.\"\"\"\n", " r = requests.get(f\"{BASE_URL}/api/access/datafile/{_row(filename, directory).id}\",\n", " headers=headers, timeout=300)\n", " r.raise_for_status()\n", " return pd.read_csv(io.BytesIO(r.content), **kwargs)\n", "\n", "\n", "# Schema of the 106 MB national station layer, without downloading it:\n", "# pyogrio.read_info(remote_url(\"wk_am_transit_frequency_station.gpkg\",\n", "# \"other/transit_frequency_station\"))\n", "print(\"streaming helpers ready: remote_url / read_remote / read_remote_csv\")" ] }, { "cell_type": "markdown", "id": "970ebe36", "metadata": {}, "source": [ "> **Note on bulk download.** The full dataset is ~2.5 GB (a 1.4 GB road network, 12 station\n", "> layers of 20–210 MB each, and 1,158 per-city grids). Fetch only what each analysis needs\n", "> — which is what `fetch()` above does." ] }, { "cell_type": "markdown", "id": "78967631", "metadata": {}, "source": [ "## 2 · Analysis 1 — Country-level point density\n", "\n", "Station-level service frequency for all of Germany. The dataset ships each time slot twice,\n", "as GeoJSON and GeoPackage; the GeoPackage is half the size and parses far faster, so use it\n", "(`wk_am_transit_frequency_station.gpkg`, 106 MB — the 211 MB GeoJSON at file id 6268 holds\n", "the same data)." ] }, { "cell_type": "code", "execution_count": null, "id": "c35ec318", "metadata": {}, "outputs": [], "source": [ "STATION_SLOT = \"wk_am\" # wk_am | wk_pm | sat_am | sat_pm | sun_am | sun_pm\n", "\n", "stations_path = fetch(f\"{STATION_SLOT}_transit_frequency_station.gpkg\",\n", " directory=\"other/transit_frequency_station\")\n", "\n", "stations = gpd.read_file(stations_path, engine=\"pyogrio\")\n", "print(f\"{len(stations):,} station records | CRS {stations.crs}\")\n", "stations.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "63962e06", "metadata": {}, "outputs": [], "source": [ "# The station layer's fields are not documented in the dataset README, so inspect them\n", "# and pick the frequency column rather than assuming a name.\n", "print(stations.dtypes, end=\"\\n\\n\")\n", "\n", "numeric_cols = stations.select_dtypes(\"number\").columns.tolist()\n", "print(\"numeric columns:\", numeric_cols)\n", "\n", "# Prefer an explicit frequency field; fall back to the first numeric column.\n", "for candidate in (\"fi\", \"perH_sum\", \"freq\", \"count\", \"Value\"):\n", " if candidate in stations.columns:\n", " WEIGHT_COL = candidate\n", " break\n", "else:\n", " WEIGHT_COL = numeric_cols[0]\n", "\n", "print(f\"\\nusing WEIGHT_COL = {WEIGHT_COL!r}\")\n", "stations[WEIGHT_COL].describe()" ] }, { "cell_type": "code", "execution_count": null, "id": "d53843f9", "metadata": {}, "outputs": [], "source": [ "# Project to ETRS89 / LAEA Europe (EPSG:3035) — the dataset's own grid CRS, and an\n", "# equal-area projection, so hexagon counts are comparable across the country.\n", "stations_3035 = stations.to_crs(3035)\n", "x = stations_3035.geometry.x.to_numpy()\n", "y = stations_3035.geometry.y.to_numpy()\n", "\n", "boundaries = gpd.read_file(\n", " fetch(\"de_city_admin_boundary.gpkg\", directory=\"other\"), engine=\"pyogrio\"\n", ").to_crs(3035)\n", "\n", "print(f\"{len(boundaries)} city boundaries | bounds {stations_3035.total_bounds.round(0)}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "855463e2", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(15, 9), constrained_layout=True)\n", "\n", "# Left: station count per hexagon (where stations are).\n", "hb = axes[0].hexbin(x, y, gridsize=140, bins=\"log\", cmap=\"magma_r\", mincnt=1, linewidths=0)\n", "axes[0].set_title(f\"Station density — {STATION_SLOT}\\n(count per hexagon, log scale)\")\n", "fig.colorbar(hb, ax=axes[0], shrink=0.6, label=\"stations per hexagon\")\n", "\n", "# Right: mean service frequency per hexagon (how good the service is).\n", "hb2 = axes[1].hexbin(\n", " x, y, C=stations_3035[WEIGHT_COL].to_numpy(),\n", " reduce_C_function=np.mean, gridsize=140, cmap=\"viridis\", mincnt=1, linewidths=0,\n", ")\n", "axes[1].set_title(f\"Mean {WEIGHT_COL} — {STATION_SLOT}\\n(per hexagon)\")\n", "fig.colorbar(hb2, ax=axes[1], shrink=0.6, label=f\"mean {WEIGHT_COL}\")\n", "\n", "for ax in axes:\n", " boundaries.boundary.plot(ax=ax, color=\"0.35\", linewidth=0.3, zorder=3)\n", " ax.set_aspect(\"equal\")\n", " ax.set_axis_off()\n", "\n", "fig.suptitle(\"Germany — transit station point density (EPSG:3035)\", fontsize=14)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "8ca5e055", "metadata": {}, "outputs": [], "source": [ "# Save for reuse outside the notebook.\n", "OUT = Path(\"outputs/figures\")\n", "OUT.mkdir(parents=True, exist_ok=True)\n", "fig.savefig(OUT / f\"station_density_{STATION_SLOT}.png\", dpi=200, bbox_inches=\"tight\")\n", "print(f\"saved {OUT / f'station_density_{STATION_SLOT}.png'}\")" ] }, { "cell_type": "markdown", "id": "c7088756", "metadata": {}, "source": [ "## 3 · Analysis 2 — City-level statistics\n", "\n", "File id **6215**, `de_transit_frequency_indicator_aggregated_statistics_all_cities.csv`:\n", "one row per city × time slot, with `min`, `max`, `mean`, `median`, `std_dev` of the\n", "frequency indicator, plus the official `ags` regional code for joining to other statistics." ] }, { "cell_type": "code", "execution_count": null, "id": "a60b02d0", "metadata": {}, "outputs": [], "source": [ "stats = pd.read_csv(fetch(\n", " \"de_transit_frequency_indicator_aggregated_statistics_all_cities.csv\",\n", " directory=\"data\",\n", "))\n", "\n", "print(f\"{stats.city.nunique()} cities × {stats.day_time.nunique()} time slots \"\n", " f\"= {len(stats)} rows\")\n", "print(\"slots:\", sorted(stats.day_time.unique()))\n", "stats.head()" ] }, { "cell_type": "code", "execution_count": null, "id": "353444e1", "metadata": {}, "outputs": [], "source": [ "SLOT_ORDER = [\"wk_am\", \"wk_pm\", \"sat_am\", \"sat_pm\", \"sun_am\", \"sun_pm\"]\n", "stats[\"day_time\"] = pd.Categorical(stats.day_time, SLOT_ORDER, ordered=True)\n", "\n", "stats.groupby(\"day_time\", observed=True)[[\"mean\", \"median\", \"std_dev\", \"max\"]].describe().T" ] }, { "cell_type": "code", "execution_count": null, "id": "d1cbf5df", "metadata": {}, "outputs": [], "source": [ "# Distribution of city mean frequency indicator, by time slot.\n", "fig, axes = plt.subplots(1, 2, figsize=(14, 5), constrained_layout=True)\n", "\n", "sns.boxplot(data=stats, x=\"mean\", y=\"day_time\", hue=\"day_time\",\n", " order=SLOT_ORDER, palette=\"viridis\", legend=False, ax=axes[0])\n", "axes[0].set(title=\"City mean frequency indicator by time slot\",\n", " xlabel=\"mean fi\", ylabel=\"\")\n", "\n", "for slot in SLOT_ORDER:\n", " sns.kdeplot(stats.loc[stats.day_time == slot, \"mean\"],\n", " label=slot, ax=axes[1], linewidth=1.8)\n", "axes[1].set(title=\"Density of city means\", xlabel=\"mean fi\")\n", "axes[1].legend(title=\"slot\")\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "7c475cbc", "metadata": {}, "outputs": [], "source": [ "# Top and bottom 15 cities on weekday morning service.\n", "wk = stats[stats.day_time == \"wk_am\"].sort_values(\"mean\", ascending=False)\n", "ends = pd.concat([wk.head(15).assign(group=\"Top 15\"),\n", " wk.tail(15).assign(group=\"Bottom 15\")])\n", "\n", "fig, ax = plt.subplots(figsize=(9, 10), constrained_layout=True)\n", "sns.barplot(data=ends, x=\"mean\", y=\"city\", hue=\"group\",\n", " palette={\"Top 15\": \"#1b7837\", \"Bottom 15\": \"#c0504d\"}, dodge=False, ax=ax)\n", "ax.set(title=\"Weekday morning (06:00–08:59) mean transit frequency indicator\",\n", " xlabel=\"mean fi\", ylabel=\"\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "8116beed", "metadata": {}, "outputs": [], "source": [ "# Heatmap: the 30 best-served cities across all six slots.\n", "top30 = wk.head(30).city\n", "grid = (stats[stats.city.isin(top30)]\n", " .pivot(index=\"city\", columns=\"day_time\", values=\"mean\")\n", " .reindex(columns=SLOT_ORDER)\n", " .sort_values(\"wk_am\", ascending=False))\n", "\n", "fig, ax = plt.subplots(figsize=(8, 11), constrained_layout=True)\n", "sns.heatmap(grid, cmap=\"viridis\", annot=True, fmt=\".1f\",\n", " cbar_kws={\"label\": \"mean fi\"}, linewidths=0.4, ax=ax)\n", "ax.set(title=\"Mean frequency indicator — 30 best-served cities\", xlabel=\"\", ylabel=\"\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "c75beb40", "metadata": {}, "outputs": [], "source": [ "# Weekend service retention: how much of weekday morning service survives to Sunday.\n", "wide = stats.pivot(index=\"city\", columns=\"day_time\", values=\"mean\").reindex(\n", " columns=SLOT_ORDER)\n", "wide[\"sunday_retention\"] = wide[\"sun_am\"] / wide[\"wk_am\"]\n", "\n", "print(wide[\"sunday_retention\"].describe().round(3), end=\"\\n\\n\")\n", "print(\"Weakest Sunday service relative to weekday:\")\n", "print(wide[\"sunday_retention\"].nsmallest(10).round(3))\n", "\n", "fig, ax = plt.subplots(figsize=(11, 5), constrained_layout=True)\n", "sns.scatterplot(data=wide.reset_index(), x=\"wk_am\", y=\"sun_am\",\n", " size=\"sunday_retention\", hue=\"sunday_retention\",\n", " palette=\"RdYlGn\", sizes=(15, 160), ax=ax)\n", "lim = [0, wide[[\"wk_am\", \"sun_am\"]].to_numpy().max() * 1.05]\n", "ax.plot(lim, lim, \"k--\", linewidth=0.8, label=\"parity\")\n", "ax.set(title=\"Sunday morning vs weekday morning service, by city\",\n", " xlabel=\"weekday AM mean fi\", ylabel=\"Sunday AM mean fi\", xlim=lim, ylim=lim)\n", "ax.legend(bbox_to_anchor=(1.02, 1), loc=\"upper left\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "9c19d5b8", "metadata": {}, "source": [ "## 3b · Service inequality within cities\n", "\n", "Two cities with the same mean can be very different: one spread evenly, one a dense core\n", "and nothing else. `std_dev`, `mean` and `median` are already in the statistics table." ] }, { "cell_type": "code", "execution_count": null, "id": "d9d4eb34", "metadata": {}, "outputs": [], "source": [ "am = stats[stats.day_time == \"wk_am\"].set_index(\"city\")\n", "\n", "ineq = pd.DataFrame({\n", " \"mean\": am[\"mean\"],\n", " \"cv\": am[\"std_dev\"] / am[\"mean\"], # spread relative to level\n", " \"skew_gap\": (am[\"mean\"] - am[\"median\"]) / am[\"mean\"], # >0 = a few hot cells\n", "}).dropna()\n", "\n", "ineq.describe().round(2)" ] }, { "cell_type": "code", "execution_count": null, "id": "a21759f2", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(8, 6), constrained_layout=True)\n", "ax.scatter(ineq[\"mean\"], ineq[\"cv\"], s=18, alpha=0.7)\n", "\n", "mx, my = ineq[\"mean\"].median(), ineq[\"cv\"].median()\n", "ax.axvline(mx, c=\"k\", lw=0.8, ls=\"--\")\n", "ax.axhline(my, c=\"k\", lw=0.8, ls=\"--\")\n", "\n", "for xf, yf, txt in [(0.97, 0.97, \"high service\\nuneven\"),\n", " (0.03, 0.97, \"low service\\nuneven\"),\n", " (0.97, 0.03, \"high service\\neven\"),\n", " (0.03, 0.03, \"low service\\neven\")]:\n", " ax.text(xf, yf, txt, transform=ax.transAxes, fontsize=8, color=\"crimson\",\n", " ha=\"right\" if xf > 0.5 else \"left\", va=\"top\" if yf > 0.5 else \"bottom\")\n", "\n", "ax.set(xlabel=\"mean fi (weekday AM)\", ylabel=\"coefficient of variation\",\n", " title=\"Service level vs evenness\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "74cb2ec4", "metadata": {}, "outputs": [], "source": [ "print(\"most uneven:\\n\", ineq.nlargest(8, \"cv\").round(2), \"\\n\")\n", "print(\"most even:\\n\", ineq.nsmallest(8, \"cv\").round(2))" ] }, { "cell_type": "markdown", "id": "6604f1ae", "metadata": {}, "source": [ "## 3c · Weekly service profiles\n", "\n", "Dividing each city by its own peak removes the level and leaves the *shape* of its week,\n", "so cities can be grouped by when they run service rather than how much." ] }, { "cell_type": "code", "execution_count": null, "id": "883dbce7", "metadata": {}, "outputs": [], "source": [ "from sklearn.cluster import KMeans\n", "\n", "prof = stats.pivot(index=\"city\", columns=\"day_time\", values=\"mean\").reindex(\n", " columns=SLOT_ORDER).dropna()\n", "\n", "shape = prof.div(prof.max(axis=1), axis=0) # level removed\n", "shape[\"cluster\"] = KMeans(4, n_init=10, random_state=0).fit_predict(shape[SLOT_ORDER])\n", "\n", "shape.cluster.value_counts().sort_index()" ] }, { "cell_type": "code", "execution_count": null, "id": "2a96049d", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(8, 5), constrained_layout=True)\n", "\n", "for c, g in shape.groupby(\"cluster\"):\n", " ax.plot(SLOT_ORDER, g[SLOT_ORDER].mean(), marker=\"o\", label=f\"type {c} (n={len(g)})\")\n", "\n", "ax.set(ylabel=\"share of the city's own peak\", title=\"Weekly service profiles\")\n", "ax.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "e50adb2b", "metadata": {}, "outputs": [], "source": [ "for c, g in shape.groupby(\"cluster\"):\n", " print(f\"type {c}: {', '.join(g.index[:6])}\")" ] }, { "cell_type": "markdown", "id": "18c7626d", "metadata": {}, "source": [ "## 4 · Analysis 3 — City-level frequency EDA\n", "\n", "Per-city 1 km grids, six time slots each. Fields (per the dataset README):\n", "`fid`, `Value` (network distance, m), `grid_name`, `net_dis` (rounded), `perH_sum`\n", "(aggregation per hour), **`fi`** (the frequency indicator), `city`, `geometry`." ] }, { "cell_type": "code", "execution_count": null, "id": "6e0cb229", "metadata": {}, "outputs": [], "source": [ "# Any of the 193 cities. The folder name is the slug; filenames follow _.gpkg,\n", "# with a few irregular spellings — so select on the directory, not the filename.\n", "CITY = \"dresden\"\n", "\n", "city_dir = f\"data/de_transit_frequency_indicator_grid_single_cities/{CITY}\"\n", "available = index[index.directory == city_dir]\n", "if available.empty:\n", " slugs = (index[index.directory.str.contains(\"single_cities\", na=False)]\n", " .directory.str.rsplit(\"/\", n=1).str[-1].unique())\n", " raise ValueError(f\"Unknown city {CITY!r}. Available: {sorted(slugs)[:20]} ...\")\n", "\n", "grids = {}\n", "for _, row in available.iterrows():\n", " slot = \"_\".join(row.filename.split(\"_\")[:2]) # e.g. \"wk_am\"\n", " grids[slot] = gpd.read_file(fetch(row.filename, directory=city_dir),\n", " engine=\"pyogrio\")\n", "\n", "grids = {s: grids[s] for s in SLOT_ORDER if s in grids}\n", "for slot, gdf in grids.items():\n", " print(f\"{slot:7s} {len(gdf):5d} cells fi: \"\n", " f\"{gdf.fi.min():.1f}–{gdf.fi.max():.1f} mean {gdf.fi.mean():.1f}\")\n", "\n", "grids[\"wk_am\"].head()" ] }, { "cell_type": "code", "execution_count": null, "id": "234e4c3e", "metadata": {}, "outputs": [], "source": [ "# Choropleth of fi across all six slots, quantile-classified so classes are comparable.\n", "fig, axes = plt.subplots(3, 2, figsize=(13, 17), constrained_layout=True)\n", "\n", "for ax, (slot, gdf) in zip(axes.flatten(), grids.items()):\n", " gdf.plot(\"fi\", cmap=\"viridis\", scheme=\"quantiles\", k=5,\n", " edgecolor=\"white\", linewidth=0.1, legend=True,\n", " legend_kwds={\"loc\": \"upper left\", \"fontsize\": 7}, ax=ax)\n", " ax.set_title(f\"{CITY.title()} — {slot}\")\n", " ax.set_axis_off()\n", "\n", "fig.suptitle(f\"Transit frequency indicator, {CITY.title()} (1 km grid)\", fontsize=15)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "4ac6ff5a", "metadata": {}, "outputs": [], "source": [ "# Distributions. fi is strongly right-skewed, which is why the spatial statistics\n", "# section below works on a transformed variable.\n", "fig, axes = plt.subplots(3, 2, figsize=(13, 11), constrained_layout=True)\n", "\n", "for ax, (slot, gdf) in zip(axes.flatten(), grids.items()):\n", " ax.hist(gdf.fi.dropna(), bins=40, color=\"#3182bd\", edgecolor=\"white\", linewidth=0.4)\n", " ax.set(title=slot, xlabel=\"fi\", ylabel=\"cells\")\n", "\n", "fig.suptitle(f\"{CITY.title()} — distribution of fi by time slot\", fontsize=14)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "4dfc48be", "metadata": {}, "outputs": [], "source": [ "# Side-by-side comparison across slots.\n", "# NOTE: drawn with seaborn rather than plt.boxplot on purpose. The matplotlib keyword\n", "# for category labels changed (`labels=` -> `tick_labels=`) and `vert=False` was\n", "# superseded by `orientation=`, so raw boxplot calls break on either side of that\n", "# version boundary. Colab and a local env rarely ship the same matplotlib.\n", "fig, axes = plt.subplots(1, 2, figsize=(14, 5), constrained_layout=True)\n", "\n", "long = pd.concat([g[[\"fi\"]].assign(slot=slot) for slot, g in grids.items()])\n", "sns.boxplot(data=long, x=\"fi\", y=\"slot\", hue=\"slot\", order=list(grids),\n", " palette=\"viridis\", legend=False, ax=axes[0])\n", "axes[0].set(title=f\"{CITY.title()} — fi by time slot\", xlabel=\"fi\", ylabel=\"\")\n", "\n", "for slot, gdf in grids.items():\n", " sns.kdeplot(gdf.fi.dropna(), label=slot, ax=axes[1], linewidth=1.8)\n", "axes[1].set(title=\"Kernel density of fi\", xlabel=\"fi\")\n", "axes[1].legend(title=\"slot\")\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "8dee0a73", "metadata": {}, "outputs": [], "source": [ "# Interactive map — pan/zoom to inspect individual cells.\n", "grids[\"wk_am\"].explore(column=\"fi\", cmap=\"viridis\", scheme=\"quantiles\", k=5,\n", " tiles=\"CartoDB positron\", tooltip=[\"grid_name\", \"fi\", \"net_dis\"])" ] }, { "cell_type": "markdown", "id": "d32cf110", "metadata": {}, "source": [ "## 5 · Analysis 3b — Spatial statistics and neighbourhood clustering\n", "\n", "Following the reference spatial-statistics notebook: transform `fi` to tame its skew,\n", "build a k-nearest-neighbour spatial weights matrix, then compute global Moran's I,\n", "local Moran (LISA) clusters, and Getis-Ord \\(G_i^*\\) hot/cold spots." ] }, { "cell_type": "code", "execution_count": null, "id": "677936c1", "metadata": {}, "outputs": [], "source": [ "from libpysal import weights\n", "from esda.moran import Moran, Moran_Local\n", "from esda.getisord import G_Local\n", "from splot import esda as esdaplot\n", "\n", "SLOT = \"wk_am\"\n", "city = grids[SLOT].copy().reset_index(drop=True)\n", "\n", "# Double square root, as in the reference notebook: fi is heavily right-skewed and the\n", "# Moran statistics assume a roughly symmetric variable.\n", "city[\"fi_sqrt\"] = np.sqrt(city[\"fi\"])\n", "city[\"fi_2ndsqrt\"] = np.sqrt(city[\"fi_sqrt\"])\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(14, 4), constrained_layout=True)\n", "for ax, col in zip(axes, [\"fi\", \"fi_sqrt\", \"fi_2ndsqrt\"]):\n", " sns.histplot(city[col].dropna(), kde=True, ax=ax, color=\"#3182bd\")\n", " ax.set_title(f\"{col} (skew {city[col].skew():.2f})\")\n", "fig.suptitle(f\"{CITY.title()} {SLOT} — effect of the variance-stabilising transform\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "f65947ea", "metadata": {}, "outputs": [], "source": [ "# k=8 nearest neighbours, row-standardised (each row of weights sums to 1).\n", "w = weights.distance.KNN.from_dataframe(city, k=8)\n", "w.transform = \"R\"\n", "\n", "VAR = \"fi_2ndsqrt\"\n", "city[\"w_\" + VAR] = weights.lag_spatial(w, city[VAR])\n", "city[VAR + \"_std\"] = city[VAR] - city[VAR].mean()\n", "city[\"w_\" + VAR + \"_std\"] = weights.lag_spatial(w, city[VAR + \"_std\"])\n", "\n", "print(f\"{w.n} cells, mean neighbours {w.mean_neighbors:.1f}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "19c65ffe", "metadata": {}, "outputs": [], "source": [ "# Global Moran's I — is high frequency spatially clustered at all?\n", "values = city[VAR].to_numpy(dtype=float) # esda's numba kernels reject int input\n", "moran = Moran(values, w)\n", "print(f\"Moran's I = {moran.I:.4f} E[I] = {moran.EI:.4f} \"\n", " f\"z = {moran.z_sim:.2f} pseudo p = {moran.p_sim:.4f}\")\n", "\n", "# Moran scatterplot: value against its spatial lag, split into the four quadrants.\n", "fig, ax = plt.subplots(figsize=(6.5, 6.5), constrained_layout=True)\n", "sns.regplot(x=VAR + \"_std\", y=\"w_\" + VAR + \"_std\", data=city, ci=None,\n", " scatter_kws={\"s\": 8, \"alpha\": 0.5}, line_kws={\"color\": \"crimson\"}, ax=ax)\n", "ax.axvline(0, c=\"k\", alpha=0.5)\n", "ax.axhline(0, c=\"k\", alpha=0.5)\n", "\n", "# NOTE: Series.ptp() was removed in pandas 3; use np.ptp on the array.\n", "xr = city[VAR + \"_std\"].to_numpy()\n", "yr = city[\"w_\" + VAR + \"_std\"].to_numpy()\n", "for label, xf, yf in [(\"HH\", 0.75, 0.80), (\"HL\", 0.75, 0.10),\n", " (\"LH\", 0.10, 0.80), (\"LL\", 0.10, 0.10)]:\n", " ax.text(xr.min() + xf * np.ptp(xr), yr.min() + yf * np.ptp(yr),\n", " label, fontsize=18, color=\"r\", alpha=0.7)\n", "\n", "ax.set(title=f\"Moran scatterplot — {CITY.title()} {SLOT} (I = {moran.I:.3f})\",\n", " xlabel=f\"{VAR} (centred)\", ylabel=\"spatial lag\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "fcf0aa0a", "metadata": {}, "outputs": [], "source": [ "# Local Moran / LISA.\n", "lisa = Moran_Local(values, w)\n", "\n", "fig, axs = plt.subplots(2, 2, figsize=(13, 13), constrained_layout=True)\n", "axs = axs.flatten()\n", "\n", "# 1 — local statistic\n", "city.assign(Is=lisa.Is).plot(column=\"Is\", cmap=\"plasma\", scheme=\"quantiles\", k=5,\n", " edgecolor=\"white\", linewidth=0.1, legend=True,\n", " legend_kwds={\"fontsize\": 7}, ax=axs[0])\n", "axs[0].set_title(\"Local Moran's $I_i$\")\n", "\n", "# 2 — quadrant of every cell (p=1 forces all cells to be coloured)\n", "esdaplot.lisa_cluster(lisa, city, p=1, ax=axs[1])\n", "axs[1].set_title(\"Quadrant (all cells)\")\n", "\n", "# 3 — which cells are significant\n", "labels = pd.Series(np.where(lisa.p_sim < 0.05, \"Significant\", \"Non-significant\"),\n", " index=city.index)\n", "city.assign(cl=labels).plot(column=\"cl\", categorical=True, cmap=\"Paired\",\n", " edgecolor=\"white\", linewidth=0.1, legend=True,\n", " legend_kwds={\"fontsize\": 8}, ax=axs[2])\n", "axs[2].set_title(\"Significance (p < 0.05)\")\n", "\n", "# 4 — significant clusters only\n", "esdaplot.lisa_cluster(lisa, city, p=0.05, ax=axs[3])\n", "axs[3].set_title(\"Significant clusters (p < 0.05)\")\n", "\n", "for ax in axs:\n", " ax.set_axis_off()\n", "fig.suptitle(f\"LISA — {CITY.title()} {SLOT}, {VAR}\", fontsize=15)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "d8198aba", "metadata": {}, "outputs": [], "source": [ "# Cluster composition.\n", "# NOTE: pandas 3 removed the top-level pandas.value_counts() the reference notebook uses.\n", "SPOT_LABELS = {0: \"Non-significant\", 1: \"HH\", 2: \"LH\", 3: \"LL\", 4: \"HL\"}\n", "city[\"lisa_label\"] = pd.Series(lisa.q * (lisa.p_sim < 0.05), index=city.index).map(\n", " SPOT_LABELS)\n", "\n", "share = city[\"lisa_label\"].value_counts()\n", "print(f\"significant cells: {(lisa.p_sim < 0.05).mean() * 100:.1f}%\\n\")\n", "print(pd.DataFrame({\"cells\": share, \"share_%\": (share / len(city) * 100).round(1)}))" ] }, { "cell_type": "code", "execution_count": null, "id": "4787ab7b", "metadata": {}, "outputs": [], "source": [ "# Getis-Ord Gi and Gi* — hot spots (clusters of high service) and cold spots.\n", "g_i = G_Local(values, w)\n", "g_i_star = G_Local(values, w, star=True)\n", "\n", "\n", "def g_map(g, gdf, ax):\n", " \"\"\"Plot significant high (hot) and low (cold) clusters from a G_Local result.\"\"\"\n", " sig = g.p_sim < 0.05\n", " gdf.loc[~sig].plot(ax=ax, color=\"lightgrey\", edgecolor=\"0.8\", linewidth=0.1)\n", " gdf.loc[sig & (g.Zs > 0)].plot(ax=ax, color=\"red\", edgecolor=\"0.8\", linewidth=0.1)\n", " gdf.loc[sig & (g.Zs < 0)].plot(ax=ax, color=\"blue\", edgecolor=\"0.8\", linewidth=0.1)\n", " ax.set_title(f\"$G_i{'^*' if g.star else ''}$ — {CITY.title()} {SLOT}\", size=13)\n", " ax.set_axis_off()\n", "\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(13, 7), constrained_layout=True)\n", "for g, ax in zip([g_i, g_i_star], axes):\n", " g_map(g, city, ax)\n", "\n", "handles = [plt.Line2D([], [], marker=\"s\", linestyle=\"\", markersize=10, color=c, label=l)\n", " for c, l in [(\"red\", \"Hot spot (p<0.05)\"), (\"blue\", \"Cold spot (p<0.05)\"),\n", " (\"lightgrey\", \"Not significant\")]]\n", "axes[1].legend(handles=handles, loc=\"lower right\", fontsize=8)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "5b497aeb", "metadata": {}, "outputs": [], "source": [ "fig.savefig(OUT / f\"getis_ord_{CITY}_{SLOT}.png\", dpi=200, bbox_inches=\"tight\")\n", "print(f\"saved {OUT / f'getis_ord_{CITY}_{SLOT}.png'}\")" ] }, { "cell_type": "markdown", "id": "1248ea7a", "metadata": {}, "source": [ "## 6 · Appendix — reading the public GTFS feed without downloading it\n", "\n", "Access to the restricted files is granted **per dataset**. If you can open one ioerDATA\n", "replication package but not this one, that is a missing grant on this DOI, not a broken\n", "token — and no transfer method works around it, because streaming and downloading hit the\n", "same 403.\n", "\n", "What *is* public here are the two raw GTFS feeds the indicator was built from:\n", "`opnv.zip` (181 MB, local/regional transit) and `sbahn.zip` (7.2 MB). Because ioerDATA\n", "serves byte ranges, a zip's central directory can be read from its tail and individual\n", "members pulled out — listing all of `opnv.zip` costs about 65 KiB instead of 181 MB." ] }, { "cell_type": "code", "execution_count": null, "id": "c4f05256", "metadata": {}, "outputs": [], "source": [ "import struct\n", "import zlib\n", "\n", "GTFS_FEEDS = {\"opnv\": 6704, \"sbahn\": 6811}\n", "\n", "\n", "def _range(file_id, start, end):\n", " r = requests.get(f\"{BASE_URL}/api/access/datafile/{file_id}\",\n", " headers={**headers, \"Range\": f\"bytes={start}-{end}\"}, timeout=120)\n", " r.raise_for_status()\n", " return r.content\n", "\n", "\n", "def _size(file_id):\n", " r = requests.get(f\"{BASE_URL}/api/access/datafile/{file_id}\",\n", " headers={**headers, \"Range\": \"bytes=0-0\"}, timeout=60)\n", " r.raise_for_status()\n", " return int(r.headers[\"Content-Range\"].split(\"/\")[1])\n", "\n", "\n", "def remote_zip_index(file_id):\n", " \"\"\"Map member name -> (method, compressed_size, uncompressed_size, local_header_offset).\n", "\n", " Reads only the zip's end-of-central-directory and central directory.\n", " \"\"\"\n", " total = _size(file_id)\n", " tail = _range(file_id, max(0, total - 65536), total - 1)\n", " i = tail.rfind(b\"PK\\x05\\x06\")\n", " if i < 0:\n", " raise ValueError(\"no end-of-central-directory found; not a zip?\")\n", " cd_size, cd_off = struct.unpack(\" `stop_times.txt` is ~977 MB uncompressed and cannot be pulled apart by range — a zip\n", "> member is a single deflate stream, so any part of it requires the whole member. Computing\n", "> the frequency indicator from scratch means downloading the full feed. The point of this\n", "> section is exploration: you can inspect the feed's structure, agencies and routes for a\n", "> few tens of KiB before committing to that." ] }, { "cell_type": "markdown", "id": "4747a1a6", "metadata": {}, "source": [ "## Notes and caveats\n", "\n", "- **Restricted access is per dataset.** 1,178 of this dataset's 1,184 files are restricted;\n", " only `README.md`, the two PDFs and the two GTFS feeds are open. A token that can download\n", " another ioerDATA package says nothing about this one — check `canDownloadFile` (section 1)\n", " rather than inferring from a token that works elsewhere. Access is granted by the\n", " custodian, Sujit Sikder (s.sikder@ioer.de), via Request Access on the file page or\n", " `POST /api/access/datafile/{id}/requestAccess`.\n", "- **`fi` vs `pti`.** Per-city grids use `fi`; the aggregated INSPIRE grid CSV uses\n", " `pti_`. They are the same indicator under different names.\n", "- **Reference-notebook code does not run as-is** on current library versions:\n", " `plt.boxplot(labels=...)` was removed in matplotlib 3.11 (renamed `tick_labels=`, and\n", " `vert=` became `orientation=`), `pandas.value_counts()` was removed in pandas 3 (use\n", " `Series.value_counts()`), and `seaborn.distplot` is deprecated (use `histplot`).\n", " Because Colab and a local environment rarely ship the same versions, the plots here\n", " avoid the keywords that moved.\n", "- **If you join to city boundaries later**, zero-pad both sides of `ags` to 8 characters.\n", " IOER boundary layers store it as `'05911000'`; the statistics CSV drops the leading zero\n", " (`'5911000'`). An unpadded join silently loses every city in states 01-09 - in a test\n", " against 165 real city polygons it matched 24 of 165 instead of all of them.\n", "- **One week only.** The indicator covers 17–23 July 2022 — school-holiday season in most\n", " German states, so service levels may sit below the annual norm.\n", "- **KNN weights on polygons.** `KNN.from_dataframe` uses polygon centroids. On a regular\n", " 1 km grid, k=8 approximates a queen-contiguity neighbourhood; at a city's edge it reaches\n", " further, since the eight nearest cells may all lie inward.\n", "\n", "**Citation** — Sikder, S. K. (2025). *Urban Public Transit Frequency Indicator in Germany.*\n", "IOER Research Data Centre. [doi:10.71830/ABPCUS](https://doi.org/10.71830/ABPCUS). CC BY 4.0." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.7" } }, "nbformat": 4, "nbformat_minor": 5 }