From bd71cbfb4fde2723b55d1b50f021ccc4f1f428f0 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 5 Mar 2021 10:54:59 +0000 Subject: [PATCH 01/57] Calculate raster grid shape --- seedpod_ground_risk/core/plot_server.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index dd5e4185..a58f89a8 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -337,3 +337,25 @@ def remove_layer(self, layer): def set_layer_order(self, layer_order): self.data_layer_order = layer_order + + def _get_raster_dimensions(self, bounds_poly: sg.Polygon, raster_resolution_m: float) -> Tuple[int, int]: + """ + Return a the (x,y) shape of a raster grid given its EPSG4326 envelope and desired raster resolution + :param bounds_poly: EPSG4326 Shapely Polygon specifying bounds + :param raster_resolution_m: raster resolution in metres + :return: 2-tuple of (width, height) + """ + + import pyproj + + if self._epsg4326_to_epsg3857_proj: + self._epsg4326_to_epsg3857_proj = pyproj.Transformer.from_crs(pyproj.CRS.from_epsg('4326'), + pyproj.CRS.from_epsg('3857'), + always_xy=True) + bounds = bounds_poly.bounds + + min_x, min_y = self._epsg4326_to_epsg3857_proj.transform(bounds[1], bounds[0]) + max_x, max_y = self._epsg4326_to_epsg3857_proj.transform(bounds[3], bounds[2]) + raster_width = abs(max_x-min_x) % raster_resolution_m + raster_height = abs(max_y - min_y) % raster_resolution_m + return raster_width, raster_height From 1c07eeb3c1ce09c9bcb2e140d15950cbf2a28f3d Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 5 Mar 2021 12:43:00 +0000 Subject: [PATCH 02/57] Modify data_layer api to include raster shape --- seedpod_ground_risk/core/plot_server.py | 21 ++++++++++++------- seedpod_ground_risk/layers/data_layer.py | 3 ++- seedpod_ground_risk/layers/osm_tag_layer.py | 4 ++-- .../layers/residential_layer.py | 4 ++-- seedpod_ground_risk/layers/roads_layer.py | 5 +++-- 5 files changed, 22 insertions(+), 15 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index a58f89a8..d7b20d4b 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -83,6 +83,9 @@ def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = No self._x_range, self._y_range = [-1.45, -1.35], [50.85, 50.95] + self.raster_resolution_m = 20 + + self._epsg4326_to_epsg3857_proj = None self._epsg3857_to_epsg4326_proj = None self._preload_started = False self._preload_complete = False @@ -196,12 +199,13 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) else: # Construct box around requested bounds bounds_poly = make_bounds_polygon(x_range, y_range) + raster_shape = self._get_raster_dimensions(bounds_poly, self.raster_resolution_m) # Ensure bounds are small enough to render without OOM or heat death of universe if bounds_poly.area < 0.2: from time import time t0 = time() - self.generate_layers(bounds_poly) + self.generate_layers(bounds_poly, raster_shape) self._progress_callback("Rendering new map...") plot = Overlay([res[0] for res in self._generated_data_layers.values()]) print("Generated all layers in ", time() - t0) @@ -270,8 +274,8 @@ def generate_layers(self, bounds_poly: sg.Polygon) -> NoReturn: layers = {} self._progress_callback('Generating layer data') with ThreadPoolExecutor() as pool: - layer_futures = [pool.submit(self.generate_layer, layer, bounds_poly, self._time_idx) for layer in - self.data_layers] + layer_futures = [pool.submit(self.generate_layer, layer, bounds_poly, raster_shape, self._time_idx) for + layer in self.data_layers] # Store generated layers as they are completed for future in as_completed(layer_futures): key, result = future.result() @@ -293,8 +297,9 @@ def generate_layers(self, bounds_poly: sg.Polygon) -> NoReturn: {k: layers[k] for k in layers.keys() if k not in self._generated_data_layers}) @staticmethod - def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, hour: int) -> Union[ + def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, raster_shape: Tuple[int, int], hour: int) -> Union[ Tuple[str, Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]], Tuple[str, None]]: + import shapely.ops as so from_cache = False @@ -308,7 +313,7 @@ def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, hour: int) -> Unio layer_bounds_poly = bounds_poly.difference(layer.cached_area) layer.cached_area = so.unary_union([layer.cached_area, bounds_poly]) try: - result = layer.key, layer.generate(layer_bounds_poly, from_cache=from_cache, hour=hour) + result = layer.key, layer.generate(layer_bounds_poly, raster_shape, from_cache=from_cache, hour=hour) return result except Exception as e: print(e) @@ -348,7 +353,7 @@ def _get_raster_dimensions(self, bounds_poly: sg.Polygon, raster_resolution_m: f import pyproj - if self._epsg4326_to_epsg3857_proj: + if self._epsg4326_to_epsg3857_proj is None: self._epsg4326_to_epsg3857_proj = pyproj.Transformer.from_crs(pyproj.CRS.from_epsg('4326'), pyproj.CRS.from_epsg('3857'), always_xy=True) @@ -356,6 +361,6 @@ def _get_raster_dimensions(self, bounds_poly: sg.Polygon, raster_resolution_m: f min_x, min_y = self._epsg4326_to_epsg3857_proj.transform(bounds[1], bounds[0]) max_x, max_y = self._epsg4326_to_epsg3857_proj.transform(bounds[3], bounds[2]) - raster_width = abs(max_x-min_x) % raster_resolution_m - raster_height = abs(max_y - min_y) % raster_resolution_m + raster_width = int(abs(max_x - min_x) // raster_resolution_m) + raster_height = int(abs(max_y - min_y) // raster_resolution_m) return raster_width, raster_height diff --git a/seedpod_ground_risk/layers/data_layer.py b/seedpod_ground_risk/layers/data_layer.py index 5265d368..d1293cf2 100644 --- a/seedpod_ground_risk/layers/data_layer.py +++ b/seedpod_ground_risk/layers/data_layer.py @@ -23,10 +23,11 @@ def __init__(self, key): self.cached_area = sg.Polygon() @abc.abstractmethod - def generate(self, bounds_polygon: sg.Polygon, from_cache: bool = False, **kwargs) -> Tuple[ + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> Tuple[ Geometry, np.ndarray, gpd.GeoDataFrame]: """ Generate the map of this layer. This is called asynchronously, so cannot access plot_server members. + :param raster_shape: :param shapely.geometry.Polygon bounds_polygon: the bounding polygon for which to generate the map :param bool from_cache: flag to indicate whether to use cached data to fulfill this request :return: an Overlay-able holoviews layer with specific options diff --git a/seedpod_ground_risk/layers/osm_tag_layer.py b/seedpod_ground_risk/layers/osm_tag_layer.py index f901c239..d1ef7a9d 100644 --- a/seedpod_ground_risk/layers/osm_tag_layer.py +++ b/seedpod_ground_risk/layers/osm_tag_layer.py @@ -26,8 +26,8 @@ def __init__(self, key, osm_tag, colour: str = None, def preload_data(self) -> NoReturn: pass - def generate(self, bounds_polygon: sg.Polygon, from_cache: bool = False, **kwargs) -> Tuple[ - Geometry, np.ndarray, gpd.GeoDataFrame]: + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> \ + Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: import geoviews as gv from holoviews.operation.datashader import rasterize diff --git a/seedpod_ground_risk/layers/residential_layer.py b/seedpod_ground_risk/layers/residential_layer.py index 9d7d460b..5330c7ad 100644 --- a/seedpod_ground_risk/layers/residential_layer.py +++ b/seedpod_ground_risk/layers/residential_layer.py @@ -24,8 +24,8 @@ def preload_data(self): print("Preloading Residential Layer") self.ingest_census_data() - def generate(self, bounds_polygon: sg.Polygon, from_cache: bool = False, **kwargs) -> Tuple[ - Geometry, np.ndarray, gpd.GeoDataFrame]: + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> \ + Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: import colorcet import datashader as ds from holoviews.operation.datashader import rasterize diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index 1738fad6..66795ba4 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -64,8 +64,9 @@ def preload_data(self) -> NoReturn: self._ingest_road_geometries() self._ingest_relative_traffic_variations() - def generate(self, bounds_polygon: sg.Polygon, from_cache: bool = False, hour: int = 0, **kwargs) -> Tuple[ - Geometry, np.ndarray, gpd.GeoDataFrame]: + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, + hour: int = 0, **kwargs) -> \ + Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: from holoviews.operation.datashader import rasterize import geoviews as gv import datashader as ds From aa84214555a4d3bb1f900a17183fe3fa383256cf Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 5 Mar 2021 12:43:31 +0000 Subject: [PATCH 03/57] Individual data layer raster shape implementations --- seedpod_ground_risk/layers/osm_tag_layer.py | 2 +- seedpod_ground_risk/layers/residential_layer.py | 7 ++++--- seedpod_ground_risk/layers/roads_layer.py | 4 ++-- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/seedpod_ground_risk/layers/osm_tag_layer.py b/seedpod_ground_risk/layers/osm_tag_layer.py index d1ef7a9d..7660f5f5 100644 --- a/seedpod_ground_risk/layers/osm_tag_layer.py +++ b/seedpod_ground_risk/layers/osm_tag_layer.py @@ -40,7 +40,7 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr self._landuse_polygons.geometry = self._landuse_polygons.to_crs('EPSG:27700') \ .buffer(self.buffer_dist).to_crs('EPSG:4326') polys = gv.Polygons(self._landuse_polygons).opts(style={'alpha': 0.8, 'color': self._colour}) - raster = rasterize(polys, + raster = rasterize(polys, width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) if self.blocking: diff --git a/seedpod_ground_risk/layers/residential_layer.py b/seedpod_ground_risk/layers/residential_layer.py index 5330c7ad..d16f8288 100644 --- a/seedpod_ground_risk/layers/residential_layer.py +++ b/seedpod_ground_risk/layers/residential_layer.py @@ -58,10 +58,11 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr buffered_df.geometry = buffered_df.to_crs('EPSG:27700') \ .buffer(self.buffer_dist).to_crs('EPSG:4326') buffered_polys = gv.Polygons(buffered_df, kdims=['Longitude', 'Latitude'], vdims=['name', 'population']) - raster = rasterize(buffered_polys, aggregator=ds.sum('population'), - x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) + raster = rasterize(buffered_polys, aggregator=ds.sum('population'), width=raster_shape[0], + height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), + dynamic=False) else: - raster = rasterize(gv_polys, aggregator=ds.sum('population'), + raster = rasterize(gv_polys, aggregator=ds.sum('population'), width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index 66795ba4..c4be86e1 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -83,7 +83,7 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr cmap=colorcet.CET_L18, color='population') bounds = bounds_polygon.bounds - raster = rasterize(points, aggregator=ds.mean('population'), + raster = rasterize(points, aggregator=ds.mean('population'), width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) return points, raster_grid, gpd.GeoDataFrame(tfc_df) @@ -145,7 +145,7 @@ def _ingest_road_geometries(self) -> NoReturn: if not self._roads_geometries.crs: self._roads_geometries = self._roads_geometries.set_crs('EPSG:27700') - def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = 25) -> List[ + def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = 20) -> List[ List[float]]: """ Interpolate traffic count values between count points along all roads. From a28673e2c7b031e733b075d3d5abfc1f417c79fb Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 5 Mar 2021 12:44:49 +0000 Subject: [PATCH 04/57] Handle variable raster shape during layer aggregation --- seedpod_ground_risk/core/plot_server.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index d7b20d4b..c9f059fc 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -213,9 +213,9 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) import matplotlib.pyplot as mpl plot = Overlay([res[0] for res in self._generated_data_layers.values()]) raw_datas = [res[2] for res in self._generated_data_layers.values()] - raster_indices = dict(Longitude=np.linspace(x_range[0], x_range[1], num=400), - Latitude=np.linspace(y_range[0], y_range[1], num=400)) - raster_grid = np.zeros((400, 400), dtype=np.float64) + raster_indices = dict(Longitude=np.linspace(x_range[0], x_range[1], num=raster_shape[0]), + Latitude=np.linspace(y_range[0], y_range[1], num=raster_shape[1])) + raster_grid = np.zeros((raster_shape[1], raster_shape[0]), dtype=np.float64) for res in self._generated_data_layers.values(): layer_raster_grid = res[1] nans = np.isnan(layer_raster_grid) @@ -260,7 +260,7 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) return plot.opts(width=self.plot_size[0], height=self.plot_size[1], tools=self.tools, active_tools=self.active_tools) - def generate_layers(self, bounds_poly: sg.Polygon) -> NoReturn: + def generate_layers(self, bounds_poly: sg.Polygon, raster_shape: Tuple[int, int]) -> NoReturn: """ Generate static layers of map From ff7cf788ff21fdbcd3540851c2b0822a1718e07c Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 5 Mar 2021 13:26:43 +0000 Subject: [PATCH 05/57] Set raster resolution as kwarg to plot_server --- seedpod_ground_risk/core/plot_server.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index c9f059fc..451a9212 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -42,6 +42,7 @@ class PlotServer: def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = None, active_tools: Optional[Iterable[str]] = None, cmap: str = 'CET_L18', + raster_resolution: float = 30, plot_size: Tuple[int, int] = (760, 735), progress_callback: Optional[Callable[[str], None]] = None, update_callback: Optional[Callable[[str], None]] = None): @@ -51,8 +52,8 @@ def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = No :param str tiles: a geoviews.tile_sources attribute string from http://geoviews.org/gallery/bokeh/tile_sources.html#bokeh-gallery-tile-sources :param List[str] tools: the bokeh tools to make available for the plot from https://docs.bokeh.org/en/latest/docs/user_guide/tools.html :param List[str] active_tools: the subset of `tools` that should be enabled by default - :param bool rasterise: Whether to opportunistically raster layers :param cmap: a colorcet attribute string for the colourmap to use from https://colorcet.holoviz.org/user_guide/Continuous.html + :param raster_resolution: resolution of a single square of the raster pixel grid in metres :param Tuple[int, int] plot_size: the plot size in (width, height) order :param progress_callback: an optional callable that takes a string updating progress :param update_callback: an optional callable that is called before an plot is rendered @@ -83,7 +84,7 @@ def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = No self._x_range, self._y_range = [-1.45, -1.35], [50.85, 50.95] - self.raster_resolution_m = 20 + self.raster_resolution_m = raster_resolution self._epsg4326_to_epsg3857_proj = None self._epsg3857_to_epsg4326_proj = None From 1dc2bc8187bed3dd622622c0fcab2108bfb38523 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 5 Mar 2021 18:13:08 +0000 Subject: [PATCH 06/57] Feature arbitrary obstacles (#13) Add arbitrary obstacle layer --- .../layers/arbitrary_obstacle_layer.py | 44 +++++++++++++++++++ .../layers/blockable_data_layer.py | 27 ++++++++++++ seedpod_ground_risk/layers/osm_tag_layer.py | 22 +++------- .../ui_resources/layer_options.py | 12 ++++- 4 files changed, 87 insertions(+), 18 deletions(-) create mode 100644 seedpod_ground_risk/layers/arbitrary_obstacle_layer.py create mode 100644 seedpod_ground_risk/layers/blockable_data_layer.py diff --git a/seedpod_ground_risk/layers/arbitrary_obstacle_layer.py b/seedpod_ground_risk/layers/arbitrary_obstacle_layer.py new file mode 100644 index 00000000..9adb382f --- /dev/null +++ b/seedpod_ground_risk/layers/arbitrary_obstacle_layer.py @@ -0,0 +1,44 @@ +from typing import NoReturn, Tuple + +import geopandas as gpd +import numpy as np +from holoviews.element import Geometry +from shapely import geometry as sg + +from seedpod_ground_risk.layers.blockable_data_layer import BlockableDataLayer + + +class ArbitraryObstacleLayer(BlockableDataLayer): + def __init__(self, key, filepath: str = '', **kwargs): + super(ArbitraryObstacleLayer, self).__init__(key, **kwargs) + self.filepath = filepath + + def preload_data(self) -> NoReturn: + + self.dataframe = gpd.read_file(self.filepath) + if self.buffer_dist: + epsg3857_geom = self.dataframe.to_crs('EPSG:3857').geometry + self.buffer_poly = gpd.GeoDataFrame( + {'geometry': epsg3857_geom.buffer(self.buffer_dist).to_crs('EPSG:4326')} + ) + + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> \ + Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: + import geoviews as gv + from holoviews.operation.datashader import rasterize + + bounds = bounds_polygon.bounds + + bounded_df = self.dataframe.cx[bounds[1]:bounds[3], bounds[0]:bounds[2]] + + polys = gv.Polygons(bounded_df).opts(style={'alpha': 0.8, 'color': self._colour}) + raster = rasterize(polys, width=raster_shape[0], height=raster_shape[1], + x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) + raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) + if self.blocking: + raster_grid[raster_grid != 0] = -1 + + return polys, raster_grid, gpd.GeoDataFrame(self.dataframe) + + def clear_cache(self): + pass diff --git a/seedpod_ground_risk/layers/blockable_data_layer.py b/seedpod_ground_risk/layers/blockable_data_layer.py new file mode 100644 index 00000000..b19c9291 --- /dev/null +++ b/seedpod_ground_risk/layers/blockable_data_layer.py @@ -0,0 +1,27 @@ +import random +from abc import ABCMeta, abstractmethod + +from seedpod_ground_risk.layers.data_layer import DataLayer + + +class BlockableDataLayer(DataLayer, metaclass=ABCMeta): + def __init__(self, key, colour: str = None, + blocking=False, buffer_dist=0): + super().__init__(key) + self.blocking = blocking + self.buffer_dist = buffer_dist + # Set a random colour + self._colour = colour if colour is not None else "#" + ''.join( + [random.choice('0123456789ABCDEF') for _ in range(6)]) + + @abstractmethod + def preload_data(self): + pass + + @abstractmethod + def generate(self, bounds_polygon, raster_shape, from_cache: bool = False, **kwargs): + pass + + @abstractmethod + def clear_cache(self): + pass diff --git a/seedpod_ground_risk/layers/osm_tag_layer.py b/seedpod_ground_risk/layers/osm_tag_layer.py index 7660f5f5..cf4a34ce 100644 --- a/seedpod_ground_risk/layers/osm_tag_layer.py +++ b/seedpod_ground_risk/layers/osm_tag_layer.py @@ -1,4 +1,3 @@ -import random from typing import NoReturn, Tuple import geopandas as gpd @@ -6,28 +5,19 @@ import shapely.geometry as sg from holoviews.element import Geometry -from seedpod_ground_risk.layers.data_layer import DataLayer +from seedpod_ground_risk.layers.blockable_data_layer import BlockableDataLayer -class OSMTagLayer(DataLayer): +class OSMTagLayer(BlockableDataLayer): _landuse_polygons: gpd.GeoDataFrame - def __init__(self, key, osm_tag, colour: str = None, - blocking=False, buffer_dist=0): - super().__init__(key) + def __init__(self, key, osm_tag, **kwargs): + super(OSMTagLayer, self).__init__(key, **kwargs) self._osm_tag = osm_tag - self.blocking = blocking - self.buffer_dist = buffer_dist - # Set a random colour - self._colour = colour if colour is not None else "#" + ''.join( - [random.choice('0123456789ABCDEF') for _ in range(6)]) self._landuse_polygons = gpd.GeoDataFrame() - def preload_data(self) -> NoReturn: - pass - def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> \ - Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: + Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: import geoviews as gv from holoviews.operation.datashader import rasterize @@ -53,7 +43,7 @@ def clear_cache(self): def query_osm_polygons(self, bound_poly: sg.Polygon) -> NoReturn: """ - Perform blocking query on OpenStreetMaps Overpass API for objects with the passed landuse. + Perform blocking query on OpenStreetMaps Overpass API for objects with the passed tag. Retain only polygons and store in GeoPandas GeoDataFrame :param shapely.Polygon bound_poly: bounding box around requested area in EPSG:4326 coordinates """ diff --git a/seedpod_ground_risk/ui_resources/layer_options.py b/seedpod_ground_risk/ui_resources/layer_options.py index fecc6bc4..3a6b8c2a 100644 --- a/seedpod_ground_risk/ui_resources/layer_options.py +++ b/seedpod_ground_risk/ui_resources/layer_options.py @@ -1,3 +1,4 @@ +from seedpod_ground_risk.layers.arbitrary_obstacle_layer import ArbitraryObstacleLayer from seedpod_ground_risk.layers.geojson_layer import GeoJSONLayer from seedpod_ground_risk.layers.osm_tag_layer import OSMTagLayer from seedpod_ground_risk.layers.pathfinding_layer import PathfindingLayer @@ -12,7 +13,8 @@ 'Generic OSM': OSMTagLayer, 'Residential - England': ResidentialLayer, 'Road Traffic Flow - England': RoadsLayer, - 'GeoJSON Geometries': GeoJSONLayer, + 'Existing Path Analysis': GeoJSONLayer, + 'External Obstacles': ArbitraryObstacleLayer, 'Pathfinding': PathfindingLayer } @@ -38,10 +40,16 @@ }, 'Road Traffic Flow - England': { }, - 'GeoJSON Geometries': { + 'Existing Path Analysis': { 'File': ('path', 'filepath', str), 'Buffer Distance [m]': (r'\d{0,3}', 'buffer_dist', int), }, + 'External Obstacles': { + 'File': ('path', 'filepath', str), + 'Colour': ('colour', 'colour', str), + 'Buffer Distance [m]': (r'\d{0,3}', 'buffer_dist', int), + 'Blocking': (bool, 'blocking', bool), + }, 'Pathfinding': { 'Buffer Distance [m]': (r'\d{0,3}', 'buffer', int), 'Start Latitude [dd]': (r'-?\d{0,3}\.\d+', 'start_lat', float), From cb19eba2f543b4119680e96b06a1fd921b10aa53 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sat, 6 Mar 2021 14:58:10 +0000 Subject: [PATCH 07/57] Include road names with roads raw data --- make_installer.iss | 4 ++-- seedpod_ground_risk/layers/roads_layer.py | 9 ++++++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/make_installer.iss b/make_installer.iss index 10f9fbc6..7890fb2c 100644 --- a/make_installer.iss +++ b/make_installer.iss @@ -22,7 +22,7 @@ DefaultDirName={autopf}\{#MyAppName} DefaultGroupName=SEEDPOD AllowNoIcons=yes ;LicenseFile=D:\PycharmProjects\seedpod_ground_risk\LICENSE -InfoBeforeFile=D:\PycharmProjects\seedpod_ground_risk\README.md +InfoBeforeFile=README.md ; Remove the following line to run in administrative install mode (install for all users.) PrivilegesRequired=lowest PrivilegesRequiredOverridesAllowed=commandline @@ -38,7 +38,7 @@ Name: "english"; MessagesFile: "compiler:Default.isl" Name: "desktopicon"; Description: "{cm:CreateDesktopIcon}"; GroupDescription: "{cm:AdditionalIcons}"; Flags: unchecked [Files] -Source: "D:\PycharmProjects\seedpod_ground_risk\dist\{#MyAppName}\*"; DestDir: "{app}"; Flags: ignoreversion recursesubdirs createallsubdirs +Source: "dist\{#MyAppName}\*"; DestDir: "{app}"; Flags: ignoreversion recursesubdirs createallsubdirs ; NOTE: Don't use "Flags: ignoreversion" on any shared system files [Icons] diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index c4be86e1..8ec84f55 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -74,10 +74,13 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr relative_variation = self.relative_variations_flat[hour] - points = np.array(self._interpolate_traffic_counts(bounds_polygon)) + road_points = self._interpolate_traffic_counts(bounds_polygon) + points = np.array(road_points)[:, :3].astype(np.float) + names = np.array(road_points)[:, 3] tfc_df = gpd.GeoDataFrame( {'geometry': [sg.Point(lon, lat) for lon, lat in zip(points[:, 0], points[:, 1])], - 'population': points[:, 2] * relative_variation}) + 'population': points[:, 2] * relative_variation, + 'road_name': names}) points = gv.Points(tfc_df, kdims=['Longitude', 'Latitude'], vdims=['population']).opts(colorbar=True, cmap=colorcet.CET_L18, @@ -206,7 +209,7 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = # Recover the true road geometry along with the interpolated values for idx, mark in enumerate(coord_spacing): c = road_ls.interpolate(mark) - all_interp_points_append([*self.proj.transform(c.y, c.x), interp_flat_proj[idx]]) + all_interp_points_append([*self.proj.transform(c.y, c.x), interp_flat_proj[idx], road_name]) print(".") return all_interp_points From 9b7cf984ff543f3978f0f2a0fcab39e246eb7fe7 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sat, 6 Mar 2021 14:59:07 +0000 Subject: [PATCH 08/57] Calculate and label each road population flow individually --- seedpod_ground_risk/layers/geojson_layer.py | 26 +++++++++++++++------ 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index 0113dee6..817310f2 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -33,7 +33,7 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np import geoviews as gv if self.buffer_poly is not None: - label_str = '' + labels = [] annotation_layers = [] for gdf in data: @@ -48,23 +48,35 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np continue proj_gdf = overlay.to_crs('epsg:27700') proj_gdf['population'] = proj_gdf.geometry.area * proj_gdf.density - label_str += f"Static Population: {proj_gdf['population'].sum():.2f}" + labels.append(gv.Text(self.endpoint[0], self.endpoint[1], + f"Static Population: {proj_gdf['population'].sum():.2f}", fontsize=20).opts( + style={'color': 'blue'})) annotation_layers.append(gv.Polygons(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'})) elif geom_type == 'Point': if 'population' in overlay: - mean_pop = overlay.population.mean() - label_str += f'Dynamic Population: {mean_pop / 60:.2f}/min' + # Group data by road name, localising the points + road_geoms = list(overlay.groupby('road_name').geometry) + road_coms = {} # store the centre of mass of points associated with a road + for name, geoms in road_geoms: + mean_lon = sum([p.x for p in geoms]) / len(geoms) + mean_lat = sum([p.y for p in geoms]) / len(geoms) + road_coms[name] = (mean_lon, mean_lat) # x,y order + # get mean population flow of points for road + road_pops = list(overlay.groupby('road_name').mean().itertuples()) + for name, pop in road_pops: # add labels at the centre of mass of each group of points + labels.append(gv.Text(road_coms[name][0], road_coms[name][1], + f'Population flow on road {name}: {pop / 60:.2f}/min', + fontsize=20).opts( + style={'color': 'blue'})) annotation_layers.append((gv.Points(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'}))) elif geom_type == 'Line': pass - label_str += '\n' - return Overlay([ gv.Contours(self.dataframe).opts(line_width=4, line_color='magenta'), gv.Polygons(self.buffer_poly).opts(style={'alpha': 0.3, 'color': 'cyan'}), # Add the path stats as a text annotation to the final path point - gv.Text(self.endpoint[0], self.endpoint[1], label_str, fontsize=20).opts(style={'color': 'blue'}), + *labels, *annotation_layers ]) else: From 20bdf4fdc476d27e190e320629b1b91ee2d14eac Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sat, 6 Mar 2021 15:12:53 +0000 Subject: [PATCH 09/57] Splitting out dynamic population footprint per road (#14) * Include road names with roads raw data * Calculate and label each road population flow individually --- make_installer.iss | 4 ++-- seedpod_ground_risk/layers/geojson_layer.py | 26 +++++++++++++++------ seedpod_ground_risk/layers/roads_layer.py | 9 ++++--- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/make_installer.iss b/make_installer.iss index 10f9fbc6..7890fb2c 100644 --- a/make_installer.iss +++ b/make_installer.iss @@ -22,7 +22,7 @@ DefaultDirName={autopf}\{#MyAppName} DefaultGroupName=SEEDPOD AllowNoIcons=yes ;LicenseFile=D:\PycharmProjects\seedpod_ground_risk\LICENSE -InfoBeforeFile=D:\PycharmProjects\seedpod_ground_risk\README.md +InfoBeforeFile=README.md ; Remove the following line to run in administrative install mode (install for all users.) PrivilegesRequired=lowest PrivilegesRequiredOverridesAllowed=commandline @@ -38,7 +38,7 @@ Name: "english"; MessagesFile: "compiler:Default.isl" Name: "desktopicon"; Description: "{cm:CreateDesktopIcon}"; GroupDescription: "{cm:AdditionalIcons}"; Flags: unchecked [Files] -Source: "D:\PycharmProjects\seedpod_ground_risk\dist\{#MyAppName}\*"; DestDir: "{app}"; Flags: ignoreversion recursesubdirs createallsubdirs +Source: "dist\{#MyAppName}\*"; DestDir: "{app}"; Flags: ignoreversion recursesubdirs createallsubdirs ; NOTE: Don't use "Flags: ignoreversion" on any shared system files [Icons] diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index 0113dee6..817310f2 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -33,7 +33,7 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np import geoviews as gv if self.buffer_poly is not None: - label_str = '' + labels = [] annotation_layers = [] for gdf in data: @@ -48,23 +48,35 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np continue proj_gdf = overlay.to_crs('epsg:27700') proj_gdf['population'] = proj_gdf.geometry.area * proj_gdf.density - label_str += f"Static Population: {proj_gdf['population'].sum():.2f}" + labels.append(gv.Text(self.endpoint[0], self.endpoint[1], + f"Static Population: {proj_gdf['population'].sum():.2f}", fontsize=20).opts( + style={'color': 'blue'})) annotation_layers.append(gv.Polygons(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'})) elif geom_type == 'Point': if 'population' in overlay: - mean_pop = overlay.population.mean() - label_str += f'Dynamic Population: {mean_pop / 60:.2f}/min' + # Group data by road name, localising the points + road_geoms = list(overlay.groupby('road_name').geometry) + road_coms = {} # store the centre of mass of points associated with a road + for name, geoms in road_geoms: + mean_lon = sum([p.x for p in geoms]) / len(geoms) + mean_lat = sum([p.y for p in geoms]) / len(geoms) + road_coms[name] = (mean_lon, mean_lat) # x,y order + # get mean population flow of points for road + road_pops = list(overlay.groupby('road_name').mean().itertuples()) + for name, pop in road_pops: # add labels at the centre of mass of each group of points + labels.append(gv.Text(road_coms[name][0], road_coms[name][1], + f'Population flow on road {name}: {pop / 60:.2f}/min', + fontsize=20).opts( + style={'color': 'blue'})) annotation_layers.append((gv.Points(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'}))) elif geom_type == 'Line': pass - label_str += '\n' - return Overlay([ gv.Contours(self.dataframe).opts(line_width=4, line_color='magenta'), gv.Polygons(self.buffer_poly).opts(style={'alpha': 0.3, 'color': 'cyan'}), # Add the path stats as a text annotation to the final path point - gv.Text(self.endpoint[0], self.endpoint[1], label_str, fontsize=20).opts(style={'color': 'blue'}), + *labels, *annotation_layers ]) else: diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index c4be86e1..8ec84f55 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -74,10 +74,13 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr relative_variation = self.relative_variations_flat[hour] - points = np.array(self._interpolate_traffic_counts(bounds_polygon)) + road_points = self._interpolate_traffic_counts(bounds_polygon) + points = np.array(road_points)[:, :3].astype(np.float) + names = np.array(road_points)[:, 3] tfc_df = gpd.GeoDataFrame( {'geometry': [sg.Point(lon, lat) for lon, lat in zip(points[:, 0], points[:, 1])], - 'population': points[:, 2] * relative_variation}) + 'population': points[:, 2] * relative_variation, + 'road_name': names}) points = gv.Points(tfc_df, kdims=['Longitude', 'Latitude'], vdims=['population']).opts(colorbar=True, cmap=colorcet.CET_L18, @@ -206,7 +209,7 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = # Recover the true road geometry along with the interpolated values for idx, mark in enumerate(coord_spacing): c = road_ls.interpolate(mark) - all_interp_points_append([*self.proj.transform(c.y, c.x), interp_flat_proj[idx]]) + all_interp_points_append([*self.proj.transform(c.y, c.x), interp_flat_proj[idx], road_name]) print(".") return all_interp_points From 65925e42a41502b23271ae6518f4cc5ecd2cb2af Mon Sep 17 00:00:00 2001 From: Zach10a <35232913+Zach10a@users.noreply.github.com> Date: Sat, 6 Mar 2021 19:59:52 +0000 Subject: [PATCH 10/57] Feature export path (#12) Exporting of path geometries --- seedpod_ground_risk/core/plot_server.py | 5 +++++ seedpod_ground_risk/core/plot_worker.py | 4 ++++ seedpod_ground_risk/main.py | 11 +++++++++++ 3 files changed, 20 insertions(+) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 451a9212..333a814c 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -344,6 +344,11 @@ def remove_layer(self, layer): def set_layer_order(self, layer_order): self.data_layer_order = layer_order + def export_path_geojson(self, layer, filepath): + import os + if layer in self.annotation_layers: + layer.dataframe.to_file(os.path.join(os.sep, f'{filepath}', 'path.geojson'), driver='GeoJSON') + def _get_raster_dimensions(self, bounds_poly: sg.Polygon, raster_resolution_m: float) -> Tuple[int, int]: """ Return a the (x,y) shape of a raster grid given its EPSG4326 envelope and desired raster resolution diff --git a/seedpod_ground_risk/core/plot_worker.py b/seedpod_ground_risk/core/plot_worker.py index c78442e0..05a9a216 100644 --- a/seedpod_ground_risk/core/plot_worker.py +++ b/seedpod_ground_risk/core/plot_worker.py @@ -75,6 +75,10 @@ def set_time(self, hour): def layers_reorder(self, layer_order): self.plot_server.set_layer_order(layer_order) + @Slot(Layer, str) + def export_path_json(self, layer, filepath): + self.plot_server.export_path_geojson(layer, filepath) + def layers_update(self, layers): self.signals.update_layers.emit(layers) diff --git a/seedpod_ground_risk/main.py b/seedpod_ground_risk/main.py index 7aa2c128..6af111d1 100644 --- a/seedpod_ground_risk/main.py +++ b/seedpod_ground_risk/main.py @@ -6,6 +6,7 @@ import PySide2 from seedpod_ground_risk.core.plot_worker import PlotWorker +from seedpod_ground_risk.layers.annotation_layer import AnnotationLayer from seedpod_ground_risk.ui_resources.add_layer_wizard import LayerWizard from seedpod_ground_risk.ui_resources.layer_options import LAYER_OBJECTS from seedpod_ground_risk.ui_resources.layerlistdelegate import Ui_delegate @@ -73,11 +74,21 @@ def change_colour(self): def delete_layer(self): self._plot_worker.remove_layer(self._layer) + def export_path_json(self): + from PySide2.QtWidgets import QFileDialog + + file_dir = QFileDialog.getExistingDirectory(self, "Save plot .json file...", os.getcwd(), + QFileDialog.ShowDirsOnly) + if file_dir: + self._plot_worker.export_path_json(self._layer, file_dir) + def mousePressEvent(self, event: PySide2.QtGui.QMouseEvent) -> None: super().mousePressEvent(event) if event.button() == Qt.RightButton: menu = QMenu() menu.addAction("Delete", self.delete_layer) + if isinstance(self._layer, AnnotationLayer): + menu.addAction("Export .GeoJSON", self.export_path_json) menu.exec_(event.globalPos()) From aef8f5a1de3f66e5ec33289f0226d1e528e819d2 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sat, 6 Mar 2021 20:41:40 +0000 Subject: [PATCH 11/57] Display overflown residential populations on hover --- seedpod_ground_risk/core/plot_server.py | 2 +- seedpod_ground_risk/layers/geojson_layer.py | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 451a9212..e8ea853d 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -58,7 +58,7 @@ def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = No :param progress_callback: an optional callable that takes a string updating progress :param update_callback: an optional callable that is called before an plot is rendered """ - self.tools = ['hover', 'crosshair'] if tools is None else tools + self.tools = ['crosshair'] if tools is None else tools self.active_tools = ['wheel_zoom'] if active_tools is None else active_tools import colorcet diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index 817310f2..f535c0a1 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -51,7 +51,9 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np labels.append(gv.Text(self.endpoint[0], self.endpoint[1], f"Static Population: {proj_gdf['population'].sum():.2f}", fontsize=20).opts( style={'color': 'blue'})) - annotation_layers.append(gv.Polygons(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'})) + proj_gdf = proj_gdf.to_crs('EPSG:4326') + annotation_layers.append( + gv.Polygons(proj_gdf, vdims=['population']).opts(alpha=0.8, color='cyan', tools=['hover'])) elif geom_type == 'Point': if 'population' in overlay: # Group data by road name, localising the points From e6b27ad15e5fba754e0f746fe01cfca36bd57b73 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 8 Mar 2021 09:34:09 +0000 Subject: [PATCH 12/57] Extract snap to grid func --- seedpod_ground_risk/layers/geojson_layer.py | 13 +++++++++++++ seedpod_ground_risk/layers/pathfinding_layer.py | 8 ++++---- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index f535c0a1..baa7c374 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -86,3 +86,16 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np def clear_cache(self) -> NoReturn: pass + + def _snap_coords_to_grid(self, grid, lon: float, lat: float) -> Tuple[int, int]: + """ + Snap coordinates to grid indices + :param grid: raster grid coordinates + :param lon: longitude to snap + :param lat: latitude to snap + :return: (x, y) tuple of grid indices + """ + lat_idx = int(np.argmin(np.abs(grid['Latitude'] - lat))) + lon_idx = int(np.argmin(np.abs(grid['Longitude'] - lon))) + + return lon_idx, lat_idx diff --git a/seedpod_ground_risk/layers/pathfinding_layer.py b/seedpod_ground_risk/layers/pathfinding_layer.py index 1832986f..0e3985ab 100644 --- a/seedpod_ground_risk/layers/pathfinding_layer.py +++ b/seedpod_ground_risk/layers/pathfinding_layer.py @@ -32,13 +32,13 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np **kwargs) -> Geometry: from seedpod_ground_risk.pathfinding import environment - snapped_start_lat_idx = np.argmin(np.abs(raster_data[0]['Latitude'] - self.start_coords[0])) - snapped_start_lon_idx = np.argmin(np.abs(raster_data[0]['Longitude'] - self.start_coords[1])) + snapped_start_lon_idx, snapped_start_lat_idx = self._snap_coords_to_grid(raster_data[0], self.start_coords[1], + self.start_coords[0]) start_node = environment.Node(snapped_start_lon_idx, snapped_start_lat_idx, raster_data[1][snapped_start_lat_idx, snapped_start_lon_idx]) - snapped_end_lat_idx = np.argmin(np.abs(raster_data[0]['Latitude'] - self.end_coords[0])) - snapped_end_lon_idx = np.argmin(np.abs(raster_data[0]['Longitude'] - self.end_coords[1])) + snapped_end_lon_idx, snapped_end_lat_idx = self._snap_coords_to_grid(raster_data[0], self.end_coords[1], + self.end_coords[0]) end_node = environment.Node(snapped_end_lon_idx, snapped_end_lat_idx, raster_data[1][snapped_end_lat_idx, snapped_end_lon_idx]) From dfbafb47ac142900871897a5f1f85818a3029536 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 8 Mar 2021 13:22:49 +0000 Subject: [PATCH 13/57] Fix bresenham at high gradients (|m|>1) --- seedpod_ground_risk/pathfinding/bresenham.py | 8 ++++---- tests/pathfinding/test_bresenham.py | 16 ++++++++++++++++ 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/seedpod_ground_risk/pathfinding/bresenham.py b/seedpod_ground_risk/pathfinding/bresenham.py index 82c48a96..2a86a626 100644 --- a/seedpod_ground_risk/pathfinding/bresenham.py +++ b/seedpod_ground_risk/pathfinding/bresenham.py @@ -36,10 +36,10 @@ def _make_line_high(x0: int, y0: int, x1: int, y1: int): dy = y1 - y0 mx = max(x0, x1) xi = 1 - if dy < 0: + if dx < 0: xi = -1 dx = -dx - D = 2 * dy - dx + D = 2 * dx - dy x = x0 n = dy + 1 line = np.zeros((n, 2), dtype=np.int32) @@ -49,9 +49,9 @@ def _make_line_high(x0: int, y0: int, x1: int, y1: int): x += xi if x > mx: x = mx - D += 2 * (dy - dx) + D += 2 * (dx - dy) else: - D += 2 * dy + D += 2 * dx line[-1] = [y1, x1] return line diff --git a/tests/pathfinding/test_bresenham.py b/tests/pathfinding/test_bresenham.py index 9a0f4854..82157c11 100644 --- a/tests/pathfinding/test_bresenham.py +++ b/tests/pathfinding/test_bresenham.py @@ -80,6 +80,22 @@ def test_bounds(self): result = (flat_points <= 4).all() and (flat_points >= 0).all() self.assertTrue(result, "Points not bounded within endpoints") + def test_high_gradient(self): + start = [0, 6] # x, y + end = [6, 0] # x, y + + path = make_line(*start, *end) + + expected = np.array([[0, 6], + [1, 5], + [2, 4], + [3, 3], + [4, 2], + [5, 1], + [6, 0]]) + result = (path == expected).all() or (path == np.flipud(expected)).all() + self.assertTrue(result, 'First Quadrant path incorrect') + if __name__ == '__main__': unittest.main() From 5fa8bde1d0b1e95a8b0134b3269cd62b21699013 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 8 Mar 2021 15:14:12 +0000 Subject: [PATCH 14/57] Calculate path PDFs --- seedpod_ground_risk/layers/geojson_layer.py | 29 ++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index baa7c374..6f93e7e7 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -5,6 +5,7 @@ from holoviews.element import Overlay from seedpod_ground_risk.layers.annotation_layer import AnnotationLayer +from seedpod_ground_risk.pathfinding import bresenham class GeoJSONLayer(AnnotationLayer): @@ -31,10 +32,36 @@ def preload_data(self) -> NoReturn: def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np.array], np.array], **kwargs) -> Overlay: import geoviews as gv + from scipy.stats import multivariate_normal as mvn + + bounds = (raster_data[0]['Longitude'].min(), raster_data[0]['Latitude'].min(), + raster_data[0]['Longitude'].max(), raster_data[0]['Latitude'].max()) if self.buffer_poly is not None: - labels = [] + # Snap the line string nodes to the raster grid + snapped_points = [self._snap_coords_to_grid(raster_data[0], *coords) for coords in + self.dataframe.iloc[0].geometry.coords] + # Generate pairs of consecutive (x,y) coords + path_pairs = list(map(list, zip(snapped_points, snapped_points[1:]))) + # Feed these pairs into the Bresenham algo to find the intermediate points + path_grid_points = [bresenham.make_line(*pair[0], *pair[1]) for pair in path_pairs] + # Bring all these points together and remove duplicate coords + # Flip left to right as bresenham spits out in (y,x) order + path_grid_points = np.fliplr(np.unique(np.concatenate(path_grid_points, axis=0), axis=0)) + + dist_mean = np.array([5, 5]) + raster_shape = raster_data[1].shape + x, y = np.mgrid[0:raster_shape[0], 0:raster_shape[1]] + eval_grid = np.dstack((x, y)) + + def point_dist(c): + pdf_mean = c + dist_mean + return mvn(pdf_mean, [[10, 0], [0, 10]]).pdf(eval_grid) + + # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? + pdf_mat = np.rot90(np.sum([point_dist(c) for c in path_grid_points], axis=0)) + labels = [] annotation_layers = [] for gdf in data: if not gdf.crs: From 699fa253e801a340a3188817a06cabb160ff6b40 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 8 Mar 2021 15:15:04 +0000 Subject: [PATCH 15/57] Remove buffer polygon display, add raster image of PDF as overlay --- seedpod_ground_risk/layers/geojson_layer.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index 6f93e7e7..e6f601fb 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -103,10 +103,13 @@ def point_dist(c): return Overlay([ gv.Contours(self.dataframe).opts(line_width=4, line_color='magenta'), - gv.Polygons(self.buffer_poly).opts(style={'alpha': 0.3, 'color': 'cyan'}), + # gv.Polygons(self.buffer_poly).opts(style={'alpha': 0.3, 'color': 'cyan'}), # Add the path stats as a text annotation to the final path point *labels, - *annotation_layers + *annotation_layers, + gv.Image(pdf_mat, bounds=bounds).options(alpha=0.4, clipping_colors={'0': 'transparent', + 'NaN': 'transparent', + '-NaN': 'transparent'}) ]) else: return gv.Contours(self.dataframe).opts(line_width=4, line_color='#000000') From 562b2aab804db1859b1696c88343f63b889beb28 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Tue, 9 Mar 2021 11:59:56 +0000 Subject: [PATCH 16/57] Chunk and parallelize path pdf accumulation --- seedpod_ground_risk/layers/geojson_layer.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/geojson_layer.py index e6f601fb..2d062771 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/geojson_layer.py @@ -49,17 +49,29 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np # Flip left to right as bresenham spits out in (y,x) order path_grid_points = np.fliplr(np.unique(np.concatenate(path_grid_points, axis=0), axis=0)) - dist_mean = np.array([5, 5]) + # Testing fixed pdf for now + dist_mean = np.array([15, 5]) + + # Create grid on which to evaluate each point of path with its pdf raster_shape = raster_data[1].shape x, y = np.mgrid[0:raster_shape[0], 0:raster_shape[1]] eval_grid = np.dstack((x, y)) + from pathos.multiprocessing import Pool + from pathos.multiprocessing import cpu_count + pool = Pool() + + # TODO Use something like Dask or OpenCV to speed this up in future as it's a simple map-reduce def point_dist(c): pdf_mean = c + dist_mean - return mvn(pdf_mean, [[10, 0], [0, 10]]).pdf(eval_grid) + return mvn(pdf_mean, [[15, 5], [5, 7]]).pdf(eval_grid) + + def map_chunk(coords): + return np.sum([point_dist(c) for c in coords], axis=0) # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? - pdf_mat = np.rot90(np.sum([point_dist(c) for c in path_grid_points], axis=0)) + pdfs = pool.map(map_chunk, np.array_split(path_grid_points, cpu_count() * 2)) + pdf_mat = np.rot90(np.sum(pdfs, axis=0)) labels = [] annotation_layers = [] From 5afc666131f4636259bd83700812d9bcf3358c9d Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Tue, 9 Mar 2021 15:01:35 +0000 Subject: [PATCH 17/57] Fix OSM layer not overriding all methods --- seedpod_ground_risk/layers/osm_tag_layer.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/seedpod_ground_risk/layers/osm_tag_layer.py b/seedpod_ground_risk/layers/osm_tag_layer.py index cf4a34ce..32a21e6d 100644 --- a/seedpod_ground_risk/layers/osm_tag_layer.py +++ b/seedpod_ground_risk/layers/osm_tag_layer.py @@ -16,6 +16,9 @@ def __init__(self, key, osm_tag, **kwargs): self._osm_tag = osm_tag self._landuse_polygons = gpd.GeoDataFrame() + def preload_data(self): + pass + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> \ Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: import geoviews as gv From 40a6efae7e1918c44b6b686a8eb6d040d048087b Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:00:48 +0000 Subject: [PATCH 18/57] Pass raster resolution to layers --- seedpod_ground_risk/core/plot_server.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index e2048cfc..34f9d8c3 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -225,7 +225,8 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) annotations = [] for layer in self.annotation_layers: - annotation = layer.annotate(raw_datas, (raster_indices, raster_grid)) + annotation = layer.annotate(raw_datas, (raster_indices, raster_grid), + resolution=self.raster_resolution_m) if annotation: annotations.append(annotation) @@ -275,8 +276,8 @@ def generate_layers(self, bounds_poly: sg.Polygon, raster_shape: Tuple[int, int] layers = {} self._progress_callback('Generating layer data') with ThreadPoolExecutor() as pool: - layer_futures = [pool.submit(self.generate_layer, layer, bounds_poly, raster_shape, self._time_idx) for - layer in self.data_layers] + layer_futures = [pool.submit(self.generate_layer, layer, bounds_poly, raster_shape, self._time_idx, + self.raster_resolution_m) for layer in self.data_layers] # Store generated layers as they are completed for future in as_completed(layer_futures): key, result = future.result() @@ -298,7 +299,8 @@ def generate_layers(self, bounds_poly: sg.Polygon, raster_shape: Tuple[int, int] {k: layers[k] for k in layers.keys() if k not in self._generated_data_layers}) @staticmethod - def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, raster_shape: Tuple[int, int], hour: int) -> Union[ + def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, raster_shape: Tuple[int, int], hour: int, + resolution: float) -> Union[ Tuple[str, Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]], Tuple[str, None]]: import shapely.ops as so @@ -314,7 +316,8 @@ def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, raster_shape: Tupl layer_bounds_poly = bounds_poly.difference(layer.cached_area) layer.cached_area = so.unary_union([layer.cached_area, bounds_poly]) try: - result = layer.key, layer.generate(layer_bounds_poly, raster_shape, from_cache=from_cache, hour=hour) + result = layer.key, layer.generate(layer_bounds_poly, raster_shape, from_cache=from_cache, hour=hour, + resolution=resolution) return result except Exception as e: print(e) From efde14cbd2bb0d5c7e3ac5980499b195a99f9bea Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:01:17 +0000 Subject: [PATCH 19/57] Fix raster grid flipping --- seedpod_ground_risk/core/plot_server.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 34f9d8c3..65996775 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -222,7 +222,7 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) nans = np.isnan(layer_raster_grid) layer_raster_grid[nans] = 0 raster_grid += res[1] - + raster_grid = np.flipud(raster_grid) annotations = [] for layer in self.annotation_layers: annotation = layer.annotate(raw_datas, (raster_indices, raster_grid), From 934c2126355b15b482300669a9f4baeb13ed9770 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:35:07 +0000 Subject: [PATCH 20/57] fixup! Pass raster resolution to layers --- seedpod_ground_risk/core/plot_server.py | 2 +- .../{geojson_layer.py => path_analysis_layer.py} | 2 +- seedpod_ground_risk/layers/pathfinding_layer.py | 14 ++++++++------ seedpod_ground_risk/layers/roads_layer.py | 2 +- 4 files changed, 11 insertions(+), 9 deletions(-) rename seedpod_ground_risk/layers/{geojson_layer.py => path_analysis_layer.py} (99%) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 65996775..8bdcfba4 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -42,7 +42,7 @@ class PlotServer: def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = None, active_tools: Optional[Iterable[str]] = None, cmap: str = 'CET_L18', - raster_resolution: float = 30, + raster_resolution: float = 40, plot_size: Tuple[int, int] = (760, 735), progress_callback: Optional[Callable[[str], None]] = None, update_callback: Optional[Callable[[str], None]] = None): diff --git a/seedpod_ground_risk/layers/geojson_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py similarity index 99% rename from seedpod_ground_risk/layers/geojson_layer.py rename to seedpod_ground_risk/layers/path_analysis_layer.py index 2d062771..791db233 100644 --- a/seedpod_ground_risk/layers/geojson_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -30,7 +30,7 @@ def preload_data(self) -> NoReturn: self.endpoint = self.dataframe.iloc[0].geometry.coords[-1] def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np.array], np.array], - **kwargs) -> Overlay: + resolution=20, **kwargs) -> Overlay: import geoviews as gv from scipy.stats import multivariate_normal as mvn diff --git a/seedpod_ground_risk/layers/pathfinding_layer.py b/seedpod_ground_risk/layers/pathfinding_layer.py index 0e3985ab..6e4cfeb0 100644 --- a/seedpod_ground_risk/layers/pathfinding_layer.py +++ b/seedpod_ground_risk/layers/pathfinding_layer.py @@ -29,27 +29,29 @@ def preload_data(self) -> NoReturn: pass def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np.array], np.array], - **kwargs) -> Geometry: + resolution=20, **kwargs) -> Geometry: from seedpod_ground_risk.pathfinding import environment + raster_grid = raster_data[1] * resolution ** 2 + snapped_start_lon_idx, snapped_start_lat_idx = self._snap_coords_to_grid(raster_data[0], self.start_coords[1], self.start_coords[0]) start_node = environment.Node(snapped_start_lon_idx, snapped_start_lat_idx, - raster_data[1][snapped_start_lat_idx, snapped_start_lon_idx]) + raster_grid[snapped_start_lat_idx, snapped_start_lon_idx]) snapped_end_lon_idx, snapped_end_lat_idx = self._snap_coords_to_grid(raster_data[0], self.end_coords[1], self.end_coords[0]) end_node = environment.Node(snapped_end_lon_idx, snapped_end_lat_idx, - raster_data[1][snapped_end_lat_idx, snapped_end_lon_idx]) + raster_grid[snapped_end_lat_idx, snapped_end_lon_idx]) - if raster_data[1][snapped_start_lon_idx, snapped_start_lat_idx] < 0: + if raster_grid[snapped_start_lon_idx, snapped_start_lat_idx] < 0: print('Start node in blocked area, path impossible') return None - elif raster_data[1][snapped_end_lon_idx, snapped_end_lat_idx] < 0: + elif raster_grid[snapped_end_lon_idx, snapped_end_lat_idx] < 0: print('End node in blocked area, path impossible') return None - env = environment.GridEnvironment(raster_data[1], diagonals=True, pruning=False) + env = environment.GridEnvironment(raster_grid, diagonals=True, pruning=False) algo = self.algo(heuristic=self.heuristic(env, risk_to_dist_ratio=self.rdr)) t0 = time() path = algo.find_path(env, start_node, end_node) diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index 8ec84f55..f9de9757 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -65,7 +65,7 @@ def preload_data(self) -> NoReturn: self._ingest_relative_traffic_variations() def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, - hour: int = 0, **kwargs) -> \ + hour: int = 0, resolution: float = 20, **kwargs) -> \ Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: from holoviews.operation.datashader import rasterize import geoviews as gv From 7d243cbdb0496b0d0c4467a994e760a2234fcbd5 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:37:20 +0000 Subject: [PATCH 21/57] Refactor Path analysis layer name --- seedpod_ground_risk/core/plot_server.py | 1 + seedpod_ground_risk/layers/path_analysis_layer.py | 4 ++-- seedpod_ground_risk/layers/pathfinding_layer.py | 4 ++-- seedpod_ground_risk/ui_resources/layer_options.py | 5 ++--- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 8bdcfba4..a95314fd 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -9,6 +9,7 @@ from seedpod_ground_risk.layers.annotation_layer import AnnotationLayer from seedpod_ground_risk.layers.data_layer import DataLayer from seedpod_ground_risk.layers.layer import Layer +from seedpod_ground_risk.layers.path_analysis_layer import PathAnalysisLayer def make_bounds_polygon(*args: Iterable[float]) -> sg.Polygon: diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 791db233..943125d1 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -8,10 +8,10 @@ from seedpod_ground_risk.pathfinding import bresenham -class GeoJSONLayer(AnnotationLayer): +class PathAnalysisLayer(AnnotationLayer): def __init__(self, key: str, filepath: str = '', buffer_dist: float = None, **kwargs): - super(GeoJSONLayer, self).__init__(key) + super(PathAnalysisLayer, self).__init__(key) self.filepath = filepath self.buffer_dist = buffer_dist diff --git a/seedpod_ground_risk/layers/pathfinding_layer.py b/seedpod_ground_risk/layers/pathfinding_layer.py index 6e4cfeb0..479d47fd 100644 --- a/seedpod_ground_risk/layers/pathfinding_layer.py +++ b/seedpod_ground_risk/layers/pathfinding_layer.py @@ -6,13 +6,13 @@ import shapely.geometry as sg from holoviews.element import Geometry -from seedpod_ground_risk.layers.geojson_layer import GeoJSONLayer +from seedpod_ground_risk.layers.path_analysis_layer import PathAnalysisLayer from seedpod_ground_risk.pathfinding.a_star import RiskGridAStar from seedpod_ground_risk.pathfinding.algorithm import Algorithm from seedpod_ground_risk.pathfinding.heuristic import ManhattanRiskHeuristic, Heuristic -class PathfindingLayer(GeoJSONLayer): +class PathfindingLayer(PathAnalysisLayer): def __init__(self, key, start_lat: float = 0, start_lon: float = 0, end_lat: float = 0, end_lon: float = 0, buffer: float = 0, algo: Algorithm = RiskGridAStar, heuristic: Heuristic = ManhattanRiskHeuristic, diff --git a/seedpod_ground_risk/ui_resources/layer_options.py b/seedpod_ground_risk/ui_resources/layer_options.py index 3a6b8c2a..2c581677 100644 --- a/seedpod_ground_risk/ui_resources/layer_options.py +++ b/seedpod_ground_risk/ui_resources/layer_options.py @@ -1,10 +1,9 @@ from seedpod_ground_risk.layers.arbitrary_obstacle_layer import ArbitraryObstacleLayer -from seedpod_ground_risk.layers.geojson_layer import GeoJSONLayer from seedpod_ground_risk.layers.osm_tag_layer import OSMTagLayer +from seedpod_ground_risk.layers.path_analysis_layer import PathAnalysisLayer from seedpod_ground_risk.layers.pathfinding_layer import PathfindingLayer from seedpod_ground_risk.layers.residential_layer import ResidentialLayer from seedpod_ground_risk.layers.roads_layer import RoadsLayer - from seedpod_ground_risk.pathfinding.a_star import * from seedpod_ground_risk.pathfinding.rjps_a_star import * @@ -13,7 +12,7 @@ 'Generic OSM': OSMTagLayer, 'Residential - England': ResidentialLayer, 'Road Traffic Flow - England': RoadsLayer, - 'Existing Path Analysis': GeoJSONLayer, + 'Existing Path Analysis': PathAnalysisLayer, 'External Obstacles': ArbitraryObstacleLayer, 'Pathfinding': PathfindingLayer } From 005b9e124e163789a5179a30054a139e41162c60 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:39:35 +0000 Subject: [PATCH 22/57] Use densities in raster grids --- seedpod_ground_risk/layers/residential_layer.py | 6 +++--- seedpod_ground_risk/layers/roads_layer.py | 13 ++++++++----- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/seedpod_ground_risk/layers/residential_layer.py b/seedpod_ground_risk/layers/residential_layer.py index d16f8288..56316dd4 100644 --- a/seedpod_ground_risk/layers/residential_layer.py +++ b/seedpod_ground_risk/layers/residential_layer.py @@ -57,12 +57,12 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr buffered_df = deepcopy(census_df) buffered_df.geometry = buffered_df.to_crs('EPSG:27700') \ .buffer(self.buffer_dist).to_crs('EPSG:4326') - buffered_polys = gv.Polygons(buffered_df, kdims=['Longitude', 'Latitude'], vdims=['name', 'population']) - raster = rasterize(buffered_polys, aggregator=ds.sum('population'), width=raster_shape[0], + buffered_polys = gv.Polygons(buffered_df, kdims=['Longitude', 'Latitude'], vdims=['name', 'density']) + raster = rasterize(buffered_polys, aggregator=ds.sum('density'), width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) else: - raster = rasterize(gv_polys, aggregator=ds.sum('population'), width=raster_shape[0], height=raster_shape[1], + raster = rasterize(gv_polys, aggregator=ds.sum('density'), width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index f9de9757..237a876b 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -79,14 +79,17 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr names = np.array(road_points)[:, 3] tfc_df = gpd.GeoDataFrame( {'geometry': [sg.Point(lon, lat) for lon, lat in zip(points[:, 0], points[:, 1])], - 'population': points[:, 2] * relative_variation, + 'population': points[:, 2] * relative_variation / 3600, 'road_name': names}) + # Assume motorways are 3 lanes each way and A roads are 2 each way + tfc_df['road_width'] = np.where(tfc_df['road_name'].str.startswith('M'), 22, 14.6) + tfc_df['density'] = tfc_df['population'] / (resolution * tfc_df['road_width']) points = gv.Points(tfc_df, - kdims=['Longitude', 'Latitude'], vdims=['population']).opts(colorbar=True, - cmap=colorcet.CET_L18, - color='population') + kdims=['Longitude', 'Latitude'], vdims=['population', 'density']).opts(colorbar=True, + cmap=colorcet.CET_L18, + color='population') bounds = bounds_polygon.bounds - raster = rasterize(points, aggregator=ds.mean('population'), width=raster_shape[0], height=raster_shape[1], + raster = rasterize(points, aggregator=ds.mean('density'), width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) return points, raster_grid, gpd.GeoDataFrame(tfc_df) From eca831b047ce690bf13998b29433cb998259b7a4 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:41:10 +0000 Subject: [PATCH 23/57] Use CasEx lib --- requirements.txt | 3 +++ seedpod_ground_risk/layers/path_analysis_layer.py | 6 ++++++ 2 files changed, 9 insertions(+) diff --git a/requirements.txt b/requirements.txt index 7164611f..79d83bf9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -91,3 +91,6 @@ typing-extensions==3.7.4.3 urllib3==1.26.2 xarray==0.16.2 zict==2.0.0 + +fastparquet~=0.5.0 +casex~=1.0.5 \ No newline at end of file diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 943125d1..b2f55f17 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -1,5 +1,6 @@ from typing import NoReturn, List, Tuple, Dict +import casex import geopandas as gpd import numpy as np from holoviews.element import Overlay @@ -49,6 +50,11 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np # Flip left to right as bresenham spits out in (y,x) order path_grid_points = np.fliplr(np.unique(np.concatenate(path_grid_points, axis=0), axis=0)) + ca_model = casex.CriticalAreaModels() # Lethal area model using default params for human dimensions + aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, 2, 2, 2) # Default aircraft + ballistic_model = casex.BallisticDescent2ndOrderDragApproximation() + ballistic_model.set_aircraft(aircraft) + # Testing fixed pdf for now dist_mean = np.array([15, 5]) From 6d70b636e45272626537ae9394f5fc1e1d82e05c Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:41:45 +0000 Subject: [PATCH 24/57] Flip pdf grid dimensions --- seedpod_ground_risk/layers/path_analysis_layer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index b2f55f17..6ce31740 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -60,7 +60,7 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np # Create grid on which to evaluate each point of path with its pdf raster_shape = raster_data[1].shape - x, y = np.mgrid[0:raster_shape[0], 0:raster_shape[1]] + x, y = np.mgrid[0:raster_shape[1], 0:raster_shape[0]] eval_grid = np.dstack((x, y)) from pathos.multiprocessing import Pool From 7a4729fd7e84b70112869bb92fa51098bafa010e Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:42:17 +0000 Subject: [PATCH 25/57] Fix multiprocess pool not being closed --- seedpod_ground_risk/layers/path_analysis_layer.py | 1 + 1 file changed, 1 insertion(+) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 6ce31740..e4cba98c 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -77,6 +77,7 @@ def map_chunk(coords): # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? pdfs = pool.map(map_chunk, np.array_split(path_grid_points, cpu_count() * 2)) + pool.close() pdf_mat = np.rot90(np.sum(pdfs, axis=0)) labels = [] From 4fed20bdbb69675ec044c2ff039e834718b8b732 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Wed, 10 Mar 2021 19:43:13 +0000 Subject: [PATCH 26/57] Fix extra columns in roads data not being unpacked --- seedpod_ground_risk/layers/path_analysis_layer.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index e4cba98c..9f7142ed 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -111,9 +111,9 @@ def map_chunk(coords): road_coms[name] = (mean_lon, mean_lat) # x,y order # get mean population flow of points for road road_pops = list(overlay.groupby('road_name').mean().itertuples()) - for name, pop in road_pops: # add labels at the centre of mass of each group of points + for name, pop, *_ in road_pops: # add labels at the centre of mass of each group of points labels.append(gv.Text(road_coms[name][0], road_coms[name][1], - f'Population flow on road {name}: {pop / 60:.2f}/min', + f'Population flow on road {name}: {pop * 60:.2f}/min', fontsize=20).opts( style={'color': 'blue'})) annotation_layers.append((gv.Points(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'}))) From 7b0643aa2adfda1a899da5a9e4df97e726dd8eee Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Thu, 11 Mar 2021 11:11:16 +0000 Subject: [PATCH 27/57] Fix overloaded methods in osm_tag_layer.py --- seedpod_ground_risk/layers/osm_tag_layer.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/seedpod_ground_risk/layers/osm_tag_layer.py b/seedpod_ground_risk/layers/osm_tag_layer.py index cf4a34ce..32a21e6d 100644 --- a/seedpod_ground_risk/layers/osm_tag_layer.py +++ b/seedpod_ground_risk/layers/osm_tag_layer.py @@ -16,6 +16,9 @@ def __init__(self, key, osm_tag, **kwargs): self._osm_tag = osm_tag self._landuse_polygons = gpd.GeoDataFrame() + def preload_data(self): + pass + def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], from_cache: bool = False, **kwargs) -> \ Tuple[Geometry, np.ndarray, gpd.GeoDataFrame]: import geoviews as gv From fd2da04459626e9627b737143a28cbef43dbdaee Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Thu, 11 Mar 2021 11:13:43 +0000 Subject: [PATCH 28/57] Generate polygons for roads instead of points --- seedpod_ground_risk/layers/roads_layer.py | 90 +++++++++++++++++------ 1 file changed, 66 insertions(+), 24 deletions(-) diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index 237a876b..c2b2cb58 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -1,4 +1,4 @@ -from typing import NoReturn, List, Tuple +from typing import NoReturn, Tuple import geopandas as gpd import numpy as np @@ -39,6 +39,37 @@ def generate_week_timesteps(): return timestep_index +# https://gis.stackexchange.com/questions/291247/interchange-y-x-to-x-y-with-geopandas-python-or-qgis +def swap_xy(geom): + if geom.is_empty: + return geom + + if geom.has_z: + def swap_xy_coords(coords): + for x, y, z in coords: + yield (y, x, z) + else: + def swap_xy_coords(coords): + for x, y in coords: + yield (y, x) + + # Process coordinates from each supported geometry type + if geom.type in ('Point', 'LineString', 'LinearRing'): + return type(geom)(list(swap_xy_coords(geom.coords))) + elif geom.type == 'Polygon': + ring = geom.exterior + shell = type(ring)(list(swap_xy_coords(ring.coords))) + holes = list(geom.interiors) + for pos, ring in enumerate(holes): + holes[pos] = type(ring)(list(swap_xy_coords(ring.coords))) + return type(geom)(shell, holes) + elif geom.type.startswith('Multi') or geom.type == 'GeometryCollection': + # Recursive call + return type(geom)([swap_xy(part) for part in geom.geoms]) + else: + raise ValueError('Type %r not recognized' % geom.type) + + class RoadsLayer(DataLayer): _traffic_counts: gpd.GeoDataFrame @@ -74,25 +105,23 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr relative_variation = self.relative_variations_flat[hour] - road_points = self._interpolate_traffic_counts(bounds_polygon) - points = np.array(road_points)[:, :3].astype(np.float) - names = np.array(road_points)[:, 3] - tfc_df = gpd.GeoDataFrame( - {'geometry': [sg.Point(lon, lat) for lon, lat in zip(points[:, 0], points[:, 1])], - 'population': points[:, 2] * relative_variation / 3600, - 'road_name': names}) - # Assume motorways are 3 lanes each way and A roads are 2 each way - tfc_df['road_width'] = np.where(tfc_df['road_name'].str.startswith('M'), 22, 14.6) - tfc_df['density'] = tfc_df['population'] / (resolution * tfc_df['road_width']) - points = gv.Points(tfc_df, - kdims=['Longitude', 'Latitude'], vdims=['population', 'density']).opts(colorbar=True, - cmap=colorcet.CET_L18, - color='population') + roads_gdf = self._interpolate_traffic_counts(bounds_polygon) + roads_gdf['population_per_hour'] = roads_gdf['population_per_hour'] * relative_variation + roads_gdf['population'] = roads_gdf['population_per_hour'] / 3600 + roads_gdf['density'] = roads_gdf['population'] / roads_gdf.geometry.area + roads_gdf = roads_gdf.set_crs('EPSG:27700').to_crs('EPSG:4326') + + points = gv.Polygons(roads_gdf, + kdims=['Longitude', 'Latitude'], vdims=['population_per_hour', 'density']).opts( + colorbar=True, + cmap=colorcet.CET_L18, + color='population', + line_color='population') bounds = bounds_polygon.bounds raster = rasterize(points, aggregator=ds.mean('density'), width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) - return points, raster_grid, gpd.GeoDataFrame(tfc_df) + return points, raster_grid, gpd.GeoDataFrame(roads_gdf) def clear_cache(self) -> NoReturn: self._traffic_counts = gpd.GeoDataFrame() @@ -151,8 +180,7 @@ def _ingest_road_geometries(self) -> NoReturn: if not self._roads_geometries.crs: self._roads_geometries = self._roads_geometries.set_crs('EPSG:27700') - def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = 20) -> List[ - List[float]]: + def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = 20) -> gpd.GeoDataFrame: """ Interpolate traffic count values between count points along all roads. Interpolation points are created at frequency depending on resolution @@ -161,6 +189,7 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = """ import numpy as np import shapely.ops as so + import pandas as pd b = bounds_poly.bounds bounds = [0] * 4 @@ -168,10 +197,13 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = bounds[2], bounds[3] = self.reverse_proj.transform(b[3], b[2]) unique_road_names = self._roads_geometries.cx[bounds[0]:bounds[2], bounds[1]:bounds[3]]['RoadNumber'].unique() - all_interp_points = [] - all_interp_points_append = all_interp_points.append + all_road_gdfs = [] for road_name in unique_road_names: print(".", end='') + if road_name.startswith('M'): + road_width = 22 + else: + road_width = 14.6 # Get the line segments that make up the road road_segments = self._roads_geometries[self._roads_geometries['RoadNumber'] == road_name].cx[ bounds[0]:bounds[2], bounds[1]:bounds[3]] @@ -201,7 +233,8 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = # If not, carry the length of this segment in the hope the next one will have a counter! carried_length += seg['geometry'].length flat_proj = np.array(flat_proj) - road_ls = so.linemerge(all_segments) # Stitch the line segments together + road_ls = swap_xy(so.linemerge(all_segments)) # Stitch the line segments together + road_width_buffer_gdf = gpd.GeoDataFrame(geometry=[road_ls.buffer(road_width)]) road_length = flat_proj[-1, 0] # Generate some intermediate points to interpolate on coord_spacing = np.linspace(0, road_length, @@ -209,10 +242,19 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = # Interpolate linearly on the 1D projection of the road to estimate the interstitial values interp_flat_proj = np.interp(coord_spacing, flat_proj[:, 0], flat_proj[:, 1]) + point_polys = [] + pops = [] # Recover the true road geometry along with the interpolated values for idx, mark in enumerate(coord_spacing): c = road_ls.interpolate(mark) - all_interp_points_append([*self.proj.transform(c.y, c.x), interp_flat_proj[idx], road_name]) + point_poly = c.buffer(resolution / 2, cap_style=sg.CAP_STYLE.square) + point_polys.append(point_poly) + pops.append(interp_flat_proj[idx]) + + road_gdf = gpd.GeoDataFrame( + {'geometry': point_polys, 'population_per_hour': pops, 'road_name': len(pops) * [road_name]}) + road_gdf = gpd.overlay(road_gdf, road_width_buffer_gdf, how='intersection') + all_road_gdfs.append(road_gdf) - print(".") - return all_interp_points + roads_gdf = pd.concat(all_road_gdfs) + return roads_gdf From 325f12d1c5b40931057616fc51eee49a62a0573e Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Thu, 11 Mar 2021 11:14:26 +0000 Subject: [PATCH 29/57] Remove edges on polygon plots --- seedpod_ground_risk/layers/osm_tag_layer.py | 2 +- seedpod_ground_risk/layers/residential_layer.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/seedpod_ground_risk/layers/osm_tag_layer.py b/seedpod_ground_risk/layers/osm_tag_layer.py index 32a21e6d..e95d74a6 100644 --- a/seedpod_ground_risk/layers/osm_tag_layer.py +++ b/seedpod_ground_risk/layers/osm_tag_layer.py @@ -32,7 +32,7 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr if self.buffer_dist > 0: self._landuse_polygons.geometry = self._landuse_polygons.to_crs('EPSG:27700') \ .buffer(self.buffer_dist).to_crs('EPSG:4326') - polys = gv.Polygons(self._landuse_polygons).opts(style={'alpha': 0.8, 'color': self._colour}) + polys = gv.Polygons(self._landuse_polygons).opts(alpha=0.8, color=self._colour, line_color=self._colour) raster = rasterize(polys, width=raster_shape[0], height=raster_shape[1], x_range=(bounds[1], bounds[3]), y_range=(bounds[0], bounds[2]), dynamic=False) raster_grid = np.copy(list(raster.data.data_vars.items())[0][1].data.astype(np.float)) diff --git a/seedpod_ground_risk/layers/residential_layer.py b/seedpod_ground_risk/layers/residential_layer.py index 56316dd4..97d351e0 100644 --- a/seedpod_ground_risk/layers/residential_layer.py +++ b/seedpod_ground_risk/layers/residential_layer.py @@ -50,8 +50,8 @@ def generate(self, bounds_polygon: sg.Polygon, raster_shape: Tuple[int, int], fr # Construct the GeoViews Polygons gv_polys = gv.Polygons(census_df, kdims=['Longitude', 'Latitude'], vdims=['name', 'population']) \ .opts(color='population', - cmap=colorcet.CET_L18, - colorbar=True, colorbar_opts={'title': 'Population'}, show_legend=False) + cmap=colorcet.CET_L18, alpha=0.8, + colorbar=True, colorbar_opts={'title': 'Population'}, show_legend=False, line_color='population') if self.buffer_dist > 0: buffered_df = deepcopy(census_df) From 66a0364d18b932bf95547582f85107eb35691fa3 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 12 Mar 2021 15:05:38 +0000 Subject: [PATCH 30/57] Fix indices flipping --- seedpod_ground_risk/core/plot_server.py | 3 ++- seedpod_ground_risk/layers/path_analysis_layer.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index a95314fd..eb586dfc 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -222,8 +222,9 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) layer_raster_grid = res[1] nans = np.isnan(layer_raster_grid) layer_raster_grid[nans] = 0 - raster_grid += res[1] + raster_grid += layer_raster_grid raster_grid = np.flipud(raster_grid) + raster_indices['Latitude'] = np.flip(raster_indices['Latitude']) annotations = [] for layer in self.annotation_layers: annotation = layer.annotate(raw_datas, (raster_indices, raster_grid), diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 9f7142ed..3b849305 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -48,7 +48,7 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np path_grid_points = [bresenham.make_line(*pair[0], *pair[1]) for pair in path_pairs] # Bring all these points together and remove duplicate coords # Flip left to right as bresenham spits out in (y,x) order - path_grid_points = np.fliplr(np.unique(np.concatenate(path_grid_points, axis=0), axis=0)) + path_grid_points = np.unique(np.concatenate(path_grid_points, axis=0), axis=0) ca_model = casex.CriticalAreaModels() # Lethal area model using default params for human dimensions aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, 2, 2, 2) # Default aircraft From 1f127b714230b1016ebde1973a2ae6502dbdef9a Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 12 Mar 2021 17:09:26 +0000 Subject: [PATCH 31/57] Add basic lethal area calculation --- .../layers/path_analysis_layer.py | 26 ++++++++++++++----- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 3b849305..ea3b2f49 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -80,13 +80,9 @@ def map_chunk(coords): pool.close() pdf_mat = np.rot90(np.sum(pdfs, axis=0)) - labels = [] - annotation_layers = [] - for gdf in data: - if not gdf.crs: - # If CRS is not set, set EPSG4326 without reprojection as it must be EPSG4326 to display properly - gdf.set_crs(epsg=4326, inplace=True) - overlay = gpd.overlay(gdf, self.buffer_poly, how='intersection') + a_exp = self.get_lethal_area(30) + # Probability * Pop. Density * Lethal Area + risk_map = pdf_mat * raster_data[1] * a_exp geom_type = overlay.geometry.geom_type.all() if geom_type == 'Polygon' or geom_type == 'MultiPolygon': @@ -148,3 +144,19 @@ def _snap_coords_to_grid(self, grid, lon: float, lat: float) -> Tuple[int, int]: lon_idx = int(np.argmin(np.abs(grid['Longitude'] - lon))) return lon_idx, lat_idx + + def get_lethal_area(self, theta: float): + """ + Calculate lethal area of UAS impact from impact angle + + Method from Smith, P.G. 2000 + + :param theta: impact angle in degrees + :return: + """ + r_person = 1 # radius of a person + h_person = 1.8 # height of a person + + r_uas = 0.5 # UAS radius/halfspan + + return ((2 * (r_person + r_uas) * h_person) / np.tan(np.deg2rad(theta))) + (np.pi * (r_uas + r_person) ** 2) From b3908b9c67e5b65995a4faf3d49425693e106385 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 12 Mar 2021 17:10:39 +0000 Subject: [PATCH 32/57] Remove basic path buffer geometries --- .../layers/path_analysis_layer.py | 192 ++++++++++-------- .../layers/pathfinding_layer.py | 4 - 2 files changed, 102 insertions(+), 94 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index ea3b2f49..ff66ef03 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -23,11 +23,6 @@ def __init__(self, key: str, filepath: str = '', buffer_dist: float = None, **kw def preload_data(self) -> NoReturn: self.dataframe = gpd.read_file(self.filepath) - if self.buffer_dist: - epsg27700_geom = self.dataframe.to_crs('EPSG:27700').geometry - self.buffer_poly = gpd.GeoDataFrame( - {'geometry': epsg27700_geom.buffer(self.buffer_dist).to_crs('EPSG:4326')} - ) self.endpoint = self.dataframe.iloc[0].geometry.coords[-1] def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np.array], np.array], @@ -38,96 +33,113 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np bounds = (raster_data[0]['Longitude'].min(), raster_data[0]['Latitude'].min(), raster_data[0]['Longitude'].max(), raster_data[0]['Latitude'].max()) - if self.buffer_poly is not None: - # Snap the line string nodes to the raster grid - snapped_points = [self._snap_coords_to_grid(raster_data[0], *coords) for coords in - self.dataframe.iloc[0].geometry.coords] - # Generate pairs of consecutive (x,y) coords - path_pairs = list(map(list, zip(snapped_points, snapped_points[1:]))) - # Feed these pairs into the Bresenham algo to find the intermediate points - path_grid_points = [bresenham.make_line(*pair[0], *pair[1]) for pair in path_pairs] - # Bring all these points together and remove duplicate coords - # Flip left to right as bresenham spits out in (y,x) order - path_grid_points = np.unique(np.concatenate(path_grid_points, axis=0), axis=0) - - ca_model = casex.CriticalAreaModels() # Lethal area model using default params for human dimensions - aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, 2, 2, 2) # Default aircraft - ballistic_model = casex.BallisticDescent2ndOrderDragApproximation() - ballistic_model.set_aircraft(aircraft) - - # Testing fixed pdf for now - dist_mean = np.array([15, 5]) - - # Create grid on which to evaluate each point of path with its pdf - raster_shape = raster_data[1].shape - x, y = np.mgrid[0:raster_shape[1], 0:raster_shape[0]] - eval_grid = np.dstack((x, y)) - - from pathos.multiprocessing import Pool - from pathos.multiprocessing import cpu_count - pool = Pool() - - # TODO Use something like Dask or OpenCV to speed this up in future as it's a simple map-reduce - def point_dist(c): - pdf_mean = c + dist_mean - return mvn(pdf_mean, [[15, 5], [5, 7]]).pdf(eval_grid) - - def map_chunk(coords): - return np.sum([point_dist(c) for c in coords], axis=0) - - # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? - pdfs = pool.map(map_chunk, np.array_split(path_grid_points, cpu_count() * 2)) - pool.close() - pdf_mat = np.rot90(np.sum(pdfs, axis=0)) + # Snap the line string nodes to the raster grid + snapped_points = [self._snap_coords_to_grid(raster_data[0], *coords) for coords in + self.dataframe.iloc[0].geometry.coords] + # Generate pairs of consecutive (x,y) coords + path_pairs = list(map(list, zip(snapped_points, snapped_points[1:]))) + # Feed these pairs into the Bresenham algo to find the intermediate points + path_grid_points = [bresenham.make_line(*pair[0], *pair[1]) for pair in path_pairs] + # Bring all these points together and remove duplicate coords + # Flip left to right as bresenham spits out in (y,x) order + path_grid_points = np.unique(np.concatenate(path_grid_points, axis=0), axis=0) + + ca_model = casex.CriticalAreaModels() # Lethal area model using default params for human dimensions + aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, 2, 2, 2) # Default aircraft + aircraft.set_ballistic_drag_coefficient(0.8) + aircraft.set_ballistic_frontal_area(3) + fc = casex.FrictionCoefficients() + aircraft.set_friction_coefficient( + fc.get_coefficient(casex.enums.AircraftMaterial.ALUMINUM, casex.enums.GroundMaterial.CONCRETE)) + ballistic_model = casex.BallisticDescent2ndOrderDragApproximation() + ballistic_model.set_aircraft(aircraft) + distance_impact, velocity_impact, angle_impact, time_impact = ballistic_model.compute_ballistic_distance( + 133, 20, 1) + ca_model.critical_area(casex.enums.CriticalAreaModel.JARUS, + aircraft, velocity_impact, + angle_impact, 0) + + # Testing fixed pdf for now + dist_mean = np.array([-5, 5]) + + # Create grid on which to evaluate each point of path with its pdf + raster_shape = raster_data[1].shape + x, y = np.mgrid[0:raster_shape[0], 0:raster_shape[1]] + eval_grid = np.dstack((x, y)) + + from pathos.multiprocessing import Pool + from pathos.multiprocessing import cpu_count + pool = Pool() + + # TODO Use something like Dask or OpenCV to speed this up in future as it's a simple map-reduce + def point_dist(c): + pdf_mean = c + dist_mean + return mvn(pdf_mean, [[15, 5], [5, 7]]).pdf(eval_grid) + + def map_chunk(coords): + return np.sum([point_dist(c) for c in coords], axis=0) + + # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? + pdfs = pool.map(map_chunk, np.array_split(path_grid_points, cpu_count() * 2)) + pool.close() + pdf_mat = np.sum(pdfs, axis=0) a_exp = self.get_lethal_area(30) # Probability * Pop. Density * Lethal Area risk_map = pdf_mat * raster_data[1] * a_exp - geom_type = overlay.geometry.geom_type.all() - if geom_type == 'Polygon' or geom_type == 'MultiPolygon': - if 'density' not in gdf: - continue - proj_gdf = overlay.to_crs('epsg:27700') - proj_gdf['population'] = proj_gdf.geometry.area * proj_gdf.density - labels.append(gv.Text(self.endpoint[0], self.endpoint[1], - f"Static Population: {proj_gdf['population'].sum():.2f}", fontsize=20).opts( - style={'color': 'blue'})) - proj_gdf = proj_gdf.to_crs('EPSG:4326') - annotation_layers.append( - gv.Polygons(proj_gdf, vdims=['population']).opts(alpha=0.8, color='cyan', tools=['hover'])) - elif geom_type == 'Point': - if 'population' in overlay: - # Group data by road name, localising the points - road_geoms = list(overlay.groupby('road_name').geometry) - road_coms = {} # store the centre of mass of points associated with a road - for name, geoms in road_geoms: - mean_lon = sum([p.x for p in geoms]) / len(geoms) - mean_lat = sum([p.y for p in geoms]) / len(geoms) - road_coms[name] = (mean_lon, mean_lat) # x,y order - # get mean population flow of points for road - road_pops = list(overlay.groupby('road_name').mean().itertuples()) - for name, pop, *_ in road_pops: # add labels at the centre of mass of each group of points - labels.append(gv.Text(road_coms[name][0], road_coms[name][1], - f'Population flow on road {name}: {pop * 60:.2f}/min', - fontsize=20).opts( - style={'color': 'blue'})) - annotation_layers.append((gv.Points(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'}))) - elif geom_type == 'Line': - pass - - return Overlay([ - gv.Contours(self.dataframe).opts(line_width=4, line_color='magenta'), - # gv.Polygons(self.buffer_poly).opts(style={'alpha': 0.3, 'color': 'cyan'}), - # Add the path stats as a text annotation to the final path point - *labels, - *annotation_layers, - gv.Image(pdf_mat, bounds=bounds).options(alpha=0.4, clipping_colors={'0': 'transparent', - 'NaN': 'transparent', - '-NaN': 'transparent'}) - ]) - else: - return gv.Contours(self.dataframe).opts(line_width=4, line_color='#000000') + risk_raster = gv.Image(risk_map, bounds=bounds).options(alpha=0.7, cmap='viridis', tools=['hover'], + clipping_colors={'min': (0, 0, 0, 0)}) + risk_raster = risk_raster.redim.range(z=(1e-11, risk_map.max())) + + # labels = [] + # annotation_layers = [] + # for gdf in data: + # if not gdf.crs: + # # If CRS is not set, set EPSG4326 without reprojection as it must be EPSG4326 to display properly + # gdf.set_crs(epsg=4326, inplace=True) + # overlay = gpd.overlay(gdf, self.buffer_poly, how='intersection') + # + # geom_type = overlay.geometry.geom_type.all() + # if geom_type == 'Polygon' or geom_type == 'MultiPolygon': + # if 'density' not in gdf: + # continue + # proj_gdf = overlay.to_crs('epsg:27700') + # proj_gdf['population'] = proj_gdf.geometry.area * proj_gdf.density + # labels.append(gv.Text(self.endpoint[0], self.endpoint[1], + # f"Static Population: {proj_gdf['population'].sum():.2f}", fontsize=20).opts( + # style={'color': 'blue'})) + # proj_gdf = proj_gdf.to_crs('EPSG:4326') + # annotation_layers.append( + # gv.Polygons(proj_gdf, vdims=['population']).opts(alpha=0.8, color='cyan', tools=['hover'])) + # elif geom_type == 'Point': + # if 'population' in overlay: + # # Group data by road name, localising the points + # road_geoms = list(overlay.groupby('road_name').geometry) + # road_coms = {} # store the centre of mass of points associated with a road + # for name, geoms in road_geoms: + # mean_lon = sum([p.x for p in geoms]) / len(geoms) + # mean_lat = sum([p.y for p in geoms]) / len(geoms) + # road_coms[name] = (mean_lon, mean_lat) # x,y order + # # get mean population flow of points for road + # road_pops = list(overlay.groupby('road_name').mean().itertuples()) + # for name, pop, *_ in road_pops: # add labels at the centre of mass of each group of points + # labels.append(gv.Text(road_coms[name][0], road_coms[name][1], + # f'Population flow on road {name}: {pop * 60:.2f}/min', + # fontsize=20).opts( + # style={'color': 'blue'})) + # annotation_layers.append((gv.Points(overlay).opts(style={'alpha': 0.8, 'color': 'cyan'}))) + # elif geom_type == 'Line': + # pass + + return Overlay([ + gv.Contours(self.dataframe).opts(line_width=4, line_color='magenta'), + # gv.Polygons(self.buffer_poly).opts(style={'alpha': 0.3, 'color': 'cyan'}), + # Add the path stats as a text annotation to the final path point + # *labels, + # *annotation_layers, + risk_raster + ]) def clear_cache(self) -> NoReturn: pass diff --git a/seedpod_ground_risk/layers/pathfinding_layer.py b/seedpod_ground_risk/layers/pathfinding_layer.py index 479d47fd..b75c8af4 100644 --- a/seedpod_ground_risk/layers/pathfinding_layer.py +++ b/seedpod_ground_risk/layers/pathfinding_layer.py @@ -77,10 +77,6 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np snapped_path = sg.LineString(snapped_path) self.dataframe = gpd.GeoDataFrame({'geometry': [snapped_path]}).set_crs('EPSG:4326') - epsg27700_geom = self.dataframe.to_crs('EPSG:27700').geometry - self.buffer_poly = gpd.GeoDataFrame( - {'geometry': epsg27700_geom.buffer(self.buffer_dist).to_crs('EPSG:4326')} - ) self.endpoint = self.dataframe.iloc[0].geometry.coords[-1] return super(PathfindingLayer, self).annotate(data, raster_data, **kwargs) From 2a64a6e0773a939b639ce68a03330957cd55a92e Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 12 Mar 2021 17:11:17 +0000 Subject: [PATCH 33/57] Add better traceback printing from worker threads --- seedpod_ground_risk/core/plot_server.py | 5 ++++- seedpod_ground_risk/layers/roads_layer.py | 1 - 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index eb586dfc..6576a0c8 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -9,7 +9,6 @@ from seedpod_ground_risk.layers.annotation_layer import AnnotationLayer from seedpod_ground_risk.layers.data_layer import DataLayer from seedpod_ground_risk.layers.layer import Layer -from seedpod_ground_risk.layers.path_analysis_layer import PathAnalysisLayer def make_bounds_polygon(*args: Iterable[float]) -> sg.Polygon: @@ -258,6 +257,8 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) except Exception as e: # Catch-all to prevent plot blanking out and/or crashing app # Just display map tiles in case this was transient + import traceback + traceback.print_exc() print(e) plot = self._base_tiles @@ -322,6 +323,8 @@ def generate_layer(layer: DataLayer, bounds_poly: sg.Polygon, raster_shape: Tupl resolution=resolution) return result except Exception as e: + import traceback + traceback.print_tb(e.__traceback__) print(e) return layer.key + ' FAILED', None diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index c2b2cb58..8da6ef8b 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -199,7 +199,6 @@ def _interpolate_traffic_counts(self, bounds_poly: sg.Polygon, resolution: int = all_road_gdfs = [] for road_name in unique_road_names: - print(".", end='') if road_name.startswith('M'): road_width = 22 else: From fed6afe0849f39e53f4cb3414d5a813cd2012fc3 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 12 Mar 2021 17:12:14 +0000 Subject: [PATCH 34/57] Handle LoS through blocked areas for path smoothing --- seedpod_ground_risk/pathfinding/a_star.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/pathfinding/a_star.py b/seedpod_ground_risk/pathfinding/a_star.py index 3b64a29c..aa23b108 100644 --- a/seedpod_ground_risk/pathfinding/a_star.py +++ b/seedpod_ground_risk/pathfinding/a_star.py @@ -56,12 +56,17 @@ def _reconstruct_path(self, end: Node, closed_list: Dict[Node, Node], grid: np.n # mpl.title(f'Full Path, length={len(path)}') # mpl.show() - return path + # return path # @jit(nopython=True, nogil=True) def get_path_sum(nx, ny, tx, ty, grid): line = bresenham.make_line(nx, ny, tx, ty) - return grid[line[:, 0], line[:, 1]].sum() + line_points = grid[line[:, 0], line[:, 1]] + # If the new line crosses any blocked areas the cost is inf + if -1 in line_points: + return np.inf + else: + return line_points.sum() # integral_val = 0 # for y, x in line: # integral_val += grid[y, x] From 6fa25156974904035835d5786eeaf2349ed1e338 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Fri, 12 Mar 2021 17:14:34 +0000 Subject: [PATCH 35/57] Cache JIT'd functions --- seedpod_ground_risk/pathfinding/bresenham.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/seedpod_ground_risk/pathfinding/bresenham.py b/seedpod_ground_risk/pathfinding/bresenham.py index 2a86a626..ec801289 100644 --- a/seedpod_ground_risk/pathfinding/bresenham.py +++ b/seedpod_ground_risk/pathfinding/bresenham.py @@ -4,7 +4,7 @@ # Reference: https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm -@jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True, cache=True) def _make_line_low(x0: int, y0: int, x1: int, y1: int): dx = x1 - x0 dy = y1 - y0 @@ -30,7 +30,7 @@ def _make_line_low(x0: int, y0: int, x1: int, y1: int): return line -@jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True, cache=True) def _make_line_high(x0: int, y0: int, x1: int, y1: int): dx = x1 - x0 dy = y1 - y0 @@ -56,7 +56,7 @@ def _make_line_high(x0: int, y0: int, x1: int, y1: int): return line -@jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True, cache=True) def make_line(x0: int, y0: int, x1: int, y1: int) -> np.array: """ Return a rasterised line between pairs of coordinates on a 2d integer grid according to the Bresenham algorithm. From 05480eed9778f13b109f863c47f059812938753e Mon Sep 17 00:00:00 2001 From: Zach10a <35232913+Zach10a@users.noreply.github.com> Date: Fri, 12 Mar 2021 19:50:58 +0000 Subject: [PATCH 36/57] Feature edit layer (#18) * Added edit layer button * Add wizard into edit button --- seedpod_ground_risk/main.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/seedpod_ground_risk/main.py b/seedpod_ground_risk/main.py index 6af111d1..8ff1e5b8 100644 --- a/seedpod_ground_risk/main.py +++ b/seedpod_ground_risk/main.py @@ -74,6 +74,15 @@ def change_colour(self): def delete_layer(self): self._plot_worker.remove_layer(self._layer) + def edit_layer(self): + self._plot_worker.remove_layer(self._layer) + wizard = LayerWizard(self, Qt.Window) + wizard.exec_() # Open wizard and block until result + if wizard.result() == QDialog.Accepted: + layerObj = list(LAYER_OBJECTS.values())[wizard.layerType] + layer = layerObj(wizard.layerKey, **wizard.opts) + self._plot_worker.add_layer(layer) + def export_path_json(self): from PySide2.QtWidgets import QFileDialog @@ -87,6 +96,7 @@ def mousePressEvent(self, event: PySide2.QtGui.QMouseEvent) -> None: if event.button() == Qt.RightButton: menu = QMenu() menu.addAction("Delete", self.delete_layer) + menu.addAction("Edit Layer", self.edit_layer) if isinstance(self._layer, AnnotationLayer): menu.addAction("Export .GeoJSON", self.export_path_json) menu.exec_(event.globalPos()) From 8a4361c21e3c4ca0c4d7262e4bcc3116077f2963 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:34:45 +0000 Subject: [PATCH 37/57] Add ballistic model multivariate gaussian ground impact distribution estimation --- .../path_analysis/ballistic_model.py | 90 +++++++++++ tests/path_analysis/test_ballistic_model.py | 147 ++++++++++++++++++ 2 files changed, 237 insertions(+) create mode 100644 seedpod_ground_risk/path_analysis/ballistic_model.py create mode 100644 tests/path_analysis/test_ballistic_model.py diff --git a/seedpod_ground_risk/path_analysis/ballistic_model.py b/seedpod_ground_risk/path_analysis/ballistic_model.py new file mode 100644 index 00000000..64e5918d --- /dev/null +++ b/seedpod_ground_risk/path_analysis/ballistic_model.py @@ -0,0 +1,90 @@ +from casex import * +from numba import njit +from sklearn.mixture import GaussianMixture + +from seedpod_ground_risk.path_analysis.descent_model import DescentModel +from seedpod_ground_risk.path_analysis.utils import rotate_2d + + +@njit(cache=True) +def paef_to_ned_with_wind(x): + """ + Transform PAE frame distances to NED frame and transform with wind. + This func is designed to be used in np apply, hence the single arg. The column ordering is very specific! + :param x: array row with ordering [paef_y (always 0), paef_x, impact_time, theta (rad), wind_vel_x, wind_vel_y] + :return: + """ + paef_c = x[0:2] + t_i = x[2] + theta = x[3] + wind_vect = x[4:6] + return rotate_2d(paef_c, theta) + wind_vect * t_i + + +class BallisticModel(DescentModel): + + def __init__(self, aircraft: AircraftSpecs, n_samples: int = 2000) -> None: + super().__init__(aircraft, n_samples) + + self.bm = BallisticDescent2ndOrderDragApproximation() + self.bm.set_aircraft(aircraft) + + def impact_distance_dist_params_ned_with_wind(self, altitude, velocity, heading, wind_vel_y, wind_vel_x, loc_x, + loc_y): + """ + Return the parameters to a multivariate normal distribution describing the ground impact probability under a + ballistic descent. + + This function takes into account wind and returns the result in the NED frame with x, y corresponding to + East and North respectively. The distribution takes in the location of the event in the existing NED frame + and transforms the Path aligned event frame (PAEF) to the NED frame at the specified location. + + If passing an array, all other arrays must be the same shape. This is usually a single dimension of samples + generated with scipy.stats..rvs + + The method is as follows: + 1. The ballistic model return one dimensional results from the specified params, if these are arrays then + a number of samples are created. + 2. The heading(s) are rotated into the NED frame + 3. A vectorised operation is performed to firstly rotate the PAEF results into the NED frame + 4. The second part of the vectorised operation then multiplies the wind vector (in NED) by the time + taken to impact the ground. This is then added to the first part. + 5. The samples are used to fit a multivariate Gaussian from which the parameters are generated. + + :param altitude: the altitude in metres + :type altitude: float or np.array + :param velocity: the velocity over the ground of the aircraft in the direction of flight in m/s + :type velocity: float or np.array + :param heading: the ground track angle of the aircraf in deg + :type heading: float or np.array + :param wind_vel_x: the x component of the wind in m/s + :type wind_vel_x: float or nd.array + :param wind_vel_y: the y component of the wind in m/s + :type wind_vel_y: float or nd.array + :param loc_x: event x location + :type loc_x: int + :param loc_y: event y location + :type loc_y: int + :return: a tuple of (means, covariances) of the distribution + :rtype: tuple of np.arrays of shape (2,) for the means and (2,2) for the covariances + """ + # Compute impact distances and times in the PAE frame + # The velocity vector is assumed to be aligned with path vector, hence v_y is 0 + d_i, _, _, t_i = self.bm.compute_ballistic_distance(altitude, velocity, 0) + + # Compensate for x,y axes being rotated compared to bearings + theta = (heading - (np.pi / 2)) * (2 * np.pi) + # Get angle distribution in between body and NED frame + # Form the array structure required and transform + arr = np.vstack((np.zeros(d_i.shape), d_i, t_i, theta, wind_vel_x, wind_vel_y)) + transformed_arr = np.apply_along_axis(paef_to_ned_with_wind, 0, arr) + gm = GaussianMixture() + gm.fit_predict(transformed_arr.T) + # If there the event and NED origins match, no need to translate + if not loc_x or not loc_y: + means = gm.means_[0] + np.array([loc_x, loc_y]) + else: + means = gm.means_[0] + # Gaussian Mixture model can deal with up to 3D distributions, but we are only dealing with 2D here, + # so take first index into the depth + return means, gm.covariances_[0] diff --git a/tests/path_analysis/test_ballistic_model.py b/tests/path_analysis/test_ballistic_model.py new file mode 100644 index 00000000..736fe57b --- /dev/null +++ b/tests/path_analysis/test_ballistic_model.py @@ -0,0 +1,147 @@ +import unittest + +import scipy.stats as ss +from casex import * + +from seedpod_ground_risk.path_analysis.ballistic_model import BallisticModel + + +class BallisticModelPAEFTestCase(unittest.TestCase): + + def setUp(self) -> None: + super().setUp() + self.ac = AircraftSpecs(enums.AircraftType.FIXED_WING, + 2, # width + 1.8, # length + 7 # mass + ) + self.ac.set_ballistic_frontal_area(2) + self.ac.set_glide_speed_ratio(15, 12) + self.ac.set_glide_drag_coefficient(0.3) + self.ac.set_ballistic_drag_coefficient(1.1) + + self.bm = BallisticDescent2ndOrderDragApproximation() + self.bm.set_aircraft(self.ac) + + def test_ballistic_dist(self): + """ + Test ballistic model impact distance distributions in the Path Aligned Event frame + """ + alt_mean = 50 + alt_std = 5 + vx_mean = 18 + vx_std = 2.5 + vy = 1 + samples = 2000 + + # Compute ballistic distances in the path aligned LTP frame with origin at the directly below the event location + # AKA PAEF + d_i, v_i, a_i, t_i = self.bm.compute_ballistic_distance(ss.norm(alt_mean, alt_std).rvs(samples), + ss.norm(vx_mean, vx_std).rvs(samples), vy) + + di_mean, di_std = ss.norm.fit(d_i) + self.assertAlmostEqual(di_mean, 13.4, delta=0.5) + self.assertAlmostEqual(di_std, 1.5, delta=0.2) + + def test_ballistic_time(self): + """ + Test ballistic model impact time distributions in the Path Aligned Event frame + """ + alt_mean = 50 + alt_std = 5 + vx_mean = 18 + vx_std = 2.5 + vy = 1 + samples = 2000 + + # Compute ballistic distances in the path aligned LTP frame with origin at the directly below the event location + # AKA PAEF + d_i, v_i, a_i, t_i = self.bm.compute_ballistic_distance(ss.norm(alt_mean, alt_std).rvs(samples), + ss.norm(vx_mean, vx_std).rvs(samples), vy) + + ti_mean, ti_std = ss.norm.fit(t_i) + self.assertAlmostEqual(ti_mean, 7.4, delta=0.5) + self.assertAlmostEqual(ti_std, 0.7, delta=0.2) + + def test_ballistic_vel(self): + """ + Test ballistic model impact velocity distributions in the Path Aligned Event frame + """ + alt_mean = 50 + alt_std = 5 + vx_mean = 18 + vx_std = 2.5 + vy = 1 + samples = 2000 + + # Compute ballistic distances in the path aligned LTP frame with origin at the directly below the event location + # AKA PAEF + d_i, v_i, a_i, t_i = self.bm.compute_ballistic_distance(ss.norm(alt_mean, alt_std).rvs(samples), + ss.norm(vx_mean, vx_std).rvs(samples), vy) + + vi_mean, vi_std = ss.norm.fit(v_i) + self.assertAlmostEqual(vi_mean, 7.1, delta=0.05) + self.assertAlmostEqual(vi_std, 0.2, delta=0.2) + + +class BallisticModelNEDWindTestCase(unittest.TestCase): + + def setUp(self) -> None: + super().setUp() + self.ac = AircraftSpecs(enums.AircraftType.FIXED_WING, + 2, # width + 1.8, # length + 7 # mass + ) + self.ac.set_ballistic_frontal_area(2) + self.ac.set_glide_speed_ratio(15, 12) + self.ac.set_glide_drag_coefficient(0.3) + self.ac.set_ballistic_drag_coefficient(1.1) + + self.bm = BallisticDescent2ndOrderDragApproximation() + self.bm.set_aircraft(self.ac) + + def test_ballistic_dist(self): + """ + Profile ballistic model impact distance distributions in the North East Down frame with wind compensation + """ + make_plot = False # Set flag to plot result + samples = 5000 + + # Conjure up our distributions for various things + alt_mean = 50 + alt_std = 5 + alt = ss.norm(alt_mean, alt_std).rvs(samples) + vx_mean = 18 + vx_std = 2.5 + vel = ss.norm(vx_mean, vx_std).rvs(samples) + heading_mean = np.deg2rad(270) + heading_std = np.deg2rad(2) + heading = ss.norm(heading_mean, heading_std).rvs(samples) + + loc_x, loc_y = 0, 50 + + wind_vel_x = ss.norm(5, 1).rvs(samples) + wind_vel_y = ss.norm(1, 1).rvs(samples) + + bm = BallisticModel(self.ac) + means, cov = bm.impact_distance_dist_params_ned_with_wind(alt, vel, heading, wind_vel_x, wind_vel_y, loc_x, + loc_y) + pdf = ss.multivariate_normal(mean=means, cov=cov).pdf + + if make_plot: + # Make a sampling grid for plotting + x, y = np.mgrid[(loc_x - 100):(loc_x + 100), (loc_y - 100):(loc_y + 100)] + pos = np.vstack([x.ravel(), y.ravel()]) + # Sample KDE PDF on these points + density = pdf(pos.T) + # Plot sampled KDE PDF + import matplotlib.pyplot as mpl + fig, ax = mpl.subplots(1, 1) + sc = ax.scatter(x, y, c=density) + fig.colorbar(sc) + fig.show() + + +if __name__ == '__main__': + unittest.main() From 154a258c41eae3e05a80b6fb110eac4f8c34e131 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:36:20 +0000 Subject: [PATCH 38/57] Add descent model utils --- .../path_analysis/descent_model.py | 8 +++ seedpod_ground_risk/path_analysis/utils.py | 34 ++++++++++ tests/path_analysis/test_grid_ops.py | 48 ++++++++++++++ tests/path_analysis/test_vector_ops.py | 66 +++++++++++++++++++ 4 files changed, 156 insertions(+) create mode 100644 seedpod_ground_risk/path_analysis/descent_model.py create mode 100644 seedpod_ground_risk/path_analysis/utils.py create mode 100644 tests/path_analysis/test_grid_ops.py create mode 100644 tests/path_analysis/test_vector_ops.py diff --git a/seedpod_ground_risk/path_analysis/descent_model.py b/seedpod_ground_risk/path_analysis/descent_model.py new file mode 100644 index 00000000..08e4ed1b --- /dev/null +++ b/seedpod_ground_risk/path_analysis/descent_model.py @@ -0,0 +1,8 @@ +from casex import AircraftSpecs + + +class DescentModel: + + def __init__(self, aircraft: AircraftSpecs, n_samples: int = 2000) -> None: + self.aircraft = aircraft + self.n_samples = n_samples diff --git a/seedpod_ground_risk/path_analysis/utils.py b/seedpod_ground_risk/path_analysis/utils.py new file mode 100644 index 00000000..326ed84b --- /dev/null +++ b/seedpod_ground_risk/path_analysis/utils.py @@ -0,0 +1,34 @@ +from typing import Tuple + +import numpy as np +from numba import njit +from numba.np.arraymath import cross2d + + +@njit(cache=True) +def rotate_2d(vec: np.array, theta: float): + """ + Rotate a 2D vector anticlockwise by a given angle + :param vec: vector to rotate + :param theta: angle in radians to rotate anticlockwise + :return: rotated vector with same shape + """ + if theta < 0: + theta = 2 * np.pi - theta + c, s = np.cos(theta), np.sin(theta) + R = np.array(((c, -s), (s, c))) + return cross2d(R, vec) + + +def snap_coords_to_grid(grid, lon: float, lat: float) -> Tuple[int, int]: + """ + Snap coordinates to grid indices + :param grid: raster grid coordinates + :param lon: longitude to snap + :param lat: latitude to snap + :return: (x, y) tuple of grid indices + """ + lat_idx = int(np.argmin(np.abs(grid['Latitude'] - lat))) + lon_idx = int(np.argmin(np.abs(grid['Longitude'] - lon))) + + return lon_idx, lat_idx diff --git a/tests/path_analysis/test_grid_ops.py b/tests/path_analysis/test_grid_ops.py new file mode 100644 index 00000000..fad3f1a6 --- /dev/null +++ b/tests/path_analysis/test_grid_ops.py @@ -0,0 +1,48 @@ +import unittest + +import numpy as np + +from seedpod_ground_risk.layers.path_analysis_layer import snap_coords_to_grid + + +class GridSnappingTestCase(unittest.TestCase): + + def test_normal(self): + grid = { + 'Latitude': np.array(list(range(11))), + 'Longitude': np.array(list(range(10, 21))) + } + lat = 5.1 + lon = 15.1 + lon_idx, lat_idx = snap_coords_to_grid(grid, lon, lat) + + self.assertEqual(lon_idx, 5) + self.assertEqual(lat_idx, 5) + + def test_out_of_bounds_max(self): + grid = { + 'Latitude': np.array(list(range(11))), + 'Longitude': np.array(list(range(10, 21))) + } + lat = 50.1 + lon = 55.1 + lon_idx, lat_idx = snap_coords_to_grid(grid, lon, lat) + + self.assertEqual(lon_idx, 10) + self.assertEqual(lat_idx, 10) + + def test_out_of_bounds_min(self): + grid = { + 'Latitude': np.array(list(range(11))), + 'Longitude': np.array(list(range(10, 21))) + } + lat = -50.1 + lon = -55.1 + lon_idx, lat_idx = snap_coords_to_grid(grid, lon, lat) + + self.assertEqual(lon_idx, 0) + self.assertEqual(lat_idx, 0) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/path_analysis/test_vector_ops.py b/tests/path_analysis/test_vector_ops.py new file mode 100644 index 00000000..06232e88 --- /dev/null +++ b/tests/path_analysis/test_vector_ops.py @@ -0,0 +1,66 @@ +import unittest + +import numpy as np + +from seedpod_ground_risk.path_analysis.utils import rotate_2d + + +class VectorRotationTestCase(unittest.TestCase): + + def test_first_quad(self): + """ + Test rotation of vectors. Used in transformation between frames + """ + theta = np.deg2rad(45) + vec = np.array([0, 1]) # y, x order + out = rotate_2d(vec, theta) + val = np.sqrt(2) / 2 + np.testing.assert_array_equal(out, np.array([val, val])) + + def test_sec_quad(self): + """ + Test rotation of vectors. Used in transformation between frames + """ + theta = np.deg2rad(45) + vec = np.array([1, 0]) # y, x order + out = rotate_2d(vec, theta) + val = np.sqrt(2) / 2 + np.testing.assert_array_equal(out, np.array([val, -val])) + + def test_third_quad(self): + """ + Test rotation of vectors. Used in transformation between frames + """ + theta = np.deg2rad(45) + vec = np.array([0, -1]) # y, x order + out = rotate_2d(vec, theta) + val = np.sqrt(2) / 2 + np.testing.assert_array_equal(out, np.array([-val, -val])) + + def test_fourth_quad(self): + """ + Test rotation of vectors. Used in transformation between frames + """ + theta = np.deg2rad(45) + vec = np.array([-1, 0]) # y, x order + out = rotate_2d(vec, theta) + val = np.sqrt(2) / 2 + np.testing.assert_array_equal(out, np.array([-val, val])) + + def test_negative_theta(self): + theta = np.deg2rad(-45) + vec = np.array([-1, 0]) # y, x order + out = rotate_2d(vec, theta) + val = np.sqrt(2) / 2 + np.testing.assert_array_almost_equal(out, np.array([-val, val])) + + def test_over_2pi(self): + theta = np.deg2rad(45) + (2 * np.pi) + vec = np.array([0, 1]) # y, x order + out = rotate_2d(vec, theta) + val = np.sqrt(2) / 2 + np.testing.assert_array_almost_equal(out, np.array([val, val])) + + +if __name__ == '__main__': + unittest.main() From 865191cc49a0a1d03e8cf3490e88675cdde07eb9 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:36:57 +0000 Subject: [PATCH 39/57] Refactor to use util functions --- seedpod_ground_risk/layers/pathfinding_layer.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/seedpod_ground_risk/layers/pathfinding_layer.py b/seedpod_ground_risk/layers/pathfinding_layer.py index b75c8af4..215f8476 100644 --- a/seedpod_ground_risk/layers/pathfinding_layer.py +++ b/seedpod_ground_risk/layers/pathfinding_layer.py @@ -7,6 +7,7 @@ from holoviews.element import Geometry from seedpod_ground_risk.layers.path_analysis_layer import PathAnalysisLayer +from seedpod_ground_risk.path_analysis.utils import snap_coords_to_grid from seedpod_ground_risk.pathfinding.a_star import RiskGridAStar from seedpod_ground_risk.pathfinding.algorithm import Algorithm from seedpod_ground_risk.pathfinding.heuristic import ManhattanRiskHeuristic, Heuristic @@ -34,13 +35,13 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np raster_grid = raster_data[1] * resolution ** 2 - snapped_start_lon_idx, snapped_start_lat_idx = self._snap_coords_to_grid(raster_data[0], self.start_coords[1], - self.start_coords[0]) + snapped_start_lon_idx, snapped_start_lat_idx = snap_coords_to_grid(raster_data[0], self.start_coords[1], + self.start_coords[0]) start_node = environment.Node(snapped_start_lon_idx, snapped_start_lat_idx, raster_grid[snapped_start_lat_idx, snapped_start_lon_idx]) - snapped_end_lon_idx, snapped_end_lat_idx = self._snap_coords_to_grid(raster_data[0], self.end_coords[1], - self.end_coords[0]) + snapped_end_lon_idx, snapped_end_lat_idx = snap_coords_to_grid(raster_data[0], self.end_coords[1], + self.end_coords[0]) end_node = environment.Node(snapped_end_lon_idx, snapped_end_lat_idx, raster_grid[snapped_end_lat_idx, snapped_end_lon_idx]) From 352c14c054cbbe05845902bd31a06c02a1c7d9e0 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:39:05 +0000 Subject: [PATCH 40/57] Fix static data import paths --- seedpod_ground_risk/layers/residential_layer.py | 5 +++-- seedpod_ground_risk/layers/roads_layer.py | 6 +++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/seedpod_ground_risk/layers/residential_layer.py b/seedpod_ground_risk/layers/residential_layer.py index 97d351e0..4c27f518 100644 --- a/seedpod_ground_risk/layers/residential_layer.py +++ b/seedpod_ground_risk/layers/residential_layer.py @@ -77,13 +77,14 @@ def ingest_census_data(self) -> NoReturn: import os # Import Census boundaries in Ordnance Survey grid and reproject - census_wards_df = gpd.read_file(os.sep.join(('static_data', 'england_wa_2011_clipped.shp'))).drop( + census_wards_df = gpd.read_file( + os.sep.join(('..', 'static_data', 'england_wa_2011_clipped.shp'))).drop( ['altname', 'oldcode'], axis=1) if not census_wards_df.crs: census_wards_df = census_wards_df.set_crs('EPSG:27700') census_wards_df = census_wards_df.to_crs('EPSG:4326') # Import census ward densities - density_df = pd.read_csv(os.sep.join(('static_data', 'density.csv')), header=0) + density_df = pd.read_csv(os.sep.join(('..', 'static_data', 'density.csv')), header=0) # Scale from hectares to m^2 density_df['area'] = density_df['area'] * 10000 density_df['density'] = density_df['density'] / 10000 diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index 8da6ef8b..0a3c3e43 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -136,7 +136,7 @@ def _ingest_traffic_counts(self) -> NoReturn: import os # Ingest raw data - counts_df = pd.read_csv(os.sep.join(('static_data', 'dft_traffic_counts_aadf.csv'))) + counts_df = pd.read_csv(os.sep.join(('..', 'static_data', 'dft_traffic_counts_aadf.csv'))) # Select only desired columns counts_df = counts_df[TRAFFIC_COUNT_COLUMNS] # Groupby year and select out only the latest year @@ -164,7 +164,7 @@ def _ingest_relative_traffic_variations(self): import os # Ingest data, ignoring header and footer info - relative_variations_df = pd.read_excel(os.sep.join(('static_data', 'tra0307.ods')), engine='odf', + relative_variations_df = pd.read_excel(os.sep.join(('..', 'static_data', 'tra0307.ods')), engine='odf', header=5, skipfooter=8) # Flatten into continuous list of hourly variations for the week self.relative_variations_flat = (relative_variations_df.iloc[:, 1:] / 100).melt()['value'] @@ -175,7 +175,7 @@ def _ingest_road_geometries(self) -> NoReturn: """ import os - self._roads_geometries = gpd.read_file(os.sep.join(('static_data', '2018-MRDB-minimal.shp'))).rename( + self._roads_geometries = gpd.read_file(os.sep.join(('..', 'static_data', '2018-MRDB-minimal.shp'))).rename( columns={'CP_Number': 'count_point_id'}) if not self._roads_geometries.crs: self._roads_geometries = self._roads_geometries.set_crs('EPSG:27700') From ca292d190ff5248301d7461c6122ae52dfa5acf0 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:45:37 +0000 Subject: [PATCH 41/57] Add lethal area func --- .../layers/path_analysis_layer.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index ff66ef03..143ecee5 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -9,6 +9,23 @@ from seedpod_ground_risk.pathfinding import bresenham +def get_lethal_area(theta: float, uas_width: float): + """ + Calculate lethal area of UAS impact from impact angle + + Method from :cite: Smith, P.G. 2000 + + :param theta: impact angle in degrees + :param uas_width: UAS width in metres + :return: + """ + r_person = 1 # radius of a person + h_person = 1.8 # height of a person + r_uas = uas_width / 2 # UAS halfspan + + return ((2 * (r_person + r_uas) * h_person) / np.tan(np.deg2rad(theta))) + (np.pi * (r_uas + r_person) ** 2) + + class PathAnalysisLayer(AnnotationLayer): def __init__(self, key: str, filepath: str = '', buffer_dist: float = None, **kwargs): From 3f0b065797e72b45a2a31536ddfc7bfca5e75529 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:46:48 +0000 Subject: [PATCH 42/57] Find headings for path segments --- seedpod_ground_risk/layers/path_analysis_layer.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 143ecee5..2a58ff3b 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -50,11 +50,19 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np bounds = (raster_data[0]['Longitude'].min(), raster_data[0]['Latitude'].min(), raster_data[0]['Longitude'].max(), raster_data[0]['Latitude'].max()) + line_coords = list(self.dataframe.iloc[0].geometry.coords) # Snap the line string nodes to the raster grid - snapped_points = [self._snap_coords_to_grid(raster_data[0], *coords) for coords in - self.dataframe.iloc[0].geometry.coords] + snapped_points = [snap_coords_to_grid(raster_data[0], *coords) for coords in line_coords] # Generate pairs of consecutive (x,y) coords path_pairs = list(map(list, zip(snapped_points, snapped_points[1:]))) + headings = [] + for i in range(1, len(line_coords)): + prev = line_coords[i - 1] + next = line_coords[i] + x = np.sin(next[0] - prev[0]) * np.cos(next[1]) + y = np.cos(prev[1]) * np.sin(next[1]) - np.sin(prev[1]) * np.cos(next[1]) * np.cos(next[0] - prev[0]) + angle = np.arctan2(x, y) % 2 * np.pi + headings.append(angle) # Feed these pairs into the Bresenham algo to find the intermediate points path_grid_points = [bresenham.make_line(*pair[0], *pair[1]) for pair in path_pairs] # Bring all these points together and remove duplicate coords From a94407b34bdc5107fe3646ca50f36c3eddced2f1 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:56:03 +0000 Subject: [PATCH 43/57] Fix ballistic model typo --- seedpod_ground_risk/path_analysis/ballistic_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/seedpod_ground_risk/path_analysis/ballistic_model.py b/seedpod_ground_risk/path_analysis/ballistic_model.py index 64e5918d..f71b0291 100644 --- a/seedpod_ground_risk/path_analysis/ballistic_model.py +++ b/seedpod_ground_risk/path_analysis/ballistic_model.py @@ -73,7 +73,7 @@ def impact_distance_dist_params_ned_with_wind(self, altitude, velocity, heading, d_i, _, _, t_i = self.bm.compute_ballistic_distance(altitude, velocity, 0) # Compensate for x,y axes being rotated compared to bearings - theta = (heading - (np.pi / 2)) * (2 * np.pi) + theta = (heading - (np.pi / 2)) % (2 * np.pi) # Get angle distribution in between body and NED frame # Form the array structure required and transform arr = np.vstack((np.zeros(d_i.shape), d_i, t_i, theta, wind_vel_x, wind_vel_y)) From 60acdefba8b95fd52b7e72a21116599670a89f1e Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 18:59:06 +0000 Subject: [PATCH 44/57] Fix arg order in ballistic model tests --- tests/path_analysis/test_ballistic_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/path_analysis/test_ballistic_model.py b/tests/path_analysis/test_ballistic_model.py index 736fe57b..e7d00b66 100644 --- a/tests/path_analysis/test_ballistic_model.py +++ b/tests/path_analysis/test_ballistic_model.py @@ -125,7 +125,7 @@ def test_ballistic_dist(self): wind_vel_y = ss.norm(1, 1).rvs(samples) bm = BallisticModel(self.ac) - means, cov = bm.impact_distance_dist_params_ned_with_wind(alt, vel, heading, wind_vel_x, wind_vel_y, loc_x, + means, cov = bm.impact_distance_dist_params_ned_with_wind(alt, vel, heading, wind_vel_y, wind_vel_x, loc_x, loc_y) pdf = ss.multivariate_normal(mean=means, cov=cov).pdf From 542307a35e41b5872049b19bfaafad6a48ddf598 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 19:00:13 +0000 Subject: [PATCH 45/57] Add aircraft and ballistic model to path analysis layer --- .../layers/path_analysis_layer.py | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 2a58ff3b..39c30cab 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -69,20 +69,16 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np # Flip left to right as bresenham spits out in (y,x) order path_grid_points = np.unique(np.concatenate(path_grid_points, axis=0), axis=0) - ca_model = casex.CriticalAreaModels() # Lethal area model using default params for human dimensions - aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, 2, 2, 2) # Default aircraft + ac_width = 2 + ac_length = 2 + ac_mass = 2 + aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, ac_width, ac_length, + ac_mass) # Default aircraft aircraft.set_ballistic_drag_coefficient(0.8) aircraft.set_ballistic_frontal_area(3) - fc = casex.FrictionCoefficients() - aircraft.set_friction_coefficient( - fc.get_coefficient(casex.enums.AircraftMaterial.ALUMINUM, casex.enums.GroundMaterial.CONCRETE)) - ballistic_model = casex.BallisticDescent2ndOrderDragApproximation() - ballistic_model.set_aircraft(aircraft) - distance_impact, velocity_impact, angle_impact, time_impact = ballistic_model.compute_ballistic_distance( - 133, 20, 1) - ca_model.critical_area(casex.enums.CriticalAreaModel.JARUS, - aircraft, velocity_impact, - angle_impact, 0) + aircraft.set_glide_speed_ratio(15, 12) + aircraft.set_glide_drag_coefficient(0.3) + bm = BallisticModel(aircraft) # Testing fixed pdf for now dist_mean = np.array([-5, 5]) From 9557fa5e6fa352544d6b2862fea7a203b018f3c7 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 19:01:22 +0000 Subject: [PATCH 46/57] Add some sample distributions --- seedpod_ground_risk/layers/path_analysis_layer.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 39c30cab..95464933 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -6,6 +6,8 @@ from holoviews.element import Overlay from seedpod_ground_risk.layers.annotation_layer import AnnotationLayer +from seedpod_ground_risk.path_analysis.ballistic_model import BallisticModel +from seedpod_ground_risk.path_analysis.utils import snap_coords_to_grid from seedpod_ground_risk.pathfinding import bresenham @@ -45,7 +47,7 @@ def preload_data(self) -> NoReturn: def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np.array], np.array], resolution=20, **kwargs) -> Overlay: import geoviews as gv - from scipy.stats import multivariate_normal as mvn + import scipy.stats as ss bounds = (raster_data[0]['Longitude'].min(), raster_data[0]['Latitude'].min(), raster_data[0]['Longitude'].max(), raster_data[0]['Latitude'].max()) @@ -80,8 +82,12 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np aircraft.set_glide_drag_coefficient(0.3) bm = BallisticModel(aircraft) - # Testing fixed pdf for now - dist_mean = np.array([-5, 5]) + samples = 3000 + # Conjure up our distributions for various things + alt = ss.norm(50, 5).rvs(samples) + vel = ss.norm(18, 2.5).rvs(samples) + wind_vel_y = ss.norm(5, 1).rvs(samples) + wind_vel_x = ss.norm(5, 1).rvs(samples) # Create grid on which to evaluate each point of path with its pdf raster_shape = raster_data[1].shape From a62223facbf8cf2b989a489b6573f0ea64753326 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 19:03:56 +0000 Subject: [PATCH 47/57] Generate and cache PDFs for each heading --- .../layers/path_analysis_layer.py | 44 +++++-------------- 1 file changed, 11 insertions(+), 33 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 95464933..6e273d7f 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -40,7 +40,6 @@ def __init__(self, key: str, filepath: str = '', buffer_dist: float = None, **kw self.buffer_poly = None def preload_data(self) -> NoReturn: - self.dataframe = gpd.read_file(self.filepath) self.endpoint = self.dataframe.iloc[0].geometry.coords[-1] @@ -67,6 +66,10 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np headings.append(angle) # Feed these pairs into the Bresenham algo to find the intermediate points path_grid_points = [bresenham.make_line(*pair[0], *pair[1]) for pair in path_pairs] + for idx, segment in enumerate(path_grid_points): + n = len(segment) + point_headings = np.full(n, headings[idx]) + path_grid_points[idx] = np.column_stack((np.array(segment), point_headings)) # Bring all these points together and remove duplicate coords # Flip left to right as bresenham spits out in (y,x) order path_grid_points = np.unique(np.concatenate(path_grid_points, axis=0), axis=0) @@ -94,9 +97,13 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np x, y = np.mgrid[0:raster_shape[0], 0:raster_shape[1]] eval_grid = np.dstack((x, y)) - from pathos.multiprocessing import Pool - from pathos.multiprocessing import cpu_count - pool = Pool() + dists_for_hdg = {} + for hdg in headings: + mean, cov = bm.impact_distance_dist_params_ned_with_wind(alt, vel, + ss.norm(hdg, np.deg2rad(2)).rvs(samples), + wind_vel_y, wind_vel_x, + 0, 0) + dists_for_hdg[hdg] = (mean / resolution, cov / resolution) # TODO Use something like Dask or OpenCV to speed this up in future as it's a simple map-reduce def point_dist(c): @@ -170,32 +177,3 @@ def map_chunk(coords): def clear_cache(self) -> NoReturn: pass - - def _snap_coords_to_grid(self, grid, lon: float, lat: float) -> Tuple[int, int]: - """ - Snap coordinates to grid indices - :param grid: raster grid coordinates - :param lon: longitude to snap - :param lat: latitude to snap - :return: (x, y) tuple of grid indices - """ - lat_idx = int(np.argmin(np.abs(grid['Latitude'] - lat))) - lon_idx = int(np.argmin(np.abs(grid['Longitude'] - lon))) - - return lon_idx, lat_idx - - def get_lethal_area(self, theta: float): - """ - Calculate lethal area of UAS impact from impact angle - - Method from Smith, P.G. 2000 - - :param theta: impact angle in degrees - :return: - """ - r_person = 1 # radius of a person - h_person = 1.8 # height of a person - - r_uas = 0.5 # UAS radius/halfspan - - return ((2 * (r_person + r_uas) * h_person) / np.tan(np.deg2rad(theta))) + (np.pi * (r_uas + r_person) ** 2) From b4176173fb21bb8c7948b8ae6f16e3efb31109a4 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 19:06:39 +0000 Subject: [PATCH 48/57] Flatten eval grid and sum shifted pdfs --- seedpod_ground_risk/layers/path_analysis_layer.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 6e273d7f..2c12acc9 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -95,7 +95,7 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np # Create grid on which to evaluate each point of path with its pdf raster_shape = raster_data[1].shape x, y = np.mgrid[0:raster_shape[0], 0:raster_shape[1]] - eval_grid = np.dstack((x, y)) + eval_grid = np.vstack((x.ravel(), y.ravel())).T dists_for_hdg = {} for hdg in headings: @@ -106,12 +106,12 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np dists_for_hdg[hdg] = (mean / resolution, cov / resolution) # TODO Use something like Dask or OpenCV to speed this up in future as it's a simple map-reduce - def point_dist(c): - pdf_mean = c + dist_mean - return mvn(pdf_mean, [[15, 5], [5, 7]]).pdf(eval_grid) + def point_distr(c): + dist_params = dists_for_hdg[c[2]] + pdf = ss.multivariate_normal(dist_params[0] + np.array([c[0], c[1]]), dist_params[1]).pdf(eval_grid) + return pdf - def map_chunk(coords): - return np.sum([point_dist(c) for c in coords], axis=0) + pdf_mat = np.sum([point_distr(c) for c in path_grid_points], axis=0).reshape(raster_shape) # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? pdfs = pool.map(map_chunk, np.array_split(path_grid_points, cpu_count() * 2)) From 4fc5d8bea71dd6cf898ab0fae029b702621dfc88 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 19:07:15 +0000 Subject: [PATCH 49/57] Use UAS width to calculate lethal area --- seedpod_ground_risk/layers/path_analysis_layer.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index 2c12acc9..a0065caf 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -113,12 +113,7 @@ def point_distr(c): pdf_mat = np.sum([point_distr(c) for c in path_grid_points], axis=0).reshape(raster_shape) - # TODO Remove all these flip and rotates; the indices must be swapping somewhere else? - pdfs = pool.map(map_chunk, np.array_split(path_grid_points, cpu_count() * 2)) - pool.close() - pdf_mat = np.sum(pdfs, axis=0) - - a_exp = self.get_lethal_area(30) + a_exp = get_lethal_area(30, aircraft.width) # Probability * Pop. Density * Lethal Area risk_map = pdf_mat * raster_data[1] * a_exp From b9c7a2da3fba0fe8d40bb717f9a00c7cd5a2baa4 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 19:09:49 +0000 Subject: [PATCH 50/57] Multiply event occurrence probability --- seedpod_ground_risk/layers/path_analysis_layer.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index a0065caf..d7757b6a 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -77,6 +77,7 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np ac_width = 2 ac_length = 2 ac_mass = 2 + event_probability = 0.05 aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, ac_width, ac_length, ac_mass) # Default aircraft aircraft.set_ballistic_drag_coefficient(0.8) @@ -115,7 +116,7 @@ def point_distr(c): a_exp = get_lethal_area(30, aircraft.width) # Probability * Pop. Density * Lethal Area - risk_map = pdf_mat * raster_data[1] * a_exp + risk_map = pdf_mat * raster_data[1] * a_exp * event_probability risk_raster = gv.Image(risk_map, bounds=bounds).options(alpha=0.7, cmap='viridis', tools=['hover'], clipping_colors={'min': (0, 0, 0, 0)}) From b526d60a61dcf3000647a408766b6958e2d9889d Mon Sep 17 00:00:00 2001 From: Zach10a <35232913+Zach10a@users.noreply.github.com> Date: Sun, 14 Mar 2021 21:32:56 +0000 Subject: [PATCH 51/57] Feature window resizing (#19) Fix window resize in ui --- seedpod_ground_risk/core/plot_worker.py | 4 + seedpod_ground_risk/main.py | 11 +- seedpod_ground_risk/ui_resources/main.ui | 137 ++++++++---------- .../ui_resources/mainwindow.py | 61 ++++---- .../ui_resources/plot_webview.py | 15 ++ 5 files changed, 122 insertions(+), 106 deletions(-) create mode 100644 seedpod_ground_risk/ui_resources/plot_webview.py diff --git a/seedpod_ground_risk/core/plot_worker.py b/seedpod_ground_risk/core/plot_worker.py index 05a9a216..20ac3ebe 100644 --- a/seedpod_ground_risk/core/plot_worker.py +++ b/seedpod_ground_risk/core/plot_worker.py @@ -79,6 +79,10 @@ def layers_reorder(self, layer_order): def export_path_json(self, layer, filepath): self.plot_server.export_path_geojson(layer, filepath) + @Slot(int, int) + def resize_plot(self, width, height): + self.plot_server.plot_size = (width, height) + def layers_update(self, layers): self.signals.update_layers.emit(layers) diff --git a/seedpod_ground_risk/main.py b/seedpod_ground_risk/main.py index 8ff1e5b8..98fece2b 100644 --- a/seedpod_ground_risk/main.py +++ b/seedpod_ground_risk/main.py @@ -138,6 +138,8 @@ def __init__(self): self.addLayerButton.clicked.connect(self.layer_add) + self.plotWebview.resize.connect(self.resize_plot) + self.actionRasterise.triggered.connect(self.menu_config_rasterise) # self.actionImport.triggered.connect(self.menu_file_import) self.actionExport.triggered.connect(self.menu_file_export) @@ -187,8 +189,8 @@ def menu_about_app(self): @Slot(str) def plot_ready(self, url): - self.webview.load(url) - self.webview.show() + self.plotWebview.load(url) + self.plotWebview.show() @Slot(str) def status_update(self, update_str: str): @@ -217,6 +219,9 @@ def layer_add(self): layer = layerObj(wizard.layerKey, **wizard.opts) self.plot_worker.add_layer(layer) + def resize_plot(self, width, height): + self.plot_worker.resize_plot(width, height) + def layer_edit(self, item): print('Editing ', item) pass @@ -230,10 +235,10 @@ def layer_reorder(self): [self.listWidget.item(n).text() for n in range(self.listWidget.count())]) def time_changed(self, value): + from seedpod_ground_risk.layers.roads_layer import generate_week_timesteps try: labels = generate_week_timesteps() except NameError: - from seedpod_ground_risk.layers.roads_layer import generate_week_timesteps labels = generate_week_timesteps() self.plot_worker.signals.set_time.emit(value) self.timeSliderLabel.setText(labels[value]) diff --git a/seedpod_ground_risk/ui_resources/main.ui b/seedpod_ground_risk/ui_resources/main.ui index d3344b55..475847f4 100644 --- a/seedpod_ground_risk/ui_resources/main.ui +++ b/seedpod_ground_risk/ui_resources/main.ui @@ -6,8 +6,8 @@ 0 0 - 1216 - 818 + 1200 + 836 @@ -19,7 +19,7 @@ 1200 - 800 + 836 @@ -40,41 +40,61 @@ - 1216 + 1200 757 - - - - 0 - 0 - 1221 - 771 - - - - - 0 - 0 - - - - - 1221 - 771 - - - - Qt::Horizontal - - - + + + + + QLayout::SetMinimumSize + - - - false + + + 0 + + + + false + + + + 0 + 0 + + + + + 405 + 490 + + + + + + + + 5 + + + 5 + + + + + Add Layer + + + + + + + + + 0 @@ -83,46 +103,15 @@ - 405 - 490 + 800 + 781 - - - - 5 - - - 5 - - - - - Add Layer - - - - - - - - - - 0 - 0 - - - - - 800 - 781 - - - - + + @@ -187,17 +176,17 @@ - - QWebEngineView - QWidget -
qwebengineview.h
- 1 - MapLayersListWidget QListWidget
maplayerslistwidget.h
 + + PlotWebview + QWidget +
