import geopandas as gpd
from shapely.geometry import Polygon, Point, LineString, MultiPolygon
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import random
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:
cell_size = 0.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_in_two(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_in_two(horizontal_section, vrt_line, -0.1)
divided_map.append(vertical_parts[0]) # Append split vertical sections to final list of shapes
# Set horizontal_section to the remaining, un-split, horizontal section for next iteration
horizontal_section = vertical_parts[1]
divided_map.append(polygon)
break
tiles = gpd.GeoDataFrame(geometry=divided_map)
# NB test plot
fig, ax = plt.subplots()
ax.set_aspect(1.5)
ax.xaxis.set_major_locator(ticker.MultipleLocator(0.2))
for i, polygon in enumerate(tiles['geometry']):
random_color = "#{:06x}".format(random.randint(0, 0xFFFFFF))
x, y = polygon.exterior.xy
ax.plot(x, y, color='blue', alpha=0.5)
ax.fill(x, y, color=random_color, alpha=1)
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, 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
if cell_size > 0: # Horizontal split
for point in exterior_coords:
point = Point(point) # Convert coordinates to Shapely Point object
if point.y < divisor: # Check if point is over or below divisor
split_shape.append(point) # Append to appropriate shape
else:
remaining_shape.append(point)
if len(split_shape) > 2:
# Get last point added to
last_point = split_shape[-1]
# Get length of the newly created edge
new_edge_len = abs(last_point.x - split_shape[0].x)
print("new_edge_len: ", new_edge_len, " cell_size: ", cell_size, " last_point.x: ", last_point.x, " split_shape[0].x: ", split_shape[0].x)
# Add points along the new edge to allow horizontal sections to be split into vertical ones
while new_edge_len > cell_size:
x_val = new_edge_len - cell_size
split_shape.insert(0, (x_val, last_point.y)) # NB may have to add/subtract small offset of 0.00001
remaining_shape.insert(0, (x_val, last_point.y)) # Prepend
new_edge_len -= cell_size
else: # Vertical split
for point in exterior_coords:
point = Point(point) # Convert coordinates to Shapely Point object
if point.x < divisor:
split_shape.append(point)
else:
remaining_shape.append(point)
# Check if the split_shape has enough coordinates to create a polygon
if len(split_shape) < 3:
# print("Not enough coordinates to create valid polygon: Split shape: ", len(split_shape))
split_shape = None
else:
split_shape.append(split_shape[0]) # Append first coord to create closed loop
# Check if the remaining_shape has enough coordinates to create a polygon
if len(remaining_shape) < 3:
# print("Not enough coordinates to create valid polygon: Remaining shape: ", len(remaining_shape))
remaining_shape = None
else:
remaining_shape.append(remaining_shape[0]) # Append first coord to create closed loop
# Return split polygons as Shapely Polygon objects
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
# Divide grid into sections of size *cell_size
x_coords = np.arange(min_x, max_x, cell_size)
y_coords = np.arange(min_y, max_y, cell_size)
# Round the coordinates to 4 decimals
x_coords_rounded = np.around(x_coords, decimals=4)
y_coords_rounded = np.around(y_coords, decimals=4)
if len(x_coords_rounded) == 0 or len(y_coords_rounded) == 0:
raise ValueError("No grid points generated")
# Return tuple of list of x coordinates and list of y coordinates
return x_coords_rounded, y_coords_rounded