Merge pull request #1708 from tilezen/zerebubuth/multiple-pass-road-m…

…erge

 Make multiple passes over road merging

integration-test/1227-improve-road-merging.py

@@ -4,7 +4,7 @@
 
 class MergeJunctionTest(FixtureTest):
 
-    def test_junction(self):
+    def test_junction_x(self):
         from tilequeue.tile import coord_to_bounds
         from shapely.geometry import LineString, asShape
         from ModestMaps.Core import Coordinate
@@ -40,17 +40,137 @@ def test_junction(self):
         )
 
         with self.features_in_tile_layer(z, x, y, 'roads') as features:
-            # should have merged into a single _feature_
-            self.assertTrue(len(features) == 1)
+            # multilinestrings which contain lines which cross (as in the X
+            # above) are "non-simple", and many geometry operations start by
+            # forcing multilinestrings to be simple. we don't want this, as
+            # it introduces an extra coordinate where the lines cross.
+            # instead, we split into features which are individually simple,
+            # which means we'll need 2 in this example.
+            self.assertTrue(len(features) == 2)
 
             # when the test suite runs in "download only mode", an empty
             # set of features is passed into this block. the assertion
             # is shorted out, so we need this additional check which is
             # trivially satisfied in the case we're doing real testing.
-            if len(features) == 1:
-                # the shape should be a multilinestring
-                shape = asShape(features[0]['geometry'])
-                self.assertTrue(shape.geom_type == 'MultiLineString')
+            if len(features) == 2:
+                for i in (0, 1):
+                    # the shapes should be single linestrings in this example.
+                    shape = asShape(features[i]['geometry'])
+                    self.assertTrue(shape.geom_type == 'LineString')
 
-                # with two internal linestrings
-                self.assertTrue(len(shape.geoms) == 2)
+                    # consisting of _only two_ points. (i.e: one didn't get
+                    # inserted into the middle)
+                    self.assertTrue(len(shape.coords) == 2)
+
+    def test_junction_hash(self):
+        from tilequeue.tile import coord_to_bounds
+        from shapely.geometry import LineString, asShape
+        from ModestMaps.Core import Coordinate
+        import dsl
+
+        z, x, y = (12, 2048, 2048)
+
+        minx, miny, maxx, maxy = coord_to_bounds(
+            Coordinate(zoom=z, column=x, row=y))
+        midl = minx + (maxx - minx) / 3
+        midr = minx + 2 * (maxx - minx) / 3
+        midd = miny + (maxy - miny) / 3
+        midu = miny + 2 * (maxy - miny) / 3
+
+        road_props = dict(
+            highway='residential',
+            source='openstreetmap.org',
+        )
+
+        # make a tile with 4 roads in a # shape, as below.
+        #
+        #      |     |
+        #      1     2
+        #      |     |
+        # --7--+--8--+--9--
+        #      |     |
+        #      3     4
+        #      |     |
+        # -10--+-11--+-12--
+        #      |     |
+        #      5     6
+        #      |     |
+        #
+        # these should get merged into two features, one with 1->3->5 and
+        # 2->4->6 and the other with 7->8->9 and 10->11->12.
+        self.generate_fixtures(
+            dsl.way(1, LineString([[midl, maxy], [midl, midu]]), road_props),
+            dsl.way(2, LineString([[midr, maxy], [midr, midu]]), road_props),
+            dsl.way(3, LineString([[midl, midu], [midl, midd]]), road_props),
+            dsl.way(4, LineString([[midr, midu], [midr, midd]]), road_props),
+            dsl.way(5, LineString([[midl, midd], [midl, miny]]), road_props),
+            dsl.way(6, LineString([[midr, midd], [midr, miny]]), road_props),
+
+            dsl.way(7, LineString([[minx, midu], [midl, midu]]), road_props),
+            dsl.way(8, LineString([[minx, midd], [midl, midd]]), road_props),
+            dsl.way(9, LineString([[midl, midu], [midr, midu]]), road_props),
+            dsl.way(10, LineString([[midl, midd], [midr, midd]]), road_props),
+            dsl.way(11, LineString([[midr, midu], [maxx, midu]]), road_props),
+            dsl.way(12, LineString([[midr, midd], [maxx, midd]]), road_props),
+        )
+
+        class ApproxCoordSet(object):
+            def __init__(self, coords, tolerance):
+                self.coords = coords
+                self.tolerance = tolerance
+
+            def check_and_remove(self, item):
+                x, y = item
+
+                for ex, ey in self.coords:
+                    if abs(x - ex) < self.tolerance and \
+                       abs(y - ey) < self.tolerance:
+                        self.coords.remove((ex, ey))
+                        return True
+
+                return False
+
+        # scale from the coordinates within the tile to the coords in the
+        # generated tile.
+        scale = 20026376.39 / 180.0
+        tolerance = scale * 1.0e-4
+
+        with self.features_in_tile_layer(z, x, y, 'roads') as features:
+            self.assertTrue(len(features) == 2,
+                            "expected 2 features, got %d" % (len(features),))
+
+            expected_coords = ApproxCoordSet([
+                (midl * scale, maxy * scale),
+                (midl * scale, miny * scale),
+                (midr * scale, maxy * scale),
+                (midr * scale, miny * scale),
+                (minx * scale, midu * scale),
+                (minx * scale, midd * scale),
+                (maxx * scale, midu * scale),
+                (maxx * scale, midd * scale),
+            ], tolerance)
+
+            if len(features) == 2:
+                for i in (0, 1):
+                    # the shapes should be multilinestrings with two lines.
+                    shape = asShape(features[i]['geometry'])
+                    self.assertTrue(shape.geom_type == 'MultiLineString')
+                    self.assertTrue(
+                        len(shape.geoms) == 2,
+                        "expected 2 geometries in the MultiLineString, but "
+                        "there are %d" % (len(shape.geoms),))
+
+                    for line_i in (0, 1):
+                        # each line should consist of only two points
+                        line = shape.geoms[line_i]
+                        self.assertTrue(
+                            len(line.coords) == 2,
+                            "expected 2 points, but line has %d" %
+                            (len(line.coords),))
+
+                        for coord_i in (0, 1):
+                            coord = line.coords[coord_i]
+                            self.assertTrue(
+                                expected_coords.check_and_remove(coord),
+                                "%r not in expected set %r" %
+                                (coord, expected_coords.coords))

