renfe: generate shapes properly and consistently

- Update OSRM container to use ALL SPAIN (sorry, Trencelta)
- Generate a shape per trip (no trying to reuse, since trains that change stop sequence got wrong shapes)
- Add more position corrections for FEVE
- Run separate generators for FEVE and Renfe, since sometimes OSRM would pick the one that shouldn't and generate a wrong shape
- Add a debug script to generate a trip's visualisation from GTFS, since I was about to lose my mind debugging this pile of crap
- Update README (before starting anything else)

Time spent: ca. 6 hours

Closes #1

10 files changed, 500 insertions, 97 deletions

diff --git a/README.md b/README.md
index 4f593e9..469ee2e 100644
--- a/README.md
+++ b/README.md
@@ -12,7 +12,10 @@ Este repositorio incluye:
   - GTFS: Urbano de A Coruña
   - GTFS: Urbano de Ourense
   - GTFS: Urbano de Vigo
-- Submódulo Git para descargar y parchear los feeds de Renfe automáticamente
+- Feeds "de la comunidad":
+  - GTFS: Santiago de Compostela (WIP)
+  - GTFS: Ourense (WIP)
+  - GTFS: Lugo
 - Proxy del GTFS RealTime de Renfe para integración con OTP
 - Configuración de OpenTripPlanner para cargar los datos descargados y el tiempo real de Renfe
 - Tareas para ejecutar OTP directamente
@@ -21,7 +24,7 @@ Este repositorio incluye:
 
 Para ejecutar los contenidos de este repositorio, es necesario tener descargado:
 