plot_webview.h
+ 1 + diff --git a/seedpod_ground_risk/ui_resources/mainwindow.py b/seedpod_ground_risk/ui_resources/mainwindow.py index 8c098492..c8c53758 100644 --- a/seedpod_ground_risk/ui_resources/mainwindow.py +++ b/seedpod_ground_risk/ui_resources/mainwindow.py @@ -10,23 +10,23 @@ from PySide2.QtCore import * from PySide2.QtGui import * -from PySide2.QtWebEngineWidgets import QWebEngineView from PySide2.QtWidgets import * from .maplayerslistwidget import MapLayersListWidget +from .plot_webview import PlotWebview class Ui_MainWindow(object): def setupUi(self, MainWindow): if not MainWindow.objectName(): MainWindow.setObjectName(u"MainWindow") - MainWindow.resize(1216, 818) + MainWindow.resize(1200, 836) sizePolicy = QSizePolicy(QSizePolicy.MinimumExpanding, QSizePolicy.MinimumExpanding) sizePolicy.setHorizontalStretch(0) sizePolicy.setVerticalStretch(0) sizePolicy.setHeightForWidth(MainWindow.sizePolicy().hasHeightForWidth()) MainWindow.setSizePolicy(sizePolicy) - MainWindow.setMinimumSize(QSize(1200, 800)) + MainWindow.setMinimumSize(QSize(1200, 836)) MainWindow.setBaseSize(QSize(1200, 800)) self.actionImport = QAction(MainWindow) self.actionImport.setObjectName(u"actionImport") @@ -45,24 +45,23 @@ def setupUi(self, MainWindow): self.centralwidget.setObjectName(u"centralwidget") sizePolicy.setHeightForWidth(self.centralwidget.sizePolicy().hasHeightForWidth()) self.centralwidget.setSizePolicy(sizePolicy) - self.centralwidget.setMinimumSize(QSize(1216, 757)) - self.splitter = QSplitter(self.centralwidget) - self.splitter.setObjectName(u"splitter") - self.splitter.setGeometry(QRect(0, 0, 1221, 771)) - sizePolicy.setHeightForWidth(self.splitter.sizePolicy().hasHeightForWidth()) - self.splitter.setSizePolicy(sizePolicy) - self.splitter.setMinimumSize(QSize(1221, 771)) - self.splitter.setOrientation(Qt.Horizontal) - self.verticalLayoutWidget = QWidget(self.splitter) - self.verticalLayoutWidget.setObjectName(u"verticalLayoutWidget") - self.verticalLayout = QVBoxLayout(self.verticalLayoutWidget) + self.centralwidget.setMinimumSize(QSize(1200, 757)) + self.gridLayout = QGridLayout(self.centralwidget) + self.gridLayout.setObjectName(u"gridLayout") + self.horizontalLayout_4 = QHBoxLayout() + self.horizontalLayout_4.setObjectName(u"horizontalLayout_4") + self.horizontalLayout_4.setSizeConstraint(QLayout.SetMinimumSize) + self.verticalLayout = QVBoxLayout() self.verticalLayout.setObjectName(u"verticalLayout") - self.verticalLayout.setContentsMargins(0, 0, 0, 0) - self.listWidget = MapLayersListWidget(self.verticalLayoutWidget) + self.verticalLayout.setContentsMargins(-1, -1, 0, -1) + self.listWidget = MapLayersListWidget(self.centralwidget) self.listWidget.setObjectName(u"listWidget") self.listWidget.setEnabled(False) - sizePolicy.setHeightForWidth(self.listWidget.sizePolicy().hasHeightForWidth()) - self.listWidget.setSizePolicy(sizePolicy) + sizePolicy1 = QSizePolicy(QSizePolicy.Fixed, QSizePolicy.MinimumExpanding) + sizePolicy1.setHorizontalStretch(0) + sizePolicy1.setVerticalStretch(0) + sizePolicy1.setHeightForWidth(self.listWidget.sizePolicy().hasHeightForWidth()) + self.listWidget.setSizePolicy(sizePolicy1) self.listWidget.setMinimumSize(QSize(405, 490)) self.verticalLayout.addWidget(self.listWidget) @@ -70,21 +69,25 @@ def setupUi(self, MainWindow): self.horizontalLayout = QHBoxLayout() self.horizontalLayout.setObjectName(u"horizontalLayout") self.horizontalLayout.setContentsMargins(-1, 5, -1, 5) - self.addLayerButton = QPushButton(self.verticalLayoutWidget) + self.addLayerButton = QPushButton(self.centralwidget) self.addLayerButton.setObjectName(u"addLayerButton") self.horizontalLayout.addWidget(self.addLayerButton) - self.verticalLayout.addLayout(self.horizontalLayout) - self.splitter.addWidget(self.verticalLayoutWidget) - self.webview = QWebEngineView(self.splitter) - self.webview.setObjectName(u"webview") - sizePolicy.setHeightForWidth(self.webview.sizePolicy().hasHeightForWidth()) - self.webview.setSizePolicy(sizePolicy) - self.webview.setMinimumSize(QSize(800, 781)) - self.splitter.addWidget(self.webview) + self.horizontalLayout_4.addLayout(self.verticalLayout) + + self.plotWebview = PlotWebview(self.centralwidget) + self.plotWebview.setObjectName(u"plotWebview") + sizePolicy.setHeightForWidth(self.plotWebview.sizePolicy().hasHeightForWidth()) + self.plotWebview.setSizePolicy(sizePolicy) + self.plotWebview.setMinimumSize(QSize(800, 781)) + + self.horizontalLayout_4.addWidget(self.plotWebview) + + self.gridLayout.addLayout(self.horizontalLayout_4, 0, 0, 1, 1) + MainWindow.setCentralWidget(self.centralwidget) self.statusBar = QStatusBar(MainWindow) self.statusBar.setObjectName(u"statusBar") @@ -116,9 +119,9 @@ def retranslateUi(self, MainWindow): self.actionAbout_App.setText(QCoreApplication.translate("MainWindow", u"About App", None)) self.actionRasterise.setText(QCoreApplication.translate("MainWindow", u"Rasterise", None)) self.actionGenerate.setText(QCoreApplication.translate("MainWindow", u"Generate", None)) - # if QT_CONFIG(tooltip) +#if QT_CONFIG(tooltip) self.actionGenerate.setToolTip(QCoreApplication.translate("MainWindow", u"Generate Map for current view", None)) - #endif // QT_CONFIG(tooltip) +#endif // QT_CONFIG(tooltip) self.addLayerButton.setText(QCoreApplication.translate("MainWindow", u"Add Layer", None)) self.toolBar.setWindowTitle(QCoreApplication.translate("MainWindow", u"toolBar", None)) # retranslateUi diff --git a/seedpod_ground_risk/ui_resources/plot_webview.py b/seedpod_ground_risk/ui_resources/plot_webview.py new file mode 100644 index 00000000..c246ec60 --- /dev/null +++ b/seedpod_ground_risk/ui_resources/plot_webview.py @@ -0,0 +1,15 @@ +import PySide2 +from PySide2.QtCore import Signal +from PySide2.QtWebEngineWidgets import QWebEngineView + + +class PlotWebview(QWebEngineView): + resize = Signal(int, int) + + def __init__(self, *args, **kwargs): + super(PlotWebview, self).__init__(*args, **kwargs) + + def resizeEvent(self, event: PySide2.QtGui.QResizeEvent) -> None: + super().resizeEvent(event) + webview_size = self.size() + self.resize.emit(webview_size.width() - 50, webview_size.height() - 30) From d0bd85c5e84a484da0700f079ae9174dee7b9599 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 23:11:36 +0000 Subject: [PATCH 52/57] Add ballistic model params to layer args --- seedpod_ground_risk/core/plot_server.py | 2 +- .../layers/path_analysis_layer.py | 46 +++++++++++-------- .../layers/pathfinding_layer.py | 7 ++- 3 files changed, 31 insertions(+), 24 deletions(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 6576a0c8..53eb8fce 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -74,7 +74,7 @@ def __init__(self, tiles: str = 'Wikipedia', tools: Optional[Iterable[str]] = No self._generated_data_layers = {} self.data_layer_order = [] self.data_layers = [ResidentialLayer('Residential Population', buffer_dist=30), - RoadsLayer('Road Traffic Population per Hour')] + RoadsLayer('Road Traffic Population/Hour')] self.annotation_layers = [] diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index d7757b6a..e65741e1 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -30,10 +30,26 @@ def get_lethal_area(theta: float, uas_width: float): class PathAnalysisLayer(AnnotationLayer): - def __init__(self, key: str, filepath: str = '', buffer_dist: float = None, **kwargs): + def __init__(self, key: str, filepath: str = '', ac_width: float = 2, ac_length: float = 2, + ac_mass: float = 2, ac_glide_ratio: float = 12, ac_glide_speed: float = 15, + ac_glide_drag_coeff: float = 0.8, ac_ballistic_drag_coeff: float = 0.8, + ac_ballistic_frontal_area: float = 3, ac_failure_prob: float = 0.05, alt: float = 50, vel: float = 18, + wind_vel: float = 5, wind_dir: float = 45, **kwargs): super(PathAnalysisLayer, self).__init__(key) self.filepath = filepath - self.buffer_dist = buffer_dist + + self.aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, ac_width, ac_length, ac_mass) + self.aircraft.set_ballistic_drag_coefficient(ac_ballistic_drag_coeff) + self.aircraft.set_ballistic_frontal_area(ac_ballistic_frontal_area) + self.aircraft.set_glide_speed_ratio(ac_glide_speed, ac_glide_ratio) + self.aircraft.set_glide_drag_coefficient(ac_glide_drag_coeff) + + self.alt = alt + self.vel = vel + self.wind_vel = wind_vel + self.wind_dir = np.deg2rad(wind_dir) + + self.event_prob = ac_failure_prob import geopandas as gpd self.dataframe = gpd.GeoDataFrame() @@ -74,24 +90,16 @@ def annotate(self, data: List[gpd.GeoDataFrame], raster_data: Tuple[Dict[str, np # Flip left to right as bresenham spits out in (y,x) order path_grid_points = np.unique(np.concatenate(path_grid_points, axis=0), axis=0) - ac_width = 2 - ac_length = 2 - ac_mass = 2 - event_probability = 0.05 - aircraft = casex.AircraftSpecs(casex.enums.AircraftType.FIXED_WING, ac_width, ac_length, - ac_mass) # Default aircraft - aircraft.set_ballistic_drag_coefficient(0.8) - aircraft.set_ballistic_frontal_area(3) - aircraft.set_glide_speed_ratio(15, 12) - aircraft.set_glide_drag_coefficient(0.3) - bm = BallisticModel(aircraft) + bm = BallisticModel(self.aircraft) samples = 3000 # Conjure up our distributions for various things - alt = ss.norm(50, 5).rvs(samples) - vel = ss.norm(18, 2.5).rvs(samples) - wind_vel_y = ss.norm(5, 1).rvs(samples) - wind_vel_x = ss.norm(5, 1).rvs(samples) + alt = ss.norm(self.alt, 5).rvs(samples) + vel = ss.norm(self.vel, 2.5).rvs(samples) + wind_vels = ss.norm(self.wind_vel, 1).rvs(samples) + wind_dirs = ss.norm(self.wind_dir, 5).rvs(samples) + wind_vel_y = wind_vels * np.sin(wind_dirs) + wind_vel_x = wind_vels * np.cos(wind_dirs) # Create grid on which to evaluate each point of path with its pdf raster_shape = raster_data[1].shape @@ -114,9 +122,9 @@ def point_distr(c): pdf_mat = np.sum([point_distr(c) for c in path_grid_points], axis=0).reshape(raster_shape) - a_exp = get_lethal_area(30, aircraft.width) + a_exp = get_lethal_area(30, self.aircraft.width) # Probability * Pop. Density * Lethal Area - risk_map = pdf_mat * raster_data[1] * a_exp * event_probability + risk_map = pdf_mat * raster_data[1] * a_exp * self.event_prob risk_raster = gv.Image(risk_map, bounds=bounds).options(alpha=0.7, cmap='viridis', tools=['hover'], clipping_colors={'min': (0, 0, 0, 0)}) diff --git a/seedpod_ground_risk/layers/pathfinding_layer.py b/seedpod_ground_risk/layers/pathfinding_layer.py index 215f8476..6fa775aa 100644 --- a/seedpod_ground_risk/layers/pathfinding_layer.py +++ b/seedpod_ground_risk/layers/pathfinding_layer.py @@ -16,12 +16,11 @@ class PathfindingLayer(PathAnalysisLayer): def __init__(self, key, start_lat: float = 0, start_lon: float = 0, end_lat: float = 0, end_lon: float = 0, - buffer: float = 0, algo: Algorithm = RiskGridAStar, heuristic: Heuristic = ManhattanRiskHeuristic, - rdr: float = 0.5): - super().__init__(key, '') + algo: Algorithm = RiskGridAStar, heuristic: Heuristic = ManhattanRiskHeuristic, + rdr: float = 0.5, **kwargs): + super().__init__(key, '', **kwargs) self.start_coords = (start_lat, start_lon) self.end_coords = (end_lat, end_lon) - self.buffer_dist = buffer self.algo = algo self.heuristic = heuristic self.rdr = rdr From 8e5296eb70dc96b62f1aba5e61ff7307522cda25 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Sun, 14 Mar 2021 23:17:28 +0000 Subject: [PATCH 53/57] Add ballistic model params to wizard --- .../ui_resources/layer_options.py | 30 +++++++++++++++++-- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/seedpod_ground_risk/ui_resources/layer_options.py b/seedpod_ground_risk/ui_resources/layer_options.py index 2c581677..197ea884 100644 --- a/seedpod_ground_risk/ui_resources/layer_options.py +++ b/seedpod_ground_risk/ui_resources/layer_options.py @@ -41,7 +41,19 @@ }, 'Existing Path Analysis': { 'File': ('path', 'filepath', str), - 'Buffer Distance [m]': (r'\d{0,3}', 'buffer_dist', int), + 'Aircraft Width [m]': (r'-?\d{0,3}\.?\d+', 'ac_width', float), + 'Aircraft Length [m]': (r'-?\d{0,3}\.?\d+', 'ac_length', float), + 'Aircraft Mass [kg]': (r'-?\d{0,3}\.?\d+', 'ac_mass', float), + 'Aircraft Glide Ratio': (r'-?\d{0,3}\.?\d+', 'ac_glide_ratio', float), + 'Aircraft Glide Speed [m/s]': (r'-?\d{0,3}\.?\d+', 'ac_glide_speed', float), + 'Aircraft Glide Drag Coeff': (r'-?\d{0,3}\.?\d+', 'ac_glide_drag_coeff', float), + 'Aircraft Ballistic Drag Coeff': (r'-?\d{0,3}\.?\d+', 'ac_ballistic_drag_coeff', float), + 'Aircraft Ballistic Frontal Area [m^2]': (r'-?\d{0,3}\.?\d+', 'ac_ballistic_frontal_area', float), + 'Aircraft Failure Probability [0-1]': (r'-?\d{0,3}\.?\d+', 'ac_failure_prob', float), + 'Flight Altitude [m]': (r'-?\d{0,3}\.?\d+', 'alt', float), + 'Flight Airspeed [m/s]': (r'-?\d{0,3}\.?\d+', 'vel', float), + 'Wind Speed [m/s]': (r'-?\d{0,3}\.?\d+', 'wind_vel', float), + 'Wind Bearing [deg]': (r'-?\d{0,3}\.?\d+', 'wind_dir', float) }, 'External Obstacles': { 'File': ('path', 'filepath', str), @@ -50,12 +62,24 @@ 'Blocking': (bool, 'blocking', bool), }, 'Pathfinding': { - 'Buffer Distance [m]': (r'\d{0,3}', 'buffer', int), 'Start Latitude [dd]': (r'-?\d{0,3}\.\d+', 'start_lat', float), 'Start Longitude [dd]': (r'-?\d{0,3}\.\d+', 'start_lon', float), 'End Latitude [dd]': (r'-?\d{0,3}\.\d+', 'end_lat', float), 'End Longitude [dd]': (r'-?\d{0,3}\.\d+', 'end_lon', float), 'Algorithm': ('algos', 'algo', eval), - 'Risk-Distance Ratio': (r'\d{0,3}(\.\d+)?', 'rdr', float) + 'Risk-Distance Ratio': (r'\d{0,3}(\.\d+)?', 'rdr', float), + 'Aircraft Width [m]': (r'-?\d{0,3}\.?\d+', 'ac_width', float), + 'Aircraft Length [m]': (r'-?\d{0,3}\.?\d+', 'ac_length', float), + 'Aircraft Mass [kg]': (r'-?\d{0,3}\.?\d+', 'ac_mass', float), + 'Aircraft Glide Ratio': (r'-?\d{0,3}\.?\d+', 'ac_glide_ratio', float), + 'Aircraft Glide Speed [m/s]': (r'-?\d{0,3}\.?\d+', 'ac_glide_speed', float), + 'Aircraft Glide Drag Coeff': (r'-?\d{0,3}\.?\d+', 'ac_glide_drag_coeff', float), + 'Aircraft Ballistic Drag Coeff': (r'-?\d{0,3}\.?\d+', 'ac_ballistic_drag_coeff', float), + 'Aircraft Ballistic Frontal Area [m^2]': (r'-?\d{0,3}\.?\d+', 'ac_ballistic_frontal_area', float), + 'Aircraft Failure Probability [0-1]': (r'-?\d{0,3}\.?\d+', 'ac_failure_prob', float), + 'Flight Altitude [m]': (r'-?\d{0,3}\.?\d+', 'alt', float), + 'Flight Airspeed [m/s]': (r'-?\d{0,3}\.?\d+', 'vel', float), + 'Wind Speed [m/s]': (r'-?\d{0,3}\.?\d+', 'wind_vel', float), + 'Wind Bearing [deg]': (r'-?\d{0,3}\.?\d+', 'wind_dir', float) } } From 5b3cf08945acdc70862ecf5f5207b966f3543cc5 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 15 Mar 2021 00:10:18 +0000 Subject: [PATCH 54/57] Add progress messages --- seedpod_ground_risk/core/plot_server.py | 1 + seedpod_ground_risk/layers/path_analysis_layer.py | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/seedpod_ground_risk/core/plot_server.py b/seedpod_ground_risk/core/plot_server.py index 53eb8fce..209b66fe 100644 --- a/seedpod_ground_risk/core/plot_server.py +++ b/seedpod_ground_risk/core/plot_server.py @@ -225,6 +225,7 @@ def compose_overlay_plot(self, x_range: Optional[Sequence[float]] = (-1.6, -1.2) raster_grid = np.flipud(raster_grid) raster_indices['Latitude'] = np.flip(raster_indices['Latitude']) annotations = [] + print('Annotating Layers...') for layer in self.annotation_layers: annotation = layer.annotate(raw_datas, (raster_indices, raster_grid), resolution=self.raster_resolution_m) diff --git a/seedpod_ground_risk/layers/path_analysis_layer.py b/seedpod_ground_risk/layers/path_analysis_layer.py index e65741e1..b1759573 100644 --- a/seedpod_ground_risk/layers/path_analysis_layer.py +++ b/seedpod_ground_risk/layers/path_analysis_layer.py @@ -128,7 +128,8 @@ def point_distr(c): risk_raster = gv.Image(risk_map, bounds=bounds).options(alpha=0.7, cmap='viridis', tools=['hover'], clipping_colors={'min': (0, 0, 0, 0)}) - risk_raster = risk_raster.redim.range(z=(1e-11, risk_map.max())) + risk_raster = risk_raster.redim.range(z=(1e-9, risk_map.max())) + print('Max probability of striking person: ', risk_map.max()) # labels = [] # annotation_layers = [] From 083af2caf22423c9e09a92c267f312f78ea2de69 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 15 Mar 2021 00:35:46 +0000 Subject: [PATCH 55/57] Fix paths --- seedpod_ground_risk/layers/residential_layer.py | 4 ++-- seedpod_ground_risk/layers/roads_layer.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/seedpod_ground_risk/layers/residential_layer.py b/seedpod_ground_risk/layers/residential_layer.py index 4c27f518..61fd23af 100644 --- a/seedpod_ground_risk/layers/residential_layer.py +++ b/seedpod_ground_risk/layers/residential_layer.py @@ -78,13 +78,13 @@ def ingest_census_data(self) -> NoReturn: # Import Census boundaries in Ordnance Survey grid and reproject census_wards_df = gpd.read_file( - os.sep.join(('..', 'static_data', 'england_wa_2011_clipped.shp'))).drop( + os.sep.join(('static_data', 'england_wa_2011_clipped.shp'))).drop( ['altname', 'oldcode'], axis=1) if not census_wards_df.crs: census_wards_df = census_wards_df.set_crs('EPSG:27700') census_wards_df = census_wards_df.to_crs('EPSG:4326') # Import census ward densities - density_df = pd.read_csv(os.sep.join(('..', 'static_data', 'density.csv')), header=0) + density_df = pd.read_csv(os.sep.join(('static_data', 'density.csv')), header=0) # Scale from hectares to m^2 density_df['area'] = density_df['area'] * 10000 density_df['density'] = density_df['density'] / 10000 diff --git a/seedpod_ground_risk/layers/roads_layer.py b/seedpod_ground_risk/layers/roads_layer.py index 0a3c3e43..8da6ef8b 100644 --- a/seedpod_ground_risk/layers/roads_layer.py +++ b/seedpod_ground_risk/layers/roads_layer.py @@ -136,7 +136,7 @@ def _ingest_traffic_counts(self) -> NoReturn: import os # Ingest raw data - counts_df = pd.read_csv(os.sep.join(('..', 'static_data', 'dft_traffic_counts_aadf.csv'))) + counts_df = pd.read_csv(os.sep.join(('static_data', 'dft_traffic_counts_aadf.csv'))) # Select only desired columns counts_df = counts_df[TRAFFIC_COUNT_COLUMNS] # Groupby year and select out only the latest year @@ -164,7 +164,7 @@ def _ingest_relative_traffic_variations(self): import os # Ingest data, ignoring header and footer info - relative_variations_df = pd.read_excel(os.sep.join(('..', 'static_data', 'tra0307.ods')), engine='odf', + relative_variations_df = pd.read_excel(os.sep.join(('static_data', 'tra0307.ods')), engine='odf', header=5, skipfooter=8) # Flatten into continuous list of hourly variations for the week self.relative_variations_flat = (relative_variations_df.iloc[:, 1:] / 100).melt()['value'] @@ -175,7 +175,7 @@ def _ingest_road_geometries(self) -> NoReturn: """ import os - self._roads_geometries = gpd.read_file(os.sep.join(('..', 'static_data', '2018-MRDB-minimal.shp'))).rename( + self._roads_geometries = gpd.read_file(os.sep.join(('static_data', '2018-MRDB-minimal.shp'))).rename( columns={'CP_Number': 'count_point_id'}) if not self._roads_geometries.crs: self._roads_geometries = self._roads_geometries.set_crs('EPSG:27700') From 033cd7390df7c826285bc8d91a081dc08c31fb9f Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 15 Mar 2021 00:51:32 +0000 Subject: [PATCH 56/57] Fix dependency packaging --- SEEDPOD Ground Risk.spec | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SEEDPOD Ground Risk.spec b/SEEDPOD Ground Risk.spec index 87a15914..05e53103 100644 --- a/SEEDPOD Ground Risk.spec +++ b/SEEDPOD Ground Risk.spec @@ -12,6 +12,7 @@ import panel import pyviz_comms import rtree import shiboken2 +import sklearn from PyInstaller.building.api import PYZ, EXE, COLLECT from PyInstaller.building.build_main import Analysis @@ -62,6 +63,7 @@ a = Analysis(['seedpod_ground_risk/main.py'], ("README.md", '.'), (os.path.join(pyviz_comms.comm_path, "notebook.js"), "pyviz_comms"), (panel.__path__[0], "panel"), + (sklearn.__path__[0], "sklearn"), (datashader.__path__[0], "datashader"), (distributed.__path__[0], "distributed"), (os.path.join(fiona.__path__[0], "*.pyd"), "fiona"), # Geospatial primitives @@ -101,7 +103,7 @@ a = Analysis(['seedpod_ground_risk/main.py'], "PyQt5", "PyQt4", "tkinter", - "pydoc", + # "pydoc", "pdb", "IPython", "jupyter", From b4ded6de2ce7e0547d2ee59239f11d7adb715d83 Mon Sep 17 00:00:00 2001 From: Aliaksei Pilko Date: Mon, 15 Mar 2021 09:32:39 +0000 Subject: [PATCH 57/57] Update installer scripts --- SEEDPOD Ground Risk.spec | 2 +- make_installer.iss | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SEEDPOD Ground Risk.spec b/SEEDPOD Ground Risk.spec index 05e53103..f16127f8 100644 --- a/SEEDPOD Ground Risk.spec +++ b/SEEDPOD Ground Risk.spec @@ -63,7 +63,7 @@ a = Analysis(['seedpod_ground_risk/main.py'], ("README.md", '.'), (os.path.join(pyviz_comms.comm_path, "notebook.js"), "pyviz_comms"), (panel.__path__[0], "panel"), - (sklearn.__path__[0], "sklearn"), + (sklearn.utils.__path__[0], "sklearn/utils"), (datashader.__path__[0], "datashader"), (distributed.__path__[0], "distributed"), (os.path.join(fiona.__path__[0], "*.pyd"), "fiona"), # Geospatial primitives diff --git a/make_installer.iss b/make_installer.iss index 7890fb2c..6a876dde 100644 --- a/make_installer.iss +++ b/make_installer.iss @@ -2,7 +2,7 @@ ; SEE THE DOCUMENTATION FOR DETAILS ON CREATING INNO SETUP SCRIPT FILES! #define MyAppName "SEEDPOD Ground Risk" -#define MyAppVersion "0.6a" +#define MyAppVersion "0.10.0" #define MyAppPublisher "CASCADE UAV" #define MyAppURL "https://cascadeuav.com/seedpod/" #define MyAppExeName "SEEDPOD Ground Risk.exe"