diff --git a/integration-test/1227-improve-road-merging.py b/integration-test/1227-improve-road-merging.py index 35230bea9..cdffa6516 100644 --- a/integration-test/1227-improve-road-merging.py +++ b/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)) diff --git a/queries.yaml b/queries.yaml index d77cf12f8..5566fc0c5 100644 --- a/queries.yaml +++ b/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! diff --git a/vectordatasource/transform.py b/vectordatasource/transform.py index d7c692e1a..bba6dca22 100644 --- a/vectordatasource/transform.py +++ b/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