-- **Java 21 LTS**: en sistemas Windows o macOS, se recomienda descargar de [Adoptium](https://adoptium.net/). En sistemas Linux, se puede instalar OpenJDK 21 desde los repositorios oficiales.
+- **Java 25 LTS**: en sistemas Windows o macOS, se recomienda descargar de [Adoptium](https://adoptium.net/). En sistemas Linux, se puede instalar OpenJDK 25 desde los repositorios oficiales.
 - **Task**: Se recomienda tener [Task](https://taskfile.dev/) instalado para gestionar las tareas definidas en el `Taskfile.yml`. Alternativamente, se pueden ejecutar directamente, pero es más sencillo con Task.
 - **Python y `uv`**: Para el proxy de tiempo real de Renfe, es necesario tener Python, con los paquetes `Flask`, `requests` y `protobuf` instalados. Utilizando [`uv`](https://docs.astral.sh/uv) se puede con solo una línea. El Taskfile asume que se está utilizando `uv` directamente.
 - **Clave de API del NAP del Ministerio de Transportes**: Para poder descargar los feeds disponibles en el Punto de Acceso Nacional (NAP), es necesario registrarse y obtener una clave de API en [https://nap.transportes.gob.es/](https://nap.transportes.gob.es/).
@@ -31,12 +34,6 @@ Para descargar los datos y ejecutar OTP, se pueden utilizar las siguientes tarea
 ```bash
 git clone https://github.com/tpgalicia/opentripplanner-galicia.git
 cd opentripplanner-galicia
-git submodule update --init --recursive
-```
-
-Descargar OpenTripPlanner, datos de OpenStreetMap y feeds:
-
-```bash
 task setup
 task download-feeds NAP_API_KEY=tu_clave_de_api_aqui
 ```
diff --git a/Taskfile.yml b/Taskfile.yml
index 9471070..14165b5 100644
--- a/Taskfile.yml
+++ b/Taskfile.yml
@@ -23,8 +23,7 @@ tasks:
     desc: "Download Renfe GTFS data (NAP MITRAMS)"
     cmds:
       - uv --directory build_renfe run ./build_static_feed.py {{.NAP_API_KEY}} --merge
-      - cp build_renfe/gtfs_renfe_galicia_general*.zip feeds/renfe.zip
-      - cp build_renfe/gtfs_renfe_galicia_cercanias*.zip feeds/feve.zip
+      - cp build_renfe/gtfs_renfe_galicia_merged.zip feeds/renfe.zip
 
   download-tranvias:
     desc: "Download A Coruña city GTFS data (NAP MITRAMS)"
diff --git a/build_renfe/Dockerfile b/build_renfe/Dockerfile
index f565320..02a8480 100644
--- a/build_renfe/Dockerfile
+++ b/build_renfe/Dockerfile
@@ -1,25 +1,42 @@
 # Use a multi-stage build to download necessary files
 FROM alpine/curl AS downloader
 
-RUN curl -L https://download.geofabrik.de/europe/spain/galicia-latest.osm.pbf -o /galicia-latest.osm.pbf
-RUN curl -L https://raw.githubusercontent.com/railnova/osrm-train-profile/refs/heads/master/basic.lua -o /opt/train.lua
+RUN curl -L https://download.geofabrik.de/europe/spain-latest.osm.pbf -o /spain-latest.osm.pbf
 
-FROM osrm/osrm-backend
+FROM osrm/osrm-backend AS builder
 
 # Copy the downloaded OSM file from the downloader stage
-COPY --from=downloader /galicia-latest.osm.pbf /data/galicia-latest.osm.pbf
-COPY --from=downloader /opt/train.lua /opt/train.lua
+RUN mkdir -p /data/standard /data/narrow
+COPY --from=downloader /spain-latest.osm.pbf /data/standard/spain-latest.osm.pbf
+COPY --from=downloader /spain-latest.osm.pbf /data/narrow/spain-latest.osm.pbf
+COPY ./train_standard.lua /opt/train_standard.lua
+COPY ./train_narrow.lua /opt/train_narrow.lua
 
 # Extract the map data using osrm-train-profile (by Railnova)
-RUN osrm-extract -p /opt/train.lua /data/galicia-latest.osm.pbf
+RUN osrm-extract -p /opt/train_standard.lua /data/standard/spain-latest.osm.pbf
+RUN osrm-partition /data/standard/spain-latest.osrm
 
 # Prepare the map data for routing
-RUN osrm-partition /data/galicia-latest.osrm
-RUN osrm-customize /data/galicia-latest.osrm
+RUN osrm-extract -p /opt/train_narrow.lua /data/narrow/spain-latest.osm.pbf
+RUN osrm-partition /data/narrow/spain-latest.osrm
+
+RUN osrm-customize /data/standard/spain-latest.osrm
+RUN osrm-customize /data/narrow/spain-latest.osrm
+
+RUN rm /data/standard/spain-latest.osm.pbf
+RUN rm /data/narrow/spain-latest.osm.pbf
+
+FROM osrm/osrm-backend
+
+RUN mkdir -p /data/standard /data/narrow
+COPY --from=builder /data/standard/spain-latest.osrm* /data/standard/
+COPY --from=builder /data/narrow/spain-latest.osrm* /data/narrow/
 
 # Expose the OSRM server port
-EXPOSE 5000
+EXPOSE 5000 5001
 
 # Start the OSRM server
-CMD ["osrm-routed", "--algorithm", "mld", "/data/galicia-latest.osrm"]
-
+COPY ./start.sh /start.sh
+RUN chmod +x /start.sh
+EXPOSE 5000 5001
+CMD ["/start.sh"]
\ No newline at end of file
diff --git a/build_renfe/build_static_feed.py b/build_renfe/build_static_feed.py
index eb247a9..6c12c1f 100644
--- a/build_renfe/build_static_feed.py
+++ b/build_renfe/build_static_feed.py
@@ -16,6 +16,7 @@ import os
 import shutil
 import tempfile
 import zipfile
+import binascii
 
 import pandas as pd
 
@@ -181,13 +182,20 @@ if __name__ == "__main__":
         help="NAP API Key (https://nap.transportes.gob.es/)"
     )
     parser.add_argument(
-        "--osrm-url",
+        "--osrm-std",
         type=str,
-        help="OSRM server URL",
+        help="OSRM standard server URL",
         default="http://localhost:5050",
         required=False,
     )
     parser.add_argument(
+        "--osrm-narrow",
+        type=str,
+        help="OSRM narrow gauge server URL",
+        default="http://localhost:5051",
+        required=False,
+    )
+    parser.add_argument(
         "--debug",
         help="Enable debug logging",
         action="store_true"
@@ -201,8 +209,9 @@ if __name__ == "__main__":
     args = parser.parse_args()
 
     try:
-        osrm_check = requests.head(args.osrm_url, timeout=5)
-        GENERATE_SHAPES = osrm_check.status_code < 500
+        osrm_check_std = requests.head(args.osrm_std, timeout=5)
+        osrm_check_narrow = requests.head(args.osrm_narrow, timeout=5)
+        GENERATE_SHAPES = osrm_check_std.status_code < 500 and osrm_check_narrow.status_code < 500
     except requests.RequestException:
         GENERATE_SHAPES = False
         logging.warning("OSRM server is not reachable. Shape generation will be skipped.")
@@ -214,6 +223,10 @@ if __name__ == "__main__":
     )
 
     for feed in FEEDS.keys():
+        def get_shape_id(trip_id: str) -> str:
+            trip_crc = binascii.crc32(trip_id.encode("utf-8"))
+            return f"Shape_{feed}_{trip_crc}_{trip_crc}"
+
         INPUT_GTFS_FD, INPUT_GTFS_ZIP = tempfile.mkstemp(suffix=".zip", prefix=f"renfe_galicia_in_{feed}_")
         INPUT_GTFS_PATH = tempfile.mkdtemp(prefix=f"renfe_galicia_in_{feed}_")
         OUTPUT_GTFS_PATH = tempfile.mkdtemp(prefix=f"renfe_galicia_out_{feed}_")
@@ -397,7 +410,7 @@ if __name__ == "__main__":
 
         for tig in trips_in_galicia:
             if GENERATE_SHAPES:
-                tig["shape_id"] = f"Shape_{tig['trip_id'][0:5]}"
+                tig["shape_id"] = get_shape_id(tig["trip_id"])
             tig["trip_headsign"] = stops_by_id[last_stop_in_trips[tig["trip_id"]]]["stop_name"]
         with open(
             os.path.join(OUTPUT_GTFS_PATH, "trips.txt"),
@@ -423,12 +436,19 @@ if __name__ == "__main__":
         logging.info("GTFS data for Galicia has been extracted successfully. Generate shapes for the trips...")
 
         if GENERATE_SHAPES:
-            shape_ids_total = len(set(f"Shape_{trip_id[0:5]}" for trip_id in trip_ids))
+            shape_ids_total = len(set(trip["shape_id"] for trip in trips_in_galicia))
             shape_ids_generated: set[str] = set()
 
             # Pre-load stops for quick lookup
             stops_dict = {stop["stop_id"]: stop for stop in stops_in_trips}
 
+            # Map trip_id to route_type for OSRM profile selection
+            trip_to_route_type = {tig["trip_id"]: tig.get("route_type", "2") for tig in trips_in_galicia}
+            # Fallback if not in trips_in_galicia (shouldn't happen)
+            route_id_to_type = {r["route_id"]: r["route_type"] for r in routes_in_trips}
+            for tig in trips_in_galicia:
+                trip_to_route_type[tig["trip_id"]] = route_id_to_type.get(tig["route_id"], "2")
+
             # Group stop times by trip_id to avoid repeated file reads
             stop_times_by_trip: dict[str, list[dict]] = {}
             for st in stop_times_in_galicia:
@@ -437,11 +457,21 @@ if __name__ == "__main__":
                     stop_times_by_trip[tid] = []
                 stop_times_by_trip[tid].append(st)
 
-            OSRM_BASE_URL = f"{args.osrm_url}/route/v1/driving/"
+            shapes_file = open("shapes_debug.txt", "w", encoding="utf-8")
+
             for trip_id in tqdm(trip_ids, total=shape_ids_total, desc="Generating shapes"):
-                shape_id = f"Shape_{trip_id[0:5]}"
+                shape_id = get_shape_id(trip_id)
                 if shape_id in shape_ids_generated:
                     continue
+                
+                route_type = trip_to_route_type.get(trip_id, "2")
+                osrm_profile = "driving" 
+                
+                # If we are on feed cercanias or the 5-digit trip ID starts with 7, it's a narrow gauge train, otherwise standard gauge
+                if feed == "cercanias" or (len(trip_id) >= 5 and trip_id[0] == "7"):
+                    OSRM_BASE_URL = f"{args.osrm_narrow}/route/v1/driving/"
+                else:
+                    OSRM_BASE_URL = f"{args.osrm_std}/route/v1/driving/"
 
                 stop_seq = stop_times_by_trip.get(trip_id, [])
                 stop_seq.sort(key=lambda x: int(x["stop_sequence"].strip()))
@@ -450,79 +480,36 @@ if __name__ == "__main__":
                     continue
 
                 final_shape_points = []
-                i = 0
-                while i < len(stop_seq) - 1:
-                    stop_a = stops_dict[stop_seq[i]["stop_id"]]
-                    lat_a, lon_a = float(stop_a["stop_lat"]), float(stop_a["stop_lon"])
+                coordinates = []
+                for st in stop_seq:
+                    s = stops_dict[st["stop_id"]]
+                    coordinates.append(f"{s['stop_lon']},{s['stop_lat']}")
 
-                    if not is_in_bounds(lat_a, lon_a):
-                        # S_i is out of bounds. Segment S_i -> S_{i+1} is straight line.
-                        stop_b = stops_dict[stop_seq[i+1]["stop_id"]]
-                        lat_b, lon_b = float(stop_b["stop_lat"]), float(stop_b["stop_lon"])
+                coords_str = ";".join(coordinates)
+                osrm_url = f"{OSRM_BASE_URL}{coords_str}?overview=full&geometries=geojson&continue_straight=false"
 
-                        segment_points = [[lon_a, lat_a], [lon_b, lat_b]]
-                        if not final_shape_points:
-                            final_shape_points.extend(segment_points)
-                        else:
-                            final_shape_points.extend(segment_points[1:])
-                        i += 1
-                    else:
-                        # S_i is in bounds. Find how many subsequent stops are also in bounds.
-                        j = i + 1
-                        while j < len(stop_seq):
-                            stop_j = stops_dict[stop_seq[j]["stop_id"]]
-                            if is_in_bounds(float(stop_j["stop_lat"]), float(stop_j["stop_lon"])):
-                                j += 1
-                            else:
-                                break
-
-                        # Stops from i to j-1 are in bounds.
-                        if j > i + 1:
-                            # We have at least two consecutive stops in bounds.
-                            in_bounds_stops = stop_seq[i:j]
-                            coordinates = []
-                            for st in in_bounds_stops:
-                                s = stops_dict[st["stop_id"]]
-                                coordinates.append(f"{s['stop_lon']},{s['stop_lat']}")
-
-                            coords_str = ";".join(coordinates)
-                            osrm_url = f"{OSRM_BASE_URL}{coords_str}?overview=full&geometries=geojson"
+                shapes_file.write(f"{trip_id} ({shape_id}): {osrm_url}\n")
 
-                            segment_points = []
-                            try:
-                                response = requests.get(osrm_url, timeout=10)
-                                if response.status_code == 200:
-                                    data = response.json()
-                                    if data.get("code") == "Ok":
-                                        segment_points = data["routes"][0]["geometry"]["coordinates"]
-                            except Exception:
-                                pass
-
-                            if not segment_points:
-                                # Fallback to straight lines for this whole sub-sequence
-                                segment_points = []
-                                for k in range(i, j):
-                                    s = stops_dict[stop_seq[k]["stop_id"]]
-                                    segment_points.append([float(s["stop_lon"]), float(s["stop_lat"])])
-
-                            if not final_shape_points:
-                                final_shape_points.extend(segment_points)
-                            else:
-                                final_shape_points.extend(segment_points[1:])
-
-                            i = j - 1 # Next iteration starts from S_{j-1}
+                try:
+                    response = requests.get(osrm_url, timeout=10)
+                    if response.status_code == 200:
+                        data = response.json()
+                        if data.get("code") == "Ok":
+                            final_shape_points = data["routes"][0]["geometry"]["coordinates"]
+                            logging.debug(f"OSRM success for {shape_id} ({len(coordinates)} stops): {len(final_shape_points)} points")
                         else:
-                            # Only S_i is in bounds, S_{i+1} is out.
-                            # Segment S_i -> S_{i+1} is straight line.
-                            stop_b = stops_dict[stop_seq[i+1]["stop_id"]]
-                            lat_b, lon_b = float(stop_b["stop_lat"]), float(stop_b["stop_lon"])
+                            logging.warning(f"OSRM returned error code {data.get('code')} for {shape_id} with {len(coordinates)} stops")
+                    else:
+                        logging.warning(f"OSRM request failed for {shape_id} with status {response.status_code}")
+                except Exception as e:
+                    logging.error(f"OSRM exception for {shape_id}: {str(e)}")
 
-                            segment_points = [[lon_a, lat_a], [lon_b, lat_b]]
-                            if not final_shape_points:
-                                final_shape_points.extend(segment_points)
-                            else:
-                                final_shape_points.extend(segment_points[1:])
-                            i += 1
+                if not final_shape_points:
+                    # Fallback to straight lines
+                    logging.info(f"Using straight-line fallback for {shape_id}")
+                    for st in stop_seq:
+                        s = stops_dict[st["stop_id"]]
+                        final_shape_points.append([float(s["stop_lon"]), float(s["stop_lat"])])
 
                 shape_ids_generated.add(shape_id)
 
diff --git a/build_renfe/compose.yml b/build_renfe/compose.yml
index ebe29cf..312d15d 100644
--- a/build_renfe/compose.yml
+++ b/build_renfe/compose.yml
@@ -5,3 +5,4 @@ services:
     restart: unless-stopped
     ports:
       - "5050:5000"
+      - "5051:5001"
diff --git a/build_renfe/start.sh b/build_renfe/start.sh
new file mode 100644
index 0000000..2e4ddbd
--- /dev/null
+++ b/build_renfe/start.sh
@@ -0,0 +1,3 @@
+#!/bin/sh
+osrm-routed --algorithm mld --port 5000 /data/standard/spain-latest.osrm &
+osrm-routed --algorithm mld --port 5001 /data/narrow/spain-latest.osrm
\ No newline at end of file
diff --git a/build_renfe/stop_overrides.json b/build_renfe/stop_overrides.json
index 5f4b5b8..4ebd4be 100644
--- a/build_renfe/stop_overrides.json
+++ b/build_renfe/stop_overrides.json
@@ -1,5 +1,11 @@
 [
   {
+    "stop_id": "08223",
+    "stop_name": "Vigo Urzáiz",
+    "stop_lat": 42.234960,
+    "stop_lon": -8.711909
+  },
+  {
     "stop_id": "31412",
     "stop_lat": 43.3504,
     "stop_lon": -8.412142,
@@ -91,5 +97,36 @@
     "stop_id": "99161",
     "stop_name": "Pontevedra Turístico",
     "_delete": true
+  },
+  {
+    "stop_id": "05156",
+    "stop_name": "Viveiro - Apeadero",
+    "stop_lat": 43.66307948733166,
+    "stop_lon": -7.592119972501223
+  },
+  {
+    "stop_id": "05299",
+    "stop_name": "Castropol",
+    "stop_lat": 43.50879091456042,
+    "stop_lon": -7.013783988194376
+  },
+  {
+    "stop_id": "05265",
+    "stop_name": "Tablizo",
+    "stop_lat": 43.548132589750224,
+    "stop_lon": -6.351521338297374
+  },
+  {
+    "stop_id": "05263",
+    "stop_name": "Ballota",
+    "stop_lat": 43.55092647470215,
+    "stop_lon": -6.322749469667542
+  },
+  {
+    "stop_id": "15211",
+    "stop_name": "Oviedo/Uviéu",
+    "stop_lat": 43.366680663195496,
+    "stop_lon": -5.855315584365599
+          
   }
 ]
diff --git a/build_renfe/train_narrow.lua b/build_renfe/train_narrow.lua
new file mode 100644
index 0000000..73032d2
--- /dev/null
+++ b/build_renfe/train_narrow.lua
@@ -0,0 +1,134 @@
+-- Copyright 2017-2019 Railnova SA <support@railnova.eu>, Nikita Marchant <nikita.marchant@gmail.com>
+-- Code under the 2-clause BSD license
+
+api_version = 4
+
+function setup()
+  return {
+    properties = {
+      max_speed_for_map_matching     = 220/3.6, -- speed conversion to m/s
+      weight_name                    = 'routability',
+      left_hand_driving              = true,
+      u_turn_penalty                 = 60 * 2, -- 2 minutes to change cabin
+      turn_duration                  = 20,
+      continue_straight_at_waypoint  = false,
+      max_angle                      = 30,
+
+      secondary_speed                = 30,
+      speed                          = 160,
+    },
+
+    default_mode              = mode.train,
+    default_speed             = 120,
+}
+
+end
+
+
+function ternary ( cond , T , F )
+    if cond then return T else return F end
+end
+
+
+function process_node(profile, node, result, relations)
+    local railway = node:get_value_by_key("railway")
+
+    -- refuse railway nodes that we cannot go through
+    result.barrier = (
+        railway == "buffer_stop" or
+        railway == "derail"
+    )
+    result.traffic_lights = false
+end
+
+function process_way(profile, way, result, relations)
+    local data = {
+        railway = way:get_value_by_key("railway"),
+        service = way:get_value_by_key("service"),
+        usage = way:get_value_by_key("usage"),
+        maxspeed = way:get_value_by_key("maxspeed"),
+        gauge = way:get_value_by_key("gauge"),
+    }
+
+    -- Remove everything that is not railway
+    if not data.railway then
+        return
+    -- Remove everything that is not a rail, a turntable, a traverser
+    elseif (
+        data.railway ~= 'rail' and
+        data.railway ~= 'turntable' and
+        data.railway ~= 'traverser' and
+        data.railway ~= 'narrow_gauge' and
+        data.railway  ~= 'ferry'
+    ) then
+        return
+    -- Remove military and tourism rails
+    elseif (
+        data.usage == "military" or
+        data.usage == "tourism"
+    ) then
+        return
+    -- Keep only most common gauges (and undefined)
+    -- uses .find() as some gauges are specified like "1668;1435"
+    elseif (
+        data.gauge ~= nil and
+        data.gauge ~= 1000 and not string.find(data.gauge, "1000")
+   ) then
+        return
+    end
+
+    local is_secondary = (
+        data.service == "siding" or
+        data.service == "spur" or
+        data.service == "yard" or
+        data.usage == "industrial"
+    )
+
+    -- by default, use 30km/h for secondary rails, else 160
+    local default_speed = ternary(is_secondary, profile.properties.secondary_speed, profile.properties.speed)
+    -- but is OSM specifies a maxspeed, use the one from OSM
+    local speed = ternary(data.maxspeed, data.maxspeed, default_speed)
+
+   -- fix speed for mph issue
+    speed = tostring(speed)
+    if speed:find(" mph") or speed:find("mph") then
+      speed = speed:gsub(" mph", "")
+      speed = speed:gsub("mph", "")
+        speed = tonumber (speed)
+        if speed == nil then speed = 20 end
+	speed = speed * 1.609344
+    else
+     speed = tonumber (speed)
+    end
+    -- fix speed for mph issue end
+
+
+    result.forward_speed = speed
+    result.backward_speed = speed
+    --
+    result.forward_mode = mode.train
+    result.backward_mode = mode.train
+    --
+    result.forward_rate = 1
+    result.backward_rate = 1
+
+end
+
+function process_turn(profile, turn)
+    -- Refuse truns that have a big angle
+    if math.abs(turn.angle) >  profile.properties.max_angle then
+        return
+    end
+
+    -- If we go backwards, add the penalty to change cabs
+    if turn.is_u_turn then
+      turn.duration = turn.duration + profile.properties.u_turn_penalty
+    end
+end
+
+return {
+    setup = setup,
+    process_way = process_way,
+    process_node = process_node,
+    process_turn = process_turn
+}
\ No newline at end of file
diff --git a/build_renfe/train_standard.lua b/build_renfe/train_standard.lua
new file mode 100644
index 0000000..8f424b1
--- /dev/null
+++ b/build_renfe/train_standard.lua
@@ -0,0 +1,136 @@
+-- Copyright 2017-2019 Railnova SA <support@railnova.eu>, Nikita Marchant <nikita.marchant@gmail.com>
+-- Code under the 2-clause BSD license
+
+api_version = 4
+
+function setup()
+  return {
+    properties = {
+      max_speed_for_map_matching     = 220/3.6, -- speed conversion to m/s
+      weight_name                    = 'routability',
+      left_hand_driving              = true,
+      u_turn_penalty                 = 60 * 2, -- 2 minutes to change cabin
+      turn_duration                  = 20,
+      continue_straight_at_waypoint  = false,
+      max_angle                      = 30,
+
+      secondary_speed                = 30,
+      speed                          = 160,
+    },
+
+    default_mode              = mode.train,
+    default_speed             = 120,
+}
+
+end
+
+
+function ternary ( cond , T , F )
+    if cond then return T else return F end
+end
+
+
+function process_node(profile, node, result, relations)
+    local railway = node:get_value_by_key("railway")
+
+    -- refuse railway nodes that we cannot go through
+    result.barrier = (
+        railway == "derail"
+    )
+    result.traffic_lights = false
+end
+
+function process_way(profile, way, result, relations)
+    local data = {
+        railway = way:get_value_by_key("railway"),
+        service = way:get_value_by_key("service"),
+        usage = way:get_value_by_key("usage"),
+        maxspeed = way:get_value_by_key("maxspeed"),
+        gauge = way:get_value_by_key("gauge"),
+    }
+
+    -- Remove everything that is not railway
+    if not data.railway then
+        return
+    -- Remove everything that is not a rail, a turntable, a traverser
+    elseif (
+        data.railway ~= 'rail' and
+        data.railway ~= 'turntable' and
+        data.railway ~= 'traverser' and
+        data.railway  ~= 'ferry'
+    ) then
+        return
+    -- Remove military and tourism rails
+    elseif (
+        data.usage == "military" or
+        data.usage == "tourism"
+    ) then
+        return
+    -- Keep only most common gauges (and undefined)
+    -- uses .find() as some gauges are specified like "1668;1435"
+    elseif (
+        data.gauge ~= nil and
+        data.gauge ~= 1435 and not string.find(data.gauge, "1435") and
+        data.gauge ~= 1520 and not string.find(data.gauge, "1520") and
+        data.gauge ~= 1524 and not string.find(data.gauge, "1524") and
+        data.gauge ~= 1600 and not string.find(data.gauge, "1600") and
+        data.gauge ~= 1668 and not string.find(data.gauge, "1668")
+   ) then
+        return
+    end
+
+    local is_secondary = (
+        data.service == "siding" or
+        data.service == "spur" or
+        data.service == "yard" or
+        data.usage == "industrial"
+    )
+
+    -- by default, use 30km/h for secondary rails, else 160
+    local default_speed = ternary(is_secondary, profile.properties.secondary_speed, profile.properties.speed)
+    -- but is OSM specifies a maxspeed, use the one from OSM
+    local speed = ternary(data.maxspeed, data.maxspeed, default_speed)
+
+   -- fix speed for mph issue
+    speed = tostring(speed)
+    if speed:find(" mph") or speed:find("mph") then
+      speed = speed:gsub(" mph", "")
+      speed = speed:gsub("mph", "")
+        speed = tonumber (speed)
+        if speed == nil then speed = 20 end
+	speed = speed * 1.609344
+    else
+     speed = tonumber (speed)
+    end
+    -- fix speed for mph issue end
+
+
+    result.forward_speed = speed
+    result.backward_speed = speed
+    --
+    result.forward_mode = mode.train
+    result.backward_mode = mode.train
+    --
+    result.forward_rate = 1
+    result.backward_rate = 1
+
+end
+
+function process_turn(profile, turn)
+    -- Refuse truns that have a big angle
+    if math.abs(turn.angle) >  profile.properties.max_angle then
+        return
+    end
+
+    -- If we go backwards, add the penalty to change cabs
+    if turn.is_u_turn then
+      turn.duration = turn.duration + profile.properties.u_turn_penalty
+    end
+end
+
+return {
+    setup = setup,
+    process_way = process_way,
+    process_node = process_node,
+    process_turn = process_turn
+}
\ No newline at end of file
diff --git a/trip_geo.py b/trip_geo.py
new file mode 100644
index 0000000..c6d88db
--- /dev/null
+++ b/trip_geo.py
@@ -0,0 +1,92 @@
+# /// script
+# requires-python = ">=3.13"
+# dependencies = [
+#     "pandas",
+# ]
+# ///
+
+import argparse
+import os
+import pandas as pd
+
+if __name__ != "__main__":
+    raise RuntimeError("This script is meant to be run as a standalone program, not imported as a module.")
+
+parser = argparse.ArgumentParser(description="Extract GeoJSON from GTFS feed and save to file")
+parser.add_argument(
+    "gtfs_path",
+    type=str,
+    help="Path to the GTFS feed directory"
+)
+
+parser.add_argument(
+    "trip_id",
+    type=str,
+    help="ID of the trip to extract from the GTFS feed"
+)
+
+args = parser.parse_args()
+
+gtfs_path = args.gtfs_path
+trip_id = args.trip_id
+
+# Load trips.txt, stop_times.txt, stops.txt and shapes.txt
+
+trips_df = pd.read_csv(os.path.join(gtfs_path, "trips.txt"))
+stop_times_df = pd.read_csv(os.path.join(gtfs_path, "stop_times.txt"))
+stops_df = pd.read_csv(os.path.join(gtfs_path, "stops.txt"))
+shapes_df = pd.read_csv(os.path.join(gtfs_path, "shapes.txt"))
+
+# Find the shape_id for the given trip_id
+trip_row = trips_df[trips_df["trip_id"] == trip_id]
+if trip_row.empty:
+    raise ValueError(f"Trip ID {trip_id} not found in trips.txt")
+
+shape_id = trip_row.iloc[0]["shape_id"]
+if pd.isna(shape_id):
+    raise ValueError(f"No shape_id found for Trip ID {trip_id}")
+
+# Extract the shape points for the shape_id
+shape_points = shapes_df[shapes_df["shape_id"] == shape_id].sort_values(by="shape_pt_sequence")
+
+# Find the stop sequence for the trip_id and get the stop coordinates
+stop_times = stop_times_df[stop_times_df["trip_id"] == trip_id].sort_values(by="stop_sequence")
+stop_ids = stop_times["stop_id"].tolist()
+stops = stops_df[stops_df["stop_id"].isin(stop_ids)]
+
+# Convert shape points to GeoJSON LineString format
+geojson = {
+    "type": "FeatureCollection",
+    "features": [
+        {
+            "type": "Feature",
+            "geometry": {
+                "type": "LineString",
+                "coordinates": shape_points.apply(lambda row: [row["shape_pt_lon"], row["shape_pt_lat"]], axis=1).tolist()
+            },
+            "properties": {
+                "headsign": trip_row.iloc[0]["trip_headsign"],
+                "shape_id": shape_id
+            }
+        },
+        *[{
+            "type": "Feature",
+            "geometry": {
+                "type": "Point",
+                "coordinates": [x.stop_lon, x.stop_lat]
+            },
+            "properties": {
+                "name": x.stop_name,
+                "stop_id": x.stop_id,
+                "code": x.stop_code
+            }
+        } for _, x in stops.iterrows()]
+    ]
+}
+
+# Save GeoJSON to file
+output_file = f"{trip_id}_shape.geojson"
+
+with open(output_file, "w") as f:
+    import json
+    json.dump(geojson, f, indent=2)