broadcast-areas: move broadcast areas into app

Signed-off-by: Toby Lorne <toby.lornewelch-richards@digital.cabinet-office.gov.uk>
This commit is contained in:
Toby Lorne
2020-08-10 12:50:20 +01:00
parent 3243faa9a9
commit 488a5440cd
12 changed files with 18716 additions and 1 deletions
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because it is too large Load Diff
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+154
View File
@@ -0,0 +1,154 @@
import geojson
import itertools
from werkzeug.utils import cached_property
from notifications_utils.serialised_model import SerialisedModelCollection
from .repo import BroadcastAreasRepository
class SortableMixin:
def __repr__(self):
return f'{self.__class__.__name__}(<{self.id}>)'
def __lt__(self, other):
# Implementing __lt__ means any classes inheriting from this
# method are sortable
return self.name < other.name
class GetItemByIdMixin:
def get(self, id):
for item in self:
if item.id == id:
return item
raise KeyError(id)
class BroadcastArea(SortableMixin):
def __init__(self, row):
id, name, feature, simple_feature = row
self.id = id
self.name = name
self._feature = feature
self._simple_feature = simple_feature
for coordinates in self.polygons:
if coordinates[0] != coordinates[-1]:
# The CAP XML format requires shapes to be closed
raise ValueError(
f'Area {self.name} is not a closed shape '
f'({coordinates[0]}, {coordinates[-1]})'
)
def __eq__(self, other):
return self.id == other.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)
@property
def unenclosed_polygons(self):
return self._unenclosed_polygons(self.feature)
@property
def simple_polygons(self):
return self._polygons(self.simple_feature)
@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):
return [
BroadcastArea(row)
for row in BroadcastAreasRepository().get_all_areas_for_group(self.id)
]
class BroadcastAreaLibrary(SerialisedModelCollection, SortableMixin, GetItemByIdMixin):
model = BroadcastArea
def __init__(self, row):
id, name, is_group = row
self.id = id
self.name = name
self.is_group = bool(is_group)
def get_examples(self):
return BroadcastAreasRepository().get_library_description(self.id)
@property
def items(self):
return BroadcastAreasRepository().get_all_areas_for_library(self.id)
class BroadcastAreaLibraries(SerialisedModelCollection, GetItemByIdMixin):
model = BroadcastAreaLibrary
def __init__(self):
self.libraries = BroadcastAreasRepository().get_libraries()
self.items = self.libraries
def get_areas(self, *area_ids):
# allow people to call `get_areas('a', 'b') or get_areas(['a', 'b'])`
if len(area_ids) == 1 and isinstance(area_ids[0], list):
area_ids = area_ids[0]
areas = BroadcastAreasRepository().get_areas(area_ids)
return [BroadcastArea(area) for area in areas]
def get_polygons_for_areas_long_lat(self, *area_ids):
return list(itertools.chain(*(
area.polygons
for area in self.get_areas(*area_ids)
)))
def get_polygons_for_areas_lat_long(self, *area_ids):
return [
[[long, lat] for lat, long in polygon]
for polygon in self.get_polygons_for_areas_long_lat(*area_ids)
]
broadcast_area_libraries = BroadcastAreaLibraries()
@@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:d92ef0b11c3586fa694f9c8afe93e64de189cb816a10264a459099f5befd51e9
size 81362944
+167
View File
@@ -0,0 +1,167 @@
#!/usr/bin/env python
from copy import deepcopy
import geojson
from pathlib import Path
import shapely.geometry as sgeom
from notifications_utils.safe_string import make_string_safe_for_id
from repo import BroadcastAreasRepository
package_path = Path(__file__).resolve().parent
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
elif feature["type"] == "MultiPolygon":
feature["coordinates"] = [
simplify_polygon(polygon)
for polygon in feature["coordinates"]
]
return feature
else:
raise Exception("Unknown type: {}".format(feature["type"]))
repo = BroadcastAreasRepository()
repo.delete_db()
repo.create_tables()
simple_datasets = [
("Countries", "ctry19cd", "ctry19nm"),
("Regions of England", "rgn18cd", "rgn18nm"),
("Counties and Unitary Authorities in England and Wales", "ctyua16cd", "ctyua16nm"),
]
for dataset_name, id_field, name_field in simple_datasets:
filepath = package_path / "{}.geojson".format(dataset_name)
dataset_id = make_string_safe_for_id(dataset_name)
dataset_geojson = geojson.loads(filepath.read_text())
repo.insert_broadcast_area_library(dataset_id, dataset_name, False)
for feature in dataset_geojson["features"]:
f_id = dataset_id + "-" + feature["properties"][id_field]
f_name = feature["properties"][name_field]
print(f_name) # noqa: T001
simple_feature = deepcopy(feature)
simple_feature["geometry"] = simplify_geometry(simple_feature["geometry"])
repo.insert_broadcast_areas([[
f_id, f_name,
dataset_id, None,
feature, simple_feature,
]])
# https://geoportal.statistics.gov.uk/datasets/wards-may-2020-boundaries-uk-bgc
# Converted to geojson manually from SHP because of GeoJSON download limits
wards_filepath = package_path / "Electoral Wards May 2020.geojson"
# http://geoportal.statistics.gov.uk/datasets/ward-to-westminster-parliamentary-constituency-to-local-authority-district-december-2019-lookup-in-the-united-kingdom/data
las_filepath = package_path / "Electoral Wards and Local Authorities 2020.geojson"
ward_code_to_la_mapping = {
f["properties"]["WD19CD"]: f["properties"]["LAD19NM"]
for f in geojson.loads(las_filepath.read_text())["features"]
}
ward_code_to_la_id_mapping = {
f["properties"]["WD19CD"]: f["properties"]["LAD19CD"]
for f in geojson.loads(las_filepath.read_text())["features"]
}
dataset_name = "Electoral Wards of the United Kingdom"
dataset_id = make_string_safe_for_id(dataset_name)
repo.insert_broadcast_area_library(dataset_id, dataset_name, True)
areas_to_add = []
for f in geojson.loads(wards_filepath.read_text())["features"]:
ward_code = f["properties"]["wd20cd"]
ward_name = f["properties"]["wd20nm"]
ward_id = dataset_id + "-" + ward_code
print(ward_name) # noqa: T001
try:
la_id = dataset_id + "-" + ward_code_to_la_id_mapping[ward_code]
la_name = ward_code_to_la_mapping[ward_code]
sf = deepcopy(f)
sf["geometry"] = simplify_geometry(sf["geometry"])
areas_to_add.append([
ward_id, ward_name,
dataset_id, la_id,
f, sf
])
except KeyError:
print("Skipping", ward_code, ward_name) # noqa: T001
repo.insert_broadcast_areas(areas_to_add)
areas_to_add = []
las_filepath = package_path / "Local Authorities May 2020.geojson"
for feature in geojson.loads(las_filepath.read_text())["features"]:
la_id = feature["properties"]["lad20cd"]
group_name = feature["properties"]["lad20nm"]
print(group_name) # noqa: T001
group_id = dataset_id + "-" + la_id
simple_feature = deepcopy(feature)
simple_feature["geometry"] = simplify_geometry(simple_feature["geometry"])
areas_to_add.append([
group_id, group_name,
dataset_id, None,
feature, simple_feature
])
repo.insert_broadcast_areas(areas_to_add)
+99
View File
@@ -0,0 +1,99 @@
#!/usr/bin/env python
import geojson
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
from random import sample
from notifications_utils.safe_string import make_string_safe_for_id
from repo import BroadcastAreasRepository
def main():
print() # noqa: T001
print("pick a library") # noqa: T001
for library in BroadcastAreasRepository().get_libraries():
print(" ", library) # noqa: T001
library = input("> ")
lid = make_string_safe_for_id(library)
features = []
simple_features = []
inp = ""
while True:
print() # noqa: T001
print("pick an area, or press enter to skip") # noqa: T001
all_areas = BroadcastAreasRepository().get_all_areas_for_library(lid)
some_areas = sample(all_areas, min(len(all_areas), 25))
for area in some_areas:
print(" ", area[0], area[1]) # noqa: T001
inp = input("> ")
if inp == "":
break
aid = inp.strip()
area = BroadcastAreasRepository().get_areas([aid])[0]
feature = area[-2]
feature_shape = sgeom.shape(geojson.loads(feature)["geometry"])
features.append(feature_shape)
simple_feature = area[-1]
simple_feature_shape = sgeom.shape(geojson.loads(simple_feature)["geometry"])
simple_features.append(simple_feature_shape)
print() # noqa: T001
print("Plotting") # noqa: T001
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.set_extent([-20, 5, 40, 60], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_geometries(
features,
ccrs.PlateCarree(),
facecolor='#00ff00',
alpha=0.25,
)
ax.add_geometries(
simple_features,
ccrs.PlateCarree(),
facecolor='#0000ff',
alpha=0.25,
)
ax.scatter(
[
p[0]
for f in simple_features
for geom in (f.geoms if hasattr(f, 'geoms') else [f])
for p in geom.exterior.coords
],
[
p[1]
for f in simple_features
for geom in (f.geoms if hasattr(f, 'geoms') else [f])
for p in geom.exterior.coords
],
transform=ccrs.PlateCarree(),
)
plt.show()
if __name__ == '__main__':
main()
+198
View File
@@ -0,0 +1,198 @@
import geojson
import os
from pathlib import Path
import sqlite3
class BroadcastAreasRepository(object):
def __init__(self):
self.database = Path(__file__).resolve().parent / 'broadcast-areas.sqlite3'
def conn(self):
return sqlite3.connect(str(self.database))
def delete_db(self):
os.remove(str(self.database))
def create_tables(self):
with self.conn() as conn:
conn.execute("""
CREATE TABLE broadcast_area_libraries (
id TEXT PRIMARY KEY,
name TEXT NOT NULL,
is_group BOOLEAN NOT NULL
)""")
conn.execute("""
CREATE TABLE broadcast_area_library_groups (
id TEXT PRIMARY KEY,
name TEXT NOT NULL,
broadcast_area_library_id TEXT NOT NULL
)""")
conn.execute("""
CREATE TABLE broadcast_areas (
id TEXT PRIMARY KEY,
name TEXT NOT NULL,
broadcast_area_library_id TEXT NOT NULL,
broadcast_area_library_group_id TEXT,
feature_geojson TEXT NOT NULL,
simple_feature_geojson TEXT NOT NULL,
FOREIGN KEY (broadcast_area_library_id)
REFERENCES broadcast_area_libraries(id),
FOREIGN KEY (broadcast_area_library_group_id)
REFERENCES broadcast_area_library_groups(id)
)""")
conn.execute("""
CREATE INDEX broadcast_areas_broadcast_area_library_id
ON broadcast_areas (broadcast_area_library_id);
""")
conn.execute("""
CREATE INDEX broadcast_areas_broadcast_area_library_group_id
ON broadcast_areas (broadcast_area_library_group_id);
""")
def insert_broadcast_area_library(self, id, name, is_group):
q = """
INSERT INTO broadcast_area_libraries (id, name, is_group)
VALUES (?, ?, ?)
"""
with self.conn() as conn:
conn.execute(q, (id, name, is_group))
def insert_broadcast_areas(self, areas):
q = """
INSERT INTO broadcast_areas (
id, name,
broadcast_area_library_id, broadcast_area_library_group_id,
feature_geojson, simple_feature_geojson
)
VALUES (?, ?, ?, ?, ?, ?)
"""
with self.conn() as conn:
for id, name, area_id, group, feature, simple_feature in areas:
conn.execute(q, (
id, name,
area_id, group,
geojson.dumps(feature), geojson.dumps(simple_feature),
))
def query(self, sql, *args):
with self.conn() as conn:
cursor = conn.cursor()
cursor.execute(sql, (*args,))
return cursor.fetchall()
def get_libraries(self):
q = "SELECT id, name, is_group FROM broadcast_area_libraries"
results = self.query(q)
libraries = [(row[0], row[1], row[2]) for row in results]
return sorted(libraries)
def get_library_description(self, library_id):
q = """
WITH
areas AS (SELECT * FROM broadcast_areas
WHERE broadcast_area_library_id = ?),
area_count AS (SELECT COUNT(*) AS c FROM areas),
subset_area_count AS (SELECT c - 4 FROM area_count),
some_area_names AS (SELECT name FROM areas LIMIT 100),
some_shuffled_area_names AS (
SELECT name FROM some_area_names ORDER BY RANDOM()
),
description_area_names AS (
SELECT name FROM some_shuffled_area_names LIMIT 4
),
description_areas_joined AS (
SELECT GROUP_CONCAT(name, ", ") FROM description_area_names
)
SELECT
CASE (SELECT * FROM subset_area_count)
WHEN 0 THEN
(SELECT * FROM description_areas_joined)
ELSE
(SELECT * FROM description_areas_joined)
|| ", and "
|| (SELECT * FROM subset_area_count)
|| " more…"
END
"""
description = self.query(q, library_id)[0][0]
return description
def get_areas(self, *area_ids):
with self.conn() as conn:
cursor = conn.cursor()
q = """
SELECT id, name, feature_geojson, simple_feature_geojson
FROM broadcast_areas
WHERE id IN ({})
""".format(("?," * len(*area_ids))[:-1])
cursor.execute(q, *area_ids)
results = cursor.fetchall()
areas = [
(row[0], row[1], row[2], row[3])
for row in results
]
return areas
def get_all_areas_for_library(self, library_id):
q = """
SELECT id, name, feature_geojson, simple_feature_geojson
FROM broadcast_areas
WHERE broadcast_area_library_id = ?
AND broadcast_area_library_group_id IS NULL
"""
results = self.query(q, library_id)
areas = [
(row[0], row[1], row[2], row[3])
for row in results
]
return areas
def get_all_areas_for_group(self, group_id):
q = """
SELECT id, name, feature_geojson, simple_feature_geojson
FROM broadcast_areas
WHERE broadcast_area_library_group_id = ?
"""
results = self.query(q, group_id)
areas = [
(row[0], row[1], row[2], row[3])
for row in results
]
return areas
def get_all_groups_for_library(self, library_id):
q = """
SELECT id, name
FROM broadcast_areas
WHERE broadcast_area_library_group_id = NULL
AND broadcast_area_library_id = ?
"""
results = self.query(q, library_id)
areas = [
(row[0], row[1])
for row in results
]
return areas