Store simplifed polygons in the SQLite database

This commit does two things:
- uses our new polygon-simplifying library to process the polygons
  before storing them, rather than processing them in real time
- stores only the polygons in the database, rather than the whole
  GeoJSON feature, because we don’t need any of the other information
  about the feature

app/broadcast_areas/__init__.py

@@ -1,6 +1,5 @@
 import itertools
 
-import geojson
 from notifications_utils.serialised_model import SerialisedModelCollection
 from werkzeug.utils import cached_property
 
@@ -34,59 +33,13 @@ class BroadcastArea(SortableMixin):
     def __eq__(self, other):
         return self.id == other.id
 
-    @property
+    @cached_property
-    def _feature(self):
-        return BroadcastAreasRepository().get_feature_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):
-        # 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)
+        return BroadcastAreasRepository().get_polygons_for_area(self.id)
 
-    @property
+    @cached_property
-    def unenclosed_polygons(self):
-        return self._unenclosed_polygons(self.feature)
-
-    @property
     def simple_polygons(self):
-        return self._polygons(self.simple_feature)
+        return 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):

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

@@ -1,68 +1,51 @@
 #!/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, simplified.as_coordinate_pairs
+
+
 repo = BroadcastAreasRepository()
 
 repo.delete_db()
@@ -88,10 +71,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 +117,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 +146,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 +162,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

@@ -1,3 +1,5 @@
+import itertools
+
 from shapely.geometry import (
     JOIN_STYLE,
     GeometryCollection,
@@ -11,31 +13,32 @@ 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
 
-    # Ratio of how much to buffer for a shape of a given perimeter. For
+    # Controls how much buffer to add for a shape of a given perimeter.
-    # example `500` means 1m of buffer for every 500m of perimeter, or
+    # Smaller number means more buffering and a smoother shape. For
-    # 40m of buffer for a 5km square. This gives us control over how
+    # 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 = 500
+    perimeter_to_buffer_ratio = 360
 
     # Ratio of how much detail a shape of a given perimeter has once
-    # simplified. Smaller number means more less detail. For example
+    # simplified. Smaller number means less detail. For example `1000`
-    # `700` means that for a shape with a perimeter of 700m, the
+    # means that for a shape with a perimeter of 1000m, the simplified
-    # simplified line will never deviate more than 1m from the original.
+    # line will never deviate more than 1m from the original.
-    # Or for a 5km square, the line won’t deviate more than 17m. This
+    # 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 = 700
+    perimeter_to_simplification_ratio = 1_750
 
-    # The absolute smallest deviation (in metres) from the original we
+    # The threshold for removing very small areas from the map. These
-    # allow no matter how big/small the shape is. Allows us to still
+    # areas are likely glitches in  the data where the shoreline hasn’t
-    # remove a bit of detail even for small shapes, for example urban
+    # been subtracted from the land properly
-    # electoral wards.
+    minimum_area_size_square_metres = 50 ** 2
-    max_resolution = 5 / approx_metres_to_degree
 
     def __init__(self, polygons):
         if not polygons:
@@ -56,30 +59,36 @@ class Polygons():
             polygon.length for polygon in self
         )
 
-    @property
+    @cached_property
     def buffer_outward_in_degrees(self):
-        return self.perimeter_length / self.perimeter_to_buffer_ratio
+        return (
-
+            # If two areas are close enough that the distance between
-    @property
+            # them is less than the minimum bleed of a cell
-    def buffer_inward_in_degrees(self):
+            # broadcast then this joins them together. The aim is to
-        return self.buffer_outward_in_degrees - (
+            # reduce the total number of polygons in areas with many
-            # We don’t want to buffer all the way back in because there
+            # small shapes like Orkney or the Isles of Scilly.
-            # needs to be a bit off wiggle room for simplifying the
+            self.approx_bleed_in_degrees / 3
-            # polygon. Theoretically we need
+        ) + (
-            # `self.simplification_tolerance_in_degrees` wiggle room,
+            self.perimeter_length / self.perimeter_to_buffer_ratio
-            # but in practice some fraction of it is enough.
-            self.simplification_tolerance_in_degrees * 2 / 3
         )
 
-    @property
-    def simplification_tolerance_in_degrees(self):
-        shape_size_adjusted_resolution = (
-            self.perimeter_length / self.perimeter_to_simplification_ratio
-        )
-        return self.max_resolution + shape_size_adjusted_resolution
-
     @cached_property
-    def buffer_and_debuffer(self):
+    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,
@@ -89,7 +98,7 @@ class Polygons():
             for polygon in self
         ]
         unioned = union_polygons(buffered)
