diff --git a/modules/inflate_geo.py b/modules/inflate_geo.py new file mode 100644 index 0000000..10c9bd4 --- /dev/null +++ b/modules/inflate_geo.py @@ -0,0 +1,193 @@ +"""Functions for manipulating polygon geometry""" + +import math +import numpy as np + + +def get_distance(vertices_org: list, vertices_scal: list) -> list: + """ + Given the cartesian vertices of two 2D polygon geometries (one original and one scaled), return a list where each element is the distance (in meters) + between a vertex on the original polygon and the same vertex on the scaled polygon. + Parameters + ----------- + vertices_org: list + List of vertices (latitude and longitude, in degrees) of the original polygon. + vertices_scal: list + List of vertices (latitude and longitude, in degrees) of the scaled polygon + + Returns + ------- + list: + Each element is the distance (in meters) between a vertex on the original polygon and the same vertex on the scaled polygon. + """ + + # Check for valid input + assert len(vertices_org) == len(vertices_scal) + assert len(vertices_org[0]) == 2 and len(vertices_scal[0]) == 2 + + distances = np.zeros((len(vertices_org), 1)) + for i, (x1, y1) in enumerate(vertices_org[:-1]): + (x2, y2) = vertices_scal[i + 1] + (dx, dy) = (x2 - x1, y2 - y1) + dist = math.sqrt(dx**2 + dy**2) + distances[i] = dist + + return distances + + +def is_valid_polygon(vertices: list) -> bool: + """ + Given the latlong vertices of two polygon geometries (one original and one scaled), return true if the there are not intersections between + the edges of the vertices. + + Parameters + ----------- + vertices: list + List of vertices (latitude and longitude, in degrees) of the polygon. + + Returns + ------- + bool: + true if the there are not intersections between the edges of the vertices. + """ + + # Check for valid input: + assert len(vertices[0]) == 2 + + # Generate a list of edges joining the vertices: + edges = np.zeros((len(vertices), 2)) + for i in range(len(edges) - 1): + (x1, y1) = vertices[i] + (x2, y2) = vertices[i + 1] + (dx, dy) = (x2 - x1, y2 - y1) + edges[i] = (dx, dy) + + # Generate the last edge (connecting last vertex in list to the first): + (x1, y1) = vertices[len(vertices) - 1] + (x2, y2) = vertices[0] + (dx, dy) = (x2 - x1, y2 - y1) + edges[len(edges) - 1] = (dx, dy) + + for i, (x1, y1) in enumerate(vertices): + for n, (x2, y2) in enumerate(vertices): + # Check that vertices are not the same. Skip iteration if they are. + if np.array_equal(vertices[i], vertices[n]): + continue + (dx1, dy1) = edges[i] + (dx2, dy2) = edges[n] + + # Solve for the intersection (if one exists): + coeffs = np.array([[dx1, -dx2], [dy1, -dy2]]) + params = np.array([x2 - x1, y2 - y1]) + try: + solution_for_intercept = np.linalg.solve(coeffs, params) + except np.linalg.LinAlgError: + continue + + # Use the solution to check if the intersection lies between the two vertices of the first edge: + if i == len(vertices) - 1: + (x3, y3) = vertices[0] + else: + (x3, y3) = vertices[i + 1] + if dx1 != 0: + solution_next_vertex = (x3 - x1) / dx1 + else: + solution_next_vertex = (y3 - y1) / dy1 + + if solution_for_intercept[0] > 0 and solution_for_intercept[0] < solution_next_vertex: + return False + + # If all checks pass, return true: + return True + + +def inflate_polygon_2d(vertices: list, scale_distance: int) -> list: + """ + Given a list of vertices (latitude and longitude) representing a 2D polygon geometry, offset the vertices such that the perpendicular + distance between the original and offset edges is equal to the scale distance. + + Parameters + ----------- + vertices: list + List of vertices (latitude and longitude, in degrees). + scale_distance: Distance (in meters) to offset the vertices by. + + Returns + ------- + list: List of the offset vertices (latitude and longitude, in degrees). + """ + # Approximate distance (in meters) for one degree changes in latitude and longitude: + deg_to_meter = 111 * (10**3) + + # Convert from latlong to cartesian: + for i, (lat, lon) in enumerate(vertices): + x = lat * deg_to_meter + y = lon * deg_to_meter + vertices[i] = (x, y) + + # Generate a list of edges joining the vertices: + edges = np.zeros((len(vertices), 2)) + for i in range(len(edges) - 1): + (x1, y1) = vertices[i] + (x2, y2) = vertices[i + 1] + (dx, dy) = (x2 - x1, y2 - y1) + edges[i] = (dx, dy) + + # Generate the last edge (connecting last vertex in list to the first): + (x1, y1) = vertices[len(vertices) - 1] + (x2, y2) = vertices[0] + (dx, dy) = (x2 - x1, y2 - y1) + edges[len(edges) - 1] = (dx, dy) + # For each edge, generate a perpendicular vector with a mangitude of the scaling distance: + perp_vecs = np.zeros((len(edges), 2)) + for i, (x, y) in enumerate(edges): + p_vec = np.array([y, -x]) + p_vec *= scale_distance / np.linalg.norm(p_vec) + perp_vecs[i] = p_vec + # Offset the vertices by the perpendicular vectors: + offset_verts = np.zeros((len(edges), 2)) + for i, (x, y) in enumerate(vertices): + (dx, dy) = perp_vecs[i] + offset_verts[i] = np.array([x + dx, y + dy]) + # Solve for the intersections of the edges originating from the offset vertices: + inf_verts = np.zeros((len(vertices), 2)) + for i in range(len(offset_verts) - 1): + coeffs = np.array([[edges[i][0], -edges[i + 1][0]], [edges[i][1], -edges[i + 1][1]]]) + params = np.array( + [ + offset_verts[i + 1][0] - offset_verts[i][0], + offset_verts[i + 1][1] - offset_verts[i][1], + ] + ) + solution = np.linalg.solve(coeffs, params) + inf_verts[i + 1] = offset_verts[i] + solution[0] * edges[i] + + # Solve for intersection of last edge to first: + coeffs = np.array( + [[edges[len(edges) - 1][0], -edges[0][0]], [edges[len(edges) - 1][1], -edges[0][1]]] + ) + params = np.array( + [ + offset_verts[0][0] - offset_verts[len(edges) - 1][0], + offset_verts[0][1] - offset_verts[len(edges) - 1][1], + ] + ) + solution = np.linalg.solve(coeffs, params) + inf_verts[0] = offset_verts[len(offset_verts) - 1] + solution[0] * edges[len(edges) - 1] + + # Perform unit tests: + distances = get_distance(vertices, inf_verts) + tolerance = 1 + for dist in distances: + assert (dist - math.sqrt(2 * (scale_distance**2))) <= tolerance + + assert is_valid_polygon(inf_verts) + + # Convert cartesian to LLA: + for i, (x, y) in enumerate(inf_verts): + lat = x / deg_to_meter + lon = y / deg_to_meter + inf_verts[i] = (lat, lon) + + # Return the offset vertices: + return inf_verts diff --git a/modules/inflate_geometry.py b/modules/inflate_geometry.py new file mode 100644 index 0000000..93e6f99 --- /dev/null +++ b/modules/inflate_geometry.py @@ -0,0 +1,104 @@ +"""Functions for manipulating polygon geometry""" + +import math +import numpy as np + + +def inflate_convex_polygon(vertices: np.ndarray, scale_distance: int) -> np.ndarray: + """ + Given a list of vertices representing a convex polygon geometry, offset the vertices such that the perpendicular + distance between the original and offset edges is equal to the scale distance. + + Parameters + ----------- + vertices: np.ndarray + List of vertices + scale_distance: Distance (in meters) to offset the vertices by + + Returns + ------- + np.ndarray: List of the offset vertices + """ + # Referance radius for altitude (Radius of earth = 6371 km): + r_ref = 6371 * (10**3) + + # Convert LLA coordinates from degrees to radians: + for i, (lat, lon, alt) in enumerate(vertices): + lat_deg, lon_deg = math.radians(lat), math.radians(lon) + vertices[i] = (lat_deg, lon_deg, alt) + + # Convert to LLA to cartesian: + for i, (lat, lon, alt) in enumerate(vertices): + x = (r_ref + alt) * math.cos(lat) * math.cos(lon) + y = (r_ref + alt) * math.cos(lat) * math.sin(lon) + z = (r_ref + alt) * math.sin(lat) + vertices[i] = (x, y, z) + + # Generate a list of edges joining the vertices: + edges = np.zeros((len(vertices), 3)) + for i in range(len(edges) - 1): + (x1, y1, z1) = vertices[i] + (x2, y2, z2) = vertices[i + 1] + (dx, dy, dz) = (x2 - x1, y2 - y1, z2 - z1) + edges[i] = (dx, dy, dz) + + # Generate the last edge (connecting last vertex in list to the first). + (x1, y1, z1) = vertices[len(vertices) - 1] + (x2, y2, z2) = vertices[0] + (dx, dy, dz) = (x2 - x1, y2 - y1, z2 - z1) + edges[len(edges) - 1] = (dx, dy, dz) + + # Calculate the direction vector of the angle bisectors: + bisectors = np.zeros((len(edges), 3)) + for i in range(len(bisectors) - 1): + edge1 = np.array(edges[i]) + edge2 = np.array(edges[i + 1]) + edge1 = edge1 / np.linalg.norm(edge1) + edge2 = edge2 / np.linalg.norm(edge2) + bisector = edge1 + edge2 + bisector /= np.linalg.norm(bisector) + + # Calculate the required norm of the bisector vector given the scaling factor + l = scale_distance / math.sqrt((1 + np.dot(edge1, edge2)) / 2) + bisector *= l + bisectors[i] = bisector + + # Calculate final bisector (final edge with first) + edge1 = np.array(edges[len(edges) - 1]) + edge2 = np.array(edges[0]) + edge1 = edge1 / np.linalg.norm(edge1) + edge2 = edge2 / np.linalg.norm(edge2) + bisector = edge1 + edge2 + bisector /= np.linalg.norm(bisector) + + # Calculate the required norm of the bisector vector given the scaling factor + l = scale_distance / math.sqrt((1 + np.dot(edge1, edge2)) / 2) + bisector *= l + bisectors[len(bisectors) - 1] = bisector + + # Offset the vertices by the bisector vectors + offset_verts = np.zeros((len(edges), 3)) + for i, ((x1, y1, z1), (x2, y2, z2)) in enumerate(zip(vertices, bisectors)): + offset_verts[i] = (x1 - x2, y1 - y2, z1 - z2) + + # Convert cartesian to LLA + for i, (x, y, z) in enumerate(offset_verts): + (x, y, z) = offset_verts[i] + lon = math.atan2(y, x) + lat = math.atan2(z, x) + alt = math.sqrt(x**2 + y**2 + z**2) - r_ref + + # Convert lat and long from radians to degrees + lon = math.degrees(lon) + lat = math.degrees(lat) + + offset_verts[i] = lat, lon, alt + + # Return the offset vertices + return offset_verts + + +triangle = np.array( + [[48.8567, 2.3508, 0], [61.4140105652, 23.7281341313, 0], [51.760597, -1.261247, 0]] +) +print(inflate_convex_polygon(triangle, 5))