""" 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 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, name_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._name_field = name_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) 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 _create_building_lod0(name, year_of_construction, function, surface_coordinates): points = igh.points_from_string(igh.remove_last_point_from_string(surface_coordinates)) points = igh.invert_points(points) polygon = Polygon(points) polygon.area = igh.ground_area(points) surface = Surface(polygon, polygon, name=f'{name}_ground') return Building(f'{name}', [surface], year_of_construction, function) @staticmethod def _create_building_lod1(name, year_of_construction, function, height, surface_coordinates): building = Geojson._create_building_lod0(name, year_of_construction, function, surface_coordinates) surfaces = [] volume = 0 for ground in building.grounds: volume += ground.solid_polygon.area * height surfaces.append(ground) 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], 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'{name}', surfaces, year_of_construction, function) building.volume = volume return building def _get_polygons(self, polygons, coordinates): if type(coordinates[0][self.X]) != float: polygons = [] for element in coordinates: polygons = self._get_polygons(polygons, element) return polygons else: transformed_coordinates = '' for coordinate in 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' polygons.append(transformed_coordinates.lstrip(' ')) return polygons @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 if percentage_height > 1: 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'] if 'id' in feature: building_name = feature['id'] if self._name_field is not None: building_name = feature['properties'][self._name_field] if str(geometry['type']).lower() == 'polygon': buildings.append(self._parse_polygon(geometry['coordinates'], building_name, function, year_of_construction, extrusion_height)) elif str(geometry['type']).lower() == 'multipolygon': buildings.append(self._parse_multi_polygon(geometry['coordinates'], building_name, 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, 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) print(surface.type, polygon.area) 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], surfaces[-1].solid_polygon.coordinates[surface_connect-1]] + 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(building_name, surfaces, year_of_construction, function) print(extrusion_height) if extrusion_height == 0: return building else: volume = 0 for ground in building.grounds: volume += ground.solid_polygon.area * extrusion_height surfaces.append(ground) 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) building.volume = volume return building """ coordinates_3d = self._polygon_coordinates_to_3d(polygon_coordinates) if extrusion_height == 0: building = Geojson._create_building_lod0(f'{building_name}', year_of_construction, function, coordinates_3d) else: if self._max_z < extrusion_height: self._max_z = extrusion_height print(building_name) building = Geojson._create_building_lod1(f'{building_name}', year_of_construction, function, extrusion_height, coordinates_3d) return building """ def _parse_multi_polygon(self, polygons_coordinates, building_name, 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(building_name, surfaces, year_of_construction, function) if extrusion_height == 0: return building else: volume = 0 for ground in building.grounds: volume += ground.solid_polygon.area * extrusion_height surfaces.append(ground) 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) building.volume = volume return building