From 71aeebabf2d31fdd6f15a12f844baa72ee7de2df Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Wed, 20 Sep 2017 12:04:42 +0100
Subject: [PATCH 01/19] Refactor fixture DataFetcher in anticipation of
 re-using this code for RAWR DataFetcher.

---
 tilequeue/query/fixture.py | 175 ++++++++++++++++++++-----------------
 1 file changed, 94 insertions(+), 81 deletions(-)

diff --git a/tilequeue/query/fixture.py b/tilequeue/query/fixture.py
index 14f9fc41..bf73e305 100644
--- a/tilequeue/query/fixture.py
+++ b/tilequeue/query/fixture.py
@@ -278,6 +278,92 @@ def mz_calculate_transit_routes_and_score(rows, rels, node_id, way_id):
                    railways=railways, trams=trams)
 
 
+# properties for a feature (fid, shape, props) in layer `layer_name` at zoom
+# level `zoom` where that feature is used in `rels` relations directly. also
+# needs `all_rows`, a list of all the features, and `all_rels` a list of all
+# the relations in the tile, even those which do not use this feature
+# directly.
+def layer_properties(fid, shape, props, layer_name, zoom, rels,
+                     all_rows, all_rels):
+    layer_props = props.copy()
+
+    # need to make sure that the name is only applied to one of
+    # the pois, landuse or buildings layers - in that order of
+    # priority.
+    #
+    # TODO: do this for all name variants & translations
+    if layer_name in ('pois', 'landuse', 'buildings'):
+        layer_props.pop('name', None)
+
+    # urgh, hack!
+    if layer_name == 'water' and shape.geom_type == 'Point':
+        layer_props['label_placement'] = True
+
+    if shape.geom_type in ('Polygon', 'MultiPolygon'):
+        layer_props['area'] = shape.area
+
+    if layer_name == 'roads' and \
+       shape.geom_type in ('LineString', 'MultiLineString'):
+        mz_networks = []
+        mz_cycling_networks = set()
+        mz_is_bus_route = False
+        for rel in rels:
+            rel_tags = deassoc(rel['tags'])
+            typ, route, network, ref = [rel_tags.get(k) for k in (
+                'type', 'route', 'network', 'ref')]
+            if route and (network or ref):
+                mz_networks.extend([route, network, ref])
+            if typ == 'route' and \
+               route in ('hiking', 'foot', 'bicycle') and \
+               network in ('icn', 'ncn', 'rcn', 'lcn'):
+                mz_cycling_networks.add(network)
+            if typ == 'route' and route in ('bus', 'trolleybus'):
+                mz_is_bus_route = True
+
+        mz_cycling_network = None
+        for cn in ('icn', 'ncn', 'rcn', 'lcn'):
+            if layer_props.get(cn) == 'yes' or \
+               ('%s_ref' % cn) in layer_props or \
+               cn in mz_cycling_networks:
+                mz_cycling_network = cn
+                break
+
+        if mz_is_bus_route and \
+           zoom >= 12 and \
+           layer_props.get('highway') in BUS_ROADS:
+            layer_props['is_bus_route'] = True
+
+        layer_props['mz_networks'] = mz_networks
+        if mz_cycling_network:
+            layer_props['mz_cycling_network'] = mz_cycling_network
+
+    is_poi = layer_name == 'pois'
+    is_railway_station = props.get('railway') == 'station'
+    is_point_or_poly = shape.geom_type in (
+        'Point', 'MultiPoint', 'Polygon', 'MultiPolygon')
+
+    if is_poi and is_railway_station and \
+       is_point_or_poly and fid >= 0:
+        node_id = None
+        way_id = None
+        if shape.geom_type in ('Point', 'MultiPoint'):
+            node_id = fid
+        else:
+            way_id = fid
+
+        transit = mz_calculate_transit_routes_and_score(
+            all_rows, all_rels, node_id, way_id)
+        layer_props['mz_transit_score'] = transit.score
+        layer_props['mz_transit_root_relation_id'] = (
+            transit.root_relation_id)
+        layer_props['train_routes'] = transit.trains
+        layer_props['subway_routes'] = transit.subways
+        layer_props['light_rail_routes'] = transit.light_rails
+        layer_props['tram_routes'] = transit.trams
+
+    return layer_props
+
+
 class DataFetcher(object):
 
     def __init__(self, layers, rows, rels, label_placement_layers):
@@ -358,88 +444,15 @@ def __call__(self, zoom, unpadded_bounds):
                 if layer_name in label_layers:
                     generate_label_placement = True
 
-                layer_props = props.copy()
-                layer_props['min_zoom'] = min_zoom
-
-                # need to make sure that the name is only applied to one of
-                # the pois, landuse or buildings layers - in that order of
-                # priority.
-                #
-                # TODO: do this for all name variants & translations
-                if layer_name in ('pois', 'landuse', 'buildings'):
-                    layer_props.pop('name', None)
-
-                # urgh, hack!
-                if layer_name == 'water' and shape.geom_type == 'Point':
-                    layer_props['label_placement'] = True
-
-                if shape.geom_type in ('Polygon', 'MultiPolygon'):
-                    layer_props['area'] = shape.area
-
-                if layer_name == 'roads' and \
-                   shape.geom_type in ('LineString', 'MultiLineString'):
-                    mz_networks = []
-                    mz_cycling_networks = set()
-                    mz_is_bus_route = False
-                    for rel in rels:
-                        rel_tags = deassoc(rel['tags'])
-                        typ, route, network, ref = [rel_tags.get(k) for k in (
-                            'type', 'route', 'network', 'ref')]
-                        if route and (network or ref):
-                            mz_networks.extend([route, network, ref])
-                        if typ == 'route' and \
-                           route in ('hiking', 'foot', 'bicycle') and \
-                           network in ('icn', 'ncn', 'rcn', 'lcn'):
-                            mz_cycling_networks.add(network)
-                        if typ == 'route' and route in ('bus', 'trolleybus'):
-                            mz_is_bus_route = True
-
-                    mz_cycling_network = None
-                    for cn in ('icn', 'ncn', 'rcn', 'lcn'):
-                        if layer_props.get(cn) == 'yes' or \
-                           ('%s_ref' % cn) in layer_props or \
-                           cn in mz_cycling_networks:
-                            mz_cycling_network = cn
-                            break
+                layer_props = layer_properties(
+                    fid, shape, props, layer_name, zoom, rels,
+                    self.rows, self.rels)
 
-                    if mz_is_bus_route and \
-                       zoom >= 12 and \
-                       layer_props.get('highway') in BUS_ROADS:
-                        layer_props['is_bus_route'] = True
-
-                    layer_props['mz_networks'] = mz_networks
-                    if mz_cycling_network:
-                        layer_props['mz_cycling_network'] = mz_cycling_network
-
-                is_poi = layer_name == 'pois'
-                is_railway_station = props.get('railway') == 'station'
-                is_point_or_poly = shape.geom_type in (
-                    'Point', 'MultiPoint', 'Polygon', 'MultiPolygon')
-
-                if is_poi and is_railway_station and \
-                   is_point_or_poly and fid >= 0:
-                    node_id = None
-                    way_id = None
-                    if shape.geom_type in ('Point', 'MultiPoint'):
-                        node_id = fid
-                    else:
-                        way_id = fid
-
-                    transit = mz_calculate_transit_routes_and_score(
-                        self.rows, self.rels, node_id, way_id)
-                    layer_props['mz_transit_score'] = transit.score
-                    layer_props['mz_transit_root_relation_id'] = (
-                        transit.root_relation_id)
-                    layer_props['train_routes'] = transit.trains
-                    layer_props['subway_routes'] = transit.subways
-                    layer_props['light_rail_routes'] = transit.light_rails
-                    layer_props['tram_routes'] = transit.trams
-
-                if layer_props:
-                    props_name = '__%s_properties__' % layer_name
-                    read_row[props_name] = layer_props
-                    if layer_name == 'water':
-                        has_water_layer = True
+                layer_props['min_zoom'] = min_zoom
+                props_name = '__%s_properties__' % layer_name
+                read_row[props_name] = layer_props
+                if layer_name == 'water':
+                    has_water_layer = True
 
             # if at least one min_zoom / properties match
             if read_row:

From b2a8e671b7eeca960eee1bd540969bca24a17f71 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Wed, 20 Sep 2017 12:05:08 +0100
Subject: [PATCH 02/19] Start adding RAWR DataFetcher.

---
 tests/test_query_rawr.py |  70 ++++++++++++++++++++++++
 tilequeue/query/rawr.py  | 111 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 181 insertions(+)
 create mode 100644 tests/test_query_rawr.py
 create mode 100644 tilequeue/query/rawr.py

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
new file mode 100644
index 00000000..ca352066
--- /dev/null
+++ b/tests/test_query_rawr.py
@@ -0,0 +1,70 @@
+import unittest
+
+
+class TestGetTable(object):
+
+    def __init__(self, tables):
+        self.tables = tables
+
+    def __call__(self, table_name):
+        return self.tables.get(table_name, [])
+
+
+class TestQueryRawr(unittest.TestCase):
+
+    def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
+              layer_name='testlayer'):
+        from tilequeue.query.fixture import LayerInfo
+        from tilequeue.query.rawr import make_rawr_data_fetcher
+        from tilequeue.tile import deg2num
+
+        layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
+        return make_rawr_data_fetcher(layers, tables, tile_pyramid)
+
+    def test_query_simple(self):
+        from shapely.geometry import Point
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        feature_min_zoom = 11
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return feature_min_zoom
+
+        shape = Point(0, 0)
+        # get_table(table_name) should return a generator of rows.
+        tables = TestGetTable({
+            'planet_osm_point': [(0, shape.wkb, {})],
+        })
+
+        zoom = 10
+        max_zoom = zoom + 5
+        coord = mercator_point_to_coord(zoom, shape.x, shape.y)
+        tile_pyramid = TilePyramid(zoom, coord.column, coord.row, max_zoom)
+
+        fetch = self._make(min_zoom_fn, None, tables, tile_pyramid)
+
+        # first, check that it can get the original item back when both the
+        # min zoom filter and geometry filter are okay.
+        feature_coord = mercator_point_to_coord(
+            feature_min_zoom, shape.x, shape.y)
+        read_rows = fetch(
+            feature_min_zoom, coord_to_mercator_bounds(feature_coord))
+
+        self.assertEquals(1, len(read_rows))
+        read_row = read_rows[0]
+        self.assertEquals(0, read_row.get('__id__'))
+        # query processing code expects WKB bytes in the __geometry__ column
+        self.assertEquals(shape.wkb, read_row.get('__geometry__'))
+        self.assertEquals({'min_zoom': 11},
+                          read_row.get('__testlayer_properties__'))
+
+        # now, check that if the min zoom or geometry filters would exclude
+        # the feature then it isn't returned.
+        read_rows = fetch(zoom, coord_to_mercator_bounds(coord))
+        self.assertEquals(0, len(read_rows))
+
+        read_rows = fetch(
+            feature_min_zoom, coord_to_mercator_bounds(feature_coord.left()))
+        self.assertEquals(0, len(read_rows))
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
new file mode 100644
index 00000000..193de3cf
--- /dev/null
+++ b/tilequeue/query/rawr.py
@@ -0,0 +1,111 @@
+from collections import namedtuple
+from shapely.geometry import box
+
+
+class TilePyramid(namedtuple('TilePyramid', 'z x y max_z')):
+
+    def tile(self):
+        from raw_tiles.tile import Tile
+        return Tile(self.z, self.x, self.y)
+
+    def bounds(self):
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+
+        coord = Coordinate(zoom=self.z, column=self.x, row=self.y)
+        bounds = coord_to_mercator_bounds(coord)
+
+        return bounds
+
+    def bbox(self):
+        return box(*self.bounds())
+
+
+class DataFetcher(object):
+
+    def __init__(self, layers, tables, tile_pyramid):
+        """
+        Expect layers to be a dict of layer name to LayerInfo (see fixture.py).
+        Tables should be a callable which returns a generator over the rows in
+        the table when called with that table's name.
+        """
+
+        from raw_tiles.index.features import FeatureTileIndex
+        from raw_tiles.index.index import index_table
+
+        self.layers = layers
+        self.tile_pyramid = tile_pyramid
+        self.layer_indexes = {}
+
+        tile = self.tile_pyramid.tile()
+        max_zoom = self.tile_pyramid.max_z
+
+        for layer_name, info in self.layers.items():
+            meta = None
+
+            def min_zoom(fid, shape, props):
+                return info.min_zoom_fn(fid, shape, props, meta)
+
+            layer_index = FeatureTileIndex(tile, max_zoom, min_zoom)
+
+            for shape_type in ('point', 'line', 'polygon'):
+                if not info.allows_shape_type(shape_type):
+                    continue
+
+                source = tables('planet_osm_' + shape_type)
+                index_table(source, 'add_feature', layer_index)
+
+            self.layer_indexes[layer_name] = layer_index
+
+    def _lookup(self, zoom, unpadded_bounds, layer_name):
+        from tilequeue.tile import mercator_point_to_coord
+        from raw_tiles.tile import Tile
+
+        minx, miny, maxx, maxy = unpadded_bounds
+        topleft = mercator_point_to_coord(zoom, minx, miny)
+        bottomright = mercator_point_to_coord(zoom, maxx, maxy)
+        index = self.layer_indexes[layer_name]
+
+        features = []
+        for x in range(int(topleft.column), int(bottomright.column) + 1):
+            for y in range(int(topleft.row), int(bottomright.row) + 1):
+                tile = Tile(zoom, x, y)
+                features.extend(index(tile))
+        return features
+
+    def __call__(self, zoom, unpadded_bounds):
+        read_rows = []
+        bbox = box(*unpadded_bounds)
+
+        # check that the call is fetching data which is actually within the
+        # bounds of the tile pyramid. we don't have data outside of that, so
+        # can't fulfil requests. if these assertions are tripping, it probably
+        # indicates a programming error - has the wrong DataFetcher been
+        # loaded?
+        assert zoom <= self.tile_pyramid.max_z
+        assert zoom >= self.tile_pyramid.z
+        assert bbox.within(self.tile_pyramid.bbox())
+
+        for layer_name, info in self.layers.items():
+
+            for (fid, shape, props) in self._lookup(
+                    zoom, unpadded_bounds, layer_name):
+                # reject any feature which doesn't intersect the given bounds
+                if bbox.disjoint(shape):
+                    continue
+
+                # place for assembing the read row as if from postgres
+                read_row = {}
+
+                read_row['__' + layer_name + '_properties__'] = props.copy()
+                read_row['__id__'] = fid
+                read_row['__geometry__'] = bytes(shape.wkb)
+                read_rows.append(read_row)
+
+        return read_rows
+
+
+# tables is a callable which should return a generator over the rows of the
+# table when called with the table name.
+def make_rawr_data_fetcher(layers, tables, tile_pyramid):
+    return DataFetcher(layers, tables, tile_pyramid)

From a62f13951658db658f558ab153b24632007976b2 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Wed, 20 Sep 2017 16:54:26 +0100
Subject: [PATCH 03/19] Add extra tests to make sure zoom filtering is working
 as expected in RAWR tiles.

---
 tests/test_query_rawr.py | 65 +++++++++++++++++++++++++++++++++++++++-
 tilequeue/query/rawr.py  |  8 +++++
 2 files changed, 72 insertions(+), 1 deletion(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index ca352066..faaf3352 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -16,7 +16,6 @@ def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
               layer_name='testlayer'):
         from tilequeue.query.fixture import LayerInfo
         from tilequeue.query.rawr import make_rawr_data_fetcher
-        from tilequeue.tile import deg2num
 
         layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
         return make_rawr_data_fetcher(layers, tables, tile_pyramid)
@@ -68,3 +67,67 @@ def min_zoom_fn(shape, props, fid, meta):
         read_rows = fetch(
             feature_min_zoom, coord_to_mercator_bounds(feature_coord.left()))
         self.assertEquals(0, len(read_rows))
+
+    def test_query_min_zoom_fraction(self):
+        from shapely.geometry import Point
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return 11.999
+
+        shape = Point(0, 0)
+        tables = TestGetTable({
+            'planet_osm_point': [(0, shape.wkb, {})]
+        })
+
+        zoom = 10
+        max_zoom = zoom + 5
+        coord = mercator_point_to_coord(zoom, shape.x, shape.y)
+        tile_pyramid = TilePyramid(zoom, coord.column, coord.row, max_zoom)
+
+        fetch = self._make(min_zoom_fn, None, tables, tile_pyramid)
+
+        # check that the fractional zoom of 11.999 means that it's included in
+        # the zoom 11 tile, but not the zoom 10 one.
+        feature_coord = mercator_point_to_coord(11, shape.x, shape.y)
+        read_rows = fetch(11, coord_to_mercator_bounds(feature_coord))
+        self.assertEquals(1, len(read_rows))
+
+        feature_coord = feature_coord.zoomBy(-1).container()
+        read_rows = fetch(10, coord_to_mercator_bounds(feature_coord))
+        self.assertEquals(0, len(read_rows))
+
+    def test_query_past_max_zoom(self):
+        from shapely.geometry import Point
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return 20
+
+        shape = Point(0, 0)
+        tables = TestGetTable({
+            'planet_osm_point': [(0, shape.wkb, {})]
+        })
+
+        zoom = 10
+        max_zoom = zoom + 6
+        coord = mercator_point_to_coord(zoom, shape.x, shape.y)
+        tile_pyramid = TilePyramid(zoom, coord.column, coord.row, max_zoom)
+
+        fetch = self._make(min_zoom_fn, None, tables, tile_pyramid)
+
+        # the min_zoom of 20 should mean that the feature is included at zoom
+        # 16, even though 16<20, because 16 is the "max zoom" at which all the
+        # data is included.
+        feature_coord = mercator_point_to_coord(16, shape.x, shape.y)
+        read_rows = fetch(16, coord_to_mercator_bounds(feature_coord))
+        self.assertEquals(1, len(read_rows))
+
+        # but it should not exist at zoom 15
+        feature_coord = feature_coord.zoomBy(-1).container()
+        read_rows = fetch(10, coord_to_mercator_bounds(feature_coord))
+        self.assertEquals(0, len(read_rows))
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 193de3cf..62ede84e 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -66,6 +66,14 @@ def _lookup(self, zoom, unpadded_bounds, layer_name):
         bottomright = mercator_point_to_coord(zoom, maxx, maxy)
         index = self.layer_indexes[layer_name]
 
