Build an LRS based on swiss opendata

build_swiss.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating an LRS for Swiss rail\n",
    "\n",
    "This notebook shows how to build an LRS for the swiss railway network using libLRS and opendata from the swiss federal railway.\n",
    "\n",
    "We need two data sources:\n",
    "- the [geometry of lines](https://data.sbb.ch/explore/dataset/linie-mit-polygon)\n",
    "- the [kilometrical reference points](https://data.sbb.ch/explore/dataset/linienkilometrierung/)\n",
    "\n",
    "We’ll make some simpflications:\n",
    "- the positionning will be at the line, not the track\n",
    "- we consider that there is alway 1km between two reference points\n",
    "\n",
    "## Building the geometry\n",
    "\n",
    "The geometry of the lines are provided by segments. We group them and order them to generate a single linestring for each line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "from typing import DefaultDict\n",
    "import json\n",
    "\n",
    "lines = defaultdict(list)\n",
    "with open(\"linie-mit-polygon.geojson\") as lines_geojson:\n",
    "    segments = json.loads(lines_geojson.read())[\"features\"]\n",
    "    for segment in segments:\n",
    "        lines[segment[\"properties\"][\"linienr\"]].append(segment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from liblrs_python import Builder, Point, SegmentOfTraversal\n",
    "\n",
    "builder = Builder()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance(a:Point, b:Point) -> float:\n",
    "    return (a.x - b.x)**2 + (a.y - b.y)**2\n",
    "\n",
    "traversals = {}\n",
    "lrm_name = {}\n",
    "for id, segments in lines.items():\n",
    "    segments.sort(key=lambda l: l[\"properties\"][\"km_agm_von\"])\n",
    "    segment_handles = []\n",
    "    last_coord = None\n",
    "    for segment in segments:\n",
    "        coords = [Point(c[0], c[1]) for c in segment[\"geometry\"][\"coordinates\"]]\n",
    "        props = segment[\"properties\"]\n",
    "        segment_id = f\"{props[\"linienr\"]}-{props[\"km_agm_von\"]}-{props[\"km_agm_bis\"]}\"\n",
    "        segment_handle = builder.add_segment(segment_id, coords, 0, 1)\n",
    "\n",
    "        reversed = (last_coord is not None) and distance(last_coord, coords[0]) > distance(last_coord, coords[-1])\n",
    "        last_coord = coords[-1] if reversed else coords[0]\n",
    "\n",
    "        segment_handles.append(SegmentOfTraversal(segment_handle, reversed))\n",
    "        \n",
    "\n",
    "\n",
    "    lrm_name[id] = segments[0][\"properties\"][\"liniename\"]\n",
    "    traversals[id] = builder.add_traversal(lrm_name[id], segment_handles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import modf\n",
    "from pathlib import Path\n",
    "from liblrs_python import AnchorOnLrm\n",
    "\n",
    "\n",
    "with open(\"linienkilometrierung.geojson\") as milestones_geojson:\n",
    "    milestones = json.loads(milestones_geojson.read())[\"features\"]\n",
    "\n",
    "    anchors = defaultdict(list)\n",
    "    for milestone in milestones:\n",
    "        km = milestone[\"properties\"][\"km\"]\n",
    "        linienr = milestone[\"properties\"][\"linienr\"]\n",
    "        offset, ref = modf(km)\n",
    "        name = str(int(ref))\n",
    "        geom = milestone[\"geometry\"][\"coordinates\"]\n",
    "        # We only create named anchors for round km\n",
    "        if offset < 1e-3:\n",
    "            anchor = builder.add_anchor(id=name, coord=Point(geom[0], geom[1]), properties={}, name=name)\n",
    "            anchors[linienr].append(AnchorOnLrm(anchor, km*1000))\n",
    "        else:\n",
    "            pass\n",
    "            # we could add unnamed anchors for better precisions, we skip for now to save time when iterationg\n",
    "            #anchor = builder.add_anchor(id=str(km), coord=Point(geom[0], geom[1]), properties={}, name=None)\n",
    "        \n",
    "    for linienr, anchors_on_lrm in anchors.items():\n",
    "        if linienr in traversals:\n",
    "            builder.add_lrm(lrm_name[linienr], traversals[linienr], anchors_on_lrm, {})\n",
    "        else:\n",
    "            print(f\"Warning: could not find line {linienr}\")\n",
    "\n",
    "    builder.save(Path(\"swiss.lrs.bin\"), {})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "osrd-dags-Z5JFnb7w-py3.12",
   "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.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}