Merge branch 'master' of ucyo/geojsoncontour into master

.gitignore

@@ -3,6 +3,9 @@ __pycache__/
 *.py[cod]
 *$py.class
 
+# C extensions
+*.so
+
 # Distribution / packaging
 .Python
 env/
@@ -21,4 +24,75 @@ var/
 .installed.cfg
 *.egg
 
-.idea
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+.hypothesis/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# IPython Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
+
+# celery beat schedule file
+celerybeat-schedule
+
+# dotenv
+.env
+
+# virtualenv
+venv/
+ENV/
+
+# Spyder project settings
+.spyderproject
+
+# Rope project settings
+.ropeproject
+
+# Macintosh
+.DS_Store
+
+
+# GeoJSON
+*.geojson
+
+.idea/

geojsoncontour/__init__.py

@@ -1 +1,2 @@
-from .contour import contour_to_geojson
+from .contour import contour_to_geojson, contourf_to_geojson
+from .contour import contourf_to_multipolygeojson, to_geojson

geojsoncontour/contour.py

@@ -1,32 +1,22 @@
-import math
+#!/usr/bin/python3.4
+# -*- encoding: utf-8 -*-
+"""Transform matplotlib.contour(f) to GeoJSON."""
 
-from matplotlib.colors import rgb2hex
-
-from geojson import Feature, LineString, FeatureCollection
 import geojson
+import numpy as np
+from matplotlib.colors import rgb2hex
+from geojson import Feature, LineString
+from geojson import Polygon, FeatureCollection
+from .utilities.multipoly import MP, keep_high_angle, set_properties
+import matplotlib
 
 
-def dotproduct(v1, v2):
-    return sum((a * b) for a, b in zip(v1, v2))
-
-
-def length(v):
-    return math.sqrt(dotproduct(v, v))
-
-
-def angle(v1, v2):
-    cos_angle = dotproduct(v1, v2) / (length(v1) * length(v2))
-    cos_angle = min(1.0, max(cos_angle, -1.0))
-    assert cos_angle <= 1.0
-    assert cos_angle >= -1.0
-    return math.acos(cos_angle)
-
-
-def contour_to_geojson(contour, geojson_filepath, contour_levels, min_angle_deg=2,
-                       ndigits=5, unit='', stroke_width=5, geojson_properties=None):
+def contour_to_geojson(contour, geojson_filepath=None, contour_levels=None, min_angle_deg=None,
+                       ndigits=5, unit='', stroke_width=1, geojson_properties=None, strdump=False):
+    """Transform matplotlib.contour to geojson."""
+    if contour_levels is None:
+        contour_levels = contour.levels
     collections = contour.collections
-    total_points = 0
-    total_points_original = 0
     contour_index = 0
     assert len(contour_levels) == len(collections)
     line_features = []
