import geopandas as gpd from shapely.geometry import Polygon, Point, LineString, MultiPolygon import matplotlib.pyplot as plt from shapely.ops import linemerge, unary_union, polygonize import matplotlib.ticker as ticker import numpy as np def get_relation(self, body_of_water: str): # Load GeoJSON data using geopandas geo_data = gpd.read_file("server/map/mjosa.geojson") # Filter only polygons, exclude points and other feature types to reduce response size polygon_data = geo_data[geo_data['geometry'].geom_type == 'Polygon'] if polygon_data.empty: raise ValueError("Failed to extract polygon data from file") # Extract coordinates from polygons and create polygon objects polygons = [Polygon(polygon.exterior) for polygon in polygon_data['geometry']] if not polygons: raise ValueError("Failed to convert to polygons") divided_map = [] # Divide all polygons into sections for polygon in polygons: # Divide the length and with of polygon into a grid grid_lines = create_grid_coords(polygon, cell_size=0.1) hrz_lines = grid_lines[0] # Horizontal coordinates vrt_lines = grid_lines[1] # Vertical coordinates # Cut polygon into horizontal sections for hrz_line in hrz_lines: # Split shape into upper and lower section as hrz_line as divider cut_poly_1 = cut_polygon_in_two(polygon, hrz_line, True) polygon_part = cut_poly_1[0] # Save upper horizontal section if not polygon_part or not cut_poly_1[0]: continue # Cut each horizontal section into vertical sections for vrt_line in vrt_lines: cut_poly_2 = cut_polygon_in_two(polygon_part, vrt_line, False) if not cut_poly_2[0]: continue divided_map.extend(cut_poly_2[0]) # Append part to final list of shapes # Set polygon_part to the remaining, un-split, horizontal section for next iteration polygon_part = cut_poly_2[1] polygon = cut_poly_1[1] # Set polygon to the remaining, un-split shape for next iteration tiles = gpd.GeoDataFrame(geometry=divided_map) # NB: test plots fig, ax = plt.subplots() ax.set_aspect(1.5) ax.xaxis.set_major_locator(ticker.MultipleLocator(0.2)) tiles.plot(ax=ax, color='blue', edgecolor='orange', alpha=0.5, label='Original Polygon') plt.show() # Convert GeoDataFrame to GeoJSON tiles_json = tiles.to_json() # Set headers self.send_response(200) self.send_header("Content-type", "application/json") self.end_headers() # Write GeoJSON to response object self.wfile.write(tiles_json.encode('utf-8')) # Takes a polygon and divides its coordinates into two shapes, where divisor is a # coordinate that defines the point of division def cut_polygon_in_two(polygon: Polygon, divisor: float, horizontal: bool): # 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 point in exterior_coords: point = Point(point) # Convert coordinates to Shapely Point object if horizontal: # Horizontal split if point.y < divisor: split_shape.append(point) else: remaining_shape.append(point) else: # Vertical split if point.x < divisor: split_shape.append(point) else: remaining_shape.append(point) # Check if polygons have enough coordinates if len(split_shape) < 3 or len(remaining_shape) < 3: print("Not enough coordinates to create valid polygons") return None, None # Append the first coordinate of the shapes to create closed loop split_shape.append(split_shape[0]) # NB: may not be necessary? remaining_shape.append(remaining_shape[0]) return Polygon(split_shape), Polygon(remaining_shape) # Return split polygons as Shapely Polygon objects # 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): min_x, min_y, max_x, max_y = polygon.bounds x_coords = np.arange(min_x, max_x, cell_size) y_coords = np.arange(min_y, max_y, cell_size) if len(x_coords) == 0 or len(y_coords) == 0: raise ValueError("No grid points generated") return x_coords, y_coords