Merge pull request #27 from magnusuMET/deprecated_collections

Replace use of deprecated collections with get_paths

geojsoncontour/contour.py

@@ -5,36 +5,43 @@
 from matplotlib.colors import rgb2hex
 from geojson import Feature, LineString
 from geojson import Polygon, FeatureCollection
-from .utilities.multipoly import MP, keep_high_angle, set_contourf_properties,get_contourf_levels
+
+from .utilities.multipoly import multi_polygon, keep_high_angle, set_contourf_properties,get_contourf_levels
+from .utilities.vertices import get_vertices_from_path
 
 
 def contour_to_geojson(contour, geojson_filepath=None, min_angle_deg=None,
                        ndigits=5, unit='', stroke_width=1, geojson_properties=None, strdump=False,
                        serialize=True):
     """Transform matplotlib.contour to geojson."""
-    collections = contour.collections
-    contour_index = 0
     line_features = []
-    for collection in collections:
-        color = collection.get_edgecolor()
-        for path in collection.get_paths():
-            v = path.vertices
-            if len(v) < 3:
+    paths = contour.get_paths()
+    colors = contour.get_edgecolors()
+    levels = contour.levels
+    for contour_index, (path, color, level) in enumerate(zip(paths, colors, levels)):
+        for coordinates in get_vertices_from_path(path):
+            if len(coordinates) < 3:
+                continue
+            if np.all(np.equal(coordinates, coordinates[0])):
+                # Matplotlib sometimes emits empty paths which
+                # can be ignored
                 continue
-            coordinates = keep_high_angle(v, min_angle_deg) if min_angle_deg else v
-            coordinates = np.around(coordinates, ndigits) if ndigits is not None else coordinates
+            if min_angle_deg:
+                coordinates = keep_high_angle(coordinates, min_angle_deg)
+            if ndigits:
+                coordinates = np.around(coordinates, ndigits)
             line = LineString(coordinates.tolist())
             properties = {
                 "stroke-width": stroke_width,
-                "stroke": rgb2hex(color[0]),
-                "title": "%.2f" % contour.levels[contour_index] + ' ' + unit,
-                "level-value": float("%.6f" % contour.levels[contour_index]),
+                "stroke": rgb2hex(color),
+                "title": f"{level:.2f} {unit}",
+                "level-value": float(f"{level:.6f}"),
                 "level-index": contour_index
             }
             if geojson_properties:
                 properties.update(geojson_properties)
             line_features.append(Feature(geometry=line, properties=properties))
-        contour_index += 1
+
     feature_collection = FeatureCollection(line_features)
     return _render_feature_collection(feature_collection, geojson_filepath, strdump, serialize)
 
@@ -44,23 +51,21 @@ def contourf_to_geojson_overlap(contourf, geojson_filepath=None, min_angle_deg=N
                                 geojson_properties=None, strdump=False, serialize=True):
     """Transform matplotlib.contourf to geojson with overlapping filled contours."""
     polygon_features = []
-    contourf_idx = 0
     contourf_levels = get_contourf_levels(contourf.levels, contourf.extend)
-    for collection in contourf.collections:
-        color = collection.get_facecolor()
-        for path in collection.get_paths():
-            for coord in path.to_polygons():
-                if min_angle_deg:
-                    coord = keep_high_angle(coord, min_angle_deg)
-                coord = np.around(coord, ndigits) if ndigits else coord
-                polygon = Polygon(coordinates=[coord.tolist()])
-                fcolor = rgb2hex(color[0])
-                properties = set_contourf_properties(stroke_width, fcolor, fill_opacity, contourf_levels[contourf_idx], unit)
-                if geojson_properties:
-                    properties.update(geojson_properties)
-                feature = Feature(geometry=polygon, properties=properties)
-                polygon_features.append(feature)
-        contourf_idx += 1
+    contourf_colors = contourf.get_facecolor()
+    for path, level, color in zip(contourf.get_paths(), contourf_levels, contourf_colors):
+        for coord in get_vertices_from_path(path):
+            if min_angle_deg:
+                coord = keep_high_angle(coord, min_angle_deg)
+            if ndigits:
+                coord = np.around(coord, ndigits)
+            polygon = Polygon(coordinates=[coord.tolist()])
+            fcolor = rgb2hex(color)
+            properties = set_contourf_properties(stroke_width, fcolor, fill_opacity, level, unit)
+            if geojson_properties:
+                properties.update(geojson_properties)
+            feature = Feature(geometry=polygon, properties=properties)
+            polygon_features.append(feature)
     feature_collection = FeatureCollection(polygon_features)
     return _render_feature_collection(feature_collection, geojson_filepath, strdump, serialize)
 
@@ -78,20 +83,20 @@ def contourf_to_geojson(contourf, geojson_filepath=None, min_angle_deg=None,
         variable_opacity = False
     polygon_features = []
     contourf_levels = get_contourf_levels(contourf.levels, contourf.extend)
-    for coll, level in zip(contourf.collections, contourf_levels):
-        color = coll.get_facecolor()
-        muli = MP(coll, min_angle_deg, ndigits)
-        polygon = muli.mpoly()
-        fcolor = rgb2hex(color[0])
-        if polygon.coordinates:
-            properties = set_contourf_properties(stroke_width, fcolor, fill_opacity, level, unit)
-            if geojson_properties:
-                properties.update(geojson_properties)
-            feature = Feature(geometry=polygon, properties=properties)
-            polygon_features.append(feature)
-            # print(len(polygon.coordinates))
-            if variable_opacity:
-                fill_opacity += opacity_increment
+    contourf_colors = contourf.get_facecolor()
+    for path, level, color in zip(contourf.get_paths(), contourf_levels, contourf_colors):
+        polygon = multi_polygon(path, min_angle_deg, ndigits)
+        if not polygon.coordinates:
+            continue
+        fcolor = rgb2hex(color)
+        properties = set_contourf_properties(stroke_width, fcolor, fill_opacity, level, unit)
+        if geojson_properties:
+            properties.update(geojson_properties)
+        feature = Feature(geometry=polygon, properties=properties)
+        polygon_features.append(feature)
+        # print(len(polygon.coordinates))
+        if variable_opacity:
+            fill_opacity += opacity_increment
     feature_collection = FeatureCollection(polygon_features)
     return _render_feature_collection(feature_collection, geojson_filepath, strdump, serialize)
 

geojsoncontour/utilities/multipoly.py

@@ -1,32 +1,61 @@
 #!/usr/bin/python3.4
 # -*- encoding: utf-8 -*-
 """Helper module for transformation of matplotlib.contour(f) to GeoJSON."""
+import enum
+
 from geojson import MultiPolygon
 import numpy as np
 
+from .vertices import get_vertices_from_path
+
 
-class MP(object):
-    """Class for easy MultiPolygon generation.
+class Orientation(enum.IntEnum):
+    CW = enum.auto()
+    CCW = enum.auto()
+
 
-    This class converts a matplotlib PathCollection into a GeoJSON MultiPolygon.
+def orientation(vertices) -> Orientation:
+    """Determine orientation for a closed polygon
+
+    See https://en.wikipedia.org/wiki/Curve_orientation for an
+    explanation of this formula
     """
+    # Ignoring last index (which is duplicate of first)
+    vertices = vertices[:-1, :]
+    lowest_x_indices = np.flatnonzero(vertices[:, 0] == np.min(vertices[:, 0]))
+    lowest_y_indices = np.flatnonzero(vertices[lowest_x_indices, 1] == np.min(vertices[lowest_x_indices, 1]))
+    index = lowest_x_indices[lowest_y_indices[0]]
+
+    prev_index = index - 1
+    next_index = (index + 1) % vertices.shape[0]
+    xa, ya = vertices[prev_index, :]
+    xb, yb = vertices[index, :]
+    xc, yc = vertices[next_index, :]
+
+    detO = (xb - xa)*(yc - ya) - (xc - xa)*(yb - ya)
+    if detO > 0:
+        return Orientation.CCW
+    else:
+        return Orientation.CW
+
+
+def multi_polygon(path, min_angle_deg, ndigits):
+    polygons = []
+    for linestring in path.to_polygons():
+        if min_angle_deg:
+            linestring = keep_high_angle(linestring, min_angle_deg)
+        if ndigits:
+            linestring = np.around(linestring, ndigits)
+
+        handedness = orientation(linestring)
+        if handedness == Orientation.CCW:
+            polygons.append([linestring.tolist()])
+        else:
+            # This is a hole, which we assume belong
+            # to the previous polygon
+            polygons[-1].extend([linestring.tolist()])
 
-    def __init__(self, path_collection, min_angle_deg, ndigits):
-        self.coords = []
-        for path in path_collection.get_paths():
-            polygon = []
-            for linestring in path.to_polygons():
-                if min_angle_deg:
-                    linestring = keep_high_angle(linestring, min_angle_deg)
-                if ndigits:
-                    linestring = np.around(linestring, ndigits)
-                polygon.append(linestring.tolist())
-            if polygon:
-                self.coords.append(polygon)
-
-    def mpoly(self):
-        """Output of GeoJSON MultiPolygon object."""
-        return MultiPolygon(coordinates=self.coords)
+    return MultiPolygon(coordinates=polygons)
 
 
 def unit_vector(vector):

geojsoncontour/utilities/vertices.py

@@ -0,0 +1,39 @@
+from collections.abc import Generator
+
+import numpy as np
+import numpy.typing as npt
+from matplotlib.path import Path
+
+def get_vertices_from_path(path: Path) -> Generator[npt.NDArray[float], None, None]:
+    """Splits the vertices from path into continous lines,
+    by taking into account path.codes
+
+    See https://matplotlib.org/stable/api/path_api.html for a
+    description of path.vertices/path.codes
+    """
+    # path = path.cleaned(curves=False, simplify=False)
+    vertices = path.vertices
+    codes = path.codes
+    if codes is None:
+        codes = [Path.MOVETO] + [Path.LINETO]*(len(vertices)-1)
+    current_vertice = []
+    for (v, c) in zip(vertices, codes):
+        if c == Path.STOP:
+            if len(current_vertice) != 0:
+                yield np.array(current_vertice)
+                current_vertice = []
+        elif c == Path.MOVETO:
+            if len(current_vertice) != 0:
+                yield np.array(current_vertice)
+            current_vertice = [v]
+        elif c == Path.LINETO:
+            current_vertice.append(v)
+        elif c == Path.CLOSEPOLY:
+            if len(current_vertice) != 0:
+                current_vertice.append(current_vertice[0])
+                yield np.array(current_vertice)
+                current_vertice = []
+        else:
+            raise Exception(f"Unknown code {c} encountered")
+    if len(current_vertice) != 0:
+        yield np.array(current_vertice)

setup.py

@@ -20,7 +20,7 @@
     install_requires=[
         'geojson',
         'numpy',
-        'matplotlib',
+        'matplotlib>=3.8',
         'xarray'
     ],
     zip_safe=False,