From 9442a623fe191048b544e150a8ec9c3527222e0b Mon Sep 17 00:00:00 2001
From: Daniel Schadt <kingdread@gmx.de>
Date: Wed, 15 Oct 2025 21:42:19 +0200
Subject: initial geo.Path implementation

Since we want to move GPX data into the database, we need to do all the
things that gpxpy currently does for us, including the length and speed
computations. This is the start.
---
 fietsboek/geo.py | 132 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 132 insertions(+)
 create mode 100644 fietsboek/geo.py

diff --git a/fietsboek/geo.py b/fietsboek/geo.py
new file mode 100644
index 0000000..c665705
--- /dev/null
+++ b/fietsboek/geo.py
@@ -0,0 +1,132 @@
+from dataclasses import dataclass
+from itertools import islice
+from math import sqrt, sin, cos, radians
+
+import gpxpy
+
+
+# WGS-84 equatorial radius, also called the semi-major axis.
+# https://en.wikipedia.org/wiki/Earth_radius
+EARTH_RADIUS = 6378137.0
+"""Radius of the earth, in meters."""
+
+# https://en.wikipedia.org/wiki/Preferred_walking_speed
+MOVING_THRESHOLD = 1.1
+"""Speed which is considered to be the moving threshold, in m/s."""
+
+
+@dataclass
+class MovementData:
+    duration: float = 0.0
+    """Duration of the path, in seconds."""
+
+    moving_duration: float = 0.0
+    """Duration spent moving, in seconds."""
+
+    stopped_duration: float = 0.0
+    """Duration spent stopped, in seconds."""
+
+    length: float = 0.0
+    """Length of the path, in meters."""
+
+    average_speed: float = 0.0
+    """Average speed, in m/s."""
+
+    maximum_speed: float = 0.0
+    """Maximum speed, in m/s."""
+
+    uphill: float = 0.0
+    """Uphill elevation, in meters."""
+
+    downhill: float = 0.0
+    """Downhill elevation, in meters."""
+
+
+@dataclass(slots=True)
+class Point:
+    longitude: float
+    latitude: float
+    elevation: float
+    time_offset: float
+
+    def distance(self, other: "Point") -> float:
+        """Returns the distance between this point and the given other point in
+        meters.
+        """
+        r_1 = EARTH_RADIUS + self.elevation
+        r_2 = EARTH_RADIUS + other.elevation
+        # The formula assumes that 0° is straight upward, but 0° in geo
+        # coordinates is actually on the equator plane.
+        t_1 = radians(90 - self.latitude)
+        t_2 = radians(90 - other.latitude)
+        p_1 = radians(self.longitude)
+        p_2 = radians(other.longitude)
+        # See
+        # https://en.wikipedia.org/wiki/Spherical_coordinate_system#Distance_in_spherical_coordinates
+        # While this is not the Haversine formula for distances along the
+        # circle curvature, it allows us to take the elevation into account,
+        # and for most GPS point differences that we encounter it should be
+        # enough.
+        radicand = (
+            r_1**2 +
+            r_2**2 -
+            2 * r_1 * r_2 * (
+                sin(t_1) * sin(t_2) * cos(p_1 - p_2) +
+                cos(t_1) * cos(t_2)
+            )
+        )
+        if radicand < 0.0:
+            return 0.0
+        return sqrt(radicand)
+
+
+class Path:
+    @classmethod
+    def from_gpx(cls, gpx: gpxpy.gpx.GPX) -> "Path":
+        points = []
+        start_time = None
+
+        for track in gpx.tracks:
+            for segment in track.segments:
+                for point in segment.points:
+                    if start_time is None:
+                        start_time = point.time
+
+                    time_offset = (point.time - start_time).total_seconds()
+                    points.append(Point(
+                        longitude=point.longitude,
+                        latitude=point.latitude,
+                        elevation=point.elevation,
+                        time_offset=time_offset,
+                    ))
+
+        return cls(points)
+
+    def __init__(self, points: list[Point]):
+        self.points = points
+
+    def _point_pairs(self):
+        return zip(self.points, islice(self.points, 1, None))
+
+    def movement_data(self) -> MovementData:
+        """Returns the movement data."""
+        movement_data = MovementData()
+        for a, b in self._point_pairs():
+            distance = a.distance(b)
+            time = b.time_offset - a.time_offset
+            speed = distance / time
+            elevation = b.elevation - a.elevation
+
+            movement_data.length += distance
+            if speed >= MOVING_THRESHOLD:
+                movement_data.moving_duration += time
+            else:
+                movement_data.stopped_duration += time
+            movement_data.maximum_speed = max(movement_data.maximum_speed, speed)
+            if elevation > 0.0:
+                movement_data.uphill += elevation
+            else:
+                movement_data.downhill += -elevation
+            movement_data.duration = b.time_offset
+        movement_data.average_speed = movement_data.length / movement_data.moving_duration
+        return movement_data
-- 
cgit v1.2.3