add_new_lake.py

import os
import json
from math import cos, sqrt, fabs
import random
import geopandas as gpd
from matplotlib import pyplot as plt
from shapely.ops import linemerge, unary_union, polygonize
from shapely.geometry import Polygon, LineString, MultiLineString

from server.consts import LAKE_RELATIONS_PATH

'''
# 0
starting coordinate (x,y)

# 1
1 deg lat = 111.32km
1 deg lng = 40075km * cos( lat ) / 360

# 2 Formulas for calculating a distance in kilometers to a distance in latitude and longitude
lat = distance_in_km/111.32km
lng = (distance_in_km × 360)/(40075km × cos(lat))

'''


# Read a json file with relation data and send to response object
def cut_map(self, cursor, lake_name: str):
    try:
        # Read relation from GeoJson file and extract all polygons
        geo_data = gpd.read_file(LAKE_RELATIONS_PATH + lake_name + ".geojson")
        polygon_data = geo_data[geo_data['geometry'].geom_type == 'Polygon']
        polygons = [Polygon(polygon.exterior) for polygon in polygon_data['geometry']]

        if len(polygons) <= 1:
            raise Exception("Failed to convert JSON object to Shapely Polygons")

        divided_map = []

        dist_in_km = 0.08  # NB placeholder
        cell_size = dist_in_km * 10
        cell_width = cell_size / 111.32
        cell_height = (cell_size * 360) / (40075 * cos(cell_width))

        for polygon in polygons:
            lines = create_grid(polygon, cell_width, cell_height)
            lines.append(polygon.boundary)
            lines = unary_union(lines)
            lines = linemerge(lines)
            lines = list(polygonize(lines))

            divided_map.extend(combine_grid_with_poly(polygon, lines))

        features = []

        sub_div_id = 0
        for tile in divided_map:

            # Calculate tile center based on bounds, and round down to two decimals
            min_x, min_y, max_x, max_y = tile.bounds
            center = round(max_y - (max_y - min_y), 6), round(max_x - (max_x - min_x), 6)

            rounded_coordinates = []
            if isinstance(tile, Polygon):
                for coords in tile.exterior.coords:
                    rounded_coords = (round(coords[0], 4), round(coords[1], 4))
                    rounded_coordinates.append(rounded_coords)
                rounded_tile = Polygon(rounded_coordinates)

            tile_feature = {
                'type': 'Feature',
                'properties': {
                    'sub_div_id': str(sub_div_id),
                    'sub_div_center': center
                },
                'geometry': rounded_tile.__geo_interface__
            }
            features.append(tile_feature)
            sub_div_id += 1

        feature_collection = {
            'type': 'FeatureCollection',
            'features': features,
            'tile_count': sub_div_id,  # Add the last subdivision ID as number of tiles
        }

        # Add lake name to database
        cursor.execute('''
            SELECT Name FROM BodyOfWater WHERE Name = ?;
        ''', (lake_name,))
        existing_lake = cursor.fetchone()

        # If lake_name doesn't exist, insert it into the database
        if existing_lake is None:
            cursor.execute('''
                INSERT INTO BodyOfWater(Name) VALUES (?);
            ''', (lake_name,))

        plot_map(divided_map)

        write_json_to_file(lake_name, feature_collection)

        self.send_response(200)
        self.send_header("Content-type", "application/json")
        self.end_headers()

        self.wfile.write(json.dumps(feature_collection).encode('utf-8'))

    except Exception as e:
        print(f"Error in adding new map: {e}")

        self.send_response(500)
        self.send_header("Content-type", "application/json")
        self.end_headers()


def create_grid(poly: Polygon, cell_width, cell_height):
    # Retrieve bounds of the entire polygon
    bounds = poly.bounds

    min_x, min_y, max_x, max_y = bounds
    grid_lines = []

    # Horizontal lines
    y = min_y
    while y <= max_y:
        line = LineString([(min_x, y), (max_x, y)])
        grid_lines.append(line)
        y += cell_height

    # Vertical lines
    x = min_x
    while x <= max_x:
        line = LineString([(x, min_y), (x, max_y)])
        grid_lines.append(line)
        x += cell_width

    return grid_lines


def combine_grid_with_poly(polygon, grid):
    intersecting_tiles = []

    for line in grid:
        if line.intersects(polygon):
            intersection = line.intersection(polygon)
            # Check if intersection is a MultiLineString
            if isinstance(intersection, MultiLineString):
                # Extend the intersecting tiles with the polygonized results
                intersecting_tiles.extend(list(polygonize(intersection)))
            else:
                intersecting_tiles.append(intersection)

    return intersecting_tiles


def write_json_to_file(lake_name: str, json_data: dict):
    # Create and write divided map to new file
    print("Writing to file...")
    if not os.path.exists(LAKE_RELATIONS_PATH):
        raise Exception("Directory from path does not exist")

    with open(LAKE_RELATIONS_PATH + '/' + lake_name + '_div.json', 'w') as f:
        json.dump(json_data, f)

    # Update all_system_lakes
    with open(LAKE_RELATIONS_PATH + 'all_lake_names.json', 'r') as file:
        data = json.load(file)

    if lake_name not in data:
        data.append(lake_name)

        with open(LAKE_RELATIONS_PATH + 'all_lake_names.json', 'w') as file:
            # json.dump(data, file, indent=2)
            json.dump(data, file, ensure_ascii=False, indent=2)


def plot_map(divided_map):
    tiles = [gpd.GeoDataFrame(geometry=[tile]) for tile in divided_map]

    print("Plotting... This may take some time...")

    fig, ax = plt.subplots()
    ax.set_aspect(1.5)

    # Plot each tile
    for tile in tiles:
        random_color = "#{:06x}".format(random.randint(0, 0xFFFFFF))
        gpd.GeoSeries(tile.geometry).plot(ax=ax, facecolor=random_color, edgecolor='none')

    plt.show()