s_ranjbar/summer_course_2024
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request 'shared_surfaces_method' (#5) from shared_surfaces_method into main

Browse Source 
Reviewed-on: https://nextgenerations-cities.encs.concordia.ca/gitea/CERC/hub/pulls/5
This commit is contained in:
Guille Gutierrez 2023-03-16 12:04:19 -04:00
commit 89254891b8
9 changed files with 212 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
hub/city_model_structure/city_object.py
View File

@ -34,6 +34,7 @@ class CityObject:
self._max_y = ConfigurationHelper().min_coordinate
self._max_z = ConfigurationHelper().min_coordinate
self._centroid = None
self._volume = None
self._external_temperature = dict()
self._global_horizontal = dict()
self._diffuse = dict()
@ -63,7 +64,13 @@ class CityObject:
Get city object volume in cubic meters
:return: float
"""
return self.simplified_polyhedron.volume
if self._volume is None:
self._volume = self.simplified_polyhedron.volume
return self._volume
@volume.setter
def volume(self, value):
self._volume = value
@property
def detailed_polyhedron(self) -> Polyhedron:

2
hub/exports/exports_factory.py
View File

@ -71,7 +71,7 @@ class ExportsFactory:
Export the city geometry to obj with grounded coordinates
:return: None
"""
return Obj(self._city, self._path).to_ground_points()
return Obj(self._city, self._path)
@property
def _sra(self):

59
hub/exports/formats/obj.py
View File

@ -11,24 +11,51 @@ from hub.exports.formats.triangular import Triangular
from hub.imports.geometry_factory import GeometryFactory
class Obj(Triangular):
class Obj:
"""
Export to obj format
"""
def __init__(self, city, path):
super().__init__(city, path, 'obj')
self._city = city
self._path = path
self._export()
def _to_vertex(self, coordinate):
x = coordinate[0] - self._city.lower_corner[0]
y = coordinate[1] - self._city.lower_corner[1]
z = coordinate[2] - self._city.lower_corner[2]
return f'v {x} {y} {z}\n'
def _export(self):
if self._city.name is None:
self._city.name = 'unknown_city'
file_name = self._city.name + '.obj'
file_path = (Path(self._path).resolve() / file_name).resolve()
vertices = {}
with open(file_path, 'w') as obj:
obj.write("# cerc-hub export\n")
vertex_index = 0
faces = []
for building in self._city.buildings:
obj.write(f'# building {building.name}\n')
for surface in building.surfaces:
obj.write(f'# surface {surface.name}\n')
face = 'f '
for coordinate in surface.perimeter_polygon.coordinates:
vertex = self._to_vertex(coordinate)
if vertex not in vertices.keys():
vertex_index += 1
vertices[vertex] = vertex_index
current = vertex_index
obj.write(vertex)
else:
current = vertices[vertex]
face = f'{face} {current}'
faces.append(f'{face} {face.split(" ")[1]}\n')
obj.writelines(faces)
faces = []
def to_ground_points(self):
"""
Move closer to the origin
"""
file_name_in = self._city.name + '.' + self._triangular_format
file_name_out = self._city.name + '_ground.' + self._triangular_format
file_path_in = (Path(self._path).resolve() / file_name_in).resolve()
file_path_out = (Path(self._path).resolve() / file_name_out).resolve()
obj = GeometryFactory('obj', path=file_path_in)
scene = obj.scene
scene.rezero()
obj_file = trimesh.exchange.obj.export_obj(scene)
with open(file_path_out, 'w') as file:
file.write(obj_file)

21
hub/helpers/geometry_helper.py
View File

@ -16,6 +16,8 @@ from hub.city_model_structure.attributes.polygon import Polygon
from hub.city_model_structure.attributes.polyhedron import Polyhedron
from hub.helpers.location import Location
from PIL import Image
class MapPoint:
def __init__(self, x, y):
@ -63,7 +65,7 @@ class GeometryHelper:
return MapPoint(((city.upper_corner[0] - coordinate[0]) * 0.5), ((city.upper_corner[1] - coordinate[1]) * 0.5))
@staticmethod
def city_mapping(city, building_names=None):
def city_mapping(city, building_names=None, plot=False):
"""
Returns a shared_information dictionary like
@ -78,6 +80,8 @@ class GeometryHelper:
y = int((city.upper_corner[1] - city.lower_corner[1]) * 0.5) + 1
city_map = [['' for _ in range(y + 1)] for _ in range(x + 1)]
map_info = [[{} for _ in range(y + 1)] for _ in range(x + 1)]
img = Image.new('RGB', (x + 1, y + 1), "black") # create a new black image
city_image = img.load() # create the pixel map
for building_name in building_names:
building = city.city_object(building_name)
line = 0
@ -103,13 +107,14 @@ class GeometryHelper:
'line_start': (coordinate[0], coordinate[1]),
'line_end': (next_coordinate[0], next_coordinate[1]),
}
city_image[x, y] = (100, 0, 0)
elif city_map[x][y] != building.name:
neighbour = city.city_object(city_map[x][y])
neighbour_info = map_info[x][y]
# prepare the keys
neighbour_start_coordinate = f'{neighbour_info["line_start"][0]}_{neighbour_info["line_start"][1]}'
building_start_coordinate = f'{coordinate[0]}_{coordinate[1]}'
neighbour_start_coordinate = f'{GeometryHelper.coordinate_to_map_point(neighbour_info["line_start"], city)}'
building_start_coordinate = f'{GeometryHelper.coordinate_to_map_point(coordinate, city)}'
neighbour_key = f'{neighbour.name}_{neighbour_start_coordinate}_{building_start_coordinate}'
building_key = f'{building.name}_{building_start_coordinate}_{neighbour_start_coordinate}'
@ -126,6 +131,10 @@ class GeometryHelper:
'line_end': (next_coordinate[0], next_coordinate[1]),
'neighbour_line_start': neighbour_info['line_start'],
'neighbour_line_end': neighbour_info['line_end'],
'coordinate_start': f"{GeometryHelper.coordinate_to_map_point(coordinate, city)}",
'coordinate_end': f"{GeometryHelper.coordinate_to_map_point(next_coordinate, city)}",
'neighbour_start': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_start'], city)}",
'neighbour_end': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_end'], city)}",
'shared_points': 1
}
@ -142,6 +151,10 @@ class GeometryHelper:
'line_end': neighbour_info['line_end'],
'neighbour_line_start': (coordinate[0], coordinate[1]),
'neighbour_line_end': (next_coordinate[0], next_coordinate[1]),
'neighbour_start': f"{GeometryHelper.coordinate_to_map_point(coordinate, city)}",
'neighbour_end': f"{GeometryHelper.coordinate_to_map_point(next_coordinate, city)}",
'coordinate_start': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_start'], city)}",
'coordinate_end': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_end'], city)}",
'shared_points': 1
}
@ -154,6 +167,8 @@ class GeometryHelper:
elif building not in neighbour.neighbours:
neighbour.neighbours.append(building)
line += 1
if plot:
img.show()
return lines_information
@staticmethod

108
hub/imports/geometry/geojson.py
View File

@ -6,13 +6,15 @@ Project Coder Guillermo Gutierrez Guillermo.GutierrezMorote@concordia.ca
"""
import json
import numpy as np
import trimesh.creation
from pyproj import Transformer
from shapely.geometry import Polygon as ShapelyPolygon
import hub.helpers.constants as cte
from hub.imports.geometry.helpers.geometry_helper import GeometryHelper
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
@ -62,9 +64,11 @@ class Geojson:
surfaces = []
buildings = []
for zone, surface_coordinates in enumerate(surfaces_coordinates):
points = GeometryHelper.points_from_string(GeometryHelper.remove_last_point_from_string(surface_coordinates))
points = igh.points_from_string(igh.remove_last_point_from_string(surface_coordinates))
# geojson provides the roofs, need to be transform into grounds
points = igh.invert_points(points)
polygon = Polygon(points)
surfaces.append(Surface(polygon, polygon, surface_type=cte.GROUND))
surfaces.append(Surface(polygon, polygon))
buildings.append(Building(f'{name}_zone_{zone}', surfaces, year_of_construction, function))
return buildings
@ -73,22 +77,41 @@ class Geojson:
lod0_buildings = Geojson._create_buildings_lod0(name, year_of_construction, function, surface_coordinates)
surfaces = []
buildings = []
for zone, lod0_building in enumerate(lod0_buildings):
for surface in lod0_building.surfaces:
shapely_polygon = ShapelyPolygon(surface.solid_polygon.coordinates)
if not shapely_polygon.is_valid:
print(surface.solid_polygon.area)
print('error?', name, surface_coordinates)
continue
mesh = trimesh.creation.extrude_polygon(shapely_polygon, height)
for face in mesh.faces:
points = []
for vertex_index in face:
points.append(mesh.vertices[vertex_index])
polygon = Polygon(points)
surface = Surface(polygon, polygon)
surfaces.append(surface)
buildings.append(Building(f'{name}_zone_{zone}', surfaces, year_of_construction, function))
for surface in lod0_building.grounds:
volume = surface.solid_polygon.area * height
surfaces.append(surface)
roof_coordinates = []
# adding a roof means invert the polygon coordinates and change the Z value
for coordinate in surface.solid_polygon.coordinates:
roof_coordinate = np.array([coordinate[0], coordinate[1], height])
# insert the roof rotated already
roof_coordinates.insert(0, roof_coordinate)
polygon = Polygon(roof_coordinates)
roof = Surface(polygon, 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}_zone_{zone}', surfaces, year_of_construction, function)
building.volume = volume
buildings.append(building)
return buildings
def _get_polygons(self, polygons, coordinates):
@ -106,6 +129,50 @@ class Geojson:
polygons.append(transformed_coordinates.lstrip(' '))
return polygons
@staticmethod
def _find_wall(line_1, line_2):
for i in range(0, 2):
point_1 = line_1[i]
point_2 = line_2[i]
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():
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)
# todo: erase when we have no triangulation
if len(ground_line) < 2:
continue
# todo: erase down to here
for entry in building_mapped:
if building_mapped[entry]['shared_points'] <= 5:
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:
"""
@ -115,6 +182,7 @@ class Geojson:
missing_functions = []
buildings = []
building_id = 0
lod = 1
for feature in self._geojson['features']:
extrusion_height = 0
if self._extrusion_height_field is not None:
@ -140,7 +208,6 @@ class Geojson:
building_name = f'building_{building_id}'
building_id += 1
polygons = []
lod = 1
for part, coordinates in enumerate(geometry['coordinates']):
polygons = self._get_polygons(polygons, coordinates)
for zone, polygon in enumerate(polygons):
@ -163,6 +230,9 @@ class Geojson:
for building in buildings:
self._city.add_city_object(building)
self._city.level_of_detail.geometry = lod
if lod == 1:
lines_information = GeometryHelper.city_mapping(self._city)
self._store_shared_percentage_to_walls(self._city, lines_information)
if len(missing_functions) > 0:
print(f'There are unknown functions {missing_functions}')
return self._city

7
hub/imports/geometry/helpers/geometry_helper.py
View File

@ -45,3 +45,10 @@ class GeometryHelper:
array = points.split(' ')
res = " "
return res.join(array[0:len(array) - 3])
@staticmethod
def invert_points(points):
res = []
for point in points:
res.insert(0,point)
return res

1
hub/requirements.txt
View File

@ -23,3 +23,4 @@ shapely
geopandas
triangle
psycopg2-binary
PIL

18
hub/unittests/test_geometry_factory.py
View File

@ -4,13 +4,11 @@ SPDX - License - Identifier: LGPL - 3.0 - or -later
Copyright © 2022 Concordia CERC group
Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
"""
import datetime
from pathlib import Path
from unittest import TestCase
from hub.helpers.geometry_helper import GeometryHelper
from numpy import inf
import hub.exports.exports_factory
from hub.imports.construction_factory import ConstructionFactory
from hub.imports.geometry_factory import GeometryFactory
@ -118,7 +116,6 @@ class TestGeometryFactory(TestCase):
city = self._get_city(file, 'rhino')
self.assertIsNotNone(city, 'city is none')
self.assertTrue(len(city.buildings) == 36)
i = 0
def test_import_obj(self):
"""
@ -142,7 +139,7 @@ class TestGeometryFactory(TestCase):
function_field='CODE_UTILI')
hub.exports.exports_factory.ExportsFactory('obj', city, self._output_path).export()
self.assertEqual(207, len(city.buildings), 'wrong number of buildings')
self.assertEqual(195, len(city.buildings), 'wrong number of buildings')
self._check_buildings(city)
def test_map_neighbours(self):
@ -150,13 +147,21 @@ class TestGeometryFactory(TestCase):
Test neighbours map creation
"""
file = 'neighbours.geojson'
city = self._get_city(file, 'geojson',
year_of_construction_field='ANNEE_CONS',
function_field='LIBELLE_UT')
info_lod0 = GeometryHelper.city_mapping(city, plot=False)
city = self._get_city(file, 'geojson',
height_field='citygml_me',
year_of_construction_field='ANNEE_CONS',
function_field='LIBELLE_UT')
print(GeometryHelper.city_mapping(city))
info_lod1 = GeometryHelper.city_mapping(city, plot=False)
hub.exports.exports_factory.ExportsFactory('obj', city, self._output_path).export()
self.assertEqual(info_lod0, info_lod1)
for building in city.buildings:
self.assertEqual(2, len(building.neighbours))
print(building.volume)
self.assertEqual('2_part_0_zone_0',city.city_object('1_part_0_zone_0').neighbours[0].name)
self.assertEqual('3_part_0_zone_0',city.city_object('1_part_0_zone_0').neighbours[1].name)
@ -164,4 +169,3 @@ class TestGeometryFactory(TestCase):
self.assertEqual('3_part_0_zone_0',city.city_object('2_part_0_zone_0').neighbours[1].name)
self.assertEqual('1_part_0_zone_0', city.city_object('3_part_0_zone_0').neighbours[0].name)
self.assertEqual('2_part_0_zone_0', city.city_object('3_part_0_zone_0').neighbours[1].name)

