Fix polygons intersections issues

Intersection of input polygons and tiles bounding box are not always polygons. Properly discard lines and points intersections.

Close #42

README.md

@@ -143,6 +143,15 @@ Here is an example originating from a "view1" source with 10m step (lon6.00_7.00
 
 This project uses ![hatch](https://hatch.pypa.io/latest/).
 
+(Mini)Conda can be used to easily setup given version of Python and GDAL:
+```bash
+conda create -n python39 -c conda-forge python=3.9
+conda activate python39
+conda install -c conda-forge gdal hatch
+hatch -e geotiff shell
+# To start VSCode using the proper python env
+code .
+```
 The one-liner for formatting and local validation:
 ```bash
 hatch run fmt && hatch run all

pyhgtmap/hgt/file.py

@@ -230,6 +230,39 @@ def calcHgtArea(
 BBOX_EXPAND_EPSILON = 0.1
 
 
+def clip_polygons(
+    polygons: List[List[Tuple[float, float]]],
+    clip_polygon: Iterable[Tuple[float, float]],
+) -> List[List[Tuple[float, float]]]:
+    """
+    Clips a list of polygons to a given clip polygon.
+
+    Args:
+        polygons: A list of polygons to be clipped.
+        clip_polygon: The clip polygon to be used for clipping.
+
+    Returns:
+        A list of clipped polygons.
+    """
+    bbox_shape = shapely.Polygon(clip_polygon)
+    clipped_polygons: List[List[Tuple[float, float]]] = []
+    for p in polygons:
+        # Intersect each input polygon with the clip one
+        clipped_p = shapely.intersection(shapely.Polygon(p), bbox_shape)
+        # Resulting intersection(s) might have several forms
+        if isinstance(clipped_p, (shapely.MultiPolygon, shapely.GeometryCollection)):
+            # Keep only polygons intersections
+            clipped_polygons += [
+                [(x, y) for x, y in poly.exterior.coords]
+                for poly in clipped_p.geoms
+                if isinstance(poly, shapely.Polygon) and not poly.is_empty
+            ]
+        elif isinstance(clipped_p, shapely.Polygon) and not clipped_p.is_empty:
+            clipped_polygons.append([(x, y) for x, y in clipped_p.exterior.coords])
+
+    return clipped_polygons
+
+
 def polygon_mask(
     x_data: numpy.ndarray,
     y_data: numpy.ndarray,
@@ -261,18 +294,8 @@ def polygon_mask(
     if transform is not None:
         xyPoints = transform(xyPoints)
         bbox_points = transform(bbox_points)
-    bbox_shape = shapely.Polygon(bbox_points)
-    clipped_polygons = []
-    for p in polygons:
-        clipped_p = shapely.intersection(shapely.Polygon(p), bbox_shape)
-        if isinstance(clipped_p, shapely.MultiPolygon):
-            clipped_polygons += [
-                [(x, y) for x, y in poly.exterior.coords]
-                for poly in clipped_p.geoms
-                if not poly.is_empty
-            ]
-        elif not clipped_p.is_empty:
-            clipped_polygons.append([(x, y) for x, y in clipped_p.exterior.coords])
+
+    clipped_polygons = clip_polygons(polygons, bbox_points)
 
     if not clipped_polygons:
         # Empty intersection: data is fully masked

tests/hgt/test_file.py

@@ -8,8 +8,7 @@
 import pytest
 
 from pyhgtmap import hgt
-from pyhgtmap.hgt.file import calcHgtArea, hgtFile, polygon_mask
-from pyhgtmap.hgt.tile import hgtTile
+from pyhgtmap.hgt.file import calcHgtArea, clip_polygons, hgtFile, hgtTile, polygon_mask
 
 from .. import TEST_DATA_PATH
 
@@ -289,3 +288,69 @@ def test_calcHgtArea(file_name: str) -> None:
             [(os.path.join(TEST_DATA_PATH, file_name), False)], 0, 0
         )
         assert bbox == (MIN_LON, MIN_LAT, MAX_LON, MAX_LAT)
+
+
+def test_clip_polygons() -> None:
+    clip_polygon: List[Tuple[float, float]] = [
+        (-0.1, 48.900000000009435),
+        (-0.1, 50.1),
+        (1.1, 50.1),
+        (1.1, 48.900000000009435),
+        (-0.1, 48.900000000009435),
+    ]
+    # Multi-polygons in input
+    polygons: List[List[Tuple[float, float]]] = [
+        # Real intersection is a polygon + a line; line must be discarded properly
+        [
+            (2.3, 51.6),
+            (2.5, 51.3),
+            (2.4, 50.9),
+            (1.3, 50.1),
+            (0.7, 50.1),
+            (0.4, 49.9),
+            (-0.5, 50.0),
+            (-0.9, 49.8),
+            (-2.2, 49.7),
+            (-2.9, 49.8),
+        ],
+        # No intersection
+        [
+            (-14.6, 57.6),
+            (-14.6, 57.9),
+            (-13.9, 58.4),
+            (-13.2, 58.3),
+            (-12.8, 57.9),
+            (-12.9, 57.1),
+            (-13.4, 56.8),
+            (-14.2, 56.9),
+            (-14.6, 57.3),
+            (-14.6, 57.6),
+        ],
+        # Single point intersection
+        [
+            (2, 52),
+            (2, 50.1),
+            (1.1, 50.1),
+            (1.1, 52),
+            (2, 52),
+        ],
+        # Single line intersection
+        [
+            (2, 48),
+            (2, 50),
+            (1.1, 50),
+            (1.1, 48),
+            (2, 48),
+        ],
+    ]
+
+    clipped_polygons = clip_polygons(polygons, clip_polygon)
+    assert clipped_polygons == [
+        [
+            (0.4, 49.9),
+            (-0.1, 49.955555555555556),
+            (-0.1, 50.1),
+            (0.7, 50.1),
+            (0.4, 49.9),
+        ]
+    ]