From 3470c5bb31c71a5744c667e572a965adcfcc6db6 Mon Sep 17 00:00:00 2001 From: Chris Hill-Scott Date: Mon, 24 Aug 2020 12:51:23 +0100 Subject: [PATCH] Make simplification of polygons more sophisticated MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Simplifying polygons means reducing the number of points used to render them. This commit implements simplification such that, for any given input polygons, the combined point count of the simplified polygons is less than 100. When simplifying the polygons we are trying to get the smallest number of points while meeting these two rules: 1. No part of the area the user has chosen can be cut off 2. The area of the simplified polygon should be as small as possible This commit introduces two techniques we weren’t using before: 1. Dilating and eroding the area to fill in concave details of the shape, like inlets and harbours[1] 2. Making the simplification threshold proportionate to the perimeter of all polygons, so bigger and crinklier polygons get more simplification applied It also shows the estimated bleed as a separate polygon. This lets us make it bigger (so it’s more closer the the approximate bleed) without having to send a bigger area to the CBC and compounding the amount of actual bleed. 1. Inspired by this blog post about ‘removing the crinkley bits’ from Vancouver Island: http://blog.cleverelephant.ca/2010/11/removing-complexities.html --- app/broadcast_areas/polygons.py | 142 ++++++++++++++++++ app/models/broadcast_message.py | 30 +--- .../views/broadcast/preview-areas.html | 20 ++- 3 files changed, 164 insertions(+), 28 deletions(-) create mode 100644 app/broadcast_areas/polygons.py diff --git a/app/broadcast_areas/polygons.py b/app/broadcast_areas/polygons.py new file mode 100644 index 000000000..276dc278c --- /dev/null +++ b/app/broadcast_areas/polygons.py @@ -0,0 +1,142 @@ +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 + + # 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 + # example `500` means 1m of buffer for every 500m of perimeter, or + # 40m 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 + + # Ratio of how much detail a shape of a given perimeter has once + # simplified. Smaller number means more less detail. For example + # `700` means that for a shape with a perimeter of 700m, the + # simplified line will never deviate more than 1m from the original. + # Or for a 5km square, the line won’t deviate more than 17m. This + # gives us approximate control over the total number of points. + perimeter_to_simplification_ratio = 700 + + # The absolute smallest deviation (in metres) from the original we + # allow no matter how big/small the shape is. Allows us to still + # remove a bit of detail even for small shapes, for example urban + # electoral wards. + max_resolution = 5 / approx_metres_to_degree + + 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] + + @cached_property + def perimeter_length(self): + return sum( + polygon.length for polygon in self + ) + + @property + def buffer_outward_in_degrees(self): + return self.perimeter_length / self.perimeter_to_buffer_ratio + + @property + def buffer_inward_in_degrees(self): + return self.buffer_outward_in_degrees - ( + # We don’t want to buffer all the way back in because there + # needs to be a bit off wiggle room for simplifying the + # polygon. Theoretically we need + # `self.simplification_tolerance_in_degrees` wiggle room, + # 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): + buffered = [ + polygon.buffer( + self.buffer_outward_in_degrees, + resolution=4, + join_style=JOIN_STYLE.round, + ) + for polygon in self + ] + unioned = union_polygons(buffered) + polygons_debuffered = [ + polygon.buffer( + -1 * self.buffer_inward_in_degrees, + resolution=1, + join_style=JOIN_STYLE.bevel, + ) + for polygon in unioned + ] + return Polygons(polygons_debuffered) + + @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 + ])) + + @property + def as_coordinate_pairs(self): + return [ + [ + [x, y] for x, y in p.exterior.coords + ] + for p in self + ] + + +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)) diff --git a/app/models/broadcast_message.py b/app/models/broadcast_message.py index 88eba53b0..efc0bc16d 100644 --- a/app/models/broadcast_message.py +++ b/app/models/broadcast_message.py @@ -2,11 +2,10 @@ 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,31 +71,14 @@ 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( - *self._dict['areas'] + return Polygons( + broadcast_area_libraries.get_polygons_for_areas_lat_long( + *self._dict['areas'] + ) ) - @property - def simple_polygons(self): - simple_polygons = broadcast_area_libraries.get_simple_polygons_for_areas_lat_long( - *self._dict['areas'] - ) - unioned_polygons = unary_union([ - Polygon(i) for i in simple_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): response = service_api_client.get_service_template( diff --git a/app/templates/views/broadcast/preview-areas.html b/app/templates/views/broadcast/preview-areas.html index 026a80aa5..faf22d948 100644 --- a/app/templates/views/broadcast/preview-areas.html +++ b/app/templates/views/broadcast/preview-areas.html @@ -36,19 +36,31 @@ attribution: '© OpenStreetMap contributors' }).addTo(mymap); - {% for polygon in broadcast_message.simple_polygons %} + {% for polygon in broadcast_message.polygons.buffer_and_debuffer.simplify.bleed.as_coordinate_pairs %} 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.polygons.buffer_and_debuffer.simplify.as_coordinate_pairs %} + 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 %} polygons.push( L.polygon({{polygon}}, { color: '#0b0b0c', // $black