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']] '''