+        # make sure that the bottom right coordinate is below and to the right
+        # of the top left coordinate. it can happen that the coordinates are
+        # mixed up due to small numerical precision artefacts being enlarged
+        # by the conversion to integer and y-coordinate flip.
+        assert topleft.zoom == bottomright.zoom
+        bottomright.column = max(bottomright.column, topleft.column)
+        bottomright.row = max(bottomright.row, topleft.row)
+
         features = []
         for x in range(int(topleft.column), int(bottomright.column) + 1):
             for y in range(int(topleft.row), int(bottomright.row) + 1):

From c9d4bc37d4003edea68d13bbbdcc93569ccc9eab Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Thu, 21 Sep 2017 18:14:46 +0100
Subject: [PATCH 04/19] Factor out common code between RAWR and fixtures.

This is now in `common.py`, making it easier to see where the shared code is. Also factored out an interface, implemented by OsmFixtureLookup, which is used to make queries about relationships between elements; mainly station relation complexes and highway routes, etc...
---
 tests/test_query_fixture.py |   2 +-
 tests/test_query_rawr.py    |  71 +++++-
 tilequeue/query/common.py   | 321 ++++++++++++++++++++++++
 tilequeue/query/fixture.py  | 471 +++++++++---------------------------
 tilequeue/query/rawr.py     |  60 ++++-
 5 files changed, 562 insertions(+), 363 deletions(-)
 create mode 100644 tilequeue/query/common.py

diff --git a/tests/test_query_fixture.py b/tests/test_query_fixture.py
index 7dbb04b8..28c7e8f3 100644
--- a/tests/test_query_fixture.py
+++ b/tests/test_query_fixture.py
@@ -5,7 +5,7 @@ class TestQueryFixture(unittest.TestCase):
 
     def _make(self, rows, min_zoom_fn, props_fn, relations=[],
               layer_name='testlayer'):
-        from tilequeue.query.fixture import LayerInfo
+        from tilequeue.query.common import LayerInfo
         from tilequeue.query.fixture import make_fixture_data_fetcher
         layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
         return make_fixture_data_fetcher(layers, rows, relations=relations)
diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index faaf3352..b7c914be 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -14,7 +14,7 @@ class TestQueryRawr(unittest.TestCase):
 
     def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
               layer_name='testlayer'):
-        from tilequeue.query.fixture import LayerInfo
+        from tilequeue.query.common import LayerInfo
         from tilequeue.query.rawr import make_rawr_data_fetcher
 
         layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
@@ -131,3 +131,72 @@ def min_zoom_fn(shape, props, fid, meta):
         feature_coord = feature_coord.zoomBy(-1).container()
         read_rows = fetch(10, coord_to_mercator_bounds(feature_coord))
         self.assertEquals(0, len(read_rows))
+
+    # TODO!
+    # this isn't ready yet! need to implement OsmRawrLookup to use the RAWR
+    # tile indexes.
+    def _test_root_relation_id(self):
+        from shapely.geometry import Point
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return 10
+
+        def _test(rels, expected_root_id):
+            shape = Point(0, 0)
+            props = {
+                'railway': 'station',
+                'name': 'Foo Station',
+            }
+            tables = TestGetTable({
+                'planet_osm_point': [(1, shape.wkb, props)],
+                'planet_osm_rels': rels,
+            })
+
+            zoom = 10
+            max_zoom = zoom + 6
+            coord = mercator_point_to_coord(zoom, shape.x, shape.y)
+            tile_pyramid = TilePyramid(zoom, coord.column, coord.row, max_zoom)
+
+            fetch = self._make(min_zoom_fn, None, tables, tile_pyramid,
+                               layer_name='pois')
+
+            feature_coord = mercator_point_to_coord(16, shape.x, shape.y)
+            read_rows = fetch(16, coord_to_mercator_bounds(feature_coord))
+            self.assertEquals(1, len(read_rows))
+
+            props = read_rows[0]['__pois_properties__']
+            self.assertEquals(expected_root_id,
+                              props.get('mz_transit_root_relation_id'))
+
+        # the fixture code expects "raw" relations as if they come straight
+        # from osm2pgsql. the structure is a little cumbersome, so this
+        # utility function constructs it from a more readable function call.
+        def _rel(id, nodes=None, ways=None, rels=None):
+            way_off = len(nodes) if nodes else 0
+            rel_off = way_off + (len(ways) if ways else 0)
+            return {
+                'id': id,
+                'tags': ['type', 'site'],
+                'way_off': way_off,
+                'rel_off': rel_off,
+                'parts': (nodes or []) + (ways or []) + (rels or []),
+            }
+
+        # one level of relations - this one directly contains the station
+        # node.
+        _test([_rel(2, nodes=[1])], 2)
+
+        # two levels of relations r3 contains r2 contains n1.
+        _test([_rel(2, nodes=[1]), _rel(3, rels=[2])], 3)
+
+        # asymmetric diamond pattern. r2 and r3 both contain n1, r4 contains
+        # r3 and r5 contains both r4 and r2, making it the "top" relation.
+        _test([
+            _rel(2, nodes=[1]),
+            _rel(3, nodes=[1]),
+            _rel(4, rels=[3]),
+            _rel(5, rels=[2, 4]),
+        ], 5)
diff --git a/tilequeue/query/common.py b/tilequeue/query/common.py
new file mode 100644
index 00000000..ebe41f46
--- /dev/null
+++ b/tilequeue/query/common.py
@@ -0,0 +1,321 @@
+from collections import namedtuple
+from collections import defaultdict
+from itertools import izip
+
+
+def namedtuple_with_defaults(name, props, defaults):
+    t = namedtuple(name, props)
+    t.__new__.__defaults__ = defaults
+    return t
+
+
+class LayerInfo(namedtuple_with_defaults(
+        'LayerInfo', 'min_zoom_fn props_fn shape_types', (None,))):
+
+    def allows_shape_type(self, shape):
+        if self.shape_types is None:
+            return True
+        typ = shape_type_lookup(shape)
+        return typ in self.shape_types
+
+
+def deassoc(x):
+    """
+    Turns an array consisting of alternating key-value pairs into a
+    dictionary.
+
+    Osm2pgsql stores the tags for ways and relations in the planet_osm_ways and
+    planet_osm_rels tables in this format. Hstore would make more sense now,
+    but this encoding pre-dates the common availability of hstore.
+
+    Example:
+    >>> from raw_tiles.index.util import deassoc
+    >>> deassoc(['a', 1, 'b', 'B', 'c', 3.14])
+    {'a': 1, 'c': 3.14, 'b': 'B'}
+    """
+
+    pairs = [iter(x)] * 2
+    return dict(izip(*pairs))
+
+
+# fixtures extend metadata to include ways and relations for the feature.
+# this is unnecessary for SQL, as the ways and relations tables are
+# "ambiently available" and do not need to be passed in arguments.
+Metadata = namedtuple('Metadata', 'source ways relations')
+
+
+def shape_type_lookup(shape):
+    typ = shape.geom_type
+    if typ.startswith('Multi'):
+        typ = typ[len('Multi'):]
+    return typ.lower()
+
+
+# list of road types which are likely to have buses on them. used to cut
+# down the number of queries the SQL used to do for relations. although this
+# isn't necessary for fixtures, we replicate the logic to keep the behaviour
+# the same.
+BUS_ROADS = set([
+    'motorway', 'motorway_link', 'trunk', 'trunk_link', 'primary',
+    'primary_link', 'secondary', 'secondary_link', 'tertiary',
+    'tertiary_link', 'residential', 'unclassified', 'road', 'living_street',
+])
+
+
+class Relation(object):
+    def __init__(self, obj):
+        self.id = obj['id']
+        self.tags = deassoc(obj['tags'])
+        way_off = obj['way_off']
+        rel_off = obj['rel_off']
+        self.node_ids = obj['parts'][0:way_off]
+        self.way_ids = obj['parts'][way_off:rel_off]
+        self.rel_ids = obj['parts'][rel_off:]
+
+
+def mz_is_interesting_transit_relation(tags):
+    public_transport = tags.get('public_transport')
+    typ = tags.get('type')
+    return public_transport in ('stop_area', 'stop_area_group') or \
+        typ in ('stop_area', 'stop_area_group', 'site')
+
+
+# starting with the IDs in seed_relations, recurse up the transit relations
+# of which they are members. returns the set of all the relation IDs seen
+# and the "root" relation ID, which was the "furthest" relation from any
+# leaf relation.
+def mz_recurse_up_transit_relations(seed_relations, osm):
+    root_relation_ids = set()
+    root_relation_level = 0
+    all_relations = set()
+
+    for rel_id in seed_relations:
+        front = set([rel_id])
+        seen = set([rel_id])
+        level = 0
+
+        if root_relation_level == 0:
+            root_relation_ids.add(rel_id)
+
+        while front:
+            new_rels = set()
+            for r in front:
+                new_rels |= osm.transit_relations(r)
+            new_rels -= seen
+            level += 1
+            if new_rels and level > root_relation_level:
+                root_relation_ids = new_rels
+                root_relation_level = level
+            elif new_rels and level == root_relation_level:
+                root_relation_ids |= new_rels
+            front = new_rels
+            seen |= front
+
+        all_relations |= seen
+
+    root_relation_id = min(root_relation_ids) if root_relation_ids else None
+    return all_relations, root_relation_id
+
+
+# extract a name for a transit route relation. this can expand comma
+# separated lists and prefers to use the ref rather than the name.
+def mz_transit_route_name(tags):
+    # prefer ref as it's less likely to contain the destination name
+    name = tags.get('ref')
+    if not name:
+        name = tags.get('name')
+    if name:
+        name = name.strip()
+    return name
+
+
+Transit = namedtuple(
+    'Transit', 'score root_relation_id '
+    'trains subways light_rails trams railways')
+
+
+def mz_calculate_transit_routes_and_score(osm, node_id, way_id):
+    candidate_relations = set()
+    if node_id:
+        candidate_relations.update(osm.relations_using_node(node_id))
+    if way_id:
+        candidate_relations.update(osm.relations_using_way(way_id))
+
+    seed_relations = set()
+    for rel_id in candidate_relations:
+        rel = osm.relation(rel_id)
+        if mz_is_interesting_transit_relation(rel.tags):
+            seed_relations.add(rel_id)
+    del candidate_relations
+
+    # TODO: if the station is also a multipolygon relation?
+
+    # this complex query does two recursive sweeps of the relations
+    # table starting from a seed set of relations which are or contain
+    # the original station.
+    #
+    # the first sweep goes "upwards" from relations to "parent" relations. if
+    # a relation R1 is a member of relation R2, then R2 will be included in
+    # this sweep as long as it has "interesting" tags, as defined by the
+    # function mz_is_interesting_transit_relation.
+    #
+    # the second sweep goes "downwards" from relations to "child" relations.
+    # if a relation R1 has a member R2 which is also a relation, then R2 will
+    # be included in this sweep as long as it also has "interesting" tags.
+    all_relations, root_relation_id = mz_recurse_up_transit_relations(
+        seed_relations, osm)
+    del seed_relations
+
+    # collect all the interesting nodes - this includes the station node (if
+    # any) and any nodes which are members of found relations which have
+    # public transport tags indicating that they're stations or stops.
+    stations_and_stops = set()
+    for rel_id in all_relations:
+        rel = osm.relation(rel_id)
+        for node_id in rel.node_ids:
+            fid, shape, props = osm.node(node_id)
+            railway = props.get('railway') in ('station', 'stop', 'tram_stop')
+            public_transport = props.get('public_transport') in \
+                ('stop', 'stop_position', 'tram_stop')
+            if railway or public_transport:
+                stations_and_stops.add(fid)
+
+    if node_id:
+        stations_and_stops.add(node_id)
+
+    # collect any physical railway which includes any of the above
+    # nodes.
+    stations_and_lines = set()
+    for node_id in stations_and_stops:
+        for way_id in osm.ways_using_node(node_id):
+            fid, shape, props = osm.way(way_id)
+            railway = props.get('railway')
+            if railway in ('subway', 'light_rail', 'tram', 'rail'):
+                stations_and_lines.add(way_id)
+
+    if way_id:
+        stations_and_lines.add(way_id)
+
+    # collect all IDs together in one array to intersect with the parts arrays
+    # of route relations which may include them.
+    all_routes = set()
+    for lookup, ids in ((osm.relations_using_node, stations_and_stops),
+                        (osm.relations_using_way, stations_and_lines),
+                        (osm.relations_using_rel, all_relations)):
+        for i in ids:
+            for rel_id in lookup(i):
+                rel = osm.relation(rel_id)
+                if rel.tags.get('type') == 'route' and \
+                   rel.tags.get('route') in ('subway', 'light_rail', 'tram',
+                                             'train', 'railway'):
+                    all_routes.add(rel_id)
+
+    routes_lookup = defaultdict(set)
+    for rel_id in all_routes:
+        rel = osm.relations(rel_id)
+        route = rel.tags.get('route')
+        if route:
+            route_name = mz_transit_route_name(rel.tags)
+            routes_lookup[route].add(route_name)
+    trains = routes_lookup['train']
+    subways = routes_lookup['subway']
+    light_rails = routes_lookup['light_rail']
+    trams = routes_lookup['tram']
+    railways = routes_lookup['railway']
+    del routes_lookup
+
+    # if a station is an interchange between mainline rail and subway or
+    # light rail, then give it a "bonus" boost of importance.
+    bonus = 2 if trains and (subways or light_rails) else 1
+
+    score = (100 * min(9, bonus * len(trains)) +
+             10 * min(9, bonus * (len(subways) + len(light_rails))) +
+             min(9, len(trams) + len(railways)))
+
+    return Transit(score=score, root_relation_id=root_relation_id,
+                   trains=trains, subways=subways, light_rails=light_rails,
+                   railways=railways, trams=trams)
+
+
+# properties for a feature (fid, shape, props) in layer `layer_name` at zoom
+# level `zoom`. also takes an `osm` parameter, which is an object which can
+# be used to look up nodes, ways and relations and the relationships between
+# them.
+def layer_properties(fid, shape, props, layer_name, zoom, osm):
+    layer_props = props.copy()
+
+    # need to make sure that the name is only applied to one of
+    # the pois, landuse or buildings layers - in that order of
+    # priority.
+    #
+    # TODO: do this for all name variants & translations
+    if layer_name in ('pois', 'landuse', 'buildings'):
+        layer_props.pop('name', None)
+
+    # urgh, hack!
+    if layer_name == 'water' and shape.geom_type == 'Point':
+        layer_props['label_placement'] = True
+
+    if shape.geom_type in ('Polygon', 'MultiPolygon'):
+        layer_props['area'] = shape.area
+
+    if layer_name == 'roads' and \
+       shape.geom_type in ('LineString', 'MultiLineString') and \
+       fid >= 0:
+        mz_networks = []
+        mz_cycling_networks = set()
+        mz_is_bus_route = False
+        for rel in osm.relations_using_way(fid):
+            typ, route, network, ref = [rel.tags.get(k) for k in (
+                'type', 'route', 'network', 'ref')]
+            if route and (network or ref):
+                mz_networks.extend([route, network, ref])
+            if typ == 'route' and \
+               route in ('hiking', 'foot', 'bicycle') and \
+               network in ('icn', 'ncn', 'rcn', 'lcn'):
+                mz_cycling_networks.add(network)
+            if typ == 'route' and route in ('bus', 'trolleybus'):
+                mz_is_bus_route = True
+
+        mz_cycling_network = None
+        for cn in ('icn', 'ncn', 'rcn', 'lcn'):
+            if layer_props.get(cn) == 'yes' or \
+               ('%s_ref' % cn) in layer_props or \
+               cn in mz_cycling_networks:
+                mz_cycling_network = cn
+                break
+
+        if mz_is_bus_route and \
+           zoom >= 12 and \
+           layer_props.get('highway') in BUS_ROADS:
+            layer_props['is_bus_route'] = True
+
+        layer_props['mz_networks'] = mz_networks
+        if mz_cycling_network:
+            layer_props['mz_cycling_network'] = mz_cycling_network
+
+    is_poi = layer_name == 'pois'
+    is_railway_station = props.get('railway') == 'station'
+    is_point_or_poly = shape.geom_type in (
+        'Point', 'MultiPoint', 'Polygon', 'MultiPolygon')
+
+    if is_poi and is_railway_station and \
+       is_point_or_poly and fid >= 0:
+        node_id = None
+        way_id = None
+        if shape.geom_type in ('Point', 'MultiPoint'):
+            node_id = fid
+        else:
+            way_id = fid
+
+        transit = mz_calculate_transit_routes_and_score(
+            osm, node_id, way_id)
+        layer_props['mz_transit_score'] = transit.score
+        layer_props['mz_transit_root_relation_id'] = (
+            transit.root_relation_id)
+        layer_props['train_routes'] = transit.trains
+        layer_props['subway_routes'] = transit.subways
+        layer_props['light_rail_routes'] = transit.light_rails
+        layer_props['tram_routes'] = transit.trams
+
+    return layer_props
diff --git a/tilequeue/query/fixture.py b/tilequeue/query/fixture.py
index bf73e305..bed3f685 100644
--- a/tilequeue/query/fixture.py
+++ b/tilequeue/query/fixture.py
@@ -1,367 +1,118 @@
-from collections import namedtuple
-from collections import defaultdict
 from shapely.geometry import box
 from tilequeue.process import lookup_source
