Don’t combine polygons with different projections

When we combine simplified polygons from different areas we try to use
the polygons that are defined in metres as a speed optimisation.

This doesn’t work when those polygons are in several different
projections because the distance in metres is relative to the edge of
the projection. Just naively choosing one projection means that some
of the polygons will shift position on the earth by 100s of kilometers
once they are converted back to degrees.

To avoid this we need to:
- check for situations where we have polygons in multiple different
  projections
- transforms these polygons back into degrees and let `Polygons` pick a
  common projection for all of them before doing any
  merging/smoothing/etc.
This commit is contained in:
Chris Hill-Scott
2021-12-07 12:35:45 +00:00
parent 787cb3ef1f
commit 99a8a4fd48
2 changed files with 67 additions and 6 deletions

View File

@@ -218,12 +218,33 @@ class BroadcastMessage(JSONModel):
)
def get_polygons_from_areas(self, area_attribute):
polygons = Polygons(
list(itertools.chain(*(
getattr(area, area_attribute) for area in self.areas
))),
utm_crs=self.areas[0].simple_polygons.utm_polygons.utm_crs,
)
areas_polygons = [
getattr(area, area_attribute) for area in self.areas
]
coordinate_reference_systems = {
polygons.utm_crs for polygons in areas_polygons
}
if len(coordinate_reference_systems) == 1:
# All our polygons are defined in the same coordinate
# reference system so we just have to flatten the list and
# say which coordinate reference system we are using
polygons = Polygons(
list(itertools.chain(*areas_polygons)),
utm_crs=next(iter(coordinate_reference_systems)),
)
else:
# Our polygons are in different coordinate reference systems
# We need to convert them back to degrees and make a new
# instance of `Polygons` which will determine a common
# coordinate reference system
polygons = Polygons(
list(itertools.chain(*(
area_polygon.as_wgs84_coordinates
for area_polygon in areas_polygons
)))
)
if area_attribute != 'polygons' and len(self.areas) > 1:
# Were combining simplified polygons from multiple areas so we
# need to re-simplify the combined polygons to keep the point

View File

@@ -91,3 +91,43 @@ def test_areas_treats_missing_ids_as_custom_broadcast(notify_admin):
assert len(list(broadcast_message.areas)) == 2
assert type(broadcast_message.areas) == CustomBroadcastAreas
@pytest.mark.parametrize('area_ids, approx_bounds', (
([
'ctry19-N92000002', # Northern Ireland (UTM zone 29N)
'ctry19-W92000004', # Wales (UTM zone 30N)
], [
-8.2, 51.5, -2.1, 55.1
]),
([
'lad20-E06000031', # Peterborough (UTM zone 30N)
'lad20-E07000146', # Kings Lyn and West Norfolk (UTM zone 31N)
], [
-0.5, 52.5, 0.8, 53.0
]),
([
'wd20-E05009372', # Hackney Central (UTM zone 30N)
'wd20-E05009374', # Hackney Wick (UTM zone 30N)
], [
-0.1, 51.5, -0.0, 51.6
]),
([
'wd20-E05009372', # Hackney Central (UTM zone 30N)
'test-santa-claus-village-rovaniemi-a', # (UTM zone 35N)
], [
-0.1, 51.5, 25.9, 66.6
]),
))
def test_combining_multiple_areas_keeps_same_bounds(area_ids, approx_bounds):
broadcast_message = BroadcastMessage(broadcast_message_json(
areas={'ids': area_ids}
))
assert [
round(coordinate, 1) for coordinate in broadcast_message.polygons.bounds
] == [
round(coordinate, 1) for coordinate in broadcast_message.simple_polygons.bounds
] == (
approx_bounds
)