get_relation.py

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:
                split_shape.append(point)
            else:
                remaining_shape.append(point)

        if len(split_shape) > 2:
            # Get last point added to
            last_point = split_shape[len(split_shape)-1]
            # Get length of the newly created edge
            new_edge_len = last_point.x - split_shape[0].x - 0.0001
            print("new_edge_len: ", new_edge_len, "    cell_size: ", cell_size)

            # Add points along the new edge to allow horizontal sections to be split into vertical ones
            while new_edge_len > cell_size:
                print("Hit")
                split_shape.append((new_edge_len-cell_size, last_point.y))
                remaining_shape.append((new_edge_len-cell_size, last_point.y))
                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