-from itertools import izip
 from tilequeue.transform import calculate_padded_bounds
+from tilequeue.query.common import Metadata
+from tilequeue.query.common import Relation
+from tilequeue.query.common import layer_properties
+from tilequeue.query.common import shape_type_lookup
+from tilequeue.query.common import mz_is_interesting_transit_relation
+from collections import defaultdict
+
+
+class OsmFixtureLookup(object):
+
+    def __init__(self, rows, rels):
+        # extract out all relations and index by ID. this is helpful when
+        # looking them up later.
+        relations = {}
+        nodes = {}
+        ways = {}
+        ways_using_node = {}
+
+        for (fid, shape, props) in rows:
+            if fid >= 0:
+                if shape.geom_type in ('Point', 'MultiPoint'):
+                    nodes[fid] = (fid, shape, props)
+                    features = props.get('__ways__', [])
+                    ways_using_node[fid] = [f[0] for f in features]
+                else:
+                    ways[fid] = (fid, shape, props)
+
+        for r in rels:
+            r = Relation(r)
+            assert r.id not in relations
+            relations[r.id] = r
+
+        relations_using_node = defaultdict(list)
+        relations_using_way = defaultdict(list)
+        relations_using_rel = defaultdict(list)
+
+        for rel_id, rel in relations.items():
+            for (ids, index) in ((rel.node_ids, relations_using_node),
+                                 (rel.way_ids, relations_using_way),
+                                 (rel.rel_ids, relations_using_rel)):
+                for osm_id in ids:
+                    index[osm_id].append(rel_id)
+
+        transit_relations = defaultdict(set)
+        for rel_id, rel in relations.items():
+            if mz_is_interesting_transit_relation(rel.tags):
+                for member in rel.rel_ids:
+                    transit_relations[member].add(rel_id)
+
+        # looks up relation IDs
+        self._relations_using_node = relations_using_node
+        self._relations_using_way = relations_using_way
+        self._relations_using_rel = relations_using_rel
+        # looks up way IDs
+        self._ways_using_node = ways_using_node
+        # looks up Relation objects
+        self._relations = relations
+        # looks up (fid, shape, props) feature objects
+        self._ways = ways
+        self._nodes = nodes
+        # returns the set of transit relation IDs that contain the given
+        # relation IDs
+        self._transit_relations = transit_relations
+
+    def relations_using_node(self, node_id):
+        "Returns a list of relation IDs which contain the node with that ID."
+
+        return self._relations_using_node.get(node_id, [])
+
+    def relations_using_way(self, way_id):
+        "Returns a list of relation IDs which contain the way with that ID."
+
+        return self._relations_using_way.get(way_id, [])
+
+    def relations_using_rel(self, rel_id):
+        """
+        Returns a list of relation IDs which contain the relation with that
+        ID.
+        """
+
+        return self._relations_using_rel.get(rel_id, [])
+
+    def ways_using_node(self, node_id):
+        "Returns a list of way IDs which contain the node with that ID."
+
+        return self._ways_using_node.get(node_id, [])
 
+    def relation(self, rel_id):
+        "Returns the Relation object with the given ID."
+
+        return self._relations[rel_id]
+
+    def way(self, way_id):
+        """
+        Returns the feature (fid, shape, props) which was generated from the
+        given way.
+        """
+
+        return self._ways[way_id]
+
+    def node(self, node_id):
+        """
+        Returns the feature (fid, shape, props) which was generated from the
+        given node.
+        """
 
-def namedtuple_with_defaults(name, props, defaults):
-    t = namedtuple(name, props)
-    t.__new__.__defaults__ = defaults
-    return t
+        return self._nodes[node_id]
 
+    def transit_relations(self, rel_id):
+        "Return transit relations containing the relation with the given ID."
 
-class LayerInfo(namedtuple_with_defaults(
-        'LayerInfo', 'min_zoom_fn props_fn shape_types', (None,))):
-
-    def allows_shape_type(self, shape):
-        if self.shape_types is None:
-            return True
-        typ = _shape_type_lookup(shape)
-        return typ in self.shape_types
-
-
-def deassoc(x):
-    """
-    Turns an array consisting of alternating key-value pairs into a
-    dictionary.
-
-    Osm2pgsql stores the tags for ways and relations in the planet_osm_ways and
-    planet_osm_rels tables in this format. Hstore would make more sense now,
-    but this encoding pre-dates the common availability of hstore.
-
-    Example:
-    >>> from raw_tiles.index.util import deassoc
-    >>> deassoc(['a', 1, 'b', 'B', 'c', 3.14])
-    {'a': 1, 'c': 3.14, 'b': 'B'}
-    """
-
-    pairs = [iter(x)] * 2
-    return dict(izip(*pairs))
-
-
-# fixtures extend metadata to include ways and relations for the feature.
-# this is unnecessary for SQL, as the ways and relations tables are
-# "ambiently available" and do not need to be passed in arguments.
-Metadata = namedtuple('Metadata', 'source ways relations')
-
-
-def _shape_type_lookup(shape):
-    typ = shape.geom_type
-    if typ.startswith('Multi'):
-        typ = typ[len('Multi'):]
-    return typ.lower()
-
-
-# list of road types which are likely to have buses on them. used to cut
-# down the number of queries the SQL used to do for relations. although this
-# isn't necessary for fixtures, we replicate the logic to keep the behaviour
-# the same.
-BUS_ROADS = set([
-    'motorway', 'motorway_link', 'trunk', 'trunk_link', 'primary',
-    'primary_link', 'secondary', 'secondary_link', 'tertiary',
-    'tertiary_link', 'residential', 'unclassified', 'road', 'living_street',
-])
-
-
-class Relation(object):
-    def __init__(self, obj):
-        self.id = obj['id']
-        self.tags = deassoc(obj['tags'])
-        way_off = obj['way_off']
-        rel_off = obj['rel_off']
-        self.node_ids = obj['parts'][0:way_off]
-        self.way_ids = obj['parts'][way_off:rel_off]
-        self.rel_ids = obj['parts'][rel_off:]
-
-
-def mz_is_interesting_transit_relation(tags):
-    public_transport = tags.get('public_transport')
-    typ = tags.get('type')
-    return public_transport in ('stop_area', 'stop_area_group') or \
-        typ in ('stop_area', 'stop_area_group', 'site')
-
-
-# starting with the IDs in seed_relations, recurse up the transit relations
-# of which they are members. returns the set of all the relation IDs seen
-# and the "root" relation ID, which was the "furthest" relation from any
-# leaf relation.
-def mz_recurse_up_transit_relations(seed_relations, relations):
-    transit_relations = defaultdict(set)
-    for rel_id, rel in relations.items():
-        if mz_is_interesting_transit_relation(rel.tags):
-            for member in rel.rel_ids:
-                transit_relations[member].add(rel_id)
-
-    root_relation_ids = set()
-    root_relation_level = 0
-    all_relations = set()
-
-    for rel_id in seed_relations:
-        front = set([rel_id])
-        seen = set([rel_id])
-        level = 0
-
-        if root_relation_level == 0:
-            root_relation_ids.add(rel_id)
-
-        while front:
-            new_rels = set()
-            for r in front:
-                new_rels |= transit_relations[r]
-            new_rels -= seen
-            level += 1
-            if new_rels and level > root_relation_level:
-                root_relation_ids = new_rels
-                root_relation_level = level
-            elif new_rels and level == root_relation_level:
-                root_relation_ids |= new_rels
-            front = new_rels
-            seen |= front
-
-        all_relations |= seen
-
-    root_relation_id = min(root_relation_ids) if root_relation_ids else None
-    return all_relations, root_relation_id
-
-
-# extract a name for a transit route relation. this can expand comma
-# separated lists and prefers to use the ref rather than the name.
-def mz_transit_route_name(tags):
-    # prefer ref as it's less likely to contain the destination name
-    name = tags.get('ref')
-    if not name:
-        name = tags.get('name')
-    if name:
-        name = name.strip()
-    return name
-
-
-Transit = namedtuple(
-    'Transit', 'score root_relation_id '
-    'trains subways light_rails trams railways')
-
-
-def mz_calculate_transit_routes_and_score(rows, rels, node_id, way_id):
-    # extract out all relations and index by ID. this is helpful when
-    # looking them up later.
-    relations = {}
-    nodes = {}
-    ways = {}
-    ways_using_node = {}
-
-    for (fid, shape, props) in rows:
-        if fid >= 0:
-            if shape.geom_type in ('Point', 'MultiPoint'):
-                nodes[fid] = (fid, shape, props)
-                features = props.get('__ways__', [])
-                ways_using_node[fid] = [f[0] for f in features]
-            else:
-                ways[fid] = (fid, shape, props)
-
-    for r in rels:
-        r = Relation(r)
-        assert r.id not in relations
-        relations[r.id] = r
-
-    relations_using_node = defaultdict(list)
-    relations_using_way = defaultdict(list)
-    relations_using_rel = defaultdict(list)
-
-    for rel_id, rel in relations.items():
-        for (ids, index) in ((rel.node_ids, relations_using_node),
-                             (rel.way_ids, relations_using_way),
-                             (rel.rel_ids, relations_using_rel)):
-            for osm_id in ids:
-                index[osm_id].append(rel_id)
-
-    candidate_relations = set()
-    if node_id:
-        candidate_relations.update(relations_using_node.get(node_id, []))
-    if way_id:
-        candidate_relations.update(relations_using_way.get(way_id, []))
-
-    seed_relations = set()
-    for rel_id in candidate_relations:
-        rel = relations[rel_id]
-        if mz_is_interesting_transit_relation(rel.tags):
-            seed_relations.add(rel_id)
-    del candidate_relations
-
-    # TODO: if the station is also a multipolygon relation?
-
-    # this complex query does two recursive sweeps of the relations
-    # table starting from a seed set of relations which are or contain
-    # the original station.
-    #
-    # the first sweep goes "upwards" from relations to "parent" relations. if
-    # a relation R1 is a member of relation R2, then R2 will be included in
-    # this sweep as long as it has "interesting" tags, as defined by the
-    # function mz_is_interesting_transit_relation.
-    #
-    # the second sweep goes "downwards" from relations to "child" relations.
-    # if a relation R1 has a member R2 which is also a relation, then R2 will
-    # be included in this sweep as long as it also has "interesting" tags.
-    all_relations, root_relation_id = mz_recurse_up_transit_relations(
-        seed_relations, relations)
-    del seed_relations
-
-    # collect all the interesting nodes - this includes the station node (if
-    # any) and any nodes which are members of found relations which have
-    # public transport tags indicating that they're stations or stops.
-    stations_and_stops = set()
-    for rel_id in all_relations:
-        rel = relations[rel_id]
-        for node_id in rel.node_ids:
-            fid, shape, props = nodes[node_id]
-            railway = props.get('railway') in ('station', 'stop', 'tram_stop')
-            public_transport = props.get('public_transport') in \
-                ('stop', 'stop_position', 'tram_stop')
-            if railway or public_transport:
-                stations_and_stops.add(fid)
-
-    if node_id:
-        stations_and_stops.add(node_id)
-
-    # collect any physical railway which includes any of the above
-    # nodes.
-    stations_and_lines = set()
-    for node_id in stations_and_stops:
-        for way_id in ways_using_node[node_id]:
-            fid, shape, props = ways[way_id]
-            railway = props.get('railway')
-            if railway in ('subway', 'light_rail', 'tram', 'rail'):
-                stations_and_lines.add(way_id)
-
-    if way_id:
-        stations_and_lines.add(way_id)
-
-    # collect all IDs together in one array to intersect with the parts arrays
-    # of route relations which may include them.
-    all_routes = set()
-    for lookup, ids in ((relations_using_node, stations_and_stops),
-                        (relations_using_way, stations_and_lines),
-                        (relations_using_rel, all_relations)):
-        for i in ids:
-            for rel_id in lookup.get(i, []):
-                rel = relations[rel_id]
-                if rel.tags.get('type') == 'route' and \
-                   rel.tags.get('route') in ('subway', 'light_rail', 'tram',
-                                             'train', 'railway'):
-                    all_routes.add(rel_id)
-
-    routes_lookup = defaultdict(set)
-    for rel_id in all_routes:
-        rel = relations[rel_id]
-        route = rel.tags.get('route')
-        if route:
-            route_name = mz_transit_route_name(rel.tags)
-            routes_lookup[route].add(route_name)
-    trains = routes_lookup['train']
-    subways = routes_lookup['subway']
-    light_rails = routes_lookup['light_rail']
-    trams = routes_lookup['tram']
-    railways = routes_lookup['railway']
-    del routes_lookup
-
-    # if a station is an interchange between mainline rail and subway or
-    # light rail, then give it a "bonus" boost of importance.
-    bonus = 2 if trains and (subways or light_rails) else 1
-
-    score = (100 * min(9, bonus * len(trains)) +
-             10 * min(9, bonus * (len(subways) + len(light_rails))) +
-             min(9, len(trams) + len(railways)))
-
-    return Transit(score=score, root_relation_id=root_relation_id,
-                   trains=trains, subways=subways, light_rails=light_rails,
-                   railways=railways, trams=trams)
-
-
-# properties for a feature (fid, shape, props) in layer `layer_name` at zoom
-# level `zoom` where that feature is used in `rels` relations directly. also
-# needs `all_rows`, a list of all the features, and `all_rels` a list of all
-# the relations in the tile, even those which do not use this feature
-# directly.
-def layer_properties(fid, shape, props, layer_name, zoom, rels,
-                     all_rows, all_rels):
-    layer_props = props.copy()
-
-    # need to make sure that the name is only applied to one of
-    # the pois, landuse or buildings layers - in that order of
-    # priority.
-    #
-    # TODO: do this for all name variants & translations
-    if layer_name in ('pois', 'landuse', 'buildings'):
-        layer_props.pop('name', None)
-
-    # urgh, hack!
-    if layer_name == 'water' and shape.geom_type == 'Point':
-        layer_props['label_placement'] = True
-
-    if shape.geom_type in ('Polygon', 'MultiPolygon'):
-        layer_props['area'] = shape.area
-
-    if layer_name == 'roads' and \
-       shape.geom_type in ('LineString', 'MultiLineString'):
-        mz_networks = []
-        mz_cycling_networks = set()
-        mz_is_bus_route = False
-        for rel in rels:
-            rel_tags = deassoc(rel['tags'])
-            typ, route, network, ref = [rel_tags.get(k) for k in (
-                'type', 'route', 'network', 'ref')]
-            if route and (network or ref):
-                mz_networks.extend([route, network, ref])
-            if typ == 'route' and \
-               route in ('hiking', 'foot', 'bicycle') and \
-               network in ('icn', 'ncn', 'rcn', 'lcn'):
-                mz_cycling_networks.add(network)
-            if typ == 'route' and route in ('bus', 'trolleybus'):
-                mz_is_bus_route = True
-
-        mz_cycling_network = None
-        for cn in ('icn', 'ncn', 'rcn', 'lcn'):
-            if layer_props.get(cn) == 'yes' or \
-               ('%s_ref' % cn) in layer_props or \
-               cn in mz_cycling_networks:
-                mz_cycling_network = cn
-                break
-
-        if mz_is_bus_route and \
-           zoom >= 12 and \
-           layer_props.get('highway') in BUS_ROADS:
-            layer_props['is_bus_route'] = True
-
-        layer_props['mz_networks'] = mz_networks
-        if mz_cycling_network:
-            layer_props['mz_cycling_network'] = mz_cycling_network
-
-    is_poi = layer_name == 'pois'
-    is_railway_station = props.get('railway') == 'station'
-    is_point_or_poly = shape.geom_type in (
-        'Point', 'MultiPoint', 'Polygon', 'MultiPolygon')
-
-    if is_poi and is_railway_station and \
-       is_point_or_poly and fid >= 0:
-        node_id = None
-        way_id = None
-        if shape.geom_type in ('Point', 'MultiPoint'):
-            node_id = fid
-        else:
-            way_id = fid
-
-        transit = mz_calculate_transit_routes_and_score(
-            all_rows, all_rels, node_id, way_id)
-        layer_props['mz_transit_score'] = transit.score
-        layer_props['mz_transit_root_relation_id'] = (
-            transit.root_relation_id)
-        layer_props['train_routes'] = transit.trains
-        layer_props['subway_routes'] = transit.subways
-        layer_props['light_rail_routes'] = transit.light_rails
-        layer_props['tram_routes'] = transit.trams
-
-    return layer_props
+        return self._transit_relations.get(rel_id, set())
 
 
 class DataFetcher(object):
