232 lines
8.8 KiB
Python
232 lines
8.8 KiB
Python
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] |