Update utils to do linear transformation of polygons

Brings in https://github.com/alphagov/notifications-utils/pull/889/files

At the moment, we are not doing any transformation of features before
applying geometric algorithms to them. This is, in effect, assuming that
the earth is flat.

This new version of utils implements the transformation of our polygons
to a Cartesian plane. In other words, it converts them from being
defined in spherical degrees to metres.

For the admin app this means we need to convert places where the code
expects things to be measured in degrees to work in metres instead.

app/__init__.py

@@ -76,6 +76,7 @@ from app.formatters import (
     recipient_count,
     recipient_count_label,
     round_to_significant_figures,
+    square_metres_to_square_miles,
     valid_phone_number,
 )
 from app.models.organisation import Organisation
@@ -572,6 +573,7 @@ def add_template_filters(application):
         message_count_noun,
         format_mobile_network,
         format_yes_no,
+        square_metres_to_square_miles,
     ]:
         application.add_template_filter(fn)
 

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

@@ -44,14 +44,30 @@ def simplify_geometry(feature):
 
 
 def clean_up_invalid_polygons(polygons, indent="    "):
+    """
+    This function expects a list of lists of coordinates defined in degrees
+    """
     for index, polygon in enumerate(polygons):
         shapely_polygon = Polygon(polygon)
 
-        if shapely_polygon.is_valid:
+        # Some of our data has points which are incredibly close
+        # together. In some cases they are close enough to be duplicates
+        # at a given precision, which makes an invalid topology. In
+        # other cases they are close enough that, when converting from
+        # one coordinate system to another, they shift about enough to
+        # create self-intersection. The fix in both cases is to reduce
+        # the precision of the coordinates and then apply simplification
+        # with a tolerance of 0.
+        simplified_polygon = wkt.loads(wkt.dumps(
+            shapely_polygon,
+            rounding_precision=Polygons.output_precision_in_decimal_places - 1
+        )).simplify(0)
+
+        if simplified_polygon.is_valid:
             print(  # noqa: T001
                 f"{indent}Polygon {index + 1}/{len(polygons)} is valid"
             )
-            yield polygon
+            yield simplified_polygon
 
         else:
             invalid_polygons.append(shapely_polygon)
@@ -59,7 +75,7 @@ def clean_up_invalid_polygons(polygons, indent="    "):
             # We’ve found polygons where all the points line up, so they
             # don’t have an area. They wouldn’t contribute to a broadcast
             # so we can ignore them.
-            if shapely_polygon.area == 0:
+            if simplified_polygon.area == 0:
                 print(  # noqa: T001
                     f"{indent}Polygon {index + 1}/{len(polygons)} has 0 area, skipping"
                 )
@@ -69,18 +85,6 @@ def clean_up_invalid_polygons(polygons, indent="    "):
                 f"{indent}Polygon {index + 1}/{len(polygons)} needs fixing..."
             )
 
-            # The simplest kind of invalid polygon is one that has two
-            # duplicate points in a row at a given precision, so the
-            # first thing to try is removing those points using
-            # simplification with a tolerance of 0
-            polygon_with_duplicate_points_removed = wkt.loads(
-                wkt.dumps(shapely_polygon, rounding_precision=5)
-            ).simplify(0)
-
-            if polygon_with_duplicate_points_removed.is_valid:
-                yield polygon_with_duplicate_points_removed
-                continue
-
             # Buffering with a size of 0 is a trick to make valid
             # geometries from polygons that self intersect
             buffered = shapely_polygon.buffer(0)
@@ -101,7 +105,7 @@ def clean_up_invalid_polygons(polygons, indent="    "):
             # Make sure the polygon is now valid, and that we haven’t
             # drastically transformed the polygon by ‘fixing’ it
             assert fixed_polygon.is_valid