@@ -377,6 +128,7 @@ def __init__(self, layers, rows, rels, label_placement_layers):
         self.rows = rows
         self.rels = rels
         self.label_placement_layers = label_placement_layers
+        self.osm = OsmFixtureLookup(self.rows, self.rels)
 
     def __call__(self, zoom, unpadded_bounds):
         read_rows = []
@@ -440,13 +192,12 @@ def __call__(self, zoom, unpadded_bounds):
                 # if the feature exists in any label placement layer, then we
                 # should consider generating a centroid
                 label_layers = self.label_placement_layers.get(
-                    _shape_type_lookup(shape), {})
+                    shape_type_lookup(shape), {})
                 if layer_name in label_layers:
                     generate_label_placement = True
 
                 layer_props = layer_properties(
-                    fid, shape, props, layer_name, zoom, rels,
-                    self.rows, self.rels)
+                    fid, shape, props, layer_name, zoom, self.osm)
 
                 layer_props['min_zoom'] = min_zoom
                 props_name = '__%s_properties__' % layer_name
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 62ede84e..9afc7359 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -1,5 +1,6 @@
 from collections import namedtuple
 from shapely.geometry import box
+from tilequeue.query.common import layer_properties
 
 
 class TilePyramid(namedtuple('TilePyramid', 'z x y max_z')):
@@ -21,6 +22,58 @@ def bbox(self):
         return box(*self.bounds())
 
 
+class OsmRawrLookup(object):
+
+    def relations_using_node(self, node_id):
+        "Returns a list of relation IDs which contain the node with that ID."
+
+        return []
+
+    def relations_using_way(self, way_id):
+        "Returns a list of relation IDs which contain the way with that ID."
+
+        return []
+
+    def relations_using_rel(self, rel_id):
+        """
+        Returns a list of relation IDs which contain the relation with that
+        ID.
+        """
+
+        return []
+
+    def ways_using_node(self, node_id):
+        "Returns a list of way IDs which contain the node with that ID."
+
+        return []
+
+    def relation(self, rel_id):
+        "Returns the Relation object with the given ID."
+
+        return None
+
+    def way(self, way_id):
+        """
+        Returns the feature (fid, shape, props) which was generated from the
+        given way.
+        """
+
+        return None
+
+    def node(self, node_id):
+        """
+        Returns the feature (fid, shape, props) which was generated from the
+        given node.
+        """
+
+        return None
+
+    def transit_relations(self, rel_id):
+        "Return transit relations containing the relation with the given ID."
+
+        return set()
+
+
 class DataFetcher(object):
 
     def __init__(self, layers, tables, tile_pyramid):
@@ -57,6 +110,8 @@ def min_zoom(fid, shape, props):
 
             self.layer_indexes[layer_name] = layer_index
 
+        self.osm = OsmRawrLookup()
+
     def _lookup(self, zoom, unpadded_bounds, layer_name):
         from tilequeue.tile import mercator_point_to_coord
         from raw_tiles.tile import Tile
@@ -105,7 +160,10 @@ def __call__(self, zoom, unpadded_bounds):
                 # place for assembing the read row as if from postgres
                 read_row = {}
 
-                read_row['__' + layer_name + '_properties__'] = props.copy()
+                layer_props = layer_properties(
+                    fid, shape, props, layer_name, zoom, self.osm)
+
+                read_row['__' + layer_name + '_properties__'] = layer_props
                 read_row['__id__'] = fid
                 read_row['__geometry__'] = bytes(shape.wkb)
                 read_rows.append(read_row)

From 7e21eb6877be006aca4c84d0c4aad76d234c6e07 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Fri, 22 Sep 2017 19:57:44 +0100
Subject: [PATCH 05/19] Added root relation ID test.

---
 tests/test_query_rawr.py  |  16 ++----
 tilequeue/query/common.py |  35 ++++++++----
 tilequeue/query/rawr.py   | 117 +++++++++++++++++++++++++++++++++-----
 3 files changed, 132 insertions(+), 36 deletions(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index b7c914be..699d3a6e 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -132,10 +132,7 @@ def min_zoom_fn(shape, props, fid, meta):
         read_rows = fetch(10, coord_to_mercator_bounds(feature_coord))
         self.assertEquals(0, len(read_rows))
 
-    # TODO!
-    # this isn't ready yet! need to implement OsmRawrLookup to use the RAWR
-    # tile indexes.
-    def _test_root_relation_id(self):
+    def test_root_relation_id(self):
         from shapely.geometry import Point
         from tilequeue.query.rawr import TilePyramid
         from tilequeue.tile import coord_to_mercator_bounds
@@ -177,13 +174,10 @@ def _test(rels, expected_root_id):
         def _rel(id, nodes=None, ways=None, rels=None):
             way_off = len(nodes) if nodes else 0
             rel_off = way_off + (len(ways) if ways else 0)
-            return {
-                'id': id,
-                'tags': ['type', 'site'],
-                'way_off': way_off,
-                'rel_off': rel_off,
-                'parts': (nodes or []) + (ways or []) + (rels or []),
-            }
+            parts = (nodes or []) + (ways or []) + (rels or [])
+            members = [""] * len(parts)
+            tags = ['type', 'site']
+            return (id, way_off, rel_off, parts, members, tags)
 
         # one level of relations - this one directly contains the station
         # node.
diff --git a/tilequeue/query/common.py b/tilequeue/query/common.py
index ebe41f46..bbbc29ab 100644
--- a/tilequeue/query/common.py
+++ b/tilequeue/query/common.py
@@ -129,6 +129,24 @@ def mz_transit_route_name(tags):
     return name
 
 
+def is_station_or_stop(fid, shape, props):
+    "Returns true if the given (point) feature is a station or stop."
+    return (
+        props.get('railway') in ('station', 'stop', 'tram_stop') or
+        props.get('public_transport') in ('stop', 'stop_position', 'tram_stop')
+    )
+
+
+def is_station_or_line(fid, shape, props):
+    """
+    Returns true if the given (line or polygon from way) feature is a station
+    or transit line.
+    """
+
+    railway = props.get('railway')
+    return railway in ('subway', 'light_rail', 'tram', 'rail')
+
+
 Transit = namedtuple(
     'Transit', 'score root_relation_id '
     'trains subways light_rails trams railways')
@@ -173,12 +191,8 @@ def mz_calculate_transit_routes_and_score(osm, node_id, way_id):
     for rel_id in all_relations:
         rel = osm.relation(rel_id)
         for node_id in rel.node_ids:
-            fid, shape, props = osm.node(node_id)
-            railway = props.get('railway') in ('station', 'stop', 'tram_stop')
-            public_transport = props.get('public_transport') in \
-                ('stop', 'stop_position', 'tram_stop')
-            if railway or public_transport:
-                stations_and_stops.add(fid)
+            if is_station_or_stop(*osm.node(node_id)):
+                stations_and_stops.add(node_id)
 
     if node_id:
         stations_and_stops.add(node_id)
@@ -188,9 +202,7 @@ def mz_calculate_transit_routes_and_score(osm, node_id, way_id):
     stations_and_lines = set()
     for node_id in stations_and_stops:
         for way_id in osm.ways_using_node(node_id):
-            fid, shape, props = osm.way(way_id)
-            railway = props.get('railway')
-            if railway in ('subway', 'light_rail', 'tram', 'rail'):
+            if is_station_or_line(*osm.way(way_id)):
                 stations_and_lines.add(way_id)
 
     if way_id:
@@ -212,7 +224,7 @@ def mz_calculate_transit_routes_and_score(osm, node_id, way_id):
 
     routes_lookup = defaultdict(set)
     for rel_id in all_routes:
-        rel = osm.relations(rel_id)
+        rel = osm.relation(rel_id)
         route = rel.tags.get('route')
         if route:
             route_name = mz_transit_route_name(rel.tags)
@@ -265,7 +277,8 @@ def layer_properties(fid, shape, props, layer_name, zoom, osm):
         mz_networks = []
         mz_cycling_networks = set()
         mz_is_bus_route = False
-        for rel in osm.relations_using_way(fid):
+        for rel_id in osm.relations_using_way(fid):
+            rel = osm.relation(rel_id)
             typ, route, network, ref = [rel.tags.get(k) for k in (
                 'type', 'route', 'network', 'ref')]
             if route and (network or ref):
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 9afc7359..7c2867a4 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -1,6 +1,19 @@
-from collections import namedtuple
+from collections import namedtuple, defaultdict
 from shapely.geometry import box
 from tilequeue.query.common import layer_properties
+from tilequeue.query.common import is_station_or_stop
+from tilequeue.query.common import is_station_or_line
+from tilequeue.query.common import deassoc
+from tilequeue.query.common import mz_is_interesting_transit_relation
+
+
+class Relation(object):
+    def __init__(self, rel_id, way_off, rel_off, parts, members, tags):
+        self.id = rel_id
+        self.tags = deassoc(tags)
+        self.node_ids = parts[0:way_off]
+        self.way_ids = parts[way_off:rel_off]
+        self.rel_ids = parts[rel_off:]
 
 
 class TilePyramid(namedtuple('TilePyramid', 'z x y max_z')):
@@ -22,17 +35,84 @@ def bbox(self):
         return box(*self.bounds())
 
 
+# weak type of enum type
+class ShapeType(object):
+    point = 1
+    line = 2
+    polygon = 3
+
+
+def _wkb_shape(wkb):
+    reverse = ord(wkb[0]) == 1
+    type_bytes = map(ord, wkb[1:5])
+    if reverse:
+        type_bytes.reverse()
+    typ = type_bytes[3]
+    if typ == 1 or typ == 4:
+        return ShapeType.point
+    elif typ == 2 or typ == 5:
+        return ShapeType.line
+    elif typ == 3 or typ == 6:
+        return ShapeType.polygon
+    else:
+        assert False, "WKB shape type %d not understood." % (typ,)
+
+
 class OsmRawrLookup(object):
 
+    def __init__(self):
+        self.nodes = {}
+        self.ways = {}
+        self.relations = {}
+
+        self._ways_using_node = defaultdict(list)
+        self._relations_using_node = defaultdict(list)
+        self._relations_using_way = defaultdict(list)
+        self._relations_using_rel = defaultdict(list)
+
+    def add_feature(self, fid, shape_wkb, props):
+        if fid < 0:
+            return
+
+        shape_type = _wkb_shape(shape_wkb)
+        if is_station_or_stop(fid, None, props) and \
+           shape_type == ShapeType.point:
+            # must be a station or stop node
+            self.nodes[fid] = (fid, shape_wkb, props)
+
+        elif (is_station_or_line(fid, None, props) and
+              shape_type != ShapeType.point):
+            # must be a station polygon or stop line
+            self.ways[fid] = (fid, shape_wkb, props)
+
+    def add_way(self, way_id, nodes, tags):
+        for node_id in nodes:
+            if node_id in self.nodes:
+                self._ways_using_node[node_id] = way_id
+                assert way_id in self.ways
+
+    def add_relation(self, rel_id, way_off, rel_off, parts, members, tags):
+        r = Relation(rel_id, way_off, rel_off, parts, members, tags)
+        if mz_is_interesting_transit_relation(r.tags):
+            self.relations[r.id] = r
+            for node_id in r.node_ids:
+                if node_id in self.nodes:
+                    self._relations_using_node[node_id].append(rel_id)
+            for way_id in r.way_ids:
+                if way_id in self.ways:
+                    self._relations_using_way[way_id].append(rel_id)
+            for member_rel_id in r.rel_ids:
+                self._relations_using_rel[member_rel_id].append(rel_id)
+
     def relations_using_node(self, node_id):
         "Returns a list of relation IDs which contain the node with that ID."
 
-        return []
+        return self._relations_using_node.get(node_id, [])
 
     def relations_using_way(self, way_id):
         "Returns a list of relation IDs which contain the way with that ID."
 
-        return []
+        return self._relations_using_way.get(way_id, [])
 
     def relations_using_rel(self, rel_id):
         """
@@ -40,17 +120,17 @@ def relations_using_rel(self, rel_id):
         ID.
         """
 
-        return []
+        return self._relations_using_rel.get(rel_id, [])
 
     def ways_using_node(self, node_id):
         "Returns a list of way IDs which contain the node with that ID."
 
-        return []
+        return self._ways_using_node.get(node_id, [])
 
     def relation(self, rel_id):
         "Returns the Relation object with the given ID."
 
-        return None
+        return self.relations[rel_id]
 
     def way(self, way_id):
         """
@@ -58,7 +138,7 @@ def way(self, way_id):
         given way.
         """
 
-        return None
+        return self.ways[way_id]
 
     def node(self, node_id):
         """
@@ -66,12 +146,12 @@ def node(self, node_id):
         given node.
         """
 
-        return None
+        return self.nodes[node_id]
 
     def transit_relations(self, rel_id):
         "Return transit relations containing the relation with the given ID."
 
-        return set()
+        return set(self.relations_using_rel(rel_id))
 
 
 class DataFetcher(object):
@@ -93,6 +173,8 @@ def __init__(self, layers, tables, tile_pyramid):
         tile = self.tile_pyramid.tile()
         max_zoom = self.tile_pyramid.max_z
 
+        table_indexes = defaultdict(list)
+
         for layer_name, info in self.layers.items():
             meta = None
 
@@ -102,15 +184,22 @@ def min_zoom(fid, shape, props):
             layer_index = FeatureTileIndex(tile, max_zoom, min_zoom)
 
             for shape_type in ('point', 'line', 'polygon'):
-                if not info.allows_shape_type(shape_type):
-                    continue
-
-                source = tables('planet_osm_' + shape_type)
-                index_table(source, 'add_feature', layer_index)
+                if info.allows_shape_type(shape_type):
+                    table_name = 'planet_osm_' + shape_type
+                    table_indexes[table_name].append(layer_index)
 
             self.layer_indexes[layer_name] = layer_index
 
         self.osm = OsmRawrLookup()
+        for fn, typ in (('add_feature', 'point'),
+                        ('add_feature', 'line'),
+                        ('add_feature', 'polygon'),
+                        ('add_way', 'ways'),
+                        ('add_relation', 'rels')):
+            table_name = 'planet_osm_' + typ
+            source = tables(table_name)
+            extra_indexes = table_indexes[table_name]
+            index_table(source, fn, self.osm, *extra_indexes)
 
     def _lookup(self, zoom, unpadded_bounds, layer_name):
         from tilequeue.tile import mercator_point_to_coord

From 36fd4928ac060b719fae848d7878f4e54258eb77 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Mon, 25 Sep 2017 17:17:49 +0100
Subject: [PATCH 06/19] Add test for label placement generation for the fixture
 data fetcher.

---
 tests/test_query_fixture.py | 53 ++++++++++++++++++++++++++++++++++---
 tilequeue/query/fixture.py  |  4 ++-
 2 files changed, 53 insertions(+), 4 deletions(-)

diff --git a/tests/test_query_fixture.py b/tests/test_query_fixture.py
index 28c7e8f3..3d9d4237 100644
--- a/tests/test_query_fixture.py
+++ b/tests/test_query_fixture.py
@@ -1,14 +1,19 @@
 import unittest
 
 
-class TestQueryFixture(unittest.TestCase):
+class FixtureTestCase(unittest.TestCase):
 
     def _make(self, rows, min_zoom_fn, props_fn, relations=[],
-              layer_name='testlayer'):
+              layer_name='testlayer', label_placement_layers={}):
         from tilequeue.query.common import LayerInfo
         from tilequeue.query.fixture import make_fixture_data_fetcher
         layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
-        return make_fixture_data_fetcher(layers, rows, relations=relations)
+        return make_fixture_data_fetcher(
+            layers, rows, label_placement_layers=label_placement_layers,
+            relations=relations)
+
+
+class TestQueryFixture(FixtureTestCase):
 
     def test_query_simple(self):
         from shapely.geometry import Point
@@ -141,3 +146,45 @@ def _rel(id, nodes=None, ways=None, rels=None):
             _rel(4, rels=[3]),
             _rel(5, rels=[2, 4]),
         ], 5)
