import json from pathlib import Path from shapely.geometry import Polygon """ CityRegionSplitter is a tool for splitting a city into a subregions with overlaps allowing for proper SRA calculations. SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2023 Concordia CERC group Project Coder Koa Wells kekoa.wells@concordia.ca """ class CityRegionSplitter: def __init__(self, city_file): self._city_file = city_file self._city = json.load(city_file) self._buildings = self._city['features'] self._regions_dimension = -1 self._regions = [] def _find_neighbors(self): regions = [] count = 1 for row in range(self._regions_dimension): regions.append([]) for col in range(self._regions_dimension): regions[row].append(count) count += 1 rows = len(regions) cols = len(regions[0]) neighbors = {} # Iterate over each element in the grid for i in range(rows): for j in range(cols): current_region = f'{regions[i][j]}_v1' while len(current_region) < 7: current_region = f'0{current_region}' upper_neighbor = 'None' lower_neighbor = 'None' left_neighbor = 'None' right_neighbor = 'None' upper_left_neighbor = 'None' upper_right_neighbor = 'None' lower_left_neighbor = 'None' lower_right_neighbor = 'None' # find upper, lower, left, and right adjacent if i > 0: upper_neighbor = f'{regions[i - 1][j]}_v1' if i < rows - 1: lower_neighbor = f'{regions[i + 1][j]}_v1' if j > 0: left_neighbor = f'{regions[i][j - 1]}_v1' if j < cols - 1: right_neighbor = f'{regions[i][j + 1]}_v1' # find diagonal neighbors # corners if i == 0 and j == 0: lower_right_neighbor = f'{regions[i + 1][j + 1]}_v1' if i == 0 and j == self._regions_dimension - 1: lower_left_neighbor = f'{regions[i + 1][j - 1]}_v1' if i == self._regions_dimension - 1 and j == 0: upper_right_neighbor = f'{regions[i - 1][j + 1]}_v1' if i == self._regions_dimension - 1 and j == self._regions_dimension - 1: upper_left_neighbor = f'{regions[i - 1][j - 1]}_v1' # middle top row if i == 0 and 0 < j < self._regions_dimension - 1: lower_left_neighbor = f'{regions[i + 1][j - 1]}_v1' lower_right_neighbor = f'{regions[i + 1][j + 1]}_v1' # middle bottom row if i == self._regions_dimension - 1 and 0 < j < self._regions_dimension - 1: upper_left_neighbor = f'{regions[i - 1][j - 1]}_v1' upper_right_neighbor = f'{regions[i - 1][j + 1]}_v1' # middle left column if 0 < i < self._regions_dimension - 1 and j == 0: upper_right_neighbor = f'{regions[i - 1][j + 1]}_v1' lower_right_neighbor = f'{regions[i + 1][j + 1]}_v1' # middle right column if 0 < i < self._regions_dimension - 1 and j == self._regions_dimension - 1: upper_left_neighbor = f'{regions[i - 1][j - 1]}_v1' lower_left_neighbor = f'{regions[i + 1][j - 1]}_v1' # everything else if 0 < i < self._regions_dimension - 1 and 0 < j < self._regions_dimension - 1: lower_right_neighbor = f'{regions[i + 1][j + 1]}_v1' lower_left_neighbor = f'{regions[i + 1][j - 1]}_v1' upper_right_neighbor = f'{regions[i - 1][j + 1]}_v1' upper_left_neighbor = f'{regions[i - 1][j - 1]}_v1' neighbors[current_region] = [f'{upper_neighbor}', f'{lower_neighbor}', f'{left_neighbor}', f'{right_neighbor}', f'{upper_left_neighbor}', f'{upper_right_neighbor}', f'{lower_left_neighbor}', f'{lower_right_neighbor}'] final_neighbors = _add_trailing_zeros(neighbors) return final_neighbors """ Generates n x n subregions with a percentage overlap using the top_left_coordinate and bottom_right_coordinate :param: top_left_coordinate: [longitude, latitude] of top left corner of city :param: bottom_right_coordinate: [longitude, latitude] of bottom right corner of city :param: region_dimensions: the number of regions in each directiory; for example 20 creates 20 x 20 subregions :param: percent_overlap: percentage of overlap between neighboring regions :param: output_file_name: name out output file :param: output_path: file path to save the region geojson file :return: None """ def generate_regions(self, top_left_coordinate, bottom_right_coordinate, regions_dimension, percent_overlap, output_file_name, output_path): # points for # top_left_coordinate = [-74.022687, 45.722604] # bottom_right_coordinate = [-73.397182, 45.384048] top_left_coordinate = top_left_coordinate bottom_right_coordinate = bottom_right_coordinate regions_dimension = regions_dimension self._regions_dimension = regions_dimension percent_overlap = percent_overlap output_file_name = output_file_name output_path = Path(f'{output_path}/{output_file_name}').resolve() x1 = top_left_coordinate[0] y1 = top_left_coordinate[1] x2 = bottom_right_coordinate[0] y2 = bottom_right_coordinate[1] regions = [[[] for i in range(regions_dimension)] for j in range(regions_dimension)] region_x_side_length = (x2 - x1) / regions_dimension region_y_side_length = (y2 - y1) / regions_dimension x_current = x1 y_current = y1 row_count = 1 for row in regions: for column in row: # append top left corner column.append([x_current, y_current]) # append bottom left corner column.append([x_current, y_current + region_y_side_length * (1 + percent_overlap / 2)]) # append bottom right corner column.append([x_current + region_x_side_length * (1 + percent_overlap / 2), y_current + region_y_side_length * (1 + percent_overlap / 2)]) # append top right corner column.append([x_current + region_x_side_length * (1 + percent_overlap / 2), y_current]) # append top left corner again to close the polygon column.append([x_current, y_current]) x_current = x_current + region_x_side_length * (1 - percent_overlap / 2) x_current = x1 y_current = y_current + region_y_side_length * (1 - percent_overlap / 2) row_count += 1 id = 1 feature_collection = { "type": "FeatureCollection", "features": [] } for row in regions: for column in row: id_with_zeros = str(id) while len(id_with_zeros) < 4: id_with_zeros = f'0{id_with_zeros}' feature_collection['features'].append( { "type": "Feature", "id": id, "geometry": { "type": "Polygon", "coordinates": [column] }, "properties": { "region_id": f'{id_with_zeros}.bz2' } }) id += 1 with open(f'{output_path}.geojson', 'w') as file: file.write(json.dumps(feature_collection, indent=2)) print(f'Saving region file to {output_path}') self._regions = feature_collection['features'] def assign_buildings_to_regions(self): for building in self._buildings: if building['geometry']['type'] == 'Polygon': building_centroid = Polygon(building["geometry"]["coordinates"][0]).centroid elif building['geometry']['type'] == 'MultiPolygon': # use the centroid of the first polygon inside of the multipolygon building_centroid = Polygon(building["geometry"]["coordinates"][0][0]).centroid building['properties']['centroid'] = [building_centroid.x, building_centroid.y] building['properties']['district_property'] = [] target_regions = [] region_assigned = False target_region = self._regions[len(self._regions) / 2] while not region_assigned: if building_centroid.within(Polygon(target_region['geometry']['coordinates'])): region_assigned = True break target_region_centroid = Polygon(target_region['geometry']['coordinates']).centroid if building_centroid.x <= target_region_centroid.x: if building_centroid.y >= target_region_centroid.y: regions = regions[0:len(regions) / 2] elif building_centroid.y < target_region_centroid.y: regions = regions[len(regions) / 2:len(regions) - 1] elif building_centroid.x > target_region_centroid.x: if building_centroid.y >= target_region_centroid.y: regions = regions[0:len(regions) / 2] elif building_centroid.y < target_region_centroid.y: regions = regions[len(regions) / 2:len(regions) - 1]