-            assert isclose(fixed_polygon.area, polygon.area, rel_tol=0.001)
+            assert isclose(fixed_polygon.area, shapely_polygon.area, rel_tol=0.001)
 
             print(  # noqa: T001
                 f"{indent}Polygon {index + 1}/{len(polygons)} fixed!"
@@ -115,14 +119,22 @@ def polygons_and_simplified_polygons(feature):
         # cheat and shortcut out
         return [], []
 
-    polygons = Polygons(simplify_geometry(feature))
+    raw_polygons = simplify_geometry(feature)
-    polygons = list(clean_up_invalid_polygons(polygons))
+    clean_raw_polygons = [
-    polygons = Polygons(polygons)
+        [[x, y] for x, y in polygon.exterior.coords]
+        for polygon in clean_up_invalid_polygons(raw_polygons)
+    ]
+    polygons = Polygons(clean_raw_polygons)
 
     full_resolution = polygons.remove_too_small
     smoothed = full_resolution.smooth
     simplified = smoothed.simplify
 
+    if not (len(full_resolution) or len(simplified)):
+        raise RuntimeError(
+            'Polygon of 0 size found'
+        )
+
     print(  # noqa: T001
         f'    Original:{full_resolution.point_count: >5} points'
         f'    Smoothed:{smoothed.point_count: >5} points'

app/broadcast_areas/models.py

@@ -7,6 +7,8 @@ from notifications_utils.serialised_model import SerialisedModelCollection
 from rtreelib import Rect
 from werkzeug.utils import cached_property
 
+from app.formatters import square_metres_to_square_miles
+
 from .populations import CITY_OF_LONDON
 from .repo import BroadcastAreasRepository, rtree_index
 
@@ -55,13 +57,13 @@ class BaseBroadcastArea(ABC):
 
     @cached_property
     def simple_polygons_with_bleed(self):
-        return self.simple_polygons.bleed_by(self.estimated_bleed_in_degrees)
+        return self.simple_polygons.bleed_by(self.estimated_bleed_in_m)
 
     @cached_property
     def phone_density(self):
         if not self.polygons.estimated_area:
             return 0
-        return self.count_of_phones / self.polygons.estimated_area
+        return self.count_of_phones / square_metres_to_square_miles(self.polygons.estimated_area)
 
     @property
     def estimated_bleed_in_m(self):
@@ -72,14 +74,10 @@ class BaseBroadcastArea(ABC):
         range masts, so the typical bleed will be high (up to 5,000m).
         '''
         if self.phone_density < 1:
-            return Polygons.approx_bleed_in_degrees * Polygons.approx_metres_to_degree
+            return Polygons.approx_bleed_in_m
         estimated_bleed = 5_900 - (math.log(self.phone_density, 10) * 1_250)
         return max(500, min(estimated_bleed, 5000))
 
-    @property
-    def estimated_bleed_in_degrees(self):
-        return self.estimated_bleed_in_m / Polygons.approx_metres_to_degree
-
 
 class BroadcastArea(BaseBroadcastArea, SortableMixin):
 
@@ -123,7 +121,7 @@ class BroadcastArea(BaseBroadcastArea, SortableMixin):
     def count_of_phones(self):
         if self.id.endswith(CITY_OF_LONDON.WARDS):
             return CITY_OF_LONDON.DAYTIME_POPULATION * (
-                self.polygons.estimated_area / CITY_OF_LONDON.AREA_SQUARE_MILES
+                self.polygons.estimated_area / CITY_OF_LONDON.AREA_SQUARE_METRES
             )
         if self.sub_areas:
             return sum(area.count_of_phones for area in self.sub_areas)
@@ -169,7 +167,10 @@ class CustomBroadcastArea(BaseBroadcastArea):
         return Polygons(
             # Polygons in the DB are stored with the coordinate pair
             # order flipped – this flips them back again
-            Polygons(self._polygons).as_coordinate_pairs_lat_long
+            [
+                [[lat, long] for long, lat in polygon]
+                for polygon in self._polygons
+            ]
         )
 
     simple_polygons = polygons

app/broadcast_areas/populations.py

@@ -31,8 +31,9 @@ class CITY_OF_LONDON:
     )
     # https://data.london.gov.uk/blog/daytime-population-of-london-2014/
     DAYTIME_POPULATION = 553_000
-    # Approx area of the polygons we’re storing, not the actual area
+    # Exact area of the polygons we’re storing, which matches the 2.9km²
-    AREA_SQUARE_MILES = 1.78
+    # given by https://en.wikipedia.org/wiki/City_of_London
+    AREA_SQUARE_METRES = 2_885_598
 
 
 class BRYHER:

app/formatters.py

@@ -524,3 +524,7 @@ def format_yes_no(value, yes='Yes', no='No', none='No'):
     if value is None:
         return none
     return yes if value else no
+
+
+def square_metres_to_square_miles(area):
+    return area * 3.86e-7

app/models/broadcast_message.py

@@ -221,7 +221,8 @@ class BroadcastMessage(JSONModel):
         polygons = Polygons(
             list(itertools.chain(*(
                 getattr(area, area_attribute) for area in self.areas
-            )))
+            ))),
+            utm_crs=self.areas[0].polygons.utm_polygons.utm_crs,
         )
         if area_attribute != 'polygons' and len(self.areas) > 1:
             # We’re combining simplified polygons from multiple areas so we

app/templates/views/broadcast/macros/area-map.html

@@ -6,7 +6,7 @@
         <polygon points="25 5, 45 25, 25 45, 5 25" stroke="#0B0B0C" stroke-width="2" fill="#96C6E2" />
       </svg>
       <span class="govuk-visually-hidden">
-        An area of {{ (broadcast_message.simple_polygons.estimated_area)|round_to_significant_figures(1)|format_thousands }} square miles&nbsp;
+        An area of {{ (broadcast_message.simple_polygons.estimated_area)|square_metres_to_square_miles|round_to_significant_figures(1)|format_thousands }} square miles&nbsp;
       </span>
       Will get
       <span class="govuk-visually-hidden">
@@ -20,7 +20,7 @@
         <polygon points="25 5, 45 25, 25 45, 5 25" stroke="#005ea5" stroke-opacity="1" stroke-width="2" stroke-linecap="square" stroke-linejoin="round" stroke-dasharray="4,7.5,5,7.5,8,8,5,8,7.5,8,5,8,7,8,5,8,4" fill="#2B8CC4" fill-opacity="0.15" />
       </svg>
       <span class="govuk-visually-hidden">
-        An extra area of {{ (broadcast_message.simple_polygons_with_bleed.estimated_area - broadcast_message.simple_polygons.estimated_area)|round_to_significant_figures(1)|format_thousands }} square miles is&nbsp;
+        An extra area of {{ (broadcast_message.simple_polygons_with_bleed.estimated_area - broadcast_message.simple_polygons.estimated_area)|square_metres_to_square_miles|round_to_significant_figures(1)|format_thousands }} square miles is&nbsp;
       </span>
       Likely to get
       <span class="govuk-visually-hidden">

requirements.in

@@ -31,7 +31,7 @@ pyproj==3.2.1
 awscli-cwlogs>=1.4,<1.5
 itsdangerous==1.1.0  # pyup: <2
 
-notifications-utils @ git+https://github.com/alphagov/notifications-utils.git@48.0.0
+notifications-utils @ git+https://github.com/alphagov/notifications-utils.git@49.0.0
 govuk-frontend-jinja @ git+https://github.com/alphagov/govuk-frontend-jinja.git@v0.5.8-alpha
 
 # cryptography 3.4+ incorporates Rust code, which isn't supported on PaaS

requirements.txt

@@ -123,7 +123,7 @@ mistune==0.8.4
     # via notifications-utils
 notifications-python-client==6.3.0
     # via -r requirements.in
-notifications-utils @ git+https://github.com/alphagov/notifications-utils.git@48.0.0
+notifications-utils @ git+https://github.com/alphagov/notifications-utils.git@49.0.0
     # via -r requirements.in
 openpyxl==3.0.7
     # via pyexcel-xlsx
@@ -165,7 +165,9 @@ pyparsing==2.4.7
 pypdf2==1.26.0
     # via notifications-utils
 pyproj==3.2.1
-    # via -r requirements.in
+    # via
+    #   -r requirements.in
+    #   notifications-utils
 python-dateutil==2.8.1
     # via
     #   awscli-cwlogs

tests/app/broadcast_areas/custom_polygons.py

@@ -15,11 +15,11 @@ SKYE = [
 ]
 
 SANTA_A = [
-    [25.8890, 66.5500],
+    [66.5500, 25.8890],
-    [25.8890, 66.551],
+    [66.551, 25.8890],
-    [25.8910, 66.551],
+    [66.551, 25.8910],
-    [25.8910, 66.5500],
+    [66.5500, 25.8910],
-    [25.889, 66.55000],
+    [66.55000, 25.889],
 ]
 
 BURFORD = [

tests/app/broadcast_areas/test_broadcast_area.py

@@ -101,7 +101,7 @@ def test_has_polygons():
     assert len(scotland.polygons) == 195
 
     assert england.polygons.as_coordinate_pairs_lat_long[0][0] == [
-        55.811085, -2.034358  # https://goo.gl/maps/wsf2LUWzYinwydMk8
+        55.81108, -2.03436  # https://goo.gl/maps/HMFHGogohXdh9ggo8
     ]
 
 
@@ -276,24 +276,24 @@ def test_estimate_number_of_smartphones_for_population(
 @pytest.mark.parametrize('area, expected_phones_per_square_mile', (
     (
         # Islington (most dense in UK)
-        'lad20-E09000019', 21_348
+        'lad20-E09000019', 34_281
     ),
     (
         # Cordwainer Ward (City of London)
         # This is higher than Islington because we inflate the
         # popualtion to account for daytime workers
-        'wd20-E05009300', 310_674
+        'wd20-E05009300', 496_480
     ),
     (
         # Crewe East
-        'wd20-E05008621', 2_078),
+        'wd20-E05008621', 3_460),
     (
         # Eden (Cumbria, least dense in England)
-        'lad20-E07000030', 25.57
+        'lad20-E07000030', 44.12
     ),
     (
         # Highland (least dense in UK)
-        'lad20-S12000017', 4.40
+        'lad20-S12000017', 8.18
     ),
 ))
 def test_phone_density(
@@ -305,44 +305,40 @@ def test_phone_density(
     )
 
 
-@pytest.mark.parametrize('area, expected_bleed_in_m, expected_bleed_in_degrees', (
+@pytest.mark.parametrize('area, expected_bleed_in_m', (
     (
         # Islington (most dense in UK)
-        'lad20-E09000019', 500, 0.00449
+        'lad20-E09000019', 500
     ),
     (
         # Cordwainer Ward (City of London)
         # Special case because of inflated daytime population
-        'wd20-E05009300', 500, 0.00449
+        'wd20-E05009300', 500
     ),
     (
         # Crewe East
-        'wd20-E05008621', 1_752, 0.01574
+        'wd20-E05008621', 1_476
     ),
     (
         # Eden (Cumbria, least dense in England)
-        'lad20-E07000030', 4_140, 0.0372
+        'lad20-E07000030', 3_844
     ),
     (
         # Highland (least dense in UK)
-        'lad20-S12000017', 5_000, 0.0449
+        'lad20-S12000017', 4_759
     ),
     (
         # No population data available
-        'test-santa-claus-village-rovaniemi-a', 1_500, 0.01347
+        'test-santa-claus-village-rovaniemi-a', 1_500
     )
 ))
 def test_estimated_bleed(
-    area, expected_bleed_in_m, expected_bleed_in_degrees,
+    area, expected_bleed_in_m
 ):
     assert close_enough(
         broadcast_area_libraries.get_areas([area])[0].estimated_bleed_in_m,
         expected_bleed_in_m,
     )
-    assert close_enough(
-        broadcast_area_libraries.get_areas([area])[0].estimated_bleed_in_degrees,
-        expected_bleed_in_degrees,
-    )
 
 
 @pytest.mark.parametrize('polygon, expected_possible_overlaps, expected_count_of_phones', (
@@ -362,7 +358,7 @@ def test_estimated_bleed(
             'Stoke Bishop',
             'Windmill Hill',
         ],
-        73_119,
+        72_817,
     ),
     (
         SKYE,
@@ -372,7 +368,7 @@ def test_estimated_bleed(
             'Na Hearadh agus Ceann a Deas nan Loch',
             'Wester Ross, Strathpeffer and Lochalsh',
         ],
-        3_534,
+        3_413,
     ),
 ))
 def test_count_of_phones_for_custom_area(

tests/app/broadcast_areas/test_models.py

@@ -5,7 +5,7 @@ from tests.app.broadcast_areas.custom_polygons import BRISTOL, SANTA_A, SKYE
 
 
 @pytest.mark.parametrize(('simple_polygon', 'expected_wards_length'), [
-    (SKYE, 2),
+    (SKYE, 1),
     (BRISTOL, 12),
     (SANTA_A, 0)  # does not overlap with UK
 ])

tests/app/main/views/test_broadcast.py

@@ -837,8 +837,8 @@ def test_broadcast_page(
         'England Remove England',
         'Scotland Remove Scotland',
     ], [
-        'An area of 200,000 square miles Will get the alert',
+        'An area of 100,000 square miles Will get the alert',
-        'An extra area of 8,000 square miles is Likely to get the alert',
+        'An extra area of 6,000 square miles is Likely to get the alert',
         '40,000,000 phones estimated',
     ]),
     ([
@@ -854,8 +854,8 @@ def test_broadcast_page(
         'Penrith South Remove Penrith South',
         'Penrith West Remove Penrith West',
     ], [
-        'An area of 6 square miles Will get the alert',
+        'An area of 4 square miles Will get the alert',
-        'An extra area of 20 square miles is Likely to get the alert',
+        'An extra area of 10 square miles is Likely to get the alert',
         '9,000 to 10,000 phones',
     ]),
     ([
@@ -863,17 +863,17 @@ def test_broadcast_page(
     ], [
         'Islington Remove Islington',
     ], [
-        'An area of 10 square miles Will get the alert',
+        'An area of 6 square miles Will get the alert',
-        'An extra area of 5 square miles is Likely to get the alert',
+        'An extra area of 4 square miles is Likely to get the alert',
-        '200,000 to 500,000 phones',
+        '200,000 to 600,000 phones',
     ]),
     ([
         'ctyua19-E10000019',
     ], [
         'Lincolnshire Remove Lincolnshire',
     ], [
-        'An area of 4,000 square miles Will get the alert',
+        'An area of 2,000 square miles Will get the alert',
-        'An extra area of 700 square miles is Likely to get the alert',
+        'An extra area of 500 square miles is Likely to get the alert',
         '500,000 to 600,000 phones',
     ]),
     ([
@@ -882,8 +882,8 @@ def test_broadcast_page(
     ], [
         'Lincolnshire Remove Lincolnshire', 'North Yorkshire Remove North Yorkshire',
     ], [
-        'An area of 10,000 square miles Will get the alert',
+        'An area of 6,000 square miles Will get the alert',
-        'An extra area of 2,000 square miles is Likely to get the alert',
+        'An extra area of 1,000 square miles is Likely to get the alert',
         '1,000,000 phones estimated',
     ]),
 ))
@@ -936,7 +936,7 @@ def test_preview_broadcast_areas_page(
             [[7, 8], [9, 10], [11, 12]],
         ],
         [
-            'An area of 700 square miles Will get the alert',
+            'An area of 1,000 square miles Will get the alert',
             'An extra area of 1,000 square miles is Likely to get the alert',
             'Unknown number of phones',
         ]
@@ -944,17 +944,17 @@ def test_preview_broadcast_areas_page(
     (
         [BRISTOL],
         [
-            'An area of 7 square miles Will get the alert',
+            'An area of 4 square miles Will get the alert',
-            'An extra area of 6 square miles is Likely to get the alert',
+            'An extra area of 3 square miles is Likely to get the alert',
             '70,000 to 100,000 phones',
         ]
     ),
     (
         [SKYE],
         [
-            'An area of 3,000 square miles Will get the alert',
+            'An area of 2,000 square miles Will get the alert',
-            'An extra area of 800 square miles is Likely to get the alert',
+            'An extra area of 600 square miles is Likely to get the alert',
-            '4,000 phones estimated',
+            '3,000 to 4,000 phones',
         ]
     ),
 ))

tests/app/models/test_broadcast_message.py

@@ -45,7 +45,7 @@ def test_simple_polygons():
         # Because the areas are close to each other, the simplification
         # and unioning process results in a single polygon with fewer
         # total coordinates
-        [55],
+        [57],
     ]