+
+
+class TestLabelPlacement(FixtureTestCase):
+
+    def _test(self, layer_name, props):
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+        from shapely.geometry import box
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return 0
+
+        tile = Coordinate(zoom=15, column=0, row=0)
+        bounds = coord_to_mercator_bounds(tile)
+        shape = box(*bounds)
+
+        rows = [
+            (1, shape, props),
+        ]
+
+        label_placement_layers = {
+            'polygon': set([layer_name]),
+        }
+        fetch = self._make(
+            rows, min_zoom_fn, None, relations=[], layer_name=layer_name,
+            label_placement_layers=label_placement_layers)
+
+        read_rows = fetch(16, bounds)
+        return read_rows
+
+    def test_named_item(self):
+        from shapely import wkb
+
+        layer_name = 'testlayer'
+        read_rows = self._test(layer_name, {'name': 'Foo'})
+
+        self.assertEquals(1, len(read_rows))
+
+        label_prop = '__label__'
+        self.assertTrue(label_prop in read_rows[0])
+        point = wkb.loads(read_rows[0][label_prop])
+        self.assertEqual(point.geom_type, 'Point')
diff --git a/tilequeue/query/fixture.py b/tilequeue/query/fixture.py
index bed3f685..0cc77c8d 100644
--- a/tilequeue/query/fixture.py
+++ b/tilequeue/query/fixture.py
@@ -121,7 +121,9 @@ def __init__(self, layers, rows, rels, label_placement_layers):
         """
         Expect layers to be a dict of layer name to LayerInfo. Expect rows to
         be a list of (fid, shape, properties). Label placement layers should
-        be a set of layer names for which to generate label placement points.
+        be a dict of geometry type ('point', 'linestring', 'polygon') to set
+        of layer names, meaning that each feature of the given type in any of
+        the named layers should additionally get a generated label point.
         """
 
         self.layers = layers

From 54236ee2ef280309c1686564fe705bb0e098800a Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Mon, 25 Sep 2017 18:16:27 +0100
Subject: [PATCH 07/19] Added water layer test to fixture data fetcher.

---
 tests/test_query_fixture.py | 69 +++++++++++++++++++++++++++++++++++++
 1 file changed, 69 insertions(+)

diff --git a/tests/test_query_fixture.py b/tests/test_query_fixture.py
index 3d9d4237..3977d305 100644
--- a/tests/test_query_fixture.py
+++ b/tests/test_query_fixture.py
@@ -188,3 +188,72 @@ def test_named_item(self):
         self.assertTrue(label_prop in read_rows[0])
         point = wkb.loads(read_rows[0][label_prop])
         self.assertEqual(point.geom_type, 'Point')
+
+
+class TestGeometryClipping(FixtureTestCase):
+
+    def _test(self, layer_name, bounds, factor):
+        from shapely.geometry import box
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return 0
+
+        boxwidth = bounds[2] - bounds[0]
+        boxheight = bounds[3] - bounds[1]
+        # make shape overlap the edges of the bounds. that way we can check to
+        # see if the shape gets clipped.
+        shape = box(bounds[0] - factor * boxwidth,
+                    bounds[1] - factor * boxheight,
+                    bounds[2] + factor * boxwidth,
+                    bounds[3] + factor * boxheight)
+
+        props = {'name': 'Foo'}
+
+        rows = [
+            (1, shape, props),
+        ]
+
+        fetch = self._make(
+            rows, min_zoom_fn, None, relations=[], layer_name=layer_name)
+
+        read_rows = fetch(16, bounds)
+        self.assertEqual(1, len(read_rows))
+        return read_rows[0]
+
+    def test_normal_layer(self):
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+        from shapely import wkb
+
+        tile = Coordinate(zoom=15, column=10, row=10)
+        bounds = coord_to_mercator_bounds(tile)
+
+        read_row = self._test('testlayer', bounds, 1.0)
+        clipped_shape = wkb.loads(read_row['__geometry__'])
+        # for normal layers, clipped shape is inside the bounds of the tile.
+        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        self.assertAlmostEqual(1.0, x_factor)
+        self.assertAlmostEqual(1.0, y_factor)
+
+    def test_water_layer(self):
+        # water layer should be expanded by 10% on each side.
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+        from shapely import wkb
+
+        tile = Coordinate(zoom=15, column=10, row=10)
+        bounds = coord_to_mercator_bounds(tile)
+
+        read_row = self._test('water', bounds, 1.0)
+        clipped_shape = wkb.loads(read_row['__geometry__'])
+        # for water layer, the geometry should be 10% larger than the tile
+        # bounds.
+        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        self.assertAlmostEqual(1.1, x_factor)
+        self.assertAlmostEqual(1.1, y_factor)

From 5a943add6d58664562629718c8c5f79fe98ac6ce Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 26 Sep 2017 15:20:05 +0100
Subject: [PATCH 08/19] Add test for name handling between pois, landuse and
 buildings layer features.

---
 tests/test_query_fixture.py | 59 +++++++++++++++++++++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/tests/test_query_fixture.py b/tests/test_query_fixture.py
index 3977d305..767d6d7f 100644
--- a/tests/test_query_fixture.py
+++ b/tests/test_query_fixture.py
@@ -257,3 +257,62 @@ def test_water_layer(self):
                     (bounds[2] - bounds[0]))
         self.assertAlmostEqual(1.1, x_factor)
         self.assertAlmostEqual(1.1, y_factor)
+
+
+class TestNameHandling(FixtureTestCase):
+
+    def _test(self, input_layer_names, expected_layer_names):
+        from shapely.geometry import Point
+        from tilequeue.query.common import LayerInfo
+        from tilequeue.query.fixture import make_fixture_data_fetcher
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return 0
+
+        def props_fn(shape, props, fid, meta):
+            return {}
+
+        shape = Point(0, 0)
+        props = {'name': 'Foo'}
+
+        rows = [
+            (1, shape, props),
+        ]
+
+        layers = {}
+        for name in input_layer_names:
+            layers[name] = LayerInfo(min_zoom_fn, props_fn)
+        fetch = make_fixture_data_fetcher(layers, rows)
+
+        feature_coord = mercator_point_to_coord(16, shape.x, shape.y)
+        read_rows = fetch(16, coord_to_mercator_bounds(feature_coord))
+        self.assertEqual(1, len(read_rows))
+
+        all_names = set(expected_layer_names) | set(input_layer_names)
+        for name in all_names:
+            properties_name = '__%s_properties__' % name
+            self.assertTrue(properties_name in read_rows[0])
+            actual_name = read_rows[0][properties_name].get('name')
+            if name in expected_layer_names:
+                expected_name = props.get('name')
+                self.assertEquals(expected_name, actual_name)
+            else:
+                # check the name doesn't appear anywhere else
+                self.assertEquals(None, actual_name)
+
+    def test_name_single_layer(self):
+        # in any oone of the pois, landuse or buildings layers, a name
+        # by itself will be output in the same layer.
+        for layer_name in ('pois', 'landuse', 'buildings'):
+            self._test([layer_name], [layer_name])
+
+    def test_precedence(self):
+        # if the feature is in the pois layer, then that should get the name
+        # and the other layers should not.
+        self._test(['pois', 'landuse'], ['pois'])
+        self._test(['pois', 'buildings'], ['pois'])
+        self._test(['pois', 'landuse', 'buildings'], ['pois'])
+        # otherwise, landuse should take precedence over buildings.
+        self._test(['landuse', 'buildings'], ['landuse'])

From 7c6b3d746f48f2341724b332e934b2942d3f4ee0 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 26 Sep 2017 15:50:11 +0100
Subject: [PATCH 09/19] Add label placement handling to RAWR tile generator.

---
 tests/test_query_rawr.py | 62 ++++++++++++++++++++++++++++++++++++++--
 tilequeue/query/rawr.py  | 18 ++++++++++--
 2 files changed, 74 insertions(+), 6 deletions(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index 699d3a6e..db52cafd 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -10,15 +10,20 @@ def __call__(self, table_name):
         return self.tables.get(table_name, [])
 
 
-class TestQueryRawr(unittest.TestCase):
+class RawrTestCase(unittest.TestCase):
 
     def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
-              layer_name='testlayer'):
+              layer_name='testlayer', label_placement_layers={}):
         from tilequeue.query.common import LayerInfo
         from tilequeue.query.rawr import make_rawr_data_fetcher
 
         layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
-        return make_rawr_data_fetcher(layers, tables, tile_pyramid)
+        return make_rawr_data_fetcher(
+            layers, tables, tile_pyramid,
+            label_placement_layers=label_placement_layers)
+
+
+class TestQueryRawr(RawrTestCase):
 
     def test_query_simple(self):
         from shapely.geometry import Point
@@ -194,3 +199,54 @@ def _rel(id, nodes=None, ways=None, rels=None):
             _rel(4, rels=[3]),
             _rel(5, rels=[2, 4]),
         ], 5)
+
+
+class TestLabelPlacement(RawrTestCase):
+
+    def _test(self, layer_name, props):
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+        from shapely.geometry import box
+        from tilequeue.query.rawr import TilePyramid
+
+        top_zoom = 10
+        max_zoom = top_zoom + 6
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return top_zoom
+
+        tile = Coordinate(zoom=15, column=0, row=0)
+        top_tile = tile.zoomTo(top_zoom).container()
+        tile_pyramid = TilePyramid(
+            top_zoom, top_tile.column, top_tile.row, max_zoom)
+
+        bounds = coord_to_mercator_bounds(tile)
+        shape = box(*bounds)
+        tables = TestGetTable({
+            'planet_osm_polygon': [
+                (1, shape.wkb, props),
+            ]
+        })
+
+        label_placement_layers = {
+            'polygon': set([layer_name]),
+        }
+        fetch = self._make(
+            min_zoom_fn, None, tables, tile_pyramid, layer_name=layer_name,
+            label_placement_layers=label_placement_layers)
+
+        read_rows = fetch(tile.zoom, bounds)
+        return read_rows
+
+    def test_named_item(self):
+        from shapely import wkb
+
+        layer_name = 'testlayer'
+        read_rows = self._test(layer_name, {'name': 'Foo'})
+
+        self.assertEquals(1, len(read_rows))
+
+        label_prop = '__label__'
+        self.assertTrue(label_prop in read_rows[0])
+        point = wkb.loads(read_rows[0][label_prop])
+        self.assertEqual(point.geom_type, 'Point')
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 7c2867a4..dd3adf5c 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -5,6 +5,7 @@
 from tilequeue.query.common import is_station_or_line
 from tilequeue.query.common import deassoc
 from tilequeue.query.common import mz_is_interesting_transit_relation
+from tilequeue.query.common import shape_type_lookup
 
 
 class Relation(object):
@@ -156,7 +157,7 @@ def transit_relations(self, rel_id):
 
 class DataFetcher(object):
 
