Merge pull request #3590 from alphagov/buffer-de-buffer

Make simplification of polygons more sophisticated

app/broadcast_areas/__init__.py

@@ -1,9 +1,7 @@
-import itertools
-
-import geojson
 from notifications_utils.serialised_model import SerialisedModelCollection
 from werkzeug.utils import cached_property
 
+from .polygons import Polygons
 from .repo import BroadcastAreasRepository
 
 
@@ -34,59 +32,17 @@ class BroadcastArea(SortableMixin):
     def __eq__(self, other):
         return self.id == other.id
 
-    @property
+    @cached_property
-    def _feature(self):
+    def polygons(self):
-        return BroadcastAreasRepository().get_feature_for_area(self.id)
+        return Polygons(
-
+            BroadcastAreasRepository().get_polygons_for_area(self.id)
-    @property
-    def _simple_feature(self):
-        return BroadcastAreasRepository().get_simple_feature_for_area(self.id)
-
-    def _polygons(self, feature):
-        if feature['geometry']['type'] == 'MultiPolygon':
-            return [
-                polygons[0]
-                for polygons in feature['geometry']['coordinates']
-            ]
-        if feature['geometry']['type'] == 'Polygon':
-            return [
-                feature['geometry']['coordinates'][0]
-            ]
-        raise TypeError(
-            f'Unknown geometry type {self.feature["geometry"]["type"]} '
-            f'in {self.__class__.__name} {self.name}'
         )
 
-    def _unenclosed_polygons(self, feature):
+    @cached_property
-        # Some mapping tools require shapes to be unenclosed, i.e. the
-        # last point joins the first point implicitly
-        return [
-            coordinates[:-1] for coordinates in self._polygons(feature)
-        ]
-
-    @property
-    def polygons(self):
-        return self._polygons(self.feature)
-
-    @property
-    def unenclosed_polygons(self):
-        return self._unenclosed_polygons(self.feature)
-
-    @property
     def simple_polygons(self):
-        return self._polygons(self.simple_feature)
+        return Polygons(
-
+            BroadcastAreasRepository().get_simple_polygons_for_area(self.id)
-    @property
+        )
-    def simple_unenclosed_polygons(self):
-        return self._unenclosed_polygons(self.simple_feature)
-
-    @cached_property
-    def feature(self):
-        return geojson.loads(self._feature)
-
-    @cached_property
-    def simple_feature(self):
-        return geojson.loads(self._simple_feature)
 
     @property
     def sub_areas(self):
@@ -127,29 +83,5 @@ class BroadcastAreaLibraries(SerialisedModelCollection, GetItemByIdMixin):
         areas = BroadcastAreasRepository().get_areas(area_ids)
         return [BroadcastArea(area) for area in areas]
 
-    def get_polygons_for_areas_long_lat(self, *area_ids):
-        return list(itertools.chain(*(
-            area.polygons
-            for area in self.get_areas(*area_ids)
-        )))
-
-    def get_polygons_for_areas_lat_long(self, *area_ids):
-        return [
-            [[long, lat] for lat, long in polygon]
-            for polygon in self.get_polygons_for_areas_long_lat(*area_ids)
-        ]
-
-    def get_simple_polygons_for_areas_long_lat(self, *area_ids):
-        return list(itertools.chain(*(
-            area.simple_polygons
-            for area in self.get_areas(*area_ids)
-        )))
-
-    def get_simple_polygons_for_areas_lat_long(self, *area_ids):
-        return [
-            [[long, lat] for lat, long in polygon]
-            for polygon in self.get_simple_polygons_for_areas_long_lat(*area_ids)
-        ]
-
 
 broadcast_area_libraries = BroadcastAreaLibraries()

app/broadcast_areas/create-broadcast-areas-db.py

@@ -1,68 +1,54 @@
 #!/usr/bin/env python
 
-from copy import deepcopy
 from pathlib import Path
 
 import geojson
-import shapely.geometry as sgeom
+from notifications_utils.formatters import formatted_list
 
+from polygons import Polygons
 from repo import BroadcastAreasRepository
 
 package_path = Path(__file__).resolve().parent
-
+point_counts = []
-
-def convert_shape_to_feature(shape):
-    return {
-        "type": "Feature",
-        "properties": {},
-        "geometry": sgeom.mapping(shape),
-    }
-
-
-def simplify_polygon(series):
-    polygon, *_holes = series  # discard holes
-
-    approx_metres_to_degree = 111320
-    # initially buffer (extend area past perimeter) by ~25m
-    buffer_degrees = 500 / approx_metres_to_degree
-    # initially simplify (snap to closest point) by ~25m
-    simplify_degrees = 50.0 / approx_metres_to_degree
-
-    starting_polygon = sgeom.LineString(polygon).buffer(buffer_degrees)
-    simplified_polygon = None
-    num_polys = len(polygon)
-    last_num_polys = []
-    while True:
-        simplified_polygon = starting_polygon.simplify(simplify_degrees)
-        simplified_polygon = [[c[0], c[1]] for c in simplified_polygon.exterior.coords]
-
-        num_polys = len(simplified_polygon)
-        simplify_degrees *= 2
-
-        if num_polys <= 99 or last_num_polys[-3:] == [num_polys, num_polys, num_polys]:
-            break
-
-        last_num_polys.append(num_polys)
-
-        print(".", end="", flush=True)  # noqa: T001
-
-    return [simplified_polygon]
 
 
 def simplify_geometry(feature):
     if feature["type"] == "Polygon":
-        feature["coordinates"] = simplify_polygon(feature["coordinates"])
+        return [feature["coordinates"][0]]
-        return feature
     elif feature["type"] == "MultiPolygon":
-        feature["coordinates"] = [
+        return [polygon for polygon, *_holes in feature["coordinates"]]
-            simplify_polygon(polygon)
-            for polygon in feature["coordinates"]
-        ]
-        return feature
     else:
         raise Exception("Unknown type: {}".format(feature["type"]))
 
 
+def polygons_and_simplified_polygons(feature):
+
+    polygons = Polygons(simplify_geometry(feature))
+    full_resolution = polygons.remove_too_small
+    smoothed = full_resolution.smooth
+    simplified = smoothed.simplify
+
+    print(  # noqa: T001
+        f'    Original:{full_resolution.point_count: >5} points'
+        f'    Smoothed:{smoothed.point_count: >5} points'
+        f'    Simplified:{simplified.point_count: >4} points'
+    )
+
+    point_counts.append(simplified.point_count)
+
+    if simplified.point_count >= 200:
+        raise RuntimeError(
+            'Too many points '
+            '(adjust Polygons.perimeter_to_simplification_ratio or '
+            'Polygons.perimeter_to_buffer_ratio)'
+        )
+
+    return (
+        full_resolution.as_coordinate_pairs_long_lat,
+        simplified.as_coordinate_pairs_long_lat,
+    )
+
+
 repo = BroadcastAreasRepository()
 
 repo.delete_db()
@@ -88,10 +74,12 @@ for dataset_name, dataset_name_singular, id_field, name_field in simple_datasets
         f_id = dataset_id + "-" + feature["properties"][id_field]
         f_name = feature["properties"][name_field]
 
+        print()  # noqa: T001
         print(f_name)  # noqa: T001
 
-        simple_feature = deepcopy(feature)
+        feature, simple_feature = (
-        simple_feature["geometry"] = simplify_geometry(simple_feature["geometry"])
+            polygons_and_simplified_polygons(feature["geometry"])
+        )
 
         repo.insert_broadcast_areas([[
             f_id, f_name,
@@ -132,19 +120,21 @@ for f in geojson.loads(wards_filepath.read_text())["features"]:
     ward_name = f["properties"]["wd20nm"]
     ward_id = "wd20-" + ward_code
 
+    print()  # noqa: T001
     print(ward_name)  # noqa: T001
 
     try:
         la_id = "lad20-" + ward_code_to_la_id_mapping[ward_code]
         la_name = ward_code_to_la_mapping[ward_code]
 
-        sf = deepcopy(f)
+        feature, simple_feature = (
-        sf["geometry"] = simplify_geometry(sf["geometry"])
+            polygons_and_simplified_polygons(f["geometry"])
+        )
 
         areas_to_add.append([
             ward_id, ward_name,
             dataset_id, la_id,
-            f, sf
+            feature, simple_feature
         ])
 
     except KeyError:
@@ -159,12 +149,14 @@ for feature in geojson.loads(las_filepath.read_text())["features"]:
     la_id = feature["properties"]["lad20cd"]
     group_name = feature["properties"]["lad20nm"]
 
+    print()  # noqa: T001
     print(group_name)  # noqa: T001
 
     group_id = "lad20-" + la_id
 
-    simple_feature = deepcopy(feature)
+    feature, simple_feature = (
-    simple_feature["geometry"] = simplify_geometry(simple_feature["geometry"])
+        polygons_and_simplified_polygons(feature["geometry"])
+    )
 
     areas_to_add.append([
         group_id, group_name,
@@ -173,3 +165,16 @@ for feature in geojson.loads(las_filepath.read_text())["features"]:
     ])
 
 repo.insert_broadcast_areas(areas_to_add)
+
+
+most_detailed_polygons = formatted_list(
+    sorted(point_counts, reverse=True)[:5],
+    before_each='',
+    after_each='',
+)
+print(  # noqa: T001
+    '\n'
+    'DONE\n'
+    f'    Processed {len(point_counts):,} polygons.\n'
+    f'    Highest point counts once simplifed: {most_detailed_polygons}\n'
+)

app/broadcast_areas/polygons.py

@@ -0,0 +1,177 @@
+import itertools
+
+from shapely.geometry import (
+    JOIN_STYLE,
+    GeometryCollection,
+    MultiPolygon,
+    Polygon,
+)
+from shapely.ops import unary_union
+from werkzeug.utils import cached_property
+
+
+class Polygons():
+
+    approx_metres_to_degree = 111_320
+    approx_square_metres_to_square_degree = approx_metres_to_degree ** 2
+
+    # Estimated amount of bleed into neigbouring areas based on typical
+    # range/separation of cell towers.
+    approx_bleed_in_degrees = 1_500 / approx_metres_to_degree
+
+    # Controls how much buffer to add for a shape of a given perimeter.
+    # Smaller number means more buffering and a smoother shape. For
+    # example `1000` means 1m of buffer for every 1km of perimeter, or
+    # 20m of buffer for a 5km square. This gives us control over how
+    # much we fill in very concave features like channels, harbours and
+    # zawns.
+    perimeter_to_buffer_ratio = 360
+
+    # Ratio of how much detail a shape of a given perimeter has once
+    # simplified. Smaller number means less detail. For example `1000`
+    # means that for a shape with a perimeter of 1000m, the simplified
+    # line will never deviate more than 1m from the original.
+    # Or for a 5km square, the line won’t deviate more than 20m. This
+    # gives us approximate control over the total number of points.
+    perimeter_to_simplification_ratio = 1_750
+
+    # The threshold for removing very small areas from the map. These
+    # areas are likely glitches in  the data where the shoreline hasn’t
+    # been subtracted from the land properly
+    minimum_area_size_square_metres = 50 ** 2
+
+    def __init__(self, polygons):
+        if not polygons:
+            self.polygons = []
+        elif isinstance(polygons[0], list):
+            self.polygons = [
+                Polygon(polygon) for polygon in polygons
+            ]
+        else:
+            self.polygons = polygons
+
+    def __getitem__(self, index):
+        return self.polygons[index]
+
+    def __len__(self):
+        return len(self.polygons)
+
+    @cached_property
+    def perimeter_length(self):
+        return sum(
+            polygon.length for polygon in self
+        )
+
+    @cached_property
+    def buffer_outward_in_degrees(self):
+        return (
+            # If two areas are close enough that the distance between
+            # them is less than the minimum bleed of a cell
+            # broadcast then this joins them together. The aim is to
+            # reduce the total number of polygons in areas with many
+            # small shapes like Orkney or the Isles of Scilly.
+            self.approx_bleed_in_degrees / 3
+        ) + (
+            self.perimeter_length / self.perimeter_to_buffer_ratio
+        )
+
+    @cached_property
+    def buffer_inward_in_degrees(self):
+        return self.buffer_outward_in_degrees - (
+            # We should leave the shape expanded by at least the
+            # simplification tolerance in all places, so the
+            # simplification never moves a point inside the original
+            # shape. In practice half of the tolerance is enough to
+            # acheive this.
+            self.simplification_tolerance_in_degrees / 2
+        )
+
+    @cached_property
+    def simplification_tolerance_in_degrees(self):
+        return self.perimeter_length / self.perimeter_to_simplification_ratio
+
+    @cached_property
+    def smooth(self):
+        buffered = [
+            polygon.buffer(
+                self.buffer_outward_in_degrees,
+                resolution=4,
+                join_style=JOIN_STYLE.round,
+            )
+            for polygon in self
+        ]
+        unioned = union_polygons(buffered)
+        debuffered = [
+            polygon.buffer(
+                -1 * self.buffer_inward_in_degrees,
+                resolution=1,
+                join_style=JOIN_STYLE.bevel,
+            )
+            for polygon in unioned
+        ]
+        flattened = list(itertools.chain(*[
+            flatten_polygons(polygon) for polygon in debuffered
+        ]))
+        return Polygons(flattened)
+
+    @cached_property
+    def simplify(self):
+        return Polygons([
+            polygon.simplify(self.simplification_tolerance_in_degrees)
+            for polygon in self
+        ])
+
+    @cached_property
+    def bleed(self):
+        return Polygons(union_polygons([
+            polygon.buffer(
+                self.approx_bleed_in_degrees,
+                resolution=4,
+                join_style=JOIN_STYLE.round,
+            )
+            for polygon in self
+        ]))
+
+    @cached_property
+    def remove_too_small(self):
+        return Polygons([
+            polygon for polygon in self
+            if (
+                polygon.area * self.approx_square_metres_to_square_degree
+            ) > (
+                self.minimum_area_size_square_metres
+            )
+        ])
+
+    @cached_property
+    def as_coordinate_pairs_long_lat(self):
+        return [
+            [[x, y] for x, y in polygon.exterior.coords]
+            for polygon in self
+        ]
+
+    @cached_property
+    def as_coordinate_pairs_lat_long(self):
+        return [
+            [[y, x] for x, y in coordinate_pairs]
+            for coordinate_pairs in self.as_coordinate_pairs_long_lat
+        ]
+
+    @cached_property
+    def point_count(self):
+        return len(list(itertools.chain(*self.as_coordinate_pairs_long_lat)))
+
+
+def flatten_polygons(polygons):
+    if isinstance(polygons, GeometryCollection):
+        return []
+    if isinstance(polygons, MultiPolygon):
+        return [
+            p for p in polygons
+        ]
+    else:
+        return [polygons]
+
+
+def union_polygons(polygons):
+    return flatten_polygons(unary_union(polygons))

app/broadcast_areas/repo.py

@@ -1,9 +1,8 @@
+import json
 import os
 import sqlite3
 from pathlib import Path
 
-import geojson
-
 
 class BroadcastAreasRepository(object):
     def __init__(self):
@@ -47,10 +46,10 @@ class BroadcastAreasRepository(object):
             )""")
 
             conn.execute("""
-            CREATE TABLE broadcast_area_features (
+            CREATE TABLE broadcast_area_polygons (
                 id TEXT PRIMARY KEY,
-                feature_geojson TEXT NOT NULL,
+                polygons TEXT NOT NULL,
-                simple_feature_geojson TEXT NOT NULL
+                simple_polygons TEXT NOT NULL
             )""")
 
             conn.execute("""
@@ -84,20 +83,20 @@ class BroadcastAreasRepository(object):
         """
 
         features_q = """
-        INSERT INTO broadcast_area_features (
+        INSERT INTO broadcast_area_polygons (
             id,
-            feature_geojson, simple_feature_geojson
+            polygons, simple_polygons
         )
         VALUES (?, ?, ?)
         """
 
         with self.conn() as conn:
-            for id, name, area_id, group, feature, simple_feature in areas:
+            for id, name, area_id, group, polygons, simple_polygons in areas:
                 conn.execute(areas_q, (
                     id, name, area_id, group,
                 ))
                 conn.execute(features_q, (
-                    id, geojson.dumps(feature), geojson.dumps(simple_feature),
+                    id, json.dumps(polygons), json.dumps(simple_polygons),
                 ))
 
     def query(self, sql, *args):
@@ -206,24 +205,24 @@ class BroadcastAreasRepository(object):
 
         return areas
 
-    def get_feature_for_area(self, area_id):
+    def get_polygons_for_area(self, area_id):
         q = """
-        SELECT feature_geojson
+        SELECT polygons
-        FROM broadcast_area_features
+        FROM broadcast_area_polygons
         WHERE id = ?
         """
 
         results = self.query(q, area_id)
 
-        return results[0][0]
+        return json.loads(results[0][0])
 
-    def get_simple_feature_for_area(self, area_id):
+    def get_simple_polygons_for_area(self, area_id):
         q = """
-        SELECT simple_feature_geojson
+        SELECT simple_polygons
-        FROM broadcast_area_features
+        FROM broadcast_area_polygons
         WHERE id = ?
         """
 
         results = self.query(q, area_id)
 
-        return results[0][0]
+        return json.loads(results[0][0])

app/models/broadcast_message.py

@@ -1,12 +1,12 @@
+import itertools
 from datetime import datetime, timedelta
 
 from notifications_utils.template import BroadcastPreviewTemplate
 from orderedset import OrderedSet
-from shapely.geometry import MultiPolygon, Polygon
-from shapely.ops import unary_union
 from werkzeug.utils import cached_property
 
 from app.broadcast_areas import broadcast_area_libraries
+from app.broadcast_areas.polygons import Polygons
 from app.models import JSONModel, ModelList
 from app.models.user import User
 from app.notify_client.broadcast_message_api_client import (
@@ -72,30 +72,24 @@ class BroadcastMessage(JSONModel):
             area.name for area in self.areas
         ][:10]
 
-    @property
+    @cached_property
     def polygons(self):
-        return broadcast_area_libraries.get_polygons_for_areas_lat_long(
+        return Polygons(
-            *self._dict['areas']
+            list(itertools.chain(*(
+                area.polygons for area in self.areas
+            )))
         )
 
-    @property
+    @cached_property
     def simple_polygons(self):
-        simple_polygons = broadcast_area_libraries.get_simple_polygons_for_areas_lat_long(
+        polygons = Polygons(
-            *self._dict['areas']
+            list(itertools.chain(*(
+                area.simple_polygons for area in self.areas
+            )))
         )
-        unioned_polygons = unary_union([
+        # If we’ve added multiple areas then we need to re-simplify the
-            Polygon(i) for i in simple_polygons
+        # combined shapes to keep the point count down
-        ])
+        return polygons.smooth.simplify if len(self.areas) > 1 else polygons
-        if isinstance(unioned_polygons, MultiPolygon):
-            return [
-                [
-                    [x, y] for x, y in p.exterior.coords
-                ]
-                for p in unioned_polygons
-            ]
-        return [[
-            [x, y] for x, y in unioned_polygons.exterior.coords
-        ]]
 
     @property
     def template(self):

app/templates/views/broadcast/preview-areas.html

@@ -36,19 +36,31 @@
       attribution: '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors'
     }).addTo(mymap);
 
-    {% for polygon in broadcast_message.simple_polygons %}
+    {% for polygon in broadcast_message.simple_polygons.bleed.as_coordinate_pairs_lat_long %}
       polygons.push(
         L.polygon({{polygon}}, {
+          opacity: 0.5,
           color: '#2B8CC4', // $light-blue
-          opacity: 0.4,
           fillColor: '#2B8CC4', // $light-blue
-          fillOpacity: 0.3,
+          fillOpacity: 0.2,
+          weight: 2
+        })
+      );
+    {% endfor %}
+
+    {% for polygon in broadcast_message.simple_polygons.as_coordinate_pairs_lat_long %}
+      polygons.push(
+        L.polygon({{polygon}}, {
+          opacity: 0.1,
+          color: '#2B8CC4', // $black
+          fillColor: '#2B8CC4', // $light-blue
+          fillOpacity: 0.4,
           weight: 1
         })
       );
     {% endfor %}
 
-    {% for polygon in broadcast_message.polygons %}
+    {% for polygon in broadcast_message.polygons.as_coordinate_pairs_lat_long %}
       polygons.push(
         L.polygon({{polygon}}, {
           color: '#0b0b0c', // $black
@@ -61,7 +73,7 @@
     var polygonGroup = L.featureGroup(polygons).addTo(mymap);
     mymap.fitBounds(
       polygonGroup.getBounds(),
-      {padding: [5, 5]}
+      {padding: [1, 1]}
     );
 
   </script>

tests/app/broadcast_areas/test_broadcast_area.py

@@ -87,52 +87,32 @@ def test_get_areas_accepts_lists():
 
 def test_has_polygons():
 
-    assert len(
+    england = broadcast_area_libraries.get_areas('ctry19-E92000001')[0]
-        broadcast_area_libraries.get_polygons_for_areas_long_lat('ctry19-E92000001')
+    scotland = broadcast_area_libraries.get_areas('ctry19-S92000003')[0]
-    ) == 35
 
-    assert len(
+    assert len(england.polygons) == 35
-        broadcast_area_libraries.get_polygons_for_areas_long_lat('ctry19-S92000003')
+    assert len(scotland.polygons) == 195
-    ) == 195
 
-    assert len(
+    assert england.polygons.as_coordinate_pairs_lat_long[0][0] == [
-        broadcast_area_libraries.get_polygons_for_areas_long_lat(
-            'ctry19-E92000001',
-            'ctry19-S92000003',
-        )
-    ) == 35 + 195 == 230
-
-    assert len(
-        broadcast_area_libraries.get_polygons_for_areas_lat_long(
-            'ctry19-E92000001',
-            'ctry19-S92000003',
-        )
-    ) == 35 + 195 == 230
-
-    assert broadcast_area_libraries.get_polygons_for_areas_lat_long('ctry19-E92000001')[0][0] == [
         55.811085, -2.034358  # https://goo.gl/maps/wsf2LUWzYinwydMk8
     ]
 
 
-def test_polygons_are_enclosed_unless_asked_not_to_be():
+def test_polygons_are_enclosed():
     england = broadcast_area_libraries.get('ctry19').get('ctry19-E92000001')
 
-    assert len(england.polygons) == len(england.unenclosed_polygons)
+    first_polygon = england.polygons.as_coordinate_pairs_lat_long[0]
-
-    first_polygon = england.polygons[0]
     assert first_polygon[0] != first_polygon[1] != first_polygon[2]
     assert first_polygon[0] == first_polygon[-1]
 
-    first_polygon_unenclosed = england.unenclosed_polygons[0]
-    assert first_polygon_unenclosed[0] == first_polygon[0]
-    assert first_polygon_unenclosed[-1] != first_polygon[-1]
-    assert first_polygon_unenclosed[-1] == first_polygon[-2]
-
 
 def test_lat_long_order():
 
-    lat_long = broadcast_area_libraries.get_polygons_for_areas_lat_long('ctry19-E92000001')
+    england = broadcast_area_libraries.get_areas('ctry19-E92000001')[0]
-    long_lat = broadcast_area_libraries.get_polygons_for_areas_long_lat('ctry19-E92000001')
+
+    lat_long = england.polygons.as_coordinate_pairs_lat_long
+    long_lat = england.polygons.as_coordinate_pairs_long_lat
+
     assert len(lat_long[0]) == len(long_lat[0]) == 2082  # Coordinates in polygon
     assert len(lat_long[0][0]) == len(long_lat[0][0]) == 2  # Axes in coordinates
     assert lat_long[0][0] == list(reversed(long_lat[0][0]))

tests/app/models/test_broadcast_message.py

@@ -18,13 +18,19 @@ def test_simple_polygons(fake_uuid):
     ))
 
     assert [
-        [len(polygon) for polygon in broadcast_message.polygons],
+        [
-        [len(polygon) for polygon in broadcast_message.simple_polygons],
+            len(polygon)
+            for polygon in broadcast_message.polygons.as_coordinate_pairs_lat_long
+        ],
+        [
+            len(polygon)
+            for polygon in broadcast_message.simple_polygons.as_coordinate_pairs_lat_long
+        ],
     ] == [
         # One polygon for each area
         [27, 31],
         # Because the areas are close to each other, the simplification
         # and unioning process results in a single polygon with fewer
         # total coordinates
-        [34],
+        [55],
     ]