queries.yaml

@@ -1145,6 +1145,12 @@ post_process:
       # roads at that point is less than 5 degrees.
       merge_junctions: true
       merge_junction_angle: 5.0
+      # setting the following will cause lines, or parts of multi-lines,
+      # shorter than 0.1px at nominal zoom to be dropped.
+      drop_short_segments: true
+      drop_length_pixels: 0.1
+      # integrated simplification step tolerance
+      simplify_tolerance: 1.0
 
   # we do want to merge at 15, but we don't want to merge junctions becase that
   # might merge across oneway information, which doesn't get dropped until
@@ -1154,30 +1160,12 @@ post_process:
       source_layer: roads
       start_zoom: 15
       end_zoom: 16
-
-  # simplify roads again, to take advantage of any opportunities opened up
-  # by merging roads with the same properties in the previous step.
-  - fn: vectordatasource.transform.simplify_layer
-    params:
-      source_layer: roads
-      start_zoom: 8
-      end_zoom: 16
-      tolerance: 1.0
-
-  # merge _again_, this time to merge any features which no longer intersect
-  # because of the simplification step and can be packed into the same
-  # MultiLineString. also, drop short segments within the MultiLineString
-  # which didn't get merged into something larger - at 0.1 pixels, these
-  # probably aren't visible anyway.
-  - fn: vectordatasource.transform.merge_line_features
-    params:
-      source_layer: roads
-      start_zoom: 8
-      end_zoom: 16
       # setting the following will cause lines, or parts of multi-lines,
       # shorter than 0.1px at nominal zoom to be dropped.
       drop_short_segments: true
       drop_length_pixels: 0.1
+      # integrated simplification step tolerance
+      simplify_tolerance: 1.0
 
   # NOTE: want to do this _before_ buildings_unify, as after that we might not
   # have a feature ID to match buildings on!

vectordatasource/transform.py

@@ -3686,25 +3686,103 @@ def _merge_junctions_in_multilinestring(mls, angle_tolerance):
         return MultiLineString(merged_geoms)
 
 
