import itertools from shapely.geometry import ( JOIN_STYLE, GeometryCollection, MultiPolygon, Polygon, ) from shapely.ops import unary_union from werkzeug.utils import cached_property class Polygons(): approx_metres_to_degree = 111_320 approx_square_metres_to_square_degree = approx_metres_to_degree ** 2 # Estimated amount of bleed into neigbouring areas based on typical # range/separation of cell towers. approx_bleed_in_degrees = 1_500 / approx_metres_to_degree # Controls how much buffer to add for a shape of a given perimeter. # Smaller number means more buffering and a smoother shape. For # example `1000` means 1m of buffer for every 1km of perimeter, or # 20m of buffer for a 5km square. This gives us control over how # much we fill in very concave features like channels, harbours and # zawns. perimeter_to_buffer_ratio = 360 # Ratio of how much detail a shape of a given perimeter has once # simplified. Smaller number means less detail. For example `1000` # means that for a shape with a perimeter of 1000m, the simplified # line will never deviate more than 1m from the original. # Or for a 5km square, the line won’t deviate more than 20m. This # gives us approximate control over the total number of points. perimeter_to_simplification_ratio = 1_750 # The threshold for removing very small areas from the map. These # areas are likely glitches in the data where the shoreline hasn’t # been subtracted from the land properly minimum_area_size_square_metres = 50 ** 2 def __init__(self, polygons): if not polygons: self.polygons = [] elif isinstance(polygons[0], list): self.polygons = [ Polygon(polygon) for polygon in polygons ] else: self.polygons = polygons def __getitem__(self, index): return self.polygons[index] def __len__(self): return len(self.polygons) @cached_property def perimeter_length(self): return sum( polygon.length for polygon in self ) @cached_property def buffer_outward_in_degrees(self): return ( # If two areas are close enough that the distance between # them is less than the minimum bleed of a cell # broadcast then this joins them together. The aim is to # reduce the total number of polygons in areas with many # small shapes like Orkney or the Isles of Scilly. self.approx_bleed_in_degrees / 3 ) + ( self.perimeter_length / self.perimeter_to_buffer_ratio ) @cached_property def buffer_inward_in_degrees(self): return self.buffer_outward_in_degrees - ( # We should leave the shape expanded by at least the # simplification tolerance in all places, so the # simplification never moves a point inside the original # shape. In practice half of the tolerance is enough to # acheive this. self.simplification_tolerance_in_degrees / 2 ) @cached_property def simplification_tolerance_in_degrees(self): return self.perimeter_length / self.perimeter_to_simplification_ratio @cached_property def smooth(self): buffered = [ polygon.buffer( self.buffer_outward_in_degrees, resolution=4, join_style=JOIN_STYLE.round, ) for polygon in self ] unioned = union_polygons(buffered) debuffered = [ polygon.buffer( -1 * self.buffer_inward_in_degrees, resolution=1, join_style=JOIN_STYLE.bevel, ) for polygon in unioned ] flattened = list(itertools.chain(*[ flatten_polygons(polygon) for polygon in debuffered ])) return Polygons(flattened) @cached_property def simplify(self): return Polygons([ polygon.simplify(self.simplification_tolerance_in_degrees) for polygon in self ]) @cached_property def bleed(self): return Polygons(union_polygons([ polygon.buffer( self.approx_bleed_in_degrees, resolution=4, join_style=JOIN_STYLE.round, ) for polygon in self ])) @cached_property def remove_too_small(self): return Polygons([ polygon for polygon in self if ( polygon.area * self.approx_square_metres_to_square_degree ) > ( self.minimum_area_size_square_metres ) ]) @cached_property def as_coordinate_pairs_long_lat(self): return [ [[x, y] for x, y in polygon.exterior.coords] for polygon in self ] @cached_property def as_coordinate_pairs_lat_long(self): return [ [[y, x] for x, y in coordinate_pairs] for coordinate_pairs in self.as_coordinate_pairs_long_lat ] @cached_property def point_count(self): return len(list(itertools.chain(*self.as_coordinate_pairs_long_lat))) def flatten_polygons(polygons): if isinstance(polygons, GeometryCollection): return [] if isinstance(polygons, MultiPolygon): return [ p for p in polygons ] else: return [polygons] def union_polygons(polygons): return flatten_polygons(unary_union(polygons))