Make simplification of polygons more sophisticated

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
This commit is contained in:
Chris Hill-Scott
2020-08-24 12:51:23 +01:00
parent 19ce1bd43a
commit 3470c5bb31
3 changed files with 164 additions and 28 deletions

View File

@@ -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 wont 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 dont 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))

View File

@@ -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(

View File

@@ -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.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