-def _merge_junctions(features, angle_tolerance):
+def _loop_merge_junctions(geom, angle_tolerance):
+    """
+    Keep applying junction merging to the MultiLineString until there are no
+    merge opportunities left.
+
+    A single merge step will only carry out one merge per LineString, which
+    means that the other endpoint might miss out on a possible merge. So we
+    loop over the merge until all opportunities are exhausted: either we end
+    up with a single LineString or we run a step and it fails to merge any
+    candidates.
+
+    For a total number of possible merges, N, we could potentially be left
+    with two thirds of these left over, depending on the order of the
+    candidates. This means we should need only O(log N) steps to merge them
+    all.
+    """
+
+    if geom.geom_type != 'MultiLineString':
+        return geom
+
+    # keep track of the number of linestrings in the multilinestring. we'll
+    # use that to figure out if we've merged as much as we possibly can.
+    mls_size = len(geom.geoms)
+
+    while True:
+        geom = _merge_junctions_in_multilinestring(geom, angle_tolerance)
+
+        # merged everything down to a single linestring
+        if geom.geom_type == 'LineString':
+            break
+
+        # made no progress
+        elif len(geom.geoms) == mls_size:
+            break
+
+        assert len(geom.geoms) < mls_size, \
+            "Number of geometries should stay the same or reduce after merge."
+
+        # otherwise, keep looping
+        mls_size = len(geom.geoms)
+
+    return geom
+
+
+def _simplify_line_collection(shape, tolerance):
+    """
+    Calling simplify on a MultiLineString doesn't always simplify if it would
+    make the MultiLineString non-simple.
+
+    However, we're trying to sort linestrings into nonoverlapping sets, and we
+    don't care whether they overlap at this point. However, we do want to make
+    sure that any colinear points in the individual LineStrings are removed.
+    """
+
+    if shape.geom_type == 'LineString':
+        shape = shape.simplify(tolerance)
+
+    elif shape.geom_type == 'MultiLineString':
+        new_geoms = []
+        for geom in shape.geoms:
+            new_geoms.append(geom.simplify(tolerance))
+        shape = MultiLineString(new_geoms)
+
+    return shape
+
+
+def _merge_junctions(features, angle_tolerance, simplify_tolerance):
     """
     Merge LineStrings within MultiLineStrings within features across junction
     boundaries where the lines appear to continue at the same angle.
 
+    If simplify_tolerance is provided, apply a simplification step. This can
+    help to remove colinear junction points left over from any merging.
+
+    Finally, group the lines into non-overlapping sets, each of which generates
+    a separate MultiLineString feature to ensure they're already simple and
+    further geometric operations won't re-introduce intersection points.
+
     Returns a new list of features.
     """
 
     new_features = []
     for shape, props, fid in features:
         if shape.geom_type == 'MultiLineString':
-            shape = _merge_junctions_in_multilinestring(
-                shape, angle_tolerance)
-            if shape.geom_type == 'MultiLineString':
-                disjoint_shapes = _linestring_nonoverlapping_partition(shape)
-                for disjoint_shape in disjoint_shapes:
-                    new_features.append((disjoint_shape, props, None))
-                continue
-        new_features.append((shape, props, fid))
+            shape = _loop_merge_junctions(shape, angle_tolerance)
+
+        if simplify_tolerance > 0.0:
+            shape = _simplify_line_collection(shape, simplify_tolerance)
+
+        if shape.geom_type == 'MultiLineString':
+            disjoint_shapes = _linestring_nonoverlapping_partition(shape)
+            for disjoint_shape in disjoint_shapes:
+                new_features.append((disjoint_shape, props, None))
+
+        else:
+            new_features.append((shape, props, fid))
+
     return new_features
 
 
@@ -3841,6 +3919,8 @@ def merge_line_features(ctx):
         'drop_short_segments', default=False, typ=bool)
     short_segment_factor = params.optional(
         'drop_length_pixels', default=0.1, typ=float)
+    simplify_tolerance = params.optional(
+        'simplify_tolerance', default=0.0, typ=float)
 
     assert source_layer, 'merge_line_features: missing source layer'
     layer = _find_layer(ctx.feature_layers, source_layer)
@@ -3862,7 +3942,7 @@ def merge_line_features(ctx):
 
     if merge_junctions:
         layer['features'] = _merge_junctions(
-            layer['features'], junction_angle_tolerance)
+            layer['features'], junction_angle_tolerance, simplify_tolerance)
 
     return layer