FixNullElevation: also take into account slope

For some reason, I have GPX tracks that have the first two points be
~100 meters apart in elevation, but only ~20 meters apart in distance.
This is quite unrealistic and produces pretty bad height plots (almost
as bad as the zero elevation).

Since the issue is very related, and the fix is pretty much the same, I
thought it would be a good idea to adapt the FixNullElevation
transformer to handle this case as well. For reference, "the internet"
says that the maximum slope for a MTB is ~15% to ~35%, depending on the
conditions - with 35% being pretty steep. I think it's fair to throw
away elevations that exceed 100% (basically a 45° angle upwards),
especially since we only discard them at the start and end.

1 files changed, 53 insertions, 13 deletions

diff --git a/fietsboek/transformers/__init__.py b/fietsboek/transformers/__init__.py
index f3afb57..f98a545 100644
--- a/fietsboek/transformers/__init__.py
+++ b/fietsboek/transformers/__init__.py
@@ -11,10 +11,10 @@ function to load and apply transformers.
 
 from abc import ABC, abstractmethod
 from collections.abc import Mapping
-from itertools import chain, islice
+from itertools import islice
 from typing import Literal, NamedTuple, TypeVar
 
-from gpxpy.gpx import GPX
+from gpxpy.gpx import GPX, GPXTrackPoint
 from pydantic import BaseModel
 from pyramid.i18n import TranslationString
 from pyramid.request import Request
@@ -164,17 +164,57 @@ class FixNullElevation(Transformer):
                 for point in segment.points
             )
 
-        points = all_points()
-        previous_points = chain([None], all_points())
-        next_points = chain(islice(all_points(), 1, None), [None])
-
-        for previous_point, point, next_point in zip(previous_points, points, next_points):
-            if point.elevation == 0.0:
-                if previous_point:
-                    point.elevation += previous_point.elevation
-                if next_point:
-                    point.elevation += next_point.elevation
-                point.elevation /= sum(1 for pt in [previous_point, next_point] if pt)
+        def rev_points():
+            return (
+                point
+                for track in reversed(gpx.tracks)
+                for segment in reversed(track.segments)
+                for point in reversed(segment.points)
+            )
+
+        max_slope = 1.0
+
+        # First, from the front, find the first point with non-zero elevation (or low enough slope)
+        bad_until = 0
+        final_elevation = 0
+        for i, (point, next_point) in enumerate(zip(all_points(), islice(all_points(), 1, None))):
+            if point.elevation != 0.0 and self.slope(point, next_point) < max_slope:
+                bad_until = i
+                final_elevation = point.elevation
+                break
+
+        for point in islice(all_points(), None, bad_until):
+            point.elevation = final_elevation
+
+        # Second, from the back
+        bad_until = 0
+        final_elevation = 0
+        for i, (point, prev_point) in enumerate(zip(rev_points(), islice(rev_points(), 1, None))):
+            if point.elevation != 0.0 and self.slope(point, prev_point) < max_slope:
+                bad_until = i
+                final_elevation = point.elevation
+                break
+
+        for point in islice(rev_points(), None, bad_until):
+            point.elevation = final_elevation
+
+    @staticmethod
+    def slope(point_a: GPXTrackPoint, point_b: GPXTrackPoint) -> float:
+        """Returns the slope between two GPX points.
+
+        This is defined as delta_h / euclid_distance.
+
+        :param point_a: First point.
+        :param point_b: Second point.
+        :return: The slope, as percentage.
+        """
+        if point_a.elevation is None or point_b.elevation is None:
+            return 0.0
+        delta_h = abs(point_a.elevation - point_b.elevation)
+        dist = point_a.distance_2d(point_b)
+        if dist == 0.0 or dist is None:
+            return 0.0
+        return delta_h / dist
 
 
 def list_transformers() -> list[type[Transformer]]: