diff --git a/server/map/__pycache__/get_relation.cpython-311.pyc b/server/map/__pycache__/get_relation.cpython-311.pyc index f13e9b873b357cf5927106718f674e4048732460..d3f0c0786e20b0ddc077d121d821cf744be96425 100644 Binary files a/server/map/__pycache__/get_relation.cpython-311.pyc and b/server/map/__pycache__/get_relation.cpython-311.pyc differ diff --git a/server/map/get_relation.py b/server/map/get_relation.py index 5003010739d56774e2fef78de1e60d02f2479123..6d7fc7ede4d0cdbb513db99a18a7e4ee40bb97d3 100644 --- a/server/map/get_relation.py +++ b/server/map/get_relation.py @@ -1,5 +1,6 @@ import geopandas as gpd -from shapely.geometry import Polygon, Point, LineString, MultiPolygon +from shapely.geometry import Polygon, LineString +from shapely.ops import linemerge, unary_union, polygonize import matplotlib.pyplot as plt import matplotlib.ticker as ticker import random @@ -18,44 +19,14 @@ def get_relation(self, body_of_water: str): # NB: implement body_of_water # Divide all polygons into sections for polygon in polygons: - - cell_size = 0.4 + cell_size = 1 # Divide the length and with of polygon into a grid of equally sized parts - grid_lines = create_grid_coords(polygon, cell_size) - - vrt_lines = grid_lines[0] # Vertical grid coordinates - hrz_lines = grid_lines[1] # Horizontal grid coordinates - - # Cut polygon into horizontal sections, bottom to top - for hrz_line in hrz_lines: - - # Split shape into upper and lower section as hrz_line as divider - divided_poly = cut_polygon_by_points(polygon, hrz_line, cell_size) - horizontal_section = divided_poly[0] # Save upper horizontal section - - polygon = divided_poly[1] # Set polygon to the remaining, un-split shape for next iteration - - if not horizontal_section or not divided_poly[0]: - continue - - # Cut each horizontal section into vertical sections, right to left - for vrt_line in vrt_lines: - if len(horizontal_section.exterior.coords) < 3: - break # Break from loop im remaining section has no coordinates - - # Split the horizontal section into two vertical parts - vertical_parts = cut_polygon_by_points(horizontal_section, vrt_line, -0.1) - - divided_map.append(vertical_parts[0]) # Append split vertical sections to final list of shapes + lines = create_grid(polygon, cell_size) - # Set horizontal_section to the remaining, un-split, horizontal section for next iteration - horizontal_section = vertical_parts[1] - - polygon = divided_poly[1] # Set polygon to the remaining, un-split shape for next iteration - - divided_map.append(polygon) - - break + lines.append(polygon.boundary) + lines = unary_union(lines) + lines = linemerge(lines) + divided_map.extend(list(polygonize(lines))) tiles = gpd.GeoDataFrame(geometry=divided_map) @@ -85,77 +56,6 @@ def get_relation(self, body_of_water: str): # NB: implement body_of_water self.wfile.write(tiles_json.encode('utf-8')) -def cut_polygon_by_points(polygon: Polygon, divisor: float, cell_size: float): - # Extract polygon exterior coordinates - exterior_coords = list(polygon.exterior.coords) - - # Initialize lists to store coordinates of new shapes after split - split_shape = [] - remaining_shape = [] - - # Loop through points and check which side of the division line they are - for i in range(len(exterior_coords) - 1): - point1 = Point(exterior_coords[i]) - point2 = Point(exterior_coords[i + 1]) - - # Check if both points are on the same side of the divisor - if (point1.y < divisor and point2.y < divisor) or (point1.y >= divisor and point2.y >= divisor): - # Both points on same side, add to appropriate shape - if point1.y < divisor: - split_shape.append(exterior_coords[i]) - else: - remaining_shape.append(exterior_coords[i]) - else: - # Points on different sides, calculate intersection with divisor - x_intersect = point1.x + (divisor - point1.y) * (point2.x - point1.x) / (point2.y - point1.y) - - # Add intersection point to both shapes - split_shape.append((x_intersect, divisor)) - remaining_shape.append((x_intersect, divisor)) - - # Determine which side each original point belongs to - if point1.y < divisor: - split_shape.append((point2.x, point2.y)) - else: - remaining_shape.append((point2.x, point2.y)) - - # Populate newly created edges with points to facilitate vertical cutting - populated_split = populate_new_edge(split_shape, cell_size, divisor) - if populated_split is not None: - split_shape = populated_split - - populated_rem = populate_new_edge(remaining_shape, cell_size, divisor) - if populated_rem is not None: - remaining_shape = populated_rem - - # Ensure that the shapes are closed loops - if len(split_shape) > 0 and split_shape[0] != split_shape[-1]: - split_shape.append(split_shape[0]) - - if len(remaining_shape) > 0 and remaining_shape[0] != remaining_shape[-1]: - remaining_shape.append(remaining_shape[0]) - - return Polygon(split_shape), Polygon(remaining_shape) - - -# Generate grid of equally spaced x and y coordinates where the grid size is determined by cell_size -def create_grid_coords(polygon: Polygon, cell_size: float): - # Define boundaries of grid - min_x, min_y, max_x, max_y = polygon.bounds - - global polygon_min_x - polygon_min_x = min_x - - # Divide grid into sections of size *cell_size - x_step = max((max_x - min_x) / 10, cell_size) - y_step = max((max_y - min_y) / 10, cell_size) - - x_coords = np.arange(min_x, max_x + x_step, x_step) - y_coords = np.arange(min_y, max_y + y_step, y_step) - - return x_coords, y_coords - - # NB: only for testing def circle_polygon(): circle_points = [] @@ -173,42 +73,25 @@ def circle_polygon(): return Polygon(circle_points) -def populate_new_edge(polygon, cell_size: float, divisor: float): - if polygon is not None and polygon_min_x is not None: - - final_shape = [] - left_edge = [] - left_most_point = None - - # Add all coordinates on the right side of the divisor first - for i in range(len(polygon)): - x, y = polygon[i] - if x > divisor: - final_shape.append((x, y)) - elif x < divisor: - left_edge.append((x, y)) - - if len(final_shape) < 1 or len(left_edge) < 1: - return None - - # Find most left point in the polygon - for point in left_edge: - if not left_most_point: - left_most_point = point[0] - elif point[0] < left_most_point: - left_most_point = point[0] +def create_grid(poly: Polygon, cell_size): + # Retrieve bounds of the entire polygon + bounds = poly.bounds - # Sort final_shape from bottom to top - if final_shape[0][1] > final_shape[-1][1]: - final_shape = sorted(final_shape, key=lambda point: final_shape[1]) + min_x, min_y, max_x, max_y = bounds + grid_lines = [] - i = 0 - while polygon[i][0] < left_most_point: # While to haven't reached left edge - final_shape.append(polygon[i][0] + cell_size) - i += 1 + # Horizontal lines + y = min_y + while y <= max_y: + line = LineString([(min_x, y), (max_x, y)]) + grid_lines.append(line) + y += cell_size - # NB VERY IMPORTANT MUST CREATE LOGIC FOR SHAPES THAT ARE SPLIT - final_shape.append(left_edge) - final_shape.append(final_shape[0]) + # Vertical lines + x = min_x + while x <= max_x: + line = LineString([(x, min_y), (x, max_y)]) + grid_lines.append(line) + x += cell_size - return final_shape + return grid_lines