-    def __init__(self, layers, tables, tile_pyramid):
+    def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
         """
         Expect layers to be a dict of layer name to LayerInfo (see fixture.py).
         Tables should be a callable which returns a generator over the rows in
@@ -168,6 +169,7 @@ def __init__(self, layers, tables, tile_pyramid):
 
         self.layers = layers
         self.tile_pyramid = tile_pyramid
+        self.label_placement_layers = label_placement_layers
         self.layer_indexes = {}
 
         tile = self.tile_pyramid.tile()
@@ -255,6 +257,15 @@ def __call__(self, zoom, unpadded_bounds):
                 read_row['__' + layer_name + '_properties__'] = layer_props
                 read_row['__id__'] = fid
                 read_row['__geometry__'] = bytes(shape.wkb)
+
+                # if the feature exists in any label placement layer, then we
+                # should consider generating a centroid
+                label_layers = self.label_placement_layers.get(
+                    shape_type_lookup(shape), {})
+                if layer_name in label_layers:
+                    read_row['__label__'] = bytes(
+                        shape.representative_point().wkb)
+
                 read_rows.append(read_row)
 
         return read_rows
@@ -262,5 +273,6 @@ def __call__(self, zoom, unpadded_bounds):
 
 # tables is a callable which should return a generator over the rows of the
 # table when called with the table name.
-def make_rawr_data_fetcher(layers, tables, tile_pyramid):
-    return DataFetcher(layers, tables, tile_pyramid)
+def make_rawr_data_fetcher(layers, tables, tile_pyramid,
+                           label_placement_layers={}):
+    return DataFetcher(layers, tables, tile_pyramid, label_placement_layers)

From 60694de4984f744b921370b2717a3a2c27ddbcb1 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 26 Sep 2017 18:31:32 +0100
Subject: [PATCH 10/19] Clip features to query bbox, except water features,
 which are clipped to an expanded box.

---
 tests/test_query_rawr.py | 81 ++++++++++++++++++++++++++++++++++++++++
 tilequeue/query/rawr.py  | 20 +++++++++-
 2 files changed, 99 insertions(+), 2 deletions(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index db52cafd..270044ff 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -250,3 +250,84 @@ def test_named_item(self):
         self.assertTrue(label_prop in read_rows[0])
         point = wkb.loads(read_rows[0][label_prop])
         self.assertEqual(point.geom_type, 'Point')
+
+
+class TestGeometryClipping(RawrTestCase):
+
+    def _test(self, layer_name, tile, factor):
+        from shapely.geometry import box
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+
+        top_zoom = 10
+        max_zoom = top_zoom + 6
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return top_zoom
+
+        top_tile = tile.zoomTo(top_zoom).container()
+        tile_pyramid = TilePyramid(
+            top_zoom, top_tile.column, top_tile.row, max_zoom)
+
+        bounds = coord_to_mercator_bounds(tile)
+        boxwidth = bounds[2] - bounds[0]
+        boxheight = bounds[3] - bounds[1]
+        # make shape overlap the edges of the bounds. that way we can check to
+        # see if the shape gets clipped.
+        shape = box(bounds[0] - factor * boxwidth,
+                    bounds[1] - factor * boxheight,
+                    bounds[2] + factor * boxwidth,
+                    bounds[3] + factor * boxheight)
+
+        props = {'name': 'Foo'}
+
+        tables = TestGetTable({
+            'planet_osm_polygon': [
+                (1, shape.wkb, props),
+            ],
+        })
+
+        fetch = self._make(
+            min_zoom_fn, None, tables, tile_pyramid, layer_name=layer_name)
+
+        read_rows = fetch(tile.zoom, bounds)
+        self.assertEqual(1, len(read_rows))
+        return read_rows[0]
+
+    def test_normal_layer(self):
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+        from shapely import wkb
+
+        tile = Coordinate(zoom=15, column=10, row=10)
+        bounds = coord_to_mercator_bounds(tile)
+
+        read_row = self._test('testlayer', tile, 1.0)
+        clipped_shape = wkb.loads(read_row['__geometry__'])
+        # for normal layers, clipped shape is inside the bounds of the tile.
+        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        self.assertAlmostEqual(1.0, x_factor)
+        self.assertAlmostEqual(1.0, y_factor)
+
+    def test_water_layer(self):
+        # water layer should be expanded by 10% on each side.
+        from ModestMaps.Core import Coordinate
+        from tilequeue.tile import coord_to_mercator_bounds
+        from shapely import wkb
+
+        tile = Coordinate(zoom=15, column=10, row=10)
+        bounds = coord_to_mercator_bounds(tile)
+
+        read_row = self._test('water', tile, 1.0)
+        clipped_shape = wkb.loads(read_row['__geometry__'])
+        # for water layer, the geometry should be 10% larger than the tile
+        # bounds.
+        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
+                    (bounds[2] - bounds[0]))
+        self.assertAlmostEqual(1.1, x_factor)
+        self.assertAlmostEqual(1.1, y_factor)
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index dd3adf5c..26eb8995 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -6,6 +6,7 @@
 from tilequeue.query.common import deassoc
 from tilequeue.query.common import mz_is_interesting_transit_relation
 from tilequeue.query.common import shape_type_lookup
+from tilequeue.transform import calculate_padded_bounds
 
 
 class Relation(object):
@@ -221,10 +222,16 @@ def _lookup(self, zoom, unpadded_bounds, layer_name):
         bottomright.row = max(bottomright.row, topleft.row)
 
         features = []
+        seen_ids = set()
         for x in range(int(topleft.column), int(bottomright.column) + 1):
             for y in range(int(topleft.row), int(bottomright.row) + 1):
                 tile = Tile(zoom, x, y)
-                features.extend(index(tile))
+                tile_features = index(tile)
+                for feature in tile_features:
+                    feature_id = id(feature)
+                    if feature_id not in seen_ids:
+                        seen_ids.add(feature_id)
+                        features.append(feature)
         return features
 
     def __call__(self, zoom, unpadded_bounds):
@@ -256,7 +263,16 @@ def __call__(self, zoom, unpadded_bounds):
 
                 read_row['__' + layer_name + '_properties__'] = layer_props
                 read_row['__id__'] = fid
-                read_row['__geometry__'] = bytes(shape.wkb)
+
+                # if this is a water layer feature, then clip to an expanded
+                # bounding box to avoid tile-edge artefacts.
+                clip_box = bbox
+                if layer_name == 'water':
+                    pad_factor = 1.1
+                    clip_box = calculate_padded_bounds(
+                        pad_factor, unpadded_bounds)
+                clip_shape = clip_box.intersection(shape)
+                read_row['__geometry__'] = bytes(clip_shape.wkb)
 
                 # if the feature exists in any label placement layer, then we
                 # should consider generating a centroid

From a50334fb32645864052af14e40162d7d7a24f9a4 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Thu, 28 Sep 2017 17:25:27 +0100
Subject: [PATCH 11/19] Update RAWR tile indexing for change in interface.

---
 tilequeue/query/rawr.py | 42 +++++++++++++++++++++++++++++++++++------
 1 file changed, 36 insertions(+), 6 deletions(-)

diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 26eb8995..233745dd 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -60,6 +60,15 @@ def _wkb_shape(wkb):
         assert False, "WKB shape type %d not understood." % (typ,)
 
 
+def _match_type(values, types):
+    if len(values) != len(types):
+        return False
+    for val, typ in zip(values, types):
+        if not isinstance(val, typ):
+            return False
+    return True
+
+
 class OsmRawrLookup(object):
 
     def __init__(self):
@@ -72,6 +81,28 @@ def __init__(self):
         self._relations_using_way = defaultdict(list)
         self._relations_using_rel = defaultdict(list)
 
+    def add_row(self, *args):
+        # there's only a single dispatch from the indexing function, which
+        # passes row data from the table. we have to figure out here what
+        # kind of row it was, and send the data on to the right function.
+
+        # IDs can be either ints or longs, and generally we don't care which,
+        # so we accept either as the type for that position in the function.
+        num = (int, long)
+
+        if _match_type(args, (num, (str, bytes), dict)):
+            self.add_feature(*args)
+
+        elif _match_type(args, (num, list, dict)):
+            self.add_way(*args)
+
+        elif _match_type(args, (num, num, num, list, list, list)):
+            self.add_relation(*args)
+
+        else:
+            raise Exception("Unknown row shape for OsmRawrLookup.add_row: %s" %
+                            (repr(map(type, args)),))
+
     def add_feature(self, fid, shape_wkb, props):
         if fid < 0:
             return
@@ -194,15 +225,14 @@ def min_zoom(fid, shape, props):
             self.layer_indexes[layer_name] = layer_index
 
         self.osm = OsmRawrLookup()
-        for fn, typ in (('add_feature', 'point'),
-                        ('add_feature', 'line'),
-                        ('add_feature', 'polygon'),
-                        ('add_way', 'ways'),
-                        ('add_relation', 'rels')):
+        # NOTE: order here is different from that in raw_tiles index()
+        # function. this is because here we want to gather up some
+        # "interesting" feature IDs before we look at the ways/rels tables.
+        for typ in ('point', 'line', 'polygon', 'ways', 'rels'):
             table_name = 'planet_osm_' + typ
             source = tables(table_name)
             extra_indexes = table_indexes[table_name]
-            index_table(source, fn, self.osm, *extra_indexes)
+            index_table(source, self.osm, *extra_indexes)
 
     def _lookup(self, zoom, unpadded_bounds, layer_name):
         from tilequeue.tile import mercator_point_to_coord

From 3299ea88d3cde8606a916208e402560a0852089a Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Thu, 28 Sep 2017 18:10:06 +0100
Subject: [PATCH 12/19] Add test for name precedence for RAWR features.

---
 tests/test_query_rawr.py | 81 ++++++++++++++++++++++++++++++++++++++++
 tilequeue/query/rawr.py  | 21 +++++++++++
 2 files changed, 102 insertions(+)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index 270044ff..eb592ca4 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -331,3 +331,84 @@ def test_water_layer(self):
                     (bounds[2] - bounds[0]))
         self.assertAlmostEqual(1.1, x_factor)
         self.assertAlmostEqual(1.1, y_factor)
+
+
+class TestNameHandling(RawrTestCase):
+
+    def _test(self, input_layer_names, expected_layer_names):
+        from shapely.geometry import Point
+        from tilequeue.query.common import LayerInfo
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.query.rawr import make_rawr_data_fetcher
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        top_zoom = 10
+        max_zoom = top_zoom + 6
+
+        def min_zoom_fn(shape, props, fid, meta):
+            return top_zoom
+
+        def props_fn(shape, props, fid, meta):
+            return {}
+
+        shape = Point(0, 0)
+        props = {'name': 'Foo'}
+
+        tables = TestGetTable({
+            'planet_osm_point': [
+                (1, shape.wkb, props),
+            ],
+        })
+
+        tile = mercator_point_to_coord(16, shape.x, shape.y)
+        top_tile = tile.zoomTo(top_zoom).container()
+        tile_pyramid = TilePyramid(
+            top_zoom, top_tile.column, top_tile.row, max_zoom)
+
+        layers = {}
+        for name in input_layer_names:
+            layers[name] = LayerInfo(min_zoom_fn, props_fn)
+        fetch = make_rawr_data_fetcher(layers, tables, tile_pyramid)
+
+        read_rows = fetch(tile.zoom, coord_to_mercator_bounds(tile))
+        # the RAWR query goes over features multiple times because of the
+        # indexing, so we can't rely on all the properties for one feature to
+        # be all together in the same place. this loops over all the features,
+        # checking that there's only really one of them and gathering together
+        # all the __%s_properties__ from all the rows for further testing.
+        all_props = {}
+        for row in read_rows:
+            self.assertEquals(1, row['__id__'])
+            self.assertEquals(shape.wkb, row['__geometry__'])
+            for key, val in row.items():
+                if key.endswith('_properties__'):
+                    self.assertFalse(key in all_props)
+                    all_props[key] = val
+
+        all_names = set(expected_layer_names) | set(input_layer_names)
+        for name in all_names:
+            properties_name = '__%s_properties__' % name
+            self.assertTrue(properties_name in all_props)
+            actual_name = all_props[properties_name].get('name')
+            if name in expected_layer_names:
+                expected_name = props.get('name')
+                self.assertEquals(expected_name, actual_name)
+            else:
+                # check the name doesn't appear anywhere else
+                self.assertEquals(None, actual_name)
+
+    def test_name_single_layer(self):
+        # in any oone of the pois, landuse or buildings layers, a name
+        # by itself will be output in the same layer.
+        for layer_name in ('pois', 'landuse', 'buildings'):
+            self._test([layer_name], [layer_name])
+
+    def test_precedence(self):
+        # if the feature is in the pois layer, then that should get the name
+        # and the other layers should not.
+        self._test(['pois', 'landuse'], ['pois'])
+        self._test(['pois', 'buildings'], ['pois'])
+        self._test(['pois', 'landuse', 'buildings'], ['pois'])
+        # otherwise, landuse should take precedence over buildings.
+        self._test(['landuse', 'buildings'], ['landuse'])
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 233745dd..8f25122d 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -210,6 +210,7 @@ def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
         table_indexes = defaultdict(list)
 
         for layer_name, info in self.layers.items():
+            # TODO! this shouldn't be none!
             meta = None
 
             def min_zoom(fid, shape, props):
@@ -234,6 +235,20 @@ def min_zoom(fid, shape, props):
             extra_indexes = table_indexes[table_name]
             index_table(source, self.osm, *extra_indexes)
 
+    def _named_layer(self, fid, shape, props):
+        # we want only one layer from ('pois', 'landuse', 'buildings') for
+        # each feature to be assigned a name. therefore, we use the presence
+        # or absence of a min zoom to check whether these features as in these
+        # layers, and therefore which should be assigned the name.
+        for layer_name in ('pois', 'landuse', 'buildings'):
+            info = self.layers.get(layer_name)
+            if info and info.min_zoom_fn:
+                # TODO! meta should not be None!
+                min_zoom = info.min_zoom_fn(shape, props, fid, None)
+                if min_zoom is not None:
+                    return layer_name
+        return None
+
     def _lookup(self, zoom, unpadded_bounds, layer_name):
         from tilequeue.tile import mercator_point_to_coord
         from raw_tiles.tile import Tile
@@ -291,6 +306,12 @@ def __call__(self, zoom, unpadded_bounds):
                 layer_props = layer_properties(
                     fid, shape, props, layer_name, zoom, self.osm)
 
+                # add name into whichever of the pois, landuse or buildings
+                # layers has claimed this feature.
+                name = props.get('name', None)
+                if name and self._named_layer(fid, shape, props) == layer_name:
+                    layer_props['name'] = name
+
                 read_row['__' + layer_name + '_properties__'] = layer_props
                 read_row['__id__'] = fid
 

From 33ec4936038f3dcfdd23f33f9f7525311511aa5c Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Fri, 29 Sep 2017 13:11:10 +0100
Subject: [PATCH 13/19] Major refactor to swap layer/feature loop order.

Previously, the code looped over layers and looked up features in a per-layer index. This caused a problem for some types of processing. For example, for a feature which might appear in `pois`, `landuse` and `buildings` layers, that only one of those features is named. This is considerably harder to do when the feature is extracted from several different indexes. Another issue was that the properties were copied for each feature and in each layer to add the `min_zoom` property, which meant that there was a lot more memory usage and GC pressure. Since the size of RAWR tiles is a concern, it would be better to keep the memory usage down where possible.

The new code has a single index, and keeps `min_zoom` for each layer in a separate `dict`. This means that the code which iterates over the layers has to do a second check for the zoom level, and the feature is included at the minimum zoom level it appears in any layer.

A side effect of this is that each `_Feature` object is unique, and so we can use `id()` to check equality, which helps when de-duplicating across tiles.
---
 tilequeue/query/rawr.py | 230 +++++++++++++++++++++++++++++-----------
 1 file changed, 167 insertions(+), 63 deletions(-)

diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 8f25122d..72e5630b 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -1,5 +1,6 @@
 from collections import namedtuple, defaultdict
 from shapely.geometry import box
+from shapely.wkb import loads as wkb_loads
 from tilequeue.query.common import layer_properties
 from tilequeue.query.common import is_station_or_stop
 from tilequeue.query.common import is_station_or_line
@@ -7,6 +8,8 @@
 from tilequeue.query.common import mz_is_interesting_transit_relation
 from tilequeue.query.common import shape_type_lookup
 from tilequeue.transform import calculate_padded_bounds
+from raw_tiles.tile import shape_tile_coverage
+from math import floor
 
 
 class Relation(object):
@@ -43,6 +46,12 @@ class ShapeType(object):
     line = 2
     polygon = 3
 
+    _LOOKUP = ['point', 'line', 'polygon']
+
+    @staticmethod
+    def lookup(typ):
+        return ShapeType._LOOKUP[typ-1]
+
 
 def _wkb_shape(wkb):
     reverse = ord(wkb[0]) == 1
@@ -187,6 +196,117 @@ def transit_relations(self, rel_id):
         return set(self.relations_using_rel(rel_id))
 
 
+def _tiles(zoom, unpadded_bounds):
+    from tilequeue.tile import mercator_point_to_coord
+    from raw_tiles.tile import Tile
+
+    minx, miny, maxx, maxy = unpadded_bounds
+    topleft = mercator_point_to_coord(zoom, minx, miny)
+    bottomright = mercator_point_to_coord(zoom, maxx, maxy)
+
+    # make sure that the bottom right coordinate is below and to the right
+    # of the top left coordinate. it can happen that the coordinates are
+    # mixed up due to small numerical precision artefacts being enlarged
+    # by the conversion to integer and y-coordinate flip.
+    assert topleft.zoom == bottomright.zoom
+    bottomright.column = max(bottomright.column, topleft.column)
+    bottomright.row = max(bottomright.row, topleft.row)
+
+    for x in range(int(topleft.column), int(bottomright.column) + 1):
+        for y in range(int(topleft.row), int(bottomright.row) + 1):
+            tile = Tile(zoom, x, y)
+            yield tile
+
+
+_Feature = namedtuple('_Feature', 'fid shape properties layer_min_zooms')
+
+
+class _LazyShape(object):
+    """
+    This proxy exists so that we can avoid parsing the WKB for a shape unless
+    it is actually needed. Parsing WKB is pretty fast, but multiplied over
+    many thousands of objects, it can become the slowest part of the indexing
+    process. Given that we reject many features on the basis of their
+    properties alone, lazily parsing the WKB can provide a significant saving.
+    """
+
+    def __init__(self, wkb):
+        self.wkb = wkb
+        self.obj = None
+
+    def __getattr__(self, name):
+        if self.obj is None:
+            self.obj = wkb_loads(self.wkb)
+        return getattr(self.obj, name)
+
+
+class _LayersIndex(object):
+
+    def __init__(self, layers, tile_pyramid):
+        self.layers = layers
+        self.tile_pyramid = tile_pyramid
+        self.tile_index = defaultdict(list)
+
+    def add_row(self, fid, shape_wkb, props):
+        shape = _LazyShape(shape_wkb)
+        # TODO! meta is not None!
+        meta = None
+
+        layer_min_zooms = {}
+        for layer_name, info in self.layers.items():
+            shape_type = _wkb_shape(shape_wkb)
+            shape_type_str = ShapeType.lookup(shape_type)
+            if info.shape_types and shape_type_str not in info.shape_types:
+                continue
+            min_zoom = info.min_zoom_fn(fid, shape, props, meta)
+            if min_zoom is not None:
+                layer_min_zooms[layer_name] = min_zoom
+
+        # quick exit if the feature didn't have a min zoom in any layer.
+        if not layer_min_zooms:
+            return
+
+        # lowest zoom that this feature appears in any layer. note that this
+        # is clamped to the max zoom, so that all features that appear at some
+        # zoom level appear at the max zoom. this is different from the min
+        # zoom in layer_min_zooms, which is a property that will be injected
+        # for each layer and is used by the _client_ to determine feature
+        # visibility.
+        min_zoom = min(self.tile_pyramid.max_z, min(layer_min_zooms.values()))
+
+        # single object (hence single id()) will be shared amongst all layers.
+        # this allows us to easily and quickly de-duplicate at later layers in
+        # the stack.
+        feature = _Feature(fid, shape, props, layer_min_zooms)
+
+        # take the minimum integer zoom - this is the min zoom tile that the
+        # feature should appear in, and a feature with min_zoom = 1.9 should
+        # appear in a tile at z=1, not 2, since the tile at z=N is used for
+        # the zoom range N to N+1.
+        #
+        # we cut this off at this index's min zoom, as we aren't interested
+        # in any tiles outside of that.
+        floor_zoom = max(self.tile_pyramid.z, int(floor(min_zoom)))
+
+        # seed initial set of tiles at maximum zoom. all features appear at
+        # least at the max zoom, even if the min_zoom function returns a
+        # value larger than the max zoom.
+        zoom = self.tile_pyramid.max_z
+        tiles = shape_tile_coverage(shape, zoom, self.tile_pyramid.tile())
+
+        while zoom >= floor_zoom:
+            parent_tiles = set()
+            for tile in tiles:
+                self.tile_index[tile].append(feature)
+                parent_tiles.add(tile.parent())
+
+            zoom -= 1
+            tiles = parent_tiles
+
+    def __call__(self, tile):
+        return self.tile_index.get(tile, [])
+
+
 class DataFetcher(object):
 
     def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
@@ -196,7 +316,6 @@ def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
         the table when called with that table's name.
         """
 
-        from raw_tiles.index.features import FeatureTileIndex
         from raw_tiles.index.index import index_table
 
         self.layers = layers
@@ -204,26 +323,12 @@ def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
         self.label_placement_layers = label_placement_layers
         self.layer_indexes = {}
 
-        tile = self.tile_pyramid.tile()
-        max_zoom = self.tile_pyramid.max_z
-
         table_indexes = defaultdict(list)
 
-        for layer_name, info in self.layers.items():
-            # TODO! this shouldn't be none!
-            meta = None
-
-            def min_zoom(fid, shape, props):
-                return info.min_zoom_fn(fid, shape, props, meta)
-
-            layer_index = FeatureTileIndex(tile, max_zoom, min_zoom)
-
-            for shape_type in ('point', 'line', 'polygon'):
-                if info.allows_shape_type(shape_type):
-                    table_name = 'planet_osm_' + shape_type
-                    table_indexes[table_name].append(layer_index)
-
-            self.layer_indexes[layer_name] = layer_index
+        self.layers_index = _LayersIndex(self.layers, self.tile_pyramid)
+        for shape_type in ('point', 'line', 'polygon'):
+            table_name = 'planet_osm_' + shape_type
+            table_indexes[table_name].append(self.layers_index)
 
         self.osm = OsmRawrLookup()
         # NOTE: order here is different from that in raw_tiles index()
@@ -249,34 +354,18 @@ def _named_layer(self, fid, shape, props):
                     return layer_name
         return None
 
-    def _lookup(self, zoom, unpadded_bounds, layer_name):
-        from tilequeue.tile import mercator_point_to_coord
-        from raw_tiles.tile import Tile
-
-        minx, miny, maxx, maxy = unpadded_bounds
-        topleft = mercator_point_to_coord(zoom, minx, miny)
-        bottomright = mercator_point_to_coord(zoom, maxx, maxy)
-        index = self.layer_indexes[layer_name]
-
-        # make sure that the bottom right coordinate is below and to the right
-        # of the top left coordinate. it can happen that the coordinates are
-        # mixed up due to small numerical precision artefacts being enlarged
-        # by the conversion to integer and y-coordinate flip.
-        assert topleft.zoom == bottomright.zoom
-        bottomright.column = max(bottomright.column, topleft.column)
-        bottomright.row = max(bottomright.row, topleft.row)
-
+    def _lookup(self, zoom, unpadded_bounds):
         features = []
         seen_ids = set()
-        for x in range(int(topleft.column), int(bottomright.column) + 1):
-            for y in range(int(topleft.row), int(bottomright.row) + 1):
-                tile = Tile(zoom, x, y)
-                tile_features = index(tile)
-                for feature in tile_features:
-                    feature_id = id(feature)
-                    if feature_id not in seen_ids:
-                        seen_ids.add(feature_id)
-                        features.append(feature)
+
+        for tile in _tiles(zoom, unpadded_bounds):
+            tile_features = self.layers_index(tile)
+            for feature in tile_features:
+                feature_id = id(feature)
+                if feature_id not in seen_ids:
+                    seen_ids.add(feature_id)
+                    features.append(feature)
+
         return features
 
     def __call__(self, zoom, unpadded_bounds):
