On invalid geometry make valid and clean

bench/bench_encode.py

@@ -2,8 +2,7 @@
 # -*- coding: utf-8 -*-
 
 import cProfile
-from itertools import islice
-from mapbox_vector_tile.encoder import VectorTile, on_invalid_geometry_ignore
+from mapbox_vector_tile.encoder import on_invalid_geometry_ignore
 from mapbox_vector_tile import encode
 from shapely.wkt import loads as loads_wkt
 from shapely.geometry import mapping
@@ -30,20 +29,21 @@ def make_layers(shapes, geom_dicts=False):
             i += 1
         except:
             pass
-    layers.append(features)
     return layers
 
+
 def run_test(layers):
     print("Running perf test")
     i = 0
     profiler = cProfile.Profile()
     for layer in layers:
         layer_description = {
-            'features' : layer,
+            'features': layer,
             'name': 'bar'
         }
         profiler.enable()
-        res = encode(layer_description, on_invalid_geometry=on_invalid_geometry_ignore, round_fn=round)
+        encode(layer_description,
+               on_invalid_geometry=on_invalid_geometry_ignore, round_fn=round)
         profiler.disable()
         if i % 100 == 0:
             print("{} tiles produced".format(i))

mapbox_vector_tile/__init__.py

@@ -16,10 +16,15 @@ def decode(tile, y_coord_down=False):
 
 
 def encode(layers, quantize_bounds=None, y_coord_down=False, extents=4096,
-           on_invalid_geometry=None, round_fn=None, check_winding_order=True):
-    vector_tile = encoder.VectorTile(extents, on_invalid_geometry,
-                                     round_fn=round_fn,
-                                     check_winding_order=check_winding_order)
+           on_invalid_geometry=None, round_fn=None,
+           check_winding_order=True,
+           max_geometry_validate_tries=5):
+    vector_tile = encoder.VectorTile(
+                     extents,
+                     on_invalid_geometry,
+                     round_fn=round_fn,
+                     check_winding_order=check_winding_order,
+                     max_geometry_validate_tries=max_geometry_validate_tries)
     if (isinstance(layers, list)):
         for layer in layers:
             vector_tile.addFeatures(layer['features'], layer['name'],

mapbox_vector_tile/encoder.py

@@ -1,4 +1,4 @@
-from mapbox_vector_tile.polygon import make_it_valid
+from mapbox_vector_tile.polygon import make_it_valid, clean_multi
 from numbers import Number
 from past.builtins import long
 from past.builtins import unicode
@@ -10,6 +10,7 @@
 from shapely.ops import transform
 from shapely.wkb import loads as load_wkb
 from shapely.wkt import loads as load_wkt
+from shapely.validation import explain_validity
 import decimal
 from .compat import PY3, vector_tile, apply_map
 from .geom_encoder import GeometryEncoder
@@ -28,6 +29,16 @@ def on_invalid_geometry_make_valid(shape):
     return make_it_valid(shape)
 
 
+def on_invalid_geometry_make_valid_and_clean(shape):
+    shape = make_it_valid(shape, asserted=False)
+    if not shape.is_valid and shape.type == 'MultiPolygon':
+        shape = clean_multi(shape)
+    assert shape.is_valid, \
+        "Not valid %s %s because %s" \
+        % (shape.type, shape.wkt, explain_validity(shape))
+    return shape
+
+
 class VectorTile:
     """
     """

mapbox_vector_tile/polygon.py

@@ -40,13 +40,14 @@ def _drop_degenerate_inners(shape):
     return Polygon(shape.exterior, new_inners)
 
 
-def _contour_to_poly(contour):
+def _contour_to_poly(contour, asserted=True):
     poly = Polygon(contour)
     if not poly.is_valid:
         poly = poly.buffer(0)
-    assert poly.is_valid, \
-        "Contour %r did not make valid polygon %s because %s" \
-        % (contour, poly.wkt, explain_validity(poly))
+    if asserted:
+        assert poly.is_valid, \
+            "Contour %r did not make valid polygon %s because %s" \
+            % (contour, poly.wkt, explain_validity(poly))
     return poly
 
 
@@ -84,7 +85,7 @@ def _generate_polys(contours):
         yield p
 
 
-def _polytree_node_to_shapely(node):
+def _polytree_node_to_shapely(node, asserted):
     """
     Recurses down a Clipper PolyTree, extracting the results as Shapely
     objects.
@@ -95,7 +96,7 @@ def _polytree_node_to_shapely(node):
     polygons = []
     children = []
     for ch in node.Childs:
-        p, c = _polytree_node_to_shapely(ch)
+        p, c = _polytree_node_to_shapely(ch, asserted)
         polygons.extend(p)
         children.extend(c)
 
@@ -110,7 +111,10 @@ def _polytree_node_to_shapely(node):
             children = []
 
     elif node.Contour:
-        poly = _contour_to_poly(node.Contour)
+        poly = _contour_to_poly(node.Contour, asserted)
+        for ch in children:
+            inner = _contour_to_poly(ch, asserted)
+            diff = poly.difference(inner)
 
         # we add each inner one-by-one so that we can reject them individually
         # if they cause the polygon to become invalid. if the shape has lots
@@ -152,7 +156,8 @@ def _polytree_node_to_shapely(node):
             if diff.is_valid:
                 poly = diff
 
-        assert poly.is_valid
+        if asserted:
+            assert poly.is_valid
         if poly.type == 'MultiPolygon':
             polygons.extend(poly.geoms)
         else:
@@ -169,19 +174,20 @@ def _polytree_node_to_shapely(node):
     return (polygons, children)
 
 
-def _polytree_to_shapely(tree):
-    polygons, children = _polytree_node_to_shapely(tree)
+def _polytree_to_shapely(tree, asserted):
+    polygons, children = _polytree_node_to_shapely(tree, asserted)
 
     # expect no left over children - should all be incorporated into polygons
     # by the time recursion returns to the root.
     assert len(children) == 0
 
     union = cascaded_union(polygons)
-    assert union.is_valid
+    if asserted:
+        assert union.is_valid
     return union
 
 
-def make_valid_pyclipper(shape):
+def make_valid_pyclipper(shape, asserted):
     """
     Use the pyclipper library to "union" a polygon on its own. This operation
     uses the even-odd rule to determine which points are in the interior of
@@ -208,10 +214,10 @@ def make_valid_pyclipper(shape):
     except pyclipper.ClipperException:
         return MultiPolygon([])
 
-    return _polytree_to_shapely(result)
+    return _polytree_to_shapely(result, asserted)
 
 
-def make_valid_polygon(shape):
+def make_valid_polygon(shape, asserted):
     """
     Make a polygon valid. Polygons can be invalid in many ways, such as
     self-intersection, self-touching and degeneracy. This process attempts to
@@ -229,19 +235,21 @@ def make_valid_polygon(shape):
 
     assert shape.geom_type == 'Polygon'
 
-    shape = make_valid_pyclipper(shape)
-    assert shape.is_valid
+    shape = make_valid_pyclipper(shape, asserted)
+
+    if asserted:
+        assert shape.is_valid
     return shape
 
 
-def make_valid_multipolygon(shape):
+def make_valid_multipolygon(shape, asserted):
     new_g = []
 
     for g in shape.geoms:
         if g.is_empty:
             continue
 
-        valid_g = make_valid_polygon(g)
+        valid_g = make_valid_polygon(g, asserted)
 
         if valid_g.type == 'MultiPolygon':
             new_g.extend(valid_g.geoms)
@@ -251,7 +259,7 @@ def make_valid_multipolygon(shape):
     return MultiPolygon(new_g)
 
 
-def make_it_valid(shape):
+def make_it_valid(shape, asserted=True):
     """
     Attempt to make any polygon or multipolygon valid.
     """
@@ -260,9 +268,44 @@ def make_it_valid(shape):
         return shape
 
     elif shape.type == 'MultiPolygon':
-        shape = make_valid_multipolygon(shape)
+        shape = make_valid_multipolygon(shape, asserted)
 
     elif shape.type == 'Polygon':
-        shape = make_valid_polygon(shape)
+        shape = make_valid_polygon(shape, asserted)
 
     return shape
+
+
+def clean_multi(shape):
+    """
+    Remove self- and ring-selfintersections from input Polygon geometries
+    """
+    polygons = []
+    exterior_lines = []
+    interior_lines = []
+    for p in shape:
+        exterior_lines = []
+        interior_lines = []
+        lnum = 0
+        boundary = p.boundary
+        if boundary.type == 'LineString':
+            if lnum == 0:
+                exterior_lines.append(boundary)
+        else:
+            for ls in boundary:
+                if lnum == 0:
+                    exterior_lines.append(ls)
+                else:
+                    interior_lines.append(ls)
+                lnum += 1
+    for el in exterior_lines:
+        if len(interior_lines) == 0:
+            polygons.append(Polygon(el).buffer(0))
+        else:
+            for il in interior_lines:
+                polygons.append(Polygon(el, Polygon(il).interiors).buffer(0))
+    poly = MultiPolygon(polygons)
+    assert poly.is_valid, \
+        "Not valid multipolygon %s because %s" \
+        % (poly.wkt, explain_validity(poly))
+    return poly

tests/test_encoder.py

@@ -913,6 +913,9 @@ def test_too_small_geometry(self):
                      on_invalid_geometry=on_invalid_geometry_make_valid)
         result = decode(pbf)
         features = result['foo']['features']
+        import sys
+        sys.stderr.write(str(features))
+        sys.stderr.write("\n")
         self.assertEqual(0, len(features))
 
 

tests/test_polygon.py

@@ -4,7 +4,7 @@
 """
 import unittest
 
-from mapbox_vector_tile.polygon import make_it_valid
+from mapbox_vector_tile.polygon import make_it_valid, clean_multi
 from shapely import wkt
 import os
 
@@ -174,6 +174,90 @@ def test_polygon_ring_of_inners_2(self):
         self.assertTrue(fixed.is_valid)
         self.assertEquals(22, fixed.area)
 
+    def test_clean_multi(self):
+        geom = wkt.loads("""MULTIPOLYGON (((1796 -189, 1794 -188, 1794 -190, 1798 -192,
+        1805 -194, 1806 -195, 1806 -194, 1807 -203, 1803 -206, 1798 -205,
+        1781 -194, 1772 -187, 1761 -185, 1748 -190, 1733 -181, 1725 -184,
+        1710 -198, 1704 -199, 1699 -195, 1695 -183, 1690 -175, 1683 -173,
+        1676 -173, 1670 -176, 1656 -190, 1642 -196, 1637 -205, 1633 -209,
+        1613 -215, 1603 -228, 1598 -233, 1564 -241, 1563 -239, 1580 -235,
+        1592 -234, 1598 -231, 1603 -227, 1611 -215, 1621 -210, 1633 -208,
+        1641 -195, 1657 -187, 1669 -174, 1675 -171, 1684 -171, 1691 -175,
+        1697 -184, 1700 -195, 1704 -197, 1710 -197, 1726 -181, 1732 -179,
+        1739 -182, 1749 -189, 1761 -184, 1773 -185, 1778 -188, 1783 -195,
+        1792 -198, 1800 -204, 1803 -204, 1804 -202, 1803 -196, 1789 -191,
+        1781 -183, 1779 -179, 1781 -174, 1786 -171, 1801 -174, 1806 -173,
+        1824 -181, 1829 -181, 1837 -174, 1838 -170, 1838 -166, 1830 -158,
+        1818 -155, 1816 -152, 1818 -149, 1822 -146, 1837 -142, 1853 -131,
+        1858 -132, 1869 -140, 1873 -140, 1875 -134, 1870 -121, 1871 -118,
+        1875 -115, 1885 -117, 1900 -117, 1915 -110, 1931 -113, 1946 -107,
+        1964 -118, 1982 -119, 1984 -122, 1987 -128, 1991 -130, 2001 -132,
+        2018 -131, 2023 -135, 2030 -144, 2037 -166, 2045 -169, 2055 -166,
+        2060 -166, 2064 -168, 2075 -182, 2095 -180, 2098 -183, 2101 -188,
+        2101 -193, 2100 -199, 2090 -209, 2085 -225, 2064 -261, 2059 -274,
+        2056 -290, 2056 -306, 2060 -321, 2066 -334, 2075 -345, 2084 -348,
+        2125 -349, 2149 -359, 2177 -389, 2193 -398, 2216 -404, 2226 -404,
+        2251 -407, 2269 -406, 2274 -407, 2277 -412, 2277 -417, 2275 -423,
+        2270 -429, 2253 -441, 2232 -443, 2224 -446, 2220 -451, 2219 -458,
+        2222 -463, 2234 -469, 2243 -471, 2250 -469, 2259 -463, 2266 -455,
+        2278 -436, 2279.176470588235 -435, 2280 -435, 2289 -427,
+        2285.5 -429.625, 2298 -419, 2302 -416, 2331 -404, 2335 -403,
+        2342 -403, 2350 -401, 2372 -384, 2381 -385, 2390.375 -390.625,
+        2393 -395, 2394 -403, 2392 -417, 2394 -423, 2398 -426, 2405 -427,
+        2414 -423, 2419 -415, 2425 -399, 2429 -392, 2435 -385, 2443 -383,
+        2450 -387, 2453 -393, 2450 -400, 2443 -410, 2442 -415, 2445 -421,
+        2453 -430, 2461 -435, 2474 -437, 2483 -435, 2486 -432, 2496 -411,
+        2497 -371, 2502 -360, 2510 -351, 2522 -340, 2531 -327, 2546 -310,
+        2552 -305, 2559 -302, 2586 -299, 2593 -296, 2582 -298, 2595 -295,
+        2601 -285, 2600 -280, 2597 -273, 2600 -271, 2603 -285, 2602 -290,
+        2597 -296, 2580 -309, 2579 -308, 2587 -302, 2574 -302, 2551 -311,
+        2538 -327, 2537 -326, 2533 -330, 2523 -345, 2513 -354, 2504 -366,
+        2500 -379, 2501 -400, 2499 -415, 2486 -435, 2480 -440, 2472 -442,
+        2459 -440, 2454 -437, 2442 -423, 2439 -416, 2439 -410, 2450 -394,
+        2449 -390, 2446 -388, 2440 -387, 2436 -388, 2430 -396, 2419 -422,
+        2412 -428, 2403 -430, 2395 -428, 2390 -423, 2389 -417, 2392 -402,
+        2390 -396, 2382 -390, 2374 -389, 2365 -392, 2351 -404, 2332 -408,
+        2314 -417, 2306 -424, 2297 -427, 2288 -432, 2278 -442, 2266 -460,
+        2252 -471, 2240 -473, 2220 -465, 2216 -459, 2216 -452, 2219 -446,
+        2227 -441, 2235 -439, 2249 -438, 2264 -429, 2274 -419, 2274 -413,
+        2272 -410, 2263 -409, 2232 -416, 2226 -416, 2184 -399, 2175 -393,
+        2162 -379, 2145 -365, 2129 -355, 2111 -352, 2085 -356, 2076 -352,
+        2069 -346, 2059 -332, 2052 -315, 2051 -302, 2053 -280, 2062 -254,
+        2059 -264, 2070 -246, 2071 -239, 2070 -238, 2071 -238, 2074 -234,
+        2087 -209, 2097 -199, 2099 -193, 2098 -187, 2094 -183, 2078 -186,
+        2073 -185, 2068 -181, 2061 -171, 2059 -169, 2042 -171, 2036 -170,
+        2031 -164, 2026 -144, 2020 -137, 2014 -134, 2000 -136, 1990 -133,
+        1986 -131, 1981 -122, 1963 -121, 1946 -109, 1942 -109, 1936 -115,
+        1932 -116, 1915 -113, 1911 -114, 1906 -118, 1900 -119, 1887 -119,
+        1880 -117, 1875 -117, 1872 -121, 1878 -136, 1876 -140, 1874 -143,
+        1869 -143, 1860 -136, 1856 -135, 1837 -146, 1821 -149, 1819 -151,
+        1819 -153, 1831 -156, 1839 -164, 1841 -170, 1839 -177, 1831 -182,
+        1826 -184, 1787 -173, 1783 -174, 1782 -177, 1785 -184, 1792 -186,
+        1796 -189),
+        (2226 -404, 2229 -406, 2231 -406, 2232 -406, 2226 -404),
+        (2563 -303, 2557 -305, 2554 -307, 2558 -307, 2563 -303),
+        (2462 -438, 2470 -441, 2478 -439, 2463 -438, 2462 -438),
+        (2336 -404, 2332 -404, 2327 -406, 2325 -407, 2336 -404),
+        (2299 -425, 2306 -423, 2316 -415, 2304 -418, 2299 -424,
+        2298 -425, 2299 -425),
+        (2206 -406, 2229 -414, 2248 -410, 2241 -409, 2206 -406),
+        (2240 -407, 2239 -408, 2242 -408, 2241 -407, 2240 -407),
+        (2204 -405, 2204 -406, 2205 -406, 2204 -405),
+        (2156 -370, 2157 -371, 2150 -362, 2144 -360, 2156 -370),
+        (2083 -353, 2089 -354, 2094 -352, 2079 -351, 2083 -353),
+        (2079 -349, 2078 -350, 2080 -350, 2080 -349, 2079 -349)),
+        ((1796 -189, 1804 -189, 1798 -190, 1796 -189)),
+        ((1804 -189, 1806 -191, 1806 -190, 1807 -191, 1806 -193, 1804 -189)),
+        ((1806 -193, 1806 -194, 1805 -194, 1806 -193)),
+        ((2390 -390, 2385 -385, 2379 -382, 2385 -384, 2390 -390)),
+        ((2062 -254, 2068 -239, 2069 -238, 2070 -238, 2062 -254)),
+        ((1786 -185, 1791 -188, 1793 -190, 1786 -185),
+        (1786 -185, 1787 -187, 1789 -189, 1793 -190, 1792 -187,
+        1786 -185)))""")
+        self.assertFalse(geom.is_valid)
+        fixed = clean_multi(geom)
+        self.assertTrue(fixed.is_valid)
+
     def test_polygon_inners_crossing_outer(self):
         geom = wkt.loads("""POLYGON (
           (2325 1015, 2329 1021, 2419 1057, 2461 944, 2369 907, 2325 1015),