""" Geojson module parses geojson files and import the geometry into the city model structure SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2022 Concordia CERC group Project Coder Guillermo Gutierrez Guillermo.GutierrezMorote@concordia.ca """ import uuid import json import numpy as np from pyproj import Transformer import hub.helpers.constants as cte from hub.helpers.geometry_helper import GeometryHelper from hub.imports.geometry.helpers.geometry_helper import GeometryHelper as igh from hub.city_model_structure.attributes.polygon import Polygon from hub.city_model_structure.building import Building from hub.city_model_structure.building_demand.surface import Surface from hub.city_model_structure.city import City class Geojson: """ Geojson class """ _X = 0 _Y = 1 def __init__(self, path, aliases_field=None, extrusion_height_field=None, year_of_construction_field=None, function_field=None, function_to_hub=None): # todo: destination epsg should change according actual the location self._transformer = Transformer.from_crs('epsg:4326', 'epsg:26911') self._min_x = cte.MAX_FLOAT self._min_y = cte.MAX_FLOAT self._max_x = cte.MIN_FLOAT self._max_y = cte.MIN_FLOAT self._max_z = 0 self._city = None self._aliases_field = aliases_field self._extrusion_height_field = extrusion_height_field self._year_of_construction_field = year_of_construction_field self._function_field = function_field self._function_to_hub = function_to_hub with open(path, 'r', encoding='utf8') as json_file: self._geojson = json.loads(json_file.read()) def _save_bounds(self, x, y): if x > self._max_x: self._max_x = x if x < self._min_x: self._min_x = x if y > self._max_y: self._max_y = y if y < self._min_y: self._min_y = y @staticmethod def _find_wall(line_1, line_2): for i in range(0, 2): j = 1 - i point_1 = line_1[i] point_2 = line_2[j] distance = GeometryHelper.distance_between_points(point_1, point_2) if distance > 1e-2: return False return True def _store_shared_percentage_to_walls(self, city, city_mapped): for building in city.buildings: if building.name not in city_mapped.keys(): for wall in building.walls: wall.percentage_shared = 0 continue building_mapped = city_mapped[building.name] for wall in building.walls: percentage = 0 ground_line = [] for point in wall.perimeter_polygon.coordinates: if point[2] < 0.5: ground_line.append(point) for entry in building_mapped: if building_mapped[entry]['shared_points'] <= 2: continue line = [building_mapped[entry]['line_start'], building_mapped[entry]['line_end']] neighbour_line = [building_mapped[entry]['neighbour_line_start'], building_mapped[entry]['neighbour_line_end']] neighbour_height = city.city_object(building_mapped[entry]['neighbour_name']).max_height if self._find_wall(line, ground_line): line_shared = (GeometryHelper.distance_between_points(line[0], line[1]) + GeometryHelper.distance_between_points(neighbour_line[0], neighbour_line[1]) - GeometryHelper.distance_between_points(line[1], neighbour_line[0]) - GeometryHelper.distance_between_points(line[0], neighbour_line[1])) / 2 percentage_ground = line_shared / GeometryHelper.distance_between_points(line[0], line[1]) percentage_height = neighbour_height / building.max_height percentage_height = min(percentage_height, 1) percentage += percentage_ground * percentage_height wall.percentage_shared = percentage @property def city(self) -> City: """ Get city out of a Geojson file """ if self._city is None: buildings = [] lod = 0 for feature in self._geojson['features']: extrusion_height = 0 if self._extrusion_height_field is not None: extrusion_height = float(feature['properties'][self._extrusion_height_field]) lod = 1 year_of_construction = None if self._year_of_construction_field is not None: year_of_construction = int(feature['properties'][self._year_of_construction_field]) function = None if self._function_field is not None: function = feature['properties'][self._function_field] if self._function_to_hub is not None: # use the transformation dictionary to retrieve the proper function if function in self._function_to_hub: function = self._function_to_hub[function] geometry = feature['geometry'] building_aliases = [] if 'id' in feature: building_name = feature['id'] else: building_name = uuid.uuid4() if self._aliases_field is not None: for alias_field in self._aliases_field: building_aliases.append(feature['properties'][alias_field]) if str(geometry['type']).lower() == 'polygon': buildings.append(self._parse_polygon(geometry['coordinates'], building_name, building_aliases, function, year_of_construction, extrusion_height)) elif str(geometry['type']).lower() == 'multipolygon': buildings.append(self._parse_multi_polygon(geometry['coordinates'], building_name, building_aliases, function, year_of_construction, extrusion_height)) else: raise NotImplementedError(f'Geojson geometry type [{geometry["type"]}] unknown') self._city = City([self._min_x, self._min_y, 0.0], [self._max_x, self._max_y, self._max_z], 'epsg:26911') for building in buildings: # Do not include "small building-like structures" to buildings if building.floor_area >= 25: self._city.add_city_object(building) self._city.level_of_detail.geometry = lod for building in self._city.buildings: building.level_of_detail.geometry = lod if lod > 0: lines_information = GeometryHelper.city_mapping(self._city, plot=False) self._store_shared_percentage_to_walls(self._city, lines_information) return self._city def _polygon_coordinates_to_3d(self, polygon_coordinates): transformed_coordinates = '' for coordinate in polygon_coordinates: transformed = self._transformer.transform(coordinate[self._Y], coordinate[self._X]) self._save_bounds(transformed[self._X], transformed[self._Y]) transformed_coordinates = f'{transformed_coordinates} {transformed[self._X]} {transformed[self._Y]} 0.0' return transformed_coordinates.lstrip(' ') def _parse_polygon(self, coordinates, building_name, building_aliases, function, year_of_construction, extrusion_height): surfaces = [] for polygon_coordinates in coordinates: points = igh.points_from_string( igh.remove_last_point_from_string( self._polygon_coordinates_to_3d(polygon_coordinates) ) ) points = igh.invert_points(points) polygon = Polygon(points) polygon.area = igh.ground_area(points) surface = Surface(polygon, polygon) if surface.type == cte.GROUND: surfaces.append(surface) else: distance = cte.MAX_FLOAT hole_connect = 0 surface_connect = 0 for hole_index, hole_coordinate in enumerate(polygon.coordinates): for surface_index, ground_coordinate in enumerate(surfaces[-1].solid_polygon.coordinates): current_distance = GeometryHelper.distance_between_points(hole_coordinate, ground_coordinate) if current_distance < distance: distance = current_distance hole_connect = hole_index surface_connect = surface_index hole = polygon.coordinates[hole_connect:] + polygon.coordinates[:hole_connect] + [polygon.coordinates[hole_connect]] prefix_coordinates = surfaces[-1].solid_polygon.coordinates[:surface_connect+1] trail_coordinates = surfaces[-1].solid_polygon.coordinates[surface_connect:] coordinates = prefix_coordinates + hole + trail_coordinates polygon = Polygon(coordinates) polygon.area = igh.ground_area(coordinates) surfaces[-1] = Surface(polygon, polygon) if len(surfaces) > 1: raise ValueError('too many surfaces!!!!') building = Building(f'{building_name}', surfaces, year_of_construction, function) for alias in building_aliases: building.add_alias(alias) if extrusion_height == 0: return building volume = 0 for ground in building.grounds: volume += ground.solid_polygon.area * extrusion_height roof_coordinates = [] # adding a roof means invert the polygon coordinates and change the Z value for coordinate in ground.solid_polygon.coordinates: roof_coordinate = np.array([coordinate[0], coordinate[1], extrusion_height]) # insert the roof rotated already roof_coordinates.insert(0, roof_coordinate) roof_polygon = Polygon(roof_coordinates) roof_polygon.area = ground.solid_polygon.area roof = Surface(roof_polygon, roof_polygon) surfaces.append(roof) # adding a wall means add the point coordinates and the next point coordinates with Z's height and 0 coordinates_length = len(roof.solid_polygon.coordinates) for i, coordinate in enumerate(roof.solid_polygon.coordinates): j = i + 1 if j == coordinates_length: j = 0 next_coordinate = roof.solid_polygon.coordinates[j] wall_coordinates = [ np.array([coordinate[0], coordinate[1], 0.0]), np.array([next_coordinate[0], next_coordinate[1], 0.0]), np.array([next_coordinate[0], next_coordinate[1], next_coordinate[2]]), np.array([coordinate[0], coordinate[1], coordinate[2]]) ] polygon = Polygon(wall_coordinates) wall = Surface(polygon, polygon) surfaces.append(wall) building = Building(f'{building_name}', surfaces, year_of_construction, function) for alias in building_aliases: building.add_alias(alias) building.volume = volume return building def _parse_multi_polygon(self, polygons_coordinates, building_name, building_aliases, function, year_of_construction, extrusion_height): surfaces = [] for coordinates in polygons_coordinates: for polygon_coordinates in coordinates: points = igh.points_from_string( igh.remove_last_point_from_string( self._polygon_coordinates_to_3d(polygon_coordinates) ) ) points = igh.invert_points(points) polygon = Polygon(points) polygon.area = igh.ground_area(points) surface = Surface(polygon, polygon) if surface.type == cte.GROUND: surfaces.append(surface) else: distance = cte.MAX_FLOAT hole_connect = 0 surface_connect = 0 for hole_index, hole_coordinate in enumerate(polygon.coordinates): for surface_index, ground_coordinate in enumerate(surfaces[-1].solid_polygon.coordinates): current_distance = GeometryHelper.distance_between_points(hole_coordinate, ground_coordinate) if current_distance < distance: distance = current_distance hole_connect = hole_index surface_connect = surface_index hole = polygon.coordinates[hole_connect:] + polygon.coordinates[:hole_connect] prefix_coordinates = surfaces[-1].solid_polygon.coordinates[:surface_connect] trail_coordinates = surfaces[-1].solid_polygon.coordinates[surface_connect:] coordinates = prefix_coordinates + hole + [hole[0]] + trail_coordinates polygon = Polygon(coordinates) polygon.area = igh.ground_area(coordinates) surfaces[-1] = Surface(polygon, polygon) building = Building(f'{building_name}', surfaces, year_of_construction, function) for alias in building_aliases: building.add_alias(alias) if extrusion_height == 0: return building volume = 0 for ground in building.grounds: volume += ground.solid_polygon.area * extrusion_height roof_coordinates = [] # adding a roof means invert the polygon coordinates and change the Z value for coordinate in ground.solid_polygon.coordinates: roof_coordinate = np.array([coordinate[0], coordinate[1], extrusion_height]) # insert the roof rotated already roof_coordinates.insert(0, roof_coordinate) roof_polygon = Polygon(roof_coordinates) roof_polygon.area = ground.solid_polygon.area roof = Surface(roof_polygon, roof_polygon) surfaces.append(roof) # adding a wall means add the point coordinates and the next point coordinates with Z's height and 0 coordinates_length = len(roof.solid_polygon.coordinates) for i, coordinate in enumerate(roof.solid_polygon.coordinates): j = i + 1 if j == coordinates_length: j = 0 next_coordinate = roof.solid_polygon.coordinates[j] wall_coordinates = [ np.array([coordinate[0], coordinate[1], 0.0]), np.array([next_coordinate[0], next_coordinate[1], 0.0]), np.array([next_coordinate[0], next_coordinate[1], next_coordinate[2]]), np.array([coordinate[0], coordinate[1], coordinate[2]]) ] polygon = Polygon(wall_coordinates) wall = Surface(polygon, polygon) surfaces.append(wall) building = Building(f'{building_name}', surfaces, year_of_construction, function) for alias in building_aliases: building.add_alias(alias) building.volume = volume return building