@@ -292,27 +381,46 @@ def __call__(self, zoom, unpadded_bounds):
         assert zoom >= self.tile_pyramid.z
         assert bbox.within(self.tile_pyramid.bbox())
 
-        for layer_name, info in self.layers.items():
-
-            for (fid, shape, props) in self._lookup(
-                    zoom, unpadded_bounds, layer_name):
-                # reject any feature which doesn't intersect the given bounds
-                if bbox.disjoint(shape):
+        for (fid, shape, props, layer_min_zooms) in self._lookup(
+                zoom, unpadded_bounds):
+            # reject any feature which doesn't intersect the given bounds
+            if bbox.disjoint(shape):
+                continue
+
+            # place for assembing the read row as if from postgres
+            read_row = {}
+            generate_label_placement = False
+
+            # add name into whichever of the pois, landuse or buildings
+            # layers has claimed this feature.
+            name = props.get('name', None)
+            named_layer = self._named_layer(fid, shape, props)
+
+            for layer_name, min_zoom in layer_min_zooms.items():
+                # we need to keep fractional zooms, e.g: 4.999 should appear
+                # in tiles at zoom level 4, but not 3. also, tiles at zooms
+                # past the max zoom should be clamped to the max zoom.
+                tile_zoom = min(self.tile_pyramid.max_z, floor(min_zoom))
+                if tile_zoom > zoom:
                     continue
 
-                # place for assembing the read row as if from postgres
-                read_row = {}
-
                 layer_props = layer_properties(
                     fid, shape, props, layer_name, zoom, self.osm)
+                layer_props['min_zoom'] = min_zoom
 
-                # add name into whichever of the pois, landuse or buildings
-                # layers has claimed this feature.
-                name = props.get('name', None)
-                if name and self._named_layer(fid, shape, props) == layer_name:
+                if name and named_layer == layer_name:
                     layer_props['name'] = name
 
                 read_row['__' + layer_name + '_properties__'] = layer_props
+
+                # if the feature exists in any label placement layer, then we
+                # should consider generating a centroid
+                label_layers = self.label_placement_layers.get(
+                    shape_type_lookup(shape), {})
+                if layer_name in label_layers:
+                    generate_label_placement = True
+
+            if read_row:
                 read_row['__id__'] = fid
 
                 # if this is a water layer feature, then clip to an expanded
@@ -325,11 +433,7 @@ def __call__(self, zoom, unpadded_bounds):
                 clip_shape = clip_box.intersection(shape)
                 read_row['__geometry__'] = bytes(clip_shape.wkb)
 
-                # if the feature exists in any label placement layer, then we
-                # should consider generating a centroid
-                label_layers = self.label_placement_layers.get(
-                    shape_type_lookup(shape), {})
-                if layer_name in label_layers:
+                if generate_label_placement:
                     read_row['__label__'] = bytes(
                         shape.representative_point().wkb)
 

From 10bb0fb4464103d24d9dd433a3c671ba5e4455f2 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Fri, 29 Sep 2017 13:22:39 +0100
Subject: [PATCH 14/19] Add raw_tiles module to requirements.

---
 requirements.txt | 1 +
 setup.py         | 1 +
 2 files changed, 2 insertions(+)

diff --git a/requirements.txt b/requirements.txt
index f31556f6..f3d9dfb3 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -25,3 +25,4 @@ wsgiref==0.1.2
 zope.dottedname==4.1.0
 edtf==0.9.3
 mapbox-vector-tile==1.2.0
+git+https://github.com/tilezen/raw_tiles@master#egg=raw_tiles
diff --git a/setup.py b/setup.py
index 6087f149..16c4a799 100644
--- a/setup.py
+++ b/setup.py
@@ -44,6 +44,7 @@
           'pyproj',
           'python-dateutil',
           'PyYAML',
+          'raw_tiles',
           'redis',
           'requests',
           'Shapely',

From c44abe73f3f8f226bd6304ba60bf72a66c13d734 Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 3 Oct 2017 10:31:11 +0100
Subject: [PATCH 15/19] Better comments.

---
 tests/test_query_rawr.py | 32 ++++++++++++++++++++++++++++-
 tilequeue/query/rawr.py  | 44 +++++++++++++++++++++++++++++++++++++++-
 2 files changed, 74 insertions(+), 2 deletions(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index eb592ca4..a3f1d830 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -2,6 +2,11 @@
 
 
 class TestGetTable(object):
+    """
+    Mocks the interface expected by raw_tiles.index.index.index_table,
+    which provides "table lookup". Here, we just return static stuff
+    previously set up in the test.
+    """
 
     def __init__(self, tables):
         self.tables = tables
@@ -11,6 +16,10 @@ def __call__(self, table_name):
 
 
 class RawrTestCase(unittest.TestCase):
+    """
+    Base layer of the tests, providing a utility function to create a data
+    fetcher with a set of mocked data.
+    """
 
     def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
               layer_name='testlayer', label_placement_layers={}):
@@ -26,6 +35,10 @@ def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
 class TestQueryRawr(RawrTestCase):
 
     def test_query_simple(self):
+        # just check that we can get back the mock data we put into a tile,
+        # and that the indexing/fetching code respects the tile boundary and
+        # min_zoom function.
+
         from shapely.geometry import Point
         from tilequeue.query.rawr import TilePyramid
         from tilequeue.tile import coord_to_mercator_bounds
@@ -74,6 +87,10 @@ def min_zoom_fn(shape, props, fid, meta):
         self.assertEquals(0, len(read_rows))
 
     def test_query_min_zoom_fraction(self):
+        # test that fractional min zooms are included in their "floor" zoom
+        # tile. this is to allow over-zooming of a zoom N tile until N+1,
+        # where the next zoom tile kicks in.
+
         from shapely.geometry import Point
         from tilequeue.query.rawr import TilePyramid
         from tilequeue.tile import coord_to_mercator_bounds
@@ -105,6 +122,10 @@ def min_zoom_fn(shape, props, fid, meta):
         self.assertEquals(0, len(read_rows))
 
     def test_query_past_max_zoom(self):
+        # check that features with a min_zoom beyond the maximum zoom are still
+        # included at the maximum zoom. since this is the last zoom level we
+        # generate, it must include everything.
+
         from shapely.geometry import Point
         from tilequeue.query.rawr import TilePyramid
         from tilequeue.tile import coord_to_mercator_bounds
@@ -138,6 +159,8 @@ def min_zoom_fn(shape, props, fid, meta):
         self.assertEquals(0, len(read_rows))
 
     def test_root_relation_id(self):
+        # check the logic for finding a root relation ID for station complexes.
+
         from shapely.geometry import Point
         from tilequeue.query.rawr import TilePyramid
         from tilequeue.tile import coord_to_mercator_bounds
@@ -239,6 +262,9 @@ def min_zoom_fn(shape, props, fid, meta):
         return read_rows
 
     def test_named_item(self):
+        # check that a label is generated for features in label placement
+        # layers.
+
         from shapely import wkb
 
         layer_name = 'testlayer'
@@ -295,6 +321,9 @@ def min_zoom_fn(shape, props, fid, meta):
         return read_rows[0]
 
     def test_normal_layer(self):
+        # check that normal layer geometries are clipped to the bounding box of
+        # the tile.
+
         from ModestMaps.Core import Coordinate
         from tilequeue.tile import coord_to_mercator_bounds
         from shapely import wkb
@@ -313,7 +342,8 @@ def test_normal_layer(self):
         self.assertAlmostEqual(1.0, y_factor)
 
     def test_water_layer(self):
-        # water layer should be expanded by 10% on each side.
+        # water layer should be clipped to the tile bounds expanded by 10%.
+
         from ModestMaps.Core import Coordinate
         from tilequeue.tile import coord_to_mercator_bounds
         from shapely import wkb
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 72e5630b..236e9aee 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -13,6 +13,13 @@
 
 
 class Relation(object):
+    """
+    Relation object holds data about a relation and provides a nicer interface
+    than the raw tuple by turning the tags array into a dict, and separating
+    out the "parts" array of IDs into separate lists for nodes, ways and other
+    relations.
+    """
+
     def __init__(self, rel_id, way_off, rel_off, parts, members, tags):
         self.id = rel_id
         self.tags = deassoc(tags)
@@ -22,6 +29,11 @@ def __init__(self, rel_id, way_off, rel_off, parts, members, tags):
 
 
 class TilePyramid(namedtuple('TilePyramid', 'z x y max_z')):
+    """
+    Represents a "tile pyramid" of all tiles which are geographically
+    contained within the tile `z/x/y` up to a maximum zoom of `max_z`. This is
+    the set of tiles corresponding to one RAWR tile.
+    """
 
     def tile(self):
         from raw_tiles.tile import Tile
@@ -40,7 +52,9 @@ def bbox(self):
         return box(*self.bounds())
 
 
-# weak type of enum type
+# weak enum type used to represent shape type, rather than use a string.
+# (or, where we have to use a string, at least use a shared single instance
+# of a string)
 class ShapeType(object):
     point = 1
     line = 2
@@ -50,9 +64,12 @@ class ShapeType(object):
 
     @staticmethod
     def lookup(typ):
+        "turn the enum into a string"
         return ShapeType._LOOKUP[typ-1]
 
 
+# determine the shape type from the raw WKB bytes. this means we don't have to
+# parse the WKB, which can be an expensive operation for large polygons.
 def _wkb_shape(wkb):
     reverse = ord(wkb[0]) == 1
     type_bytes = map(ord, wkb[1:5])
@@ -69,6 +86,9 @@ def _wkb_shape(wkb):
         assert False, "WKB shape type %d not understood." % (typ,)
 
 
+# return true if the tuple of values corresponds to, and each is an instance
+# of, the tuple of types. this is used to make sure that argument lists are
+# the right "shape" before destructuring (splatting?) them in a function call.
 def _match_type(values, types):
     if len(values) != len(types):
         return False
@@ -79,6 +99,15 @@ def _match_type(values, types):
 
 
 class OsmRawrLookup(object):
+    """
+    Implements the interface needed by the common code (e.g: layer_properties)
+    to look up information about node, way and relation IDs. For database
+    lookups, we previously did this with a JOIN, and the fixture data source
+    just iterates over the (small) number of items.
+
+    For RAWR tiles, we index the data to provide faster lookup, and are more
+    selective about what goes into the index.
+    """
 
     def __init__(self):
         self.nodes = {}
@@ -196,6 +225,7 @@ def transit_relations(self, rel_id):
         return set(self.relations_using_rel(rel_id))
 
 
+# yield all the tiles at the given zoom level which intersect the given bounds.
 def _tiles(zoom, unpadded_bounds):
     from tilequeue.tile import mercator_point_to_coord
     from raw_tiles.tile import Tile
@@ -218,6 +248,10 @@ def _tiles(zoom, unpadded_bounds):
             yield tile
 
 
+# the object which gets indexed. this is a normal (fid, shape, props) tuple
+# expanded to include a dict of layer name to min zoom in `layer_min_zooms`.
+# this means that the properties don't have to be copied and altered to
+# include the min zoom for each layer, reducing the memory footprint.
 _Feature = namedtuple('_Feature', 'fid shape properties layer_min_zooms')
 
 
@@ -241,6 +275,14 @@ def __getattr__(self, name):
 
 
 class _LayersIndex(object):
+    """
+    Index features by the tile(s) that they appear in.
+
+    This is done by calculating a min-min-zoom, the lowest min_zoom for that
+    feature across all layers, and then adding that feature to a list for each
+    tile it appears in from the min-min-zoom up to the max zoom for the tile
+    pyramid.
+    """
 
     def __init__(self, layers, tile_pyramid):
         self.layers = layers

From e552397d05a22f135de4c70eb054876488cff7fd Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 3 Oct 2017 12:46:58 +0100
Subject: [PATCH 16/19] Start threading through way and relation data via meta
 for min zoom calculations.

