get_relation.py

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

# Read a json file with relation data and send to response object
def get_relation(self, body_of_water: str):  # NB: implement body_of_water
    # Read relation from GeoJson file and extract all polygons
    geo_data = gpd.read_file("server/map/mjosa.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 = []

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

        divided_map.extend(combine_grid_with_poly(polygon, lines))

    '''
    ####################### PLOTTING ############################
    tiles = [gpd.GeoDataFrame(geometry=[tile]) for tile in divided_map]

    print("Plotting... This may take some time...")
    # NB test plot
    fig, ax = plt.subplots()
    ax.set_aspect(1.5)

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


    plt.show()
    ##################### PLOTTIND END ###########################
    '''

    features = []

    sub_div_id = 0
    for tile in divided_map: # NB temporarily limited to 5 tiles

        # Round coordinates to 4 decimals
        center = round(tile.centroid.coords[0][0], 4), round(tile.centroid.coords[0][1], 4)
        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
    }

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

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


def create_grid(poly: Polygon, cell_size):
    # 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_size

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

    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(path: str, file_name: str, json_data: dict):
    print("Writing to file...")
    if not os.path.exists(path):
        raise Exception("Directory from path does not exist")

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


'''
def get_divided_map(file_name):
    geo_data = gpd.read_file("server/map/mjosa.geojson")
    polygon_data = geo_data[geo_data['geometry'].geom_type == 'Polygon']
    polygons = [Polygon(polygon.exterior) for polygon in polygon_data['geometry']]
'''