-        polygons_debuffered = [
+        debuffered = [
             polygon.buffer(
                 -1 * self.buffer_inward_in_degrees,
                 resolution=1,
@@ -97,7 +106,10 @@ class Polygons():
             )
             for polygon in unioned
         ]
-        return Polygons(polygons_debuffered)
+        flattened = list(itertools.chain(*[
+            flatten_polygons(polygon) for polygon in debuffered
+        ]))
+        return Polygons(flattened)
 
     @cached_property
     def simplify(self):
@@ -117,7 +129,18 @@ class Polygons():
             for polygon in self
         ]))
 
-    @property
+    @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(self):
         return [
             [
@@ -126,6 +149,18 @@ class Polygons():
             for p in self
         ]
 
+    @cached_property
+    def as_unenclosed_coordinate_pairs(self):
+        # 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.as_coordinate_pairs
+        ]
+
+    @cached_property
+    def point_count(self):
+        return len(list(itertools.chain(*self.as_coordinate_pairs)))
+
 
 def flatten_polygons(polygons):
     if isinstance(polygons, GeometryCollection):

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

@@ -79,6 +79,17 @@ class BroadcastMessage(JSONModel):
             )
         )
 
+    @cached_property
+    def simple_polygons(self):
+        polygons = Polygons(
+            broadcast_area_libraries.get_simple_polygons_for_areas_lat_long(
+                *self._dict['areas']
+            )
+        )
+        # If we’ve added multiple areas then we need to re-simplify the
+        # combined shapes to keep the point count down
+        return polygons.smooth.simplify if len(self.areas) > 1 else polygons
+
     @property
     def template(self):
         response = service_api_client.get_service_template(

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

@@ -36,7 +36,7 @@
       attribution: '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors'
     }).addTo(mymap);
 
-    {% for polygon in broadcast_message.polygons.buffer_and_debuffer.simplify.bleed.as_coordinate_pairs %}
+    {% for polygon in broadcast_message.simple_polygons.bleed.as_coordinate_pairs %}
       polygons.push(
         L.polygon({{polygon}}, {
           opacity: 0.5,
@@ -48,7 +48,7 @@
       );
     {% endfor %}
 
-    {% for polygon in broadcast_message.polygons.buffer_and_debuffer.simplify.as_coordinate_pairs %}
+    {% for polygon in broadcast_message.simple_polygons.as_coordinate_pairs %}
       polygons.push(
         L.polygon({{polygon}}, {
           opacity: 0.1,

tests/app/broadcast_areas/test_broadcast_area.py

@@ -117,13 +117,13 @@ def test_has_polygons():
 def test_polygons_are_enclosed_unless_asked_not_to_be():
     england = broadcast_area_libraries.get('ctry19').get('ctry19-E92000001')
 
-    assert len(england.polygons) == len(england.unenclosed_polygons)
+    assert len(england.polygons) == len(england.polygons.as_unenclosed_coordinate_pairs)
 
-    first_polygon = england.polygons[0]
+    first_polygon = england.polygons[0].as_coordinate_pairs
     assert first_polygon[0] != first_polygon[1] != first_polygon[2]
     assert first_polygon[0] == first_polygon[-1]
 
-    first_polygon_unenclosed = england.unenclosed_polygons[0]
+    first_polygon_unenclosed = england.polygons[0].as_unenclosed_coordinate_pairs
     assert first_polygon_unenclosed[0] == first_polygon[0]
     assert first_polygon_unenclosed[-1] != first_polygon[-1]
     assert first_polygon_unenclosed[-1] == first_polygon[-2]