Add experimental "distance to alpine ridge" feature (#23)

* Add new experimental feature

* Fix tests

mlpp_features/experimental.py

@@ -0,0 +1,63 @@
+"""This module implements experimental features"""
+from typing import List, Tuple
+
+import numpy as np
+from pyproj import CRS, Transformer
+
+
+def reproject_points(
+    latlon_wgs84: List[Tuple[float, float]], dst_epsg: str
+) -> List[Tuple[float, float]]:
+    transformer = Transformer.from_crs(CRS("epsg:4326"), CRS(dst_epsg), always_xy=True)
+    lon_src = [p[1] for p in latlon_wgs84]
+    lat_src = [p[0] for p in latlon_wgs84]
+    x_dst, y_dst = transformer.transform(lon_src, lat_src)
+    return [(x, y) for x, y in zip(x_dst, y_dst)]
+
+
+def distance_point_to_segment(
+    point: Tuple[float, float],
+    segment_start: Tuple[float, float],
+    segment_end: Tuple[float, float],
+) -> float:
+    """Calculate the distance from a point to a line segment."""
+    point = np.array(point)
+    segment_start = np.array(segment_start)
+    segment_end = np.array(segment_end)
+
+    # Vector from start to point and start to end
+    start_to_point = point - segment_start
+    start_to_end = segment_end - segment_start
+
+    # Project start_to_point onto start_to_end
+    projection = np.dot(start_to_point, start_to_end) / np.dot(
+        start_to_end, start_to_end
+    )
+
+    # Check if projection is in the segment
+    if projection < 0:
+        closest_point = segment_start
+    elif projection > 1:
+        closest_point = segment_end
+    else:
+        closest_point = segment_start + projection * start_to_end
+
+    # Compute the distance from the point to the closest point
+    return np.linalg.norm(point - closest_point)
+
+
+def distances_points_to_line(
+    points: List[Tuple[float, float]], line_points: List[Tuple[float, float]]
+) -> List[float]:
+    """For a list of points, find the minimum distance a polyline."""
+    min_distances = []
+    for point in points:
+        min_distance = float("inf")
+        for i in range(len(line_points) - 1):
+            segment_start = line_points[i]
+            segment_end = line_points[i + 1]
+            distance = distance_point_to_segment(point, segment_start, segment_end)
+            if distance < min_distance:
+                min_distance = distance
+        min_distances.append(min_distance)
+    return min_distances

mlpp_features/terrain.py

@@ -5,6 +5,7 @@
 import xarray as xr
 
 from mlpp_features.decorators import out_format
+from mlpp_features import experimental as exp
 
 
 LOGGER = logging.getLogger(__name__)
@@ -100,6 +101,55 @@ def cos_valley_index_10000m(
     return cos_valley.preproc.interp(stations).astype("float32")
 
 
+@out_format()
+def distance_to_alpine_ridge(
+    data: Dict[str, xr.Dataset], stations, *args, **kwargs
+) -> xr.Dataset:
+    """
+    Compute horizontal distance to the main Alpine ridge
+
+    **Experimental feature, use with caution!**
+    """
+    if all([len(ds) == 0 for ds in data.values()]):
+        raise KeyError()
+    alpine_crest_wgs84 = [
+        [45.67975, 6.88306],
+        [45.75149, 6.80643],
+        [45.88912, 7.07724],
+        [45.86909, 7.17029],
+        [46.25074, 8.03064],
+        [46.47280, 8.38946],
+        [46.55972, 8.55968],
+        [46.56318, 8.80080],
+        [46.61256, 8.96059],
+        [46.49712, 9.17104],
+        [46.50524, 9.33031],
+        [46.39905, 9.69325],
+        [46.40885, 10.01963],
+        [46.63982, 10.29218],
+        [46.83630, 10.50783],
+        [46.90567, 11.09742],
+    ]
+    points = [(lat, lon) for lat, lon in zip(stations.latitude, stations.longitude)]
+    points_proj = exp.reproject_points(points, "epsg:2056")
+    line_proj = exp.reproject_points(alpine_crest_wgs84, "epsg:2056")
+    distances = exp.distances_points_to_line(points_proj, line_proj)
+    return xr.Dataset(
+        coords={
+            "station": stations.index,
+            "longitude": ("station", stations.longitude),
+            "latitude": ("station", stations.latitude),
+            "elevation": ("station", stations.elevation),
+        },
+        data_vars={
+            "distance_to_alpine_ridge": (
+                "station",
+                distances,
+            ),
+        },
+    ).astype("float32")
+
+
 @out_format()
 def elevation_50m(data: Dict[str, xr.Dataset], stations, *args, **kwargs) -> xr.Dataset:
     """