Coverage for src/sensai/geoanalytics/local_coords.py: 0%
104 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-08-13 22:17 +0000
1"""
2Local coordinate systems (for geographic data)
3"""
4import math
5from functools import reduce
6from typing import Tuple, Union, List
8import numpy as np
9import utm
10from shapely.geometry import polygon, multipolygon, point, LineString, mapping
11from shapely.ops import polygonize, unary_union
14class LocalCoordinateSystem(object):
15 """
16 Represents a local coordinate system for the conversion of geo-coordinates
17 (latitude, longitude) to a local Cartesian coordinate system (unit=metre) and vice versa
18 using the UTM transform
19 """
21 def __init__(self, lat, lon):
22 """
23 Parameters:
24 lat: the latitude of the origin of the coordinate system
25 lon: the longitude of the origin of the coordinate system
26 """
27 self.uRef = utm.from_latlon(lat, lon)
28 self.uRefE = self.uRef[0]
29 self.uRefN = self.uRef[1]
30 self.uRefPseudoN = self._pseudo_northing(self.uRefN)
32 def get_local_coords(self, lat, lon) -> Tuple[float, float]:
33 uE, uN, zM, zL = utm.from_latlon(lat, lon)
34 x = uE - self.uRefE
35 y = self._pseudo_northing(uN) - self.uRefPseudoN
36 return x, y
38 def get_lat_lon(self, local_x, local_y) -> Tuple[float, float]:
39 easting = local_x + self.uRefE
40 pseudo_northing = local_y + self.uRefPseudoN
41 return utm.to_latlon(easting, self._real_northing(pseudo_northing), self.uRef[2], self.uRef[3])
43 @staticmethod
44 def _pseudo_northing(real_northing):
45 if real_northing >= 10000000:
46 return real_northing - 10000000
47 else:
48 return real_northing
50 @staticmethod
51 def _real_northing(pseudo_northing):
52 if pseudo_northing < 0:
53 return pseudo_northing + 10000000
54 else:
55 return pseudo_northing
58class LocalHexagonalGrid:
59 """
60 A local hexagonal grid, where hex cells can be referenced by two integer coordinates relative to
61 the central grid cell, whose centre is at local coordinate (0, 0) and where positive x-coordinates/columns
62 are towards the east and positive y-coordinates/rows are towards the north.
63 Every odd row of cells is shifted half a hexagon to the right, i.e. column x for row 1 is half a grid cell
64 further to the right than column x for row 0.
66 For visualisation purposes, see https://www.redblobgames.com/grids/hexagons/
67 """
68 def __init__(self, radius_m):
69 """
70 :param radius_m: the radius, in metres, of each hex cell
71 """
72 self.radius_m = radius_m
73 start_angle = math.pi / 6
74 step_angle = math.pi / 3
75 self.offset_vectors = []
76 for i in range(6):
77 angle = start_angle + i * step_angle
78 x = math.cos(angle) * radius_m
79 y = math.sin(angle) * radius_m
80 self.offset_vectors.append(np.array([x, y]))
81 self.hexagon_width = 2 * self.offset_vectors[0][0]
82 self.hexagon_height = 2 * self.offset_vectors[1][1]
83 self.row_step = 0.75 * self.hexagon_height
84 self.polygon_area = 6 * self.hexagon_height * self.hexagon_width / 8
86 def get_hexagon(self, x_column: int, y_row: int) -> polygon.Polygon:
87 """
88 Gets the hexagon (polygon) for the given integer hex cell coordinates
89 :param x_column: the column coordinate
90 :param y_row: the row coordinate
91 :return: the hexagon
92 """
93 centre_x = x_column * self.hexagon_width
94 centre_y = y_row * self.row_step
95 if y_row % 2 == 1:
96 centre_x += 0.5 * self.hexagon_width
97 centre = np.array([centre_x, centre_y])
98 return polygon.Polygon([centre + o for o in self.offset_vectors])
100 def get_min_hexagon_column(self, x):
101 lowest_x_definitely_in_column0 = 0
102 return math.floor((x - lowest_x_definitely_in_column0) / self.hexagon_width)
104 def get_max_hexagon_column(self, x):
105 highest_x_definitely_in_column0 = self.hexagon_width / 2
106 return math.ceil((x - highest_x_definitely_in_column0) / self.hexagon_width)
108 def get_min_hexagon_row(self, y):
109 lowest_y_definitely_in_row0 = -self.hexagon_height / 4
110 return math.floor((y - lowest_y_definitely_in_row0) / self.row_step)
112 def get_max_hexagon_row(self, y):
113 highest_y_definitely_in_row0 = self.hexagon_height / 4
114 return math.ceil((y - highest_y_definitely_in_row0) / self.row_step)
116 def get_hexagon_coord_span_for_bounding_box(self, min_x, min_y, max_x, max_y) -> Tuple[Tuple[int, int], Tuple[int, int]]:
117 """
118 Gets the range of hex-cell coordinates that cover the given bounding box
120 :param min_x: minimum x-coordinate of bounding box
121 :param min_y: minimum y-coordinate of bounding box
122 :param max_x: maximum x-coordinate of bounding box
123 :param max_y: maximum y-coordinate of bounding box
124 :return: a pair of pairs ((minCol, min_row), (maxCol, max_row)) indicating the span of cell coordinates
125 """
126 if min_x > max_x or min_y > max_y:
127 raise ValueError()
128 min_column = self.get_min_hexagon_column(min_x)
129 max_column = self.get_max_hexagon_column(max_x)
130 min_row = self.get_min_hexagon_row(min_y)
131 max_row = self.get_max_hexagon_row(max_y)
132 return ((min_column, min_row), (max_column, max_row))
134 def get_hexagon_coords_for_point(self, x, y):
135 ((minColumn, minRow), (maxColumn, maxRow)) = self.get_hexagon_coord_span_for_bounding_box(x, y, x, y)
136 for xCol in range(minColumn, maxColumn+1):
137 for yRow in range(minRow, maxRow+1):
138 if self.get_hexagon(xCol, yRow).contains(point.Point(x, y)):
139 return xCol, yRow
140 raise Exception("No Hexagon matched; possible edge case (point on hexagon boundary)")
143def fix_polygon(poly: Union[polygon.Polygon, multipolygon.MultiPolygon], maxAreaDiff=1e-2) \
144 -> Union[polygon.Polygon, multipolygon.MultiPolygon]:
145 """
146 Fix invalid shapely polygons or multipolygons.
148 Reference:
149 https://stackoverflow.com/questions/35110632/splitting-self-intersecting-polygon-only-returned-one-polygon-in-shapely
151 :param poly: the polygon to fix
152 :param maxAreaDiff: the maximum change in area
153 :return: the fixed polygon or None if it cannot be fixed given the area change constraint
154 """
155 def _fix_polygon_component(coords: List[Tuple[float, float]]):
156 res = list(polygonize(unary_union(LineString(list(coords) + [coords[0]]))))
157 return reduce(lambda p1, p2: p1.union(p2), res)
159 if poly.is_valid:
160 return poly
161 else:
162 if isinstance(poly, polygon.Polygon):
163 exterior_coords = poly.exterior.coords[:]
164 fixed_exterior = _fix_polygon_component(exterior_coords)
165 fixed_interior = polygon.Polygon()
166 for interior in poly.interiors:
167 coords = interior.coords[:]
168 fixed_interior = fixed_interior.union(_fix_polygon_component(coords))
169 fixed_polygon = fixed_exterior.difference(fixed_interior)
170 elif isinstance(poly, multipolygon.MultiPolygon):
171 polys = list(poly)
172 fixed_polys = [fix_polygon(p, maxAreaDiff=maxAreaDiff) for p in polys]
173 fixed_polygon = reduce(lambda p1, p2: p1.union(p2), fixed_polys)
174 else:
175 raise Exception(f"Unsupported type {type(poly)}")
176 area_diff = float('Inf') if poly.area == 0 else abs(poly.area - fixed_polygon.area) / poly.area
177 #log.info(f"Invalid polygon\n{poly}\nComputed fix:\n{fixed_polygon}.\nArea error: {area_diff}")
178 if area_diff > maxAreaDiff:
179 return None
180 else:
181 return fixed_polygon