Add workaround for point segments on extracting the contours

shakyground2/shakemap.py

@@ -9,7 +9,6 @@ import rasterio
 import geopandas as gpd
 import matplotlib.pyplot as plt
 from typing import Dict, Optional, Tuple, List, Union
-from shapely.geometry import LineString, MultiLineString
 from rasterio.warp import reproject
 from openquake.hazardlib import const, imt
 from openquake.hazardlib.contexts import ContextMaker
@@ -17,6 +16,7 @@ from shakyground2 import valid
 from shakyground2.earthquake import Earthquake
 from shakyground2.site_model import SiteModel
 from shakyground2.synthetic_rupture_generator import FiniteRuptureSampler
+from shakyground2.utils import plt_contours_to_gpd_dataframe
 
 
 class Shakemap(object):
@@ -530,31 +530,15 @@ class Shakemap(object):
         ax = fig.add_subplot(111)
         # Get the contours
         contours = ax.contour(lons, lats, results, levels)
-        contour_levels = []
-        contour_geometries = []
-        # Retrieve the x, y coordinates of each contour
-        for level, segs in zip(contours.levels, contours.allsegs):
-            if not len(segs):
-                continue
-            else:
-                geometries = []
-                for seg in segs:
-                    geometries.append(LineString([(row[0], row[1]) for row in seg]))
-                if len(geometries) > 1:
-                    # Has multiple contour lines for this level
-                    contour_geometries.append(MultiLineString(geometries))
-                else:
-                    # Only a single contour for this level
-                    contour_geometries.append(geometries[0])
-                contour_levels.append(level)
+        dframe = plt_contours_to_gpd_dataframe(
+            contours,
+            name=imt,
+            crs={"init": "epsg:4326"},
+        )
         plt.close(fig)
         # Re-allow plotting
         plt.ion()
         # Create geodataframe
-        dframe = gpd.GeoDataFrame(
-            {imt: contour_levels, "geometry": gpd.GeoSeries(contour_geometries)}
-        )
-        dframe.crs = {"init": "epsg:4326"}
         return dframe
 
     def contours_to_file(

shakyground2/utils.py

+"""
+Utility classes & functions.
+"""
+from typing import Dict, Optional
+
+from geopandas import GeoDataFrame, GeoSeries
+from matplotlib.contour import QuadContourSet
+from shapely.geometry import LineString, MultiLineString
+
+
+def plt_contours_to_gpd_dataframe(
+    contours: QuadContourSet,
+    name: str,
+    crs: Optional[Dict[str, str]] = None,
+    epsilon: float = 0.000_000_001,
+) -> GeoDataFrame:
+    """
+    Transform the contours from pyplot to a geopandas dataframe.
+
+    Args:
+        contours: The contours object from plt.contours(x, y, z, level).
+        name: The name that the level series in the dataframe should get.
+        crs: The coordinate reference system for the result
+             ({"init": "egsg:4326"} for example).
+        epsilon: A difference value that can be added/substracted
+                 to a point location to convert them to linestring.
+    """
+    contour_levels = []
+    contour_geometries = []
+    # Retrieve the x, y coordinates of each contour
+    for level, segs in zip(contours.levels, contours.allsegs):
+        if not len(segs):
+            continue
+        else:
+            geometries = []
+            for seg in segs:
+                # Here is the problem
+                if len(seg) > 1:
+                    geometries.append(LineString([(row[0], row[1]) for row in seg]))
+                else:
+                    # We have a point wiht coordinates like [(1.0, 1.0)]
+                    # We can't use LineString directly as it needs at
+                    # least 2 points.
+                    # We can't add it as a Point, as the overall result
+                    # for the level should be a MultiLineString.
+                    # Just adding a LineString with two times the very
+                    # same point doesn't give us a valid geometry
+                    # (so that the MultiLineString isn't valid either).
+                    # So the idea here is to make it a small box
+                    # around that single point.
+                    # We discussed adding or substracting the smallest
+                    # value that python supports (sys.float_info.min),
+                    # but then it could happen that python handles
+                    # x + sys.float_info.min the very same way as x,
+                    # as it is not able to handle this level of precision.
+                    # For shakyground2 here, we are operate
+                    # on lat & lon values (so -180 to 180 degree maximum).
+                    # So we think adding or removing 0.000_000_001 should be ok.
+                    # For other cases the user can give another value
+                    # in the epsilon parameter.
+                    row = seg[0]
+                    x = row[0]
+                    y = row[1]
+                    x_minimal_smaller = x - epsilon
+                    x_minimal_larger = x + epsilon
+                    y_minimal_smaller = y - epsilon
+                    y_mininal_larger = y + epsilon
+
+                    points = [
+                        (x_minimal_smaller, y_minimal_smaller),
+                        (x_minimal_larger, y_minimal_smaller),
+                        (x_minimal_larger, y_mininal_larger),
+                        (x_minimal_smaller, y_mininal_larger),
+                    ]
+                    geometries.append(LineString(points))
+
+            if len(geometries) > 1:
+                # Has multiple contour lines for this level
+                contour_geometries.append(MultiLineString(geometries))
+            else:
+                # Only a single contour for this level
+                contour_geometries.append(geometries[0])
+            contour_levels.append(level)
+
+    dframe = GeoDataFrame(
+        {name: contour_levels, "geometry": GeoSeries(contour_geometries)}
+    )
+    if crs:
+        dframe.crs = crs
+    return dframe

tests/utils_test.py

+"""
+Utility function test cases.
+"""
+
+import unittest
+
+import matplotlib.pyplot as plt
+import numpy
+from shakyground2.utils import plt_contours_to_gpd_dataframe
+
+
+class PltContoursToGpdDataframeTestCase(unittest.TestCase):
+    """Test case for the conversion of plt contour object to gpd dataframe."""
+
+    def test_conversion(self):
+        """Test the conversion."""
+        # Create some test data
+        x = numpy.arange(1, 10)
+        y = x.reshape(-1, 1)
+        h = x * y
+
+        h[5, 5] = -12
+        levels = [0, 1, 5, 10, 81]
+        plt.ioff()
+        fig = plt.figure()
+        ax = fig.add_subplot(111)
+        contours = ax.contour(x, y.T[0], h, levels)
+        # First important point: Don't raise any exception here.
+        dframe = plt_contours_to_gpd_dataframe(
+            contours, name="h", crs={"init": "epsg:4326"}
+        )
+        plt.close()
+        plt.ion()
+        # And also make sure the geometry is valid
+        self.assertTrue(all(dframe.geometry.apply(lambda x: x.is_valid)))