66
hub/unittests/tests_data/neighbours.geojson
View File

@ -12,18 +12,18 @@
-73.580414175680588,
45.497641136608358
],
[
-73.581414175680588,
45.497641136608358
],
[
-73.581414175680588,
45.498641136608358
],
[
-73.580414175680588,
45.498641136608358
],
[
-73.581414175680588,
45.498641136608358
],
[
-73.581414175680588,
45.497641136608358
],
[
-73.580414175680588,
45.497641136608358
@ -204,19 +204,20 @@
[
-73.581414175680588,
45.497641136608358
]
,
[
-73.581414175680588,
45.498441136608358
],
[
-73.582214175680588,
45.498441136608358
],
[
-73.582214175680588,
45.497641136608358
],
[
-73.582214175680588,
45.498441136608358
],
[
-73.581414175680588,
45.498441136608358
],
[
-73.581414175680588,
45.497641136608358
@ -399,31 +400,30 @@
-73.581914175680588,
45.498441136608358
],
[
-73.581914175680588,
45.499641136608358
],
[
-73.580914175680588,
45.499641136608358
],
[
-73.580914175680588,
45.498641136608358
],
[
-73.581414175680588,
45.498641136608358
],
[
-73.581414175680588,
45.498441136608358
],
[
-73.581414175680588,
45.498641136608358
],
[
-73.580914175680588,
45.498641136608358
],
[
-73.580914175680588,
45.499641136608358
],
[
-73.581914175680588,
45.499641136608358
],
[
-73.581914175680588,
45.498441136608358
]
]
]
},