---
 tests/test_query_rawr.py |  78 +++++++++++++++++++++++++-
 tilequeue/query/rawr.py  | 117 +++++++++++++++++++++++++++++++++------
 2 files changed, 175 insertions(+), 20 deletions(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index a3f1d830..13933553 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -22,13 +22,14 @@ class RawrTestCase(unittest.TestCase):
     """
 
     def _make(self, min_zoom_fn, props_fn, tables, tile_pyramid,
-              layer_name='testlayer', label_placement_layers={}):
+              layer_name='testlayer', label_placement_layers={},
+              source='test'):
         from tilequeue.query.common import LayerInfo
         from tilequeue.query.rawr import make_rawr_data_fetcher
 
         layers = {layer_name: LayerInfo(min_zoom_fn, props_fn)}
         return make_rawr_data_fetcher(
-            layers, tables, tile_pyramid,
+            layers, tables, tile_pyramid, source,
             label_placement_layers=label_placement_layers)
 
 
@@ -399,7 +400,8 @@ def props_fn(shape, props, fid, meta):
         layers = {}
         for name in input_layer_names:
             layers[name] = LayerInfo(min_zoom_fn, props_fn)
-        fetch = make_rawr_data_fetcher(layers, tables, tile_pyramid)
+        source = 'test'
+        fetch = make_rawr_data_fetcher(layers, tables, tile_pyramid, source)
 
         read_rows = fetch(tile.zoom, coord_to_mercator_bounds(tile))
         # the RAWR query goes over features multiple times because of the
@@ -442,3 +444,73 @@ def test_precedence(self):
         self._test(['pois', 'landuse', 'buildings'], ['pois'])
         # otherwise, landuse should take precedence over buildings.
         self._test(['landuse', 'buildings'], ['landuse'])
+
+
+class TestMeta(RawrTestCase):
+
+    def test_meta_gate(self):
+        # test that the meta is passed to the min zoom function, and that we
+        # can use it to get information about the highways that a gate is
+        # part of.
+
+        from shapely.geometry import Point, LineString
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+
+        feature_min_zoom = 11
+
+        def min_zoom_fn(shape, props, fid, meta):
+            self.assertIsNotNone(meta)
+
+            # expect meta to have a source, which is a string name for the
+            # source of the data.
+            self.assertEquals('test', meta.source)
+
+            # expect meta to have a list of relations, which is empty for this
+            # test.
+            self.assertEquals(0, len(meta.relations))
+
+            # only do this for the node
+            if fid == 0:
+                # expect meta to have a list of ways, each of which is a (fid,
+                # shape, props) tuple, of which only props is used.
+                self.assertEquals(1, len(meta.ways))
+                way_fid, way_shape, way_props = meta.ways[0]
+                self.assertEquals(1, way_fid)
+                self.assertEquals({'highway': 'secondary'}, way_props)
+
+            # only set a min zoom for the node - this just simplifies the
+            # checking later, as there'll only be one feature.
+            return feature_min_zoom if fid == 0 else None
+
+        shape = Point(0, 0)
+        way_shape = LineString([[0, 0], [1, 1]])
+        # get_table(table_name) should return a generator of rows.
+        tables = TestGetTable({
+            'planet_osm_point': [(0, shape.wkb, {'barrier': 'gate'})],
+            'planet_osm_line': [(1, way_shape.wkb, {'highway': 'secondary'})],
+            'planet_osm_ways': [(1, [0], ['highway', 'secondary'])],
+        })
+
+        zoom = 10
+        max_zoom = zoom + 5
+        coord = mercator_point_to_coord(zoom, shape.x, shape.y)
+        tile_pyramid = TilePyramid(zoom, coord.column, coord.row, max_zoom)
+
+        fetch = self._make(min_zoom_fn, None, tables, tile_pyramid)
+
+        # first, check that it can get the original item back when both the
+        # min zoom filter and geometry filter are okay.
+        feature_coord = mercator_point_to_coord(
+            feature_min_zoom, shape.x, shape.y)
+        read_rows = fetch(
+            feature_min_zoom, coord_to_mercator_bounds(feature_coord))
+
+        self.assertEquals(1, len(read_rows))
+        read_row = read_rows[0]
+        self.assertEquals(0, read_row.get('__id__'))
+        # query processing code expects WKB bytes in the __geometry__ column
+        self.assertEquals(shape.wkb, read_row.get('__geometry__'))
+        self.assertEquals({'min_zoom': 11, 'barrier': 'gate'},
+                          read_row.get('__testlayer_properties__'))
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 236e9aee..3c625a1e 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -98,6 +98,19 @@ def _match_type(values, types):
     return True
 
 
+# return true if the tags indicate that this is a gate
+def _is_gate(props):
+    return props.get('barrier') == 'gate'
+
+
+# return true if the tags indicate that this is a highway, cycleway or footway
+# which might be part of a route relation. note that this is pretty loose, and
+# might return true for things we don't eventually render as roads, but is just
+# aimed at cutting down the number of items we need in our index.
+def _is_routeable(props):
+    return props.get('whitewater') == 'portage_way' or 'highway' in props
+
+
 class OsmRawrLookup(object):
     """
     Implements the interface needed by the common code (e.g: layer_properties)
@@ -131,7 +144,7 @@ def add_row(self, *args):
         if _match_type(args, (num, (str, bytes), dict)):
             self.add_feature(*args)
 
-        elif _match_type(args, (num, list, dict)):
+        elif _match_type(args, (num, list, list)):
             self.add_way(*args)
 
         elif _match_type(args, (num, num, num, list, list, list)):
@@ -151,16 +164,25 @@ def add_feature(self, fid, shape_wkb, props):
             # must be a station or stop node
             self.nodes[fid] = (fid, shape_wkb, props)
 
+        elif _is_gate(props) and shape_type == ShapeType.point:
+            # index the highways that use gates to influence min zoom
+            self.nodes[fid] = (fid, shape_wkb, props)
+
         elif (is_station_or_line(fid, None, props) and
               shape_type != ShapeType.point):
             # must be a station polygon or stop line
             self.ways[fid] = (fid, shape_wkb, props)
 
+        elif _is_routeable(props) and shape_type == ShapeType.line:
+            # index routable items (highways, cycleways, footpaths) to
+            # get the relations using them.
+            self.ways[fid] = (fid, shape_wkb, props)
+
     def add_way(self, way_id, nodes, tags):
         for node_id in nodes:
             if node_id in self.nodes:
-                self._ways_using_node[node_id] = way_id
                 assert way_id in self.ways
+                self._ways_using_node[node_id].append(way_id)
 
     def add_relation(self, rel_id, way_off, rel_off, parts, members, tags):
         r = Relation(rel_id, way_off, rel_off, parts, members, tags)
@@ -274,6 +296,9 @@ def __getattr__(self, name):
         return getattr(self.obj, name)
 
 
+_Metadata = namedtuple('_Metadata', 'source ways relations')
+
+
 class _LayersIndex(object):
     """
     Index features by the tile(s) that they appear in.
@@ -284,23 +309,71 @@ class _LayersIndex(object):
     pyramid.
     """
 
-    def __init__(self, layers, tile_pyramid):
+    def __init__(self, layers, tile_pyramid, source):
         self.layers = layers
         self.tile_pyramid = tile_pyramid
         self.tile_index = defaultdict(list)
+        self.source = source
+        self.delayed_features = []
 
     def add_row(self, fid, shape_wkb, props):
         shape = _LazyShape(shape_wkb)
-        # TODO! meta is not None!
-        meta = None
+        # single object (hence single id()) will be shared amongst all layers.
+        # this allows us to easily and quickly de-duplicate at later layers in
+        # the stack.
+        feature = _Feature(fid, shape, props, {})
+
+        # delay min zoom calculation in order to collect more information about
+        # the ways and relations using a particular feature.
+        self.delayed_features.append(feature)
+
+    def index(self, osm):
+        for feature in self.delayed_features:
+            self._index_feature(feature, osm)
+        del self.delayed_features
+
+    def _index_feature(self, feature, osm):
+        fid = feature.fid
+        shape = feature.shape
+        props = feature.properties
+        layer_min_zooms = feature.layer_min_zooms
+
+        # grab the shape type without decoding the WKB to save time.
+        shape_type = _wkb_shape(shape.wkb)
+
+        ways = []
+        rels = []
+
+        # fetch ways and relations for any node
+        if fid >= 0 and shape_type == ShapeType.point:
+            for way_id in osm.ways_using_node(fid):
+                ways.append(osm.way(way_id))
+            for rel_id in osm.relations_using_node(fid):
+                rels.append(osm.relation(rel_id))
+
+        # and relations for any way
+        if fid >= 0 and shape_type == ShapeType.line:
+            for rel_id in osm.relations_using_way(fid):
+                rels.append(osm.relation(rel_id))
+
+        # have to transform the Relation object into a dict, which is
+        # what the functions called on this data expect.
+        # TODO: reusing the Relation object would be better.
+        rel_dicts = []
+        for r in rels:
+            tags = []
+            for k, v in r.tags.items():
+                tags.append(k)
+                tags.append(v)
+            rel_dicts.append(dict(tags=tags))
+
+        meta = _Metadata(self.source, ways, rel_dicts)
 
-        layer_min_zooms = {}
         for layer_name, info in self.layers.items():
-            shape_type = _wkb_shape(shape_wkb)
             shape_type_str = ShapeType.lookup(shape_type)
             if info.shape_types and shape_type_str not in info.shape_types:
                 continue
-            min_zoom = info.min_zoom_fn(fid, shape, props, meta)
+            min_zoom = info.min_zoom_fn(shape, props, fid, meta)
             if min_zoom is not None:
                 layer_min_zooms[layer_name] = min_zoom
 
@@ -316,11 +389,6 @@ def add_row(self, fid, shape_wkb, props):
         # visibility.
         min_zoom = min(self.tile_pyramid.max_z, min(layer_min_zooms.values()))
 
-        # single object (hence single id()) will be shared amongst all layers.
-        # this allows us to easily and quickly de-duplicate at later layers in
-        # the stack.
-        feature = _Feature(fid, shape, props, layer_min_zooms)
-
         # take the minimum integer zoom - this is the min zoom tile that the
         # feature should appear in, and a feature with min_zoom = 1.9 should
         # appear in a tile at z=1, not 2, since the tile at z=N is used for
@@ -351,7 +419,8 @@ def __call__(self, tile):
 
 class DataFetcher(object):
 
-    def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
+    def __init__(self, layers, tables, tile_pyramid, label_placement_layers,
+                 source):
         """
         Expect layers to be a dict of layer name to LayerInfo (see fixture.py).
         Tables should be a callable which returns a generator over the rows in
@@ -363,11 +432,13 @@ def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
         self.layers = layers
         self.tile_pyramid = tile_pyramid
         self.label_placement_layers = label_placement_layers
+        self.source = source
         self.layer_indexes = {}
 
         table_indexes = defaultdict(list)
 
-        self.layers_index = _LayersIndex(self.layers, self.tile_pyramid)
+        self.layers_index = _LayersIndex(
+            self.layers, self.tile_pyramid, self.source)
         for shape_type in ('point', 'line', 'polygon'):
             table_name = 'planet_osm_' + shape_type
             table_indexes[table_name].append(self.layers_index)
@@ -382,6 +453,17 @@ def __init__(self, layers, tables, tile_pyramid, label_placement_layers):
             extra_indexes = table_indexes[table_name]
             index_table(source, self.osm, *extra_indexes)
 
+        # there's a chicken and egg problem with the indexes: we want to know
+        # which features to index, but also calculate the feature's min zoom,
+        # which might depend on ways and relations not seen yet. one solution
+        # would be to do this in two passes, but that might mean paying a cost
+        # to decompress or deserialize the data twice. instead, the index
+        # buffers the features and indexes them in the following step. this
+        # might mean we buffer more information in memory than we technically
+        # need if many of the features are not visible, but means we get one
+        # single set of _Feature objects.
+        self.layers_index.index(self.osm)
+
     def _named_layer(self, fid, shape, props):
         # we want only one layer from ('pois', 'landuse', 'buildings') for
         # each feature to be assigned a name. therefore, we use the presence
@@ -486,6 +568,7 @@ def __call__(self, zoom, unpadded_bounds):
 
 # tables is a callable which should return a generator over the rows of the
 # table when called with the table name.
-def make_rawr_data_fetcher(layers, tables, tile_pyramid,
+def make_rawr_data_fetcher(layers, tables, tile_pyramid, source,
                            label_placement_layers={}):
-    return DataFetcher(layers, tables, tile_pyramid, label_placement_layers)
+    return DataFetcher(layers, tables, tile_pyramid, label_placement_layers,
+                       source)

From 60421729641d425b43deb426051d568fb3df818d Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 3 Oct 2017 19:22:38 +0100
Subject: [PATCH 17/19] Add test for route relations via meta.

---
 tests/test_query_rawr.py | 66 ++++++++++++++++++++++++++++++++++++++++
 tilequeue/query/rawr.py  |  5 ++-
 2 files changed, 70 insertions(+), 1 deletion(-)

diff --git a/tests/test_query_rawr.py b/tests/test_query_rawr.py
index 13933553..14fdaa69 100644
--- a/tests/test_query_rawr.py
+++ b/tests/test_query_rawr.py
@@ -514,3 +514,69 @@ def min_zoom_fn(shape, props, fid, meta):
         self.assertEquals(shape.wkb, read_row.get('__geometry__'))
         self.assertEquals({'min_zoom': 11, 'barrier': 'gate'},
                           read_row.get('__testlayer_properties__'))
+
+    def test_meta_route(self):
+        # test that we can use meta in the min zoom function to find out which
+        # route(s) a road is part of.
+
+        from shapely.geometry import LineString
+        from tilequeue.query.rawr import TilePyramid
+        from tilequeue.tile import coord_to_mercator_bounds
+        from tilequeue.tile import mercator_point_to_coord
+        from tilequeue.query.common import deassoc
+
+        feature_min_zoom = 11
+
+        rel_tags = [
+            'type', 'route',
+            'route', 'road',
+            'ref', '101',
+        ]
+
+        def min_zoom_fn(shape, props, fid, meta):
+            self.assertIsNotNone(meta)
+
+            # expect meta to have a source, which is a string name for the
+            # source of the data.
+            self.assertEquals('test', meta.source)
+
+            # expect meta to have a list of ways, but empty for this test.
+            self.assertEquals(0, len(meta.ways))
+
+            # expect meta to have a list of relations, each of which is a dict
+            # containing at least the key 'tags' mapped to a list of
+            # alternating k, v suitable for passing into deassoc().
+            self.assertEquals(1, len(meta.relations))
+            rel = meta.relations[0]
+            self.assertIsInstance(rel, dict)
+            self.assertIn('tags', rel)
+            self.assertEquals(deassoc(rel_tags), deassoc(rel['tags']))
+
+            return feature_min_zoom
+
+        shape = LineString([[0, 0], [1, 1]])
+        # get_table(table_name) should return a generator of rows.
+        tables = TestGetTable({
+            'planet_osm_line': [(1, shape.wkb, {'highway': 'secondary'})],
+            'planet_osm_rels': [(2, 0, 1, [1], [''], rel_tags)],
+        })
+
+        zoom = 10
+        max_zoom = zoom + 5
+        coord = mercator_point_to_coord(zoom, *shape.coords[0])
+        tile_pyramid = TilePyramid(zoom, coord.column, coord.row, max_zoom)
+
+        fetch = self._make(min_zoom_fn, None, tables, tile_pyramid)
+
+        # first, check that it can get the original item back when both the
+        # min zoom filter and geometry filter are okay.
+        feature_coord = mercator_point_to_coord(
+            feature_min_zoom, *shape.coords[0])
+        read_rows = fetch(
+            feature_min_zoom, coord_to_mercator_bounds(feature_coord))
+
+        self.assertEquals(1, len(read_rows))
+        read_row = read_rows[0]
+        self.assertEquals(1, read_row.get('__id__'))
+        self.assertEquals({'min_zoom': 11, 'highway': 'secondary'},
+                          read_row.get('__testlayer_properties__'))
diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 3c625a1e..693ac991 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -186,7 +186,10 @@ def add_way(self, way_id, nodes, tags):
 
     def add_relation(self, rel_id, way_off, rel_off, parts, members, tags):
         r = Relation(rel_id, way_off, rel_off, parts, members, tags)
-        if mz_is_interesting_transit_relation(r.tags):
+        is_transit_relation = mz_is_interesting_transit_relation(r.tags)
+        is_route = 'route' in r.tags and \
+                   ('network' in r.tags or 'ref' in r.tags)
+        if is_route or is_transit_relation:
             self.relations[r.id] = r
             for node_id in r.node_ids:
                 if node_id in self.nodes:

From a38add4df14122d00d842f8eee7a4eef2833938b Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 3 Oct 2017 19:29:54 +0100
Subject: [PATCH 18/19] Pass meta into min zoom function when determining which
 layers are labelled.

---
 tilequeue/query/rawr.py | 64 ++++++++++++++++++++++-------------------
 1 file changed, 34 insertions(+), 30 deletions(-)

diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 693ac991..701f51cd 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -302,6 +302,36 @@ def __getattr__(self, name):
 _Metadata = namedtuple('_Metadata', 'source ways relations')
 
 
+def _make_meta(source, fid, shape_type, osm):
+    ways = []
+    rels = []
+
+    # fetch ways and relations for any node
+    if fid >= 0 and shape_type == ShapeType.point:
+        for way_id in osm.ways_using_node(fid):
+            ways.append(osm.way(way_id))
+        for rel_id in osm.relations_using_node(fid):
+            rels.append(osm.relation(rel_id))
+
+    # and relations for any way
+    if fid >= 0 and shape_type == ShapeType.line:
+        for rel_id in osm.relations_using_way(fid):
+            rels.append(osm.relation(rel_id))
+
+    # have to transform the Relation object into a dict, which is
+    # what the functions called on this data expect.
+    # TODO: reusing the Relation object would be better.
+    rel_dicts = []
+    for r in rels:
+        tags = []
+        for k, v in r.tags.items():
+            tags.append(k)
+            tags.append(v)
+        rel_dicts.append(dict(tags=tags))
+
+    return _Metadata(source, ways, rel_dicts)
+
+
 class _LayersIndex(object):
     """
     Index features by the tile(s) that they appear in.
@@ -344,34 +374,7 @@ def _index_feature(self, feature, osm):
         # grab the shape type without decoding the WKB to save time.
         shape_type = _wkb_shape(shape.wkb)
 
-        ways = []
-        rels = []
-
-        # fetch ways and relations for any node
-        if fid >= 0 and shape_type == ShapeType.point:
-            for way_id in osm.ways_using_node(fid):
-                ways.append(osm.way(way_id))
-            for rel_id in osm.relations_using_node(fid):
-                rels.append(osm.relation(rel_id))
-
-        # and relations for any way
-        if fid >= 0 and shape_type == ShapeType.line:
-            for rel_id in osm.relations_using_way(fid):
-                rels.append(osm.relation(rel_id))
-
-        # have to transform the Relation object into a dict, which is
-        # what the functions called on this data expect.
-        # TODO: reusing the Relation object would be better.
-        rel_dicts = []
-        for r in rels:
-            tags = []
-            for k, v in r.tags.items():
-                tags.append(k)
-                tags.append(v)
-            rel_dicts.append(dict(tags=tags))
-
-        meta = _Metadata(self.source, ways, rel_dicts)
-
+        meta = _make_meta(self.source, fid, shape_type, osm)
         for layer_name, info in self.layers.items():
             shape_type_str = ShapeType.lookup(shape_type)
             if info.shape_types and shape_type_str not in info.shape_types:
@@ -475,8 +478,9 @@ def _named_layer(self, fid, shape, props):
         for layer_name in ('pois', 'landuse', 'buildings'):
             info = self.layers.get(layer_name)
             if info and info.min_zoom_fn:
-                # TODO! meta should not be None!
-                min_zoom = info.min_zoom_fn(shape, props, fid, None)
+                shape_type = _wkb_shape(shape.wkb)
+                meta = _make_meta(self.source, fid, shape_type, self.osm)
+                min_zoom = info.min_zoom_fn(shape, props, fid, meta)
                 if min_zoom is not None:
                     return layer_name
         return None

From 8a73421f4e1a7d4b9a7cf7da2ce150c3bd5754bc Mon Sep 17 00:00:00 2001
From: Matt Amos <matt.amos@mapzen.com>
Date: Tue, 3 Oct 2017 19:32:19 +0100
Subject: [PATCH 19/19] Rewrote named layer selection to be much simpler, now
 that the min zooms are precomputed for all layers.

---
 tilequeue/query/rawr.py | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/tilequeue/query/rawr.py b/tilequeue/query/rawr.py
index 701f51cd..927315ca 100644
--- a/tilequeue/query/rawr.py
+++ b/tilequeue/query/rawr.py
@@ -470,19 +470,16 @@ def __init__(self, layers, tables, tile_pyramid, label_placement_layers,
         # single set of _Feature objects.
         self.layers_index.index(self.osm)
 
-    def _named_layer(self, fid, shape, props):
+    def _named_layer(self, layer_min_zooms):
         # we want only one layer from ('pois', 'landuse', 'buildings') for
         # each feature to be assigned a name. therefore, we use the presence
         # or absence of a min zoom to check whether these features as in these
-        # layers, and therefore which should be assigned the name.
+        # layers, and therefore which should be assigned the name. handily,
+        # the min zooms are already pre-calculated as layer_min_zooms from the
+        # index.
         for layer_name in ('pois', 'landuse', 'buildings'):
-            info = self.layers.get(layer_name)
-            if info and info.min_zoom_fn:
-                shape_type = _wkb_shape(shape.wkb)
-                meta = _make_meta(self.source, fid, shape_type, self.osm)
-                min_zoom = info.min_zoom_fn(shape, props, fid, meta)
-                if min_zoom is not None:
-                    return layer_name
+            if layer_name in layer_min_zooms:
+                return layer_name
         return None
 
     def _lookup(self, zoom, unpadded_bounds):
@@ -525,7 +522,7 @@ def __call__(self, zoom, unpadded_bounds):
             # add name into whichever of the pois, landuse or buildings
             # layers has claimed this feature.
             name = props.get('name', None)
-            named_layer = self._named_layer(fid, shape, props)
+            named_layer = self._named_layer(layer_min_zooms)
 
             for layer_name, min_zoom in layer_min_zooms.items():
                 # we need to keep fractional zooms, e.g: 4.999 should appear