implement Savitzky-Golay for elevation smoothing

4 files changed, 85 insertions, 4 deletions

diff --git a/fietsboek/transformers/__init__.py b/fietsboek/transformers/__init__.py
index 097fbaf..8f96af1 100644
--- a/fietsboek/transformers/__init__.py
+++ b/fietsboek/transformers/__init__.py
@@ -135,12 +135,13 @@ def list_transformers() -> list[type[Transformer]]:
     """
     # pylint: disable=import-outside-toplevel,cyclic-import
     from .breaks import RemoveBreaks
-    from .elevation import FixElevationJumps, FixNullElevation
+    from .elevation import FixElevationJumps, FixNullElevation, SmoothElevation
 
     return [
         FixNullElevation,
         FixElevationJumps,
         RemoveBreaks,
+        SmoothElevation,
     ]
 
 
diff --git a/fietsboek/transformers/elevation.py b/fietsboek/transformers/elevation.py
index 27683bb..c8874be 100644
--- a/fietsboek/transformers/elevation.py
+++ b/fietsboek/transformers/elevation.py
@@ -1,16 +1,25 @@
 """Transformers that deal with elevation changes in the track."""
 
+import logging
 from collections.abc import Callable, Iterable
 from itertools import islice, zip_longest
 
 from pyramid.i18n import TranslationString
+from SGCC import savgol  # type: ignore
 
 from .. import geo
 from . import Parameters, Transformer
 
 _ = TranslationString
 
+LOGGER = logging.getLogger(__name__)
 MAX_ORGANIC_SLOPE: float = 1.0
+ELEVATION_SMOOTHING_WINDOW_LENGTH: float = 500.0
+"""One-sided window length (in m) of the elevation smoothing filter.
+
+The filtering will estimate how many points need to be included to reach this
+length.
+"""
 
 
 def slope(point_a: geo.Point, point_b: geo.Point) -> float:
@@ -137,4 +146,62 @@ class FixElevationJumps(Transformer):
             current_point.elevation += point_adjustment
 
 
-__all__ = ["FixNullElevation", "FixElevationJumps"]
+class SmoothElevation(Transformer):
+    """A transformer that smooths the elevation."""
+
+    @classmethod
+    def identifier(cls) -> str:
+        return "smooth-elevation"
+
+    @classmethod
+    def name(cls) -> TranslationString:
+        return _("transformers.smooth-elevation")
+
+    @classmethod
+    def description(cls) -> TranslationString:
+        return _("transformers.smooth-elevation.description")
+
+    @classmethod
+    def parameter_model(cls) -> type[Parameters]:
+        return Parameters
+
+    @property
+    def parameters(self) -> Parameters:
+        return Parameters()
+
+    @parameters.setter
+    def parameters(self, value):
+        pass
+
+    def execute(self, path: geo.Path):
+        radius = self._estimate_window_size(path)
+        LOGGER.debug("Savitzky-Golay radius: %s", radius)
+        self._savitzky_golay(path, radius)
+
+    def _estimate_window_size(self, path: geo.Path) -> int:
+        length = 0.0
+        for i, (a, b) in enumerate(zip(path.points, path.points[1:])):
+            length += a.flat_distance(b)
+            if length >= ELEVATION_SMOOTHING_WINDOW_LENGTH:
+                return i
+        return 0
+
+    def _savitzky_golay(self, path: geo.Path, radius: int):
+        elevations = [0.0] * len(path.points)
+        coeffs = savgol.get_coefficients(0, 2, 2 * radius + 1, 0)
+
+        for i in range(0, radius + 1):
+            elevations[i] = path.points[i].elevation
+            elevations[-i - 1] = path.points[-i - 1].elevation
+
+        for i in range(radius, len(path.points) - radius - 1):
+            points = path.points[i - radius : i + radius + 1]
+            assert len(points) == len(coeffs)
+            smoothed = sum(coef * point.elevation for coef, point in zip(coeffs, points))
+            elevations[i] = smoothed
+
+        for i, elevation in enumerate(elevations):
+            path.points[i].elevation = elevation
+
+
+__all__ = ["FixNullElevation", "FixElevationJumps", "SmoothElevation"]
diff --git a/poetry.lock b/poetry.lock
index 3fc3a7a..251bc18 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1,4 +1,4 @@
-# This file is automatically @generated by Poetry 2.3.3 and should not be changed by hand.
+# This file is automatically @generated by Poetry 2.3.1 and should not be changed by hand.
 
 [[package]]
 name = "alabaster"
@@ -2706,6 +2706,18 @@ files = [
 ]
 
 [[package]]
+name = "savgol-calculator"
+version = "1.0.5"
+description = "Python implementation of a coefficients calculator for the Savitsky-Golay filter with any parameters."
+optional = false
+python-versions = ">=3.6.0"
+groups = ["main"]
+files = [
+    {file = "savgol-calculator-1.0.5.tar.gz", hash = "sha256:f0e49fdc1944f1aab3fda644d0b1d2fa9f6cd8f2e204599c0cf00c59234fbed3"},
+    {file = "savgol_calculator-1.0.5-py2.py3-none-any.whl", hash = "sha256:dd0db764ea9322b593f04e8bafbed383a205ae125ef4330f795e30aceaa69c7e"},
+]
+
+[[package]]
 name = "setuptools"
 version = "82.0.1"
 description = "Most extensible Python build backend with support for C/C++ extension modules"
@@ -3489,4 +3501,4 @@ test = ["zope.testing"]
 [metadata]
 lock-version = "2.1"
 python-versions = ">=3.11, <4"
-content-hash = "b4d90215cb83febb48f381f7eb9d8d2aa4fe159a34bfb93b96fbef54ea943881"
+content-hash = "1fd62aecf611f739cc79de3ea6a6f69665ca822960e9bde75954c8d255764ac3"
diff --git a/pyproject.toml b/pyproject.toml
index d3d22a1..e4c4347 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -46,6 +46,7 @@ dependencies = [
     "matplotlib (>=3.10.7,<4.0.0)",
 
     "hittekaart-py @ git+https://gitlab.com/dunj3/hittekaart@9a576a8d3a8cb1d143be0ade0310f71b490cacb9#subdirectory=hittekaart-py",
+    "savgol-calculator (==1.0.5)",
 ]
 
 [project.urls]