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