@@ -37,25 +27,10 @@ def contour_to_geojson(contour, geojson_filepath, contour_levels, min_angle_deg=
             v = path.vertices
             if len(v) < 3:
                 continue
-            coordinates = []
-            v1 = v[1] - v[0]
-            lat = round(v[0][0], ndigits)
-            long = round(v[0][1], ndigits)
-            coordinates.append((lat, long))
-            for i in range(1, len(v) - 2):
-                v2 = v[i + 1] - v[i - 1]
-                diff_angle = math.fabs(angle(v1, v2) * 180.0 / math.pi)
-                if diff_angle > min_angle_deg:
-                    lat = round(v[i][0], ndigits)
-                    long = round(v[i][1], ndigits)
-                    coordinates.append((lat, long))
-                    v1 = v[i] - v[i - 1]
-            lat = round(v[-1][0], ndigits)
-            long = round(v[-1][1], ndigits)
-            coordinates.append((lat, long))
-            total_points += len(coordinates)
-            total_points_original += len(v)
-            line = LineString(coordinates)
+            coordinates = keep_high_angle(v, min_angle_deg)
+            if ndigits:
+                coordinates = np.around(coordinates, ndigits)
+            line = LineString(coordinates.tolist())
             properties = {
                 "stroke-width": stroke_width,
                 "stroke": rgb2hex(color[0]),
@@ -67,8 +42,85 @@ def contour_to_geojson(contour, geojson_filepath, contour_levels, min_angle_deg=
                 properties.update(geojson_properties)
             line_features.append(Feature(geometry=line, properties=properties))
         contour_index += 1
-
     feature_collection = FeatureCollection(line_features)
-    dump = geojson.dumps(feature_collection, sort_keys=True)
+    if strdump or not geojson_filepath:
+        return geojson.dumps(feature_collection, sort_keys=True, separators=(',', ':'))
     with open(geojson_filepath, 'w') as fileout:
-        fileout.write(dump)
+        geojson.dump(feature_collection, fileout, sort_keys=True, separators=(',', ':'))
+
+
+def contourf_to_geojson(contourf, geojson_filepath=None, contour_levels=None, min_angle_deg=None,
+                        ndigits=5, unit='', stroke_width=1, fill_opacity=.9, strdump=False):
+    """Transform matplotlib.contourf to geojson."""
+    if contour_levels is None:
+        contour_levels = contourf.levels
+    polygon_features = []
+    contourf_idx = 0
+    for coll in contourf.collections:
+        color = coll.get_facecolor()
+        for path in coll.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_properties(stroke_width, fcolor, fill_opacity, contour_levels, contourf_idx, unit)
+                feature = Feature(geometry=polygon, properties=properties)
+                polygon_features.append(feature)
+        contourf_idx += 1
+    collection = FeatureCollection(polygon_features)
+    if strdump or not geojson_filepath:
+        return geojson.dumps(collection, sort_keys=True, separators=(',', ':'))
+    with open(geojson_filepath, 'w') as fileout:
+        geojson.dump(collection, fileout, sort_keys=True, separators=(',', ':'))
+
+
+def contourf_to_multipolygeojson(contourf, geojson_filepath=None, contour_levels=None, min_angle_deg=None,
+                                 ndigits=5, unit='', stroke_width=1, fill_opacity=.9, strdump=False):
+    """Transform matplotlib.contourf to geojson with MultiPolygons."""
+    if contour_levels is None:
+        contour_levels = contourf.levels
+    polygon_features = []
+    mps = []
+    contourf_idx = 0
+    for coll in contourf.collections:
+        color = coll.get_facecolor()
+        for path in coll.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
+                op = MP(contour_levels[contourf_idx], rgb2hex(color[0]))
+                if op in mps:
+                    for i, k in enumerate(mps):
+                        if k == op:
+                            mps[i].add_coords(coord.tolist())
+                else:
+                    op.add_coords(coord.tolist())
+                    mps.append(op)
+        contourf_idx += 1
+    # starting here the multipolys will be extracted
+    for muli in mps:
+        polygon = muli.mpoly()
+        fcolor = muli.color
+        properties = set_properties(stroke_width, fcolor, fill_opacity, contour_levels, contourf_idx, unit)
+        feature = Feature(geometry=polygon, properties=properties)
+        polygon_features.append(feature)
+    collection = FeatureCollection(polygon_features)
+    if strdump or not geojson_filepath:
+        return geojson.dumps(collection, sort_keys=True, separators=(',', ':'))
+    with open(geojson_filepath, 'w') as fileout:
+        geojson.dump(collection, fileout, sort_keys=True, separators=(',', ':'))
+
+
+def to_geojson(contour, multipolys=False, strdump=False, *args, **kwargs):
+    message = 'Expected QuadContourSet, got {}'.format(type(contour))
+    assert isinstance(contour, matplotlib.contour.QuadContourSet), message
+    unit = contour.da_unit
+    if not contour.filled:
+        return contour_to_geojson(contour, strdump=strdump, unit=unit, *args, **kwargs)
+    elif multipolys:
+        return contourf_to_multipolygeojson(contour, strdump=strdump, unit=unit, *args, **kwargs)
+    else:
+        return contourf_to_geojson(contour, unit=unit, strdump=strdump, *args, **kwargs)

geojsoncontour/utilities/multipoly.py

@@ -0,0 +1,71 @@
+#!/usr/bin/python3.4
+# -*- encoding: utf-8 -*-
+"""Helper module for transformation of matplotlib.contour(f) to GeoJSON."""
+from geojson import MultiPolygon
+import numpy as np
+
+
+class MP(object):
+    """Class for easy MultiPolygon generation.
+
+    Just a helper class for easy identification of similar matplotlib.collections.
+    """
+
+    def __init__(self, title, color):
+        """Distinction based on title and color."""
+        self.title = title
+        self.color = color
+        self.coords = []
+
+    def add_coords(self, coords):
+        """Add new coordinate set for MultiPolygon."""
+        self.coords.append(coords)
+
+    def __eq__(self, other):
+        """Comparison of two MP instances."""
+        return (self.title == getattr(other, 'title', False) and
+                self.color == getattr(other, 'color', False))
+
+    def mpoly(self):
+        """Output of GeoJSON MultiPolygon object."""
+        return MultiPolygon(coordinates=[self.coords])
+
+
+def unit_vector(vector):
+    """Return the unit vector of the vector."""
+    return vector / np.linalg.norm(vector)
+
+
+def angle(v1, v2):
+    """Return the angle in radians between vectors 'v1' and 'v2'."""
+    v1_u = unit_vector(v1)
+    v2_u = unit_vector(v2)
+    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
+
+
+def keep_high_angle(vertices, min_angle_deg):
+    """Keep vertices with angles higher then given minimum."""
+    accepted = []
+    v = vertices
+    v1 = v[1] - v[0]
+    accepted.append((v[0][0], v[0][1]))
+    for i in range(1, len(v) - 2):
+        v2 = v[i + 1] - v[i - 1]
+        diff_angle = np.fabs(angle(v1, v2) * 180.0 / np.pi)
+        if diff_angle > min_angle_deg:
+            accepted.append((v[i][0], v[i][1]))
+            v1 = v[i] - v[i - 1]
+    accepted.append((v[-1][0], v[-1][1]))
+    return np.array(accepted, dtype=vertices.dtype)
+
+
+def set_properties(stroke_width, fcolor, fill_opacity, contour_levels, contourf_idx, unit):
+    """Set property values for Polygon."""
+    return {
+        "stroke": fcolor,
+        "stroke-width": stroke_width,
+        "stroke-opacity": 1,
+        "fill": fcolor,
+        "fill-opacity": fill_opacity,
+        "title": "%.2f" % contour_levels[contourf_idx] + ' ' + unit
+    }