Merge branch 'changing_triangulate' into 'master'

Changing triangulate

See merge request Guille/hub!53
This commit is contained in:
Guillermo Gutierrez Morote 2023-02-13 04:36:52 +00:00
commit 794b9e1b85
12 changed files with 102 additions and 844 deletions

View File

@ -9,9 +9,13 @@ from __future__ import annotations
import math
import sys
from typing import List
from hub.hub_logger import logger
import numpy as np
from trimesh import Trimesh
import trimesh.intersections
import trimesh.creation
import trimesh.geometry
from shapely.geometry.polygon import Polygon as shapley_polygon
from hub.city_model_structure.attributes.plane import Plane
from hub.city_model_structure.attributes.point import Point
@ -22,6 +26,7 @@ class Polygon:
"""
Polygon class
"""
# todo: review with @Guille: Points, Coordinates, Vertices, Faces
def __init__(self, coordinates):
self._area = None
@ -66,20 +71,6 @@ class Polygon:
"""
return self._coordinates
@staticmethod
def _module(vector):
x2 = vector[0] ** 2
y2 = vector[1] ** 2
z2 = vector[2] ** 2
return math.sqrt(x2+y2+z2)
@staticmethod
def _scalar_product(vector_0, vector_1):
x = vector_0[0] * vector_1[0]
y = vector_0[1] * vector_1[1]
z = vector_0[2] * vector_1[2]
return x+y+z
def contains_point(self, point):
"""
Determines if the given point is contained by the current polygon
@ -98,9 +89,9 @@ class Polygon:
vector_1[0] = vector_1[0] - point.coordinates[0]
vector_1[1] = vector_1[1] - point.coordinates[1]
vector_1[2] = vector_1[2] - point.coordinates[2]
module = Polygon._module(vector_0) * Polygon._module(vector_1)
module = np.linalg.norm(vector_0) * np.linalg.norm(vector_1)
scalar_product = Polygon._scalar_product(vector_0, vector_1)
scalar_product = np.dot(vector_0, vector_1)
angle = np.pi/2
if module != 0:
angle = abs(np.arcsin(scalar_product / module))
@ -150,69 +141,17 @@ class Polygon:
Get surface area in square meters
:return: float
"""
# New method to calculate area
if self._area is None:
if len(self.points) < 3:
sys.stderr.write('Warning: the area of a line or point cannot be calculated 1. Area = 0\n')
return 0
alpha = 0
vec_1 = self.points[1].coordinates - self.points[0].coordinates
for i in range(2, len(self.points)):
vec_2 = self.points[i].coordinates - self.points[0].coordinates
alpha += self._angle_between_vectors(vec_1, vec_2)
if alpha == 0:
sys.stderr.write('Warning: the area of a line or point cannot be calculated 2. Area = 0\n')
return 0
horizontal_points = self._points_rotated_to_horizontal
area = 0
for i in range(0, len(horizontal_points) - 1):
point = horizontal_points[i]
next_point = horizontal_points[i + 1]
area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
next_point = horizontal_points[0]
point = horizontal_points[len(horizontal_points) - 1]
area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
self._area = abs(area)
self._area = 0
for triangle in self.triangles:
ab = np.zeros(3)
ac = np.zeros(3)
for i in range(0, 3):
ab[i] = triangle.coordinates[1][i] - triangle.coordinates[0][i]
ac[i] = triangle.coordinates[2][i] - triangle.coordinates[0][i]
self._area += np.linalg.norm(np.cross(ab, ac)) / 2
return self._area
@property
def _points_rotated_to_horizontal(self):
"""
polygon points rotated to horizontal
:return: [float]
"""
z_vector = [0, 0, 1]
normal_vector = self.normal
horizontal_points = []
x = normal_vector[0]
y = normal_vector[1]
if x == 0 and y == 0:
# Already horizontal
for point in self.points:
horizontal_points.append([point.coordinates[0], point.coordinates[1], 0])
else:
alpha = self._angle_between_vectors(normal_vector, z_vector)
rotation_line = np.cross(normal_vector, z_vector)
third_axis = np.cross(normal_vector, rotation_line)
w_1 = rotation_line / np.linalg.norm(rotation_line)
w_2 = normal_vector
w_3 = third_axis / np.linalg.norm(third_axis)
rotation_matrix = np.array([[1, 0, 0],
[0, np.cos(alpha), -np.sin(alpha)],
[0, np.sin(alpha), np.cos(alpha)]])
base_matrix = np.array([w_1, w_2, w_3])
rotation_base_matrix = np.matmul(base_matrix.transpose(), rotation_matrix.transpose())
rotation_base_matrix = np.matmul(rotation_base_matrix, base_matrix)
if rotation_base_matrix is None:
sys.stderr.write('Warning: rotation base matrix returned None\n')
else:
for point in self.points:
new_point = np.matmul(rotation_base_matrix, point.coordinates)
horizontal_points.append(new_point)
return horizontal_points
@property
def normal(self) -> np.ndarray:
"""
@ -275,284 +214,75 @@ class Polygon:
return alpha
return -alpha
def triangulate(self) -> List[Polygon]:
"""
Triangulates a polygon following the ear clipping methodology
:return: list[triangles]
"""
# todo: review triangulate_polygon in
# https://github.com/mikedh/trimesh/blob/dad11126742e140ef46ba12f8cb8643c83356467/trimesh/creation.py#L415,
# it had a problem with a class called 'triangle', but, if solved,
# it could be a very good substitute of this method
# this method is very dirty and has an infinite loop solved with a counter!!
@staticmethod
def triangle_mesh(vertices, normal):
min_x = 1e16
min_y = 1e16
min_z = 1e16
for vertex in vertices:
if vertex[0] < min_x:
min_x = vertex[0]
if vertex[1] < min_y:
min_y = vertex[1]
if vertex[2] < min_z:
min_z = vertex[2]
new_vertices = []
for vertex in vertices:
vertex = [vertex[0]-min_x, vertex[1]-min_y, vertex[2]-min_z]
new_vertices.append(vertex)
transformation_matrix = trimesh.geometry.plane_transform(origin=new_vertices[0], normal=normal)
coordinates = []
for vertex in vertices:
transformed_vertex = [vertex[0]-min_x, vertex[1]-min_y, vertex[2]-min_z, 1]
transformed_vertex = np.dot(transformation_matrix, transformed_vertex)
coordinate = [transformed_vertex[0], transformed_vertex[1]]
coordinates.append(coordinate)
polygon = shapley_polygon(coordinates)
try:
vertices_2d, faces = trimesh.creation.triangulate_polygon(polygon, engine='triangle')
mesh = Trimesh(vertices=vertices, faces=faces)
# check orientation
normal_sum = 0
for i in range(0, 3):
normal_sum += normal[i] + mesh.face_normals[0][i]
if abs(normal_sum) <= 1E-10:
new_faces = []
for face in faces:
new_face = []
for i in range(0, len(face)):
new_face.append(face[len(face)-i-1])
new_faces.append(new_face)
mesh = Trimesh(vertices=vertices, faces=new_faces)
return mesh
except ValueError:
logger.error(f'Not able to triangulate polygon\n')
sys.stderr.write(f'Not able to triangulate polygon\n')
_vertices = [[0, 0, 0], [0, 0, 1], [0, 1, 0]]
_faces = [[0, 1, 2]]
return Trimesh(vertices=_vertices, faces=_faces)
@property
def triangles(self) -> List[Polygon]:
if self._triangles is None:
points_list = self.points_list
normal = self.normal
if np.linalg.norm(normal) == 0:
sys.stderr.write('Not able to triangulate polygon\n')
return [self]
# are points concave or convex?
total_points_list, concave_points, convex_points = self._starting_lists(points_list, normal)
# list of ears
ears = []
j = 0
while (len(concave_points) > 3 or len(convex_points) != 0) and j < 100:
j += 1
for i in range(0, len(concave_points)):
ear = self._triangle(points_list, total_points_list, concave_points[i])
rest_points = []
for points in total_points_list:
rest_points.append(list(self.coordinates[points]))
if self._is_ear(ear, rest_points):
ears.append(ear)
point_to_remove = concave_points[i]
previous_point_in_list, next_point_in_list = self._enveloping_points(point_to_remove,
total_points_list)
total_points_list.remove(point_to_remove)
concave_points.remove(point_to_remove)
# Was any of the adjacent points convex? -> check if changed status to concave
for convex_point in convex_points:
if convex_point == previous_point_in_list:
concave_points, convex_points, end_loop = self._if_concave_change_status(normal,
points_list,
convex_point,
total_points_list,
concave_points,
convex_points,
previous_point_in_list)
if end_loop:
break
continue
if convex_point == next_point_in_list:
concave_points, convex_points, end_loop = self._if_concave_change_status(normal,
points_list,
convex_point,
total_points_list,
concave_points,
convex_points,
next_point_in_list)
if end_loop:
break
continue
break
if len(total_points_list) <= 3 and len(convex_points) > 0:
sys.stderr.write('Not able to triangulate polygon\n')
return [self]
if j >= 100:
sys.stderr.write('Not able to triangulate polygon\n')
return [self]
last_ear = self._triangle(points_list, total_points_list, concave_points[1])
ears.append(last_ear)
self._triangles = ears
self._triangles = []
_mesh = self.triangle_mesh(self.coordinates, self.normal)
for face in _mesh.faces:
points = []
for vertex in face:
points.append(self.coordinates[vertex])
polygon = Polygon(points)
self._triangles.append(polygon)
return self._triangles
@staticmethod
def _starting_lists(points_list, normal) -> [List[float], List[float], List[float]]:
"""
creates the list of vertices (points) that define the polygon (total_points_list), together with other two lists
separating points between convex and concave
:param points_list: points_list
:param normal: normal
:return: list[point], list[point], list[point]
"""
concave_points = []
convex_points = []
# lists of concave and convex points
# case 1: first point
point = points_list[0:3]
previous_point = points_list[len(points_list) - 3:]
next_point = points_list[3:6]
index = 0
total_points_list = [index]
if Polygon._point_is_concave(normal, point, previous_point, next_point):
concave_points.append(index)
else:
convex_points.append(index)
# case 2: all points except first and last
for i in range(0, int((len(points_list) - 6) / 3)):
point = points_list[(i + 1) * 3:(i + 2) * 3]
previous_point = points_list[i * 3:(i + 1) * 3]
next_point = points_list[(i + 2) * 3:(i + 3) * 3]
index = i + 1
total_points_list.append(index)
if Polygon._point_is_concave(normal, point, previous_point, next_point):
concave_points.append(index)
else:
convex_points.append(index)
# case 3: last point
point = points_list[len(points_list) - 3:]
previous_point = points_list[len(points_list) - 6:len(points_list) - 3]
next_point = points_list[0:3]
index = int(len(points_list) / 3) - 1
total_points_list.append(index)
if Polygon._point_is_concave(normal, point, previous_point, next_point):
concave_points.append(index)
else:
convex_points.append(index)
return total_points_list, concave_points, convex_points
@staticmethod
def _triangle(points_list, total_points_list, point_position) -> Polygon:
"""
creates a triangular polygon out of three points
:param points_list: points_list
:param total_points_list: [point]
:param point_position: int
:return: polygon
"""
index = point_position * 3
previous_point_index, next_point_index = Polygon._enveloping_points_indices(point_position, total_points_list)
points = points_list[previous_point_index:previous_point_index + 3]
points = np.append(points, points_list[index:index + 3])
points = np.append(points, points_list[next_point_index:next_point_index + 3])
rows = points.size // 3
points = points.reshape(rows, 3)
triangle = Polygon(points)
return triangle
@staticmethod
def _enveloping_points_indices(point_position, total_points_list):
"""
due to the fact that the lists are not circular, a method to find the previous and next points
of an specific one is needed
:param point_position: int
:param total_points_list: [point]
:return: int, int
"""
previous_point_index = None
next_point_index = None
if point_position == total_points_list[0]:
previous_point_index = total_points_list[len(total_points_list) - 1] * 3
next_point_index = total_points_list[1] * 3
if point_position == total_points_list[len(total_points_list) - 1]:
previous_point_index = total_points_list[len(total_points_list) - 2] * 3
next_point_index = total_points_list[0] * 3
for i in range(1, len(total_points_list) - 1):
if point_position == total_points_list[i]:
previous_point_index = total_points_list[i - 1] * 3
next_point_index = total_points_list[i + 1] * 3
return previous_point_index, next_point_index
@staticmethod
def _enveloping_points(point_to_remove, total_points_list):
"""
due to the fact that the lists are not circular, a method to find the previous and next points
of an specific one is needed
:param point_to_remove: point
:param total_points_list: [point]
:return: point, point
"""
index = total_points_list.index(point_to_remove)
if index == 0:
previous_point_in_list = total_points_list[len(total_points_list) - 1]
next_point_in_list = total_points_list[1]
elif index == len(total_points_list) - 1:
previous_point_in_list = total_points_list[len(total_points_list) - 2]
next_point_in_list = total_points_list[0]
else:
previous_point_in_list = total_points_list[index - 1]
next_point_in_list = total_points_list[index + 1]
return previous_point_in_list, next_point_in_list
@staticmethod
def _is_ear(ear, points) -> bool:
"""
finds whether a triangle is an ear of the polygon
:param ear: polygon
:param points: [point]
:return: boolean
"""
area_ear = ear.area
for point in points:
area_points = 0
point_is_not_vertex = True
for i in range(0, 3):
if abs(np.linalg.norm(point) - np.linalg.norm(ear.coordinates[i])) < 0.0001:
point_is_not_vertex = False
break
if point_is_not_vertex:
for i in range(0, 3):
if i != 2:
new_points = ear.coordinates[i][:]
new_points = np.append(new_points, ear.coordinates[i + 1][:])
new_points = np.append(new_points, point[:])
else:
new_points = ear.coordinates[i][:]
new_points = np.append(new_points, point[:])
new_points = np.append(new_points, ear.coordinates[0][:])
rows = new_points.size // 3
new_points = new_points.reshape(rows, 3)
new_triangle = Polygon(new_points)
area_points += new_triangle.area
if abs(area_points - area_ear) < 1e-6:
# point_inside_ear = True
return False
return True
@staticmethod
def _if_concave_change_status(normal, points_list, convex_point, total_points_list,
concave_points, convex_points, point_in_list) -> [List[float], List[float], bool]:
"""
checks whether an convex specific point change its status to concave after removing one ear in the polygon
returning the new convex and concave points lists together with a flag advising that the list of total points
already 3 and, therefore, the triangulation must be finished.
:param normal: normal
:param points_list: points_list
:param convex_point: int
:param total_points_list: [point]
:param concave_points: [point]
:param convex_points: [point]
:param point_in_list: int
:return: list[points], list[points], boolean
"""
end_loop = False
point = points_list[point_in_list * 3:(point_in_list + 1) * 3]
pointer = total_points_list.index(point_in_list) - 1
if pointer < 0:
pointer = len(total_points_list) - 1
previous_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3]
pointer = total_points_list.index(point_in_list) + 1
if pointer >= len(total_points_list):
pointer = 0
next_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3]
if Polygon._point_is_concave(normal, point, previous_point, next_point):
if concave_points[0] > convex_point:
concave_points.insert(0, convex_point)
elif concave_points[len(concave_points) - 1] < convex_point:
concave_points.append(convex_point)
else:
for point_index in range(0, len(concave_points) - 1):
if concave_points[point_index] < convex_point < concave_points[point_index + 1]:
concave_points.insert(point_index + 1, convex_point)
convex_points.remove(convex_point)
end_loop = True
return concave_points, convex_points, end_loop
@staticmethod
def _point_is_concave(normal, point, previous_point, next_point) -> bool:
"""
returns whether a point is concave
:param normal: normal
:param point: point
:param previous_point: point
:param next_point: point
:return: boolean
"""
is_concave = False
accepted_error = 0.1
points = np.append(previous_point, point)
points = np.append(points, next_point)
rows = points.size // 3
points = points.reshape(rows, 3)
triangle = Polygon(points)
error_sum = 0
for i in range(0, len(normal)):
error_sum += triangle.normal[i] - normal[i]
if np.abs(error_sum) < accepted_error:
is_concave = True
return is_concave
@staticmethod
def _angle_between_vectors(vec_1, vec_2):
"""
@ -652,12 +382,12 @@ class Polygon:
@property
def vertices(self) -> np.ndarray:
"""
Get polyhedron vertices
Get polygon vertices
:return: np.ndarray(int)
"""
if self._vertices is None:
vertices, self._vertices = [], []
_ = [vertices.extend(s.coordinates) for s in self.triangulate()]
_ = [vertices.extend(s.coordinates) for s in self.triangles]
for vertex_1 in vertices:
found = False
for vertex_2 in self._vertices:
@ -677,17 +407,17 @@ class Polygon:
@property
def faces(self) -> List[List[int]]:
"""
Get polyhedron triangular faces
Get polygon triangular faces
:return: [face]
"""
if self._faces is None:
self._faces = []
for polygon in self.triangulate():
for polygon in self.triangles:
face = []
points = polygon.coordinates
if len(points) != 3:
sub_polygons = polygon.triangulate()
sub_polygons = polygon.triangles
# todo: I modified this! To be checked @Guille
if len(sub_polygons) >= 1:
for sub_polygon in sub_polygons:

View File

@ -91,7 +91,7 @@ class Polyhedron:
face = []
points = polygon.coordinates
if len(points) != 3:
sub_polygons = polygon.triangulate()
sub_polygons = polygon.triangles
# todo: I modified this! To be checked @Guille
if len(sub_polygons) >= 1:
for sub_polygon in sub_polygons:

View File

@ -6,8 +6,10 @@ Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
Code contributors: Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
"""
import sys
from typing import List, Union
import numpy as np
from hub.hub_logger import logger
import hub.helpers.constants as cte
from hub.city_model_structure.building_demand.surface import Surface
from hub.city_model_structure.city_object import CityObject
@ -46,20 +48,31 @@ class Building(CityObject):
self._roofs = []
self._walls = []
self._internal_walls = []
self._ground_walls = []
self._attic_floors = []
self._interior_slabs = []
for surface_id, surface in enumerate(self.surfaces):
self._min_x = min(self._min_x, surface.lower_corner[0])
self._min_y = min(self._min_y, surface.lower_corner[1])
self._min_z = min(self._min_z, surface.lower_corner[2])
surface.id = surface_id
# todo: consider all type of surfaces, not only these four
if surface.type == cte.GROUND:
self._grounds.append(surface)
elif surface.type == cte.WALL:
self._walls.append(surface)
elif surface.type == cte.ROOF:
self._roofs.append(surface)
else:
elif surface.type == cte.INTERIOR_WALL:
self._internal_walls.append(surface)
elif surface.type == cte.GROUND_WALL:
self._ground_walls.append(surface)
elif surface.type == cte.ATTIC_FLOOR:
self._attic_floors.append(surface)
elif surface.type == cte.INTERIOR_SLAB:
self._interior_slabs.append(surface)
else:
logger.error(f'Building {self.name} [alias {self.alias}] has an unexpected surface type {surface.type}.\n')
sys.stderr.write(f'Building {self.name} [alias {self.alias}] has an unexpected surface type {surface.type}.\n')
@property
def shell(self) -> Polyhedron:

View File

@ -6,7 +6,6 @@ Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
"""
from typing import Union, List
from hub.city_model_structure.attributes.schedule import Schedule
from hub.city_model_structure.building_demand.occupant import Occupant
class Occupancy:
@ -106,19 +105,3 @@ class Occupancy:
:param value: [Schedule]
"""
self._occupancy_schedules = value
@property
def occupants(self) -> Union[None, List[Occupant]]:
"""
Get list of occupants
:return: None or List of Occupant
"""
return self._occupants
@occupants.setter
def occupants(self, value):
"""
Set list of occupants
:param value: [Occupant]
"""
self._occupants = value

View File

@ -1,145 +0,0 @@
"""
Occupant module
SPDX - License - Identifier: LGPL - 3.0 - or -later
Copyright © 2022 Concordia CERC group
Project Coder Sanam Dabirian sanam.dabirian@mail.concordia.ca
Code contributors: Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
"""
class Occupant:
"""
Occupant class
"""
def __init__(self):
"""
Constructor
"""
self._heat_dissipation = None
self._occupancy_rate = None
self._occupant_type = None
self._arrival_time = None
self._departure_time = None
self._break_time = None
self._day_of_week = None
self._pd_of_meetings_duration = None
@property
def heat_dissipation(self):
"""
Get heat dissipation of occupants in W/person
:return: float
"""
return self._heat_dissipation
@heat_dissipation.setter
def heat_dissipation(self, value):
"""
Set heat dissipation of occupants in W/person
:param value: float
"""
self._heat_dissipation = float(value)
@property
def occupancy_rate(self):
"""
Get rate of schedules
:return: float
"""
return self._occupancy_rate
@occupancy_rate.setter
def occupancy_rate(self, value):
"""
Set rate of schedules
:param value: float
"""
self._occupancy_rate = float(value)
@property
def occupant_type(self):
"""
Get type of schedules
:return: str
"""
return self._occupant_type
@occupant_type.setter
def occupant_type(self, value):
"""
Set type of schedules
:param value: float
"""
self._occupant_type = float(value)
@property
def arrival_time(self):
"""
Get the arrival time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
:return: time
"""
return self._arrival_time
@arrival_time.setter
def arrival_time(self, value):
"""
Set the arrival time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
:param value: time
"""
self._arrival_time = value
@property
def departure_time(self):
"""
Get the departure time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
:return: time
"""
return self._departure_time
@departure_time.setter
def departure_time(self, value):
"""
Set the departure time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
:param value: str
"""
self._departure_time = value
@property
def break_time(self):
"""
Get the lunch or break time of the occupant (for office building) in UTC with format ????
:return: break time
"""
# todo @Sanam: define this format, is it the starting time? is it a list with both, starting and ending time?
return self._break_time
@property
def day_of_week(self):
"""
Get the day of the week (MON, TUE, WED, THU, FRI, SAT, SUN)
:return: str
"""
# todo @Sanam: is this a property or should it be a function
# to get the day of the week of an specific day of the year?
return self._day_of_week
@property
def pd_of_meetings_duration(self):
"""
Get the probability distribution of the meeting duration
:return: ??
"""
# todo @Sanam: what format are you expecting here??
return self._pd_of_meetings_duration
@pd_of_meetings_duration.setter
def pd_of_meetings_duration(self, value):
"""
Get the probability distribution of the meeting duration
:param value: ??
:return:
"""
# todo @Sanam: what format are you expecting here??
self._pd_of_meetings_duration = value

View File

@ -72,14 +72,6 @@ class Surface:
if value is not None:
self._id = str(value)
# todo: implement share surfaces
@property
def share_surfaces(self):
"""
Raises not implemented error
"""
raise NotImplementedError
def _max_coord(self, axis):
if axis == 'x':
axis = 0
@ -163,11 +155,11 @@ class Surface:
@property
def type(self):
"""
Get surface type Ground, Wall or Roof
Get surface type Ground, Ground wall, Wall, Attic floor, Interior slab, Interior wall, Roof or Virtual internal
If the geometrical LoD is lower than 4,
the surfaces' types are not defined in the importer and can only be Ground, Wall or Roof
:return: str
"""
# todo: there are more types: internal wall, internal floor... this method must be redefined
if self._type is None:
grad = np.rad2deg(self.inclination)
if grad >= 170:
@ -304,7 +296,7 @@ class Surface:
:return: Surface, Surface, Any
"""
# todo: check return types
# todo: recheck this method for LoD3 (windows)
# recheck this method for LoD3 (windows)
origin = Point([0, 0, z])
normal = np.array([0, 0, 1])
plane = Plane(normal=normal, origin=origin)

View File

@ -23,7 +23,6 @@ from hub.city_model_structure.iot.station import Station
from hub.city_model_structure.level_of_detail import LevelOfDetail
from hub.city_model_structure.machine import Machine
from hub.city_model_structure.parts_consisting_building import PartsConsistingBuilding
from hub.city_model_structure.subway_entrance import SubwayEntrance
from hub.helpers.geometry_helper import GeometryHelper
from hub.helpers.location import Location
from hub.city_model_structure.energy_system import EnergySystem
@ -40,9 +39,7 @@ class City:
self._lower_corner = lower_corner
self._upper_corner = upper_corner
self._buildings = None
self._subway_entrances = None
self._srs_name = srs_name
# todo: right now extracted at city level, in the future should be extracted also at building level if exist
self._location = None
self._country_code = None
self._climate_reference_city = None
@ -163,9 +160,6 @@ class City:
if self.buildings is not None:
for building in self.buildings:
self._city_objects.append(building)
if self.subway_entrances is not None:
for subway_entrance in self.subway_entrances:
self._city_objects.append(subway_entrance)
if self.energy_systems is not None:
for energy_system in self.energy_systems:
self._city_objects.append(energy_system)
@ -179,14 +173,6 @@ class City:
"""
return self._buildings
@property
def subway_entrances(self) -> Union[List[SubwayEntrance], None]:
"""
Get the subway entrances belonging to the city
:return: a list of subway entrances objects or none
"""
return self._subway_entrances
@property
def lower_corner(self) -> List[float]:
"""
@ -224,10 +210,6 @@ class City:
if self._buildings is None:
self._buildings = []
self._buildings.append(new_city_object)
elif new_city_object.type == 'subway_entrance':
if self._subway_entrances is None:
self._subway_entrances = []
self._subway_entrances.append(new_city_object)
elif new_city_object.type == 'energy_system':
if self._energy_systems is None:
self._energy_systems = []

View File

@ -1,45 +0,0 @@
"""
Subway entrance module
SPDX - License - Identifier: LGPL - 3.0 - or -later
Copyright © 2022 Concordia CERC group
Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
"""
from hub.city_model_structure.city_object import CityObject
class SubwayEntrance(CityObject):
"""
SubwayEntrance(CityObject) class
"""
def __init__(self, name, latitude, longitude):
super().__init__(name, 0)
self._name = name
self._latitude = latitude
self._longitude = longitude
self._type = 'subway_entrance'
@property
def latitude(self):
# todo: to be defined the spacial point and the units
"""
Get latitude
:return: float
"""
return self._latitude
@property
def longitude(self):
# todo: to be defined the spacial point and the units
"""
Get longitude
:return: float
"""
return self._longitude
@property
def name(self):
"""
Get name
:return: str
"""
return self._name

View File

@ -39,59 +39,12 @@ class GeometryHelper:
max_distance = ConfigurationHelper().max_location_distance_for_shared_walls
return GeometryHelper.distance_between_points(location1, location2) < max_distance
def almost_same_area(self, area_1, area_2):
"""
Compare two areas and decides if they are almost equal (absolute error under delta)
:param area_1
:param area_2
:return: Boolean
"""
if area_1 == 0 or area_2 == 0:
return False
delta = math.fabs(area_1 - area_2)
return delta <= self._area_delta
def is_almost_same_surface(self, surface_1, surface_2):
"""
Compare two surfaces and decides if they are almost equal (quadratic error under delta)
:param surface_1: Surface
:param surface_2: Surface
:return: Boolean
"""
# delta is grads an need to be converted into radians
delta = np.rad2deg(self._delta)
difference = (surface_1.inclination - surface_2.inclination) % math.pi
if abs(difference) > delta:
return False
# s1 and s2 are at least almost parallel surfaces
# calculate distance point to plane using all the vertex
# select surface1 value for the point (X,Y,Z) where two of the values are 0
minimum_distance = self._delta + 1
parametric = surface_2.polygon.get_parametric()
normal_2 = surface_2.normal
for point in surface_1.points:
distance = abs(
(point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3])
normal_module = math.sqrt(pow(normal_2[0], 2) + pow(normal_2[1], 2) + pow(normal_2[2], 2))
if normal_module == 0:
continue
distance = distance / normal_module
if distance < minimum_distance:
minimum_distance = distance
if minimum_distance <= self._delta:
break
if minimum_distance > self._delta or surface_1.intersect(surface_2) is None:
return False
return True
@staticmethod
def segment_list_to_trimesh(lines) -> Trimesh:
"""
Transform a list of segments into a Trimesh
"""
# todo: trimesh has a method for this
line_points = [lines[0][0], lines[0][1]]
lines.remove(lines[0])
while len(lines) > 1:
@ -106,7 +59,7 @@ class GeometryHelper:
line_points.append(line[0])
lines.pop(i - 1)
break
polyhedron = Polyhedron(Polygon(line_points).triangulate())
polyhedron = Polyhedron(Polygon(line_points).triangles)
trimesh = Trimesh(polyhedron.vertices, polyhedron.faces)
return trimesh

View File

@ -1,138 +0,0 @@
"""
monthly_to_hourly_demand module
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 calendar as cal
import pandas as pd
from hub.city_model_structure.building_demand.occupant import Occupant
import hub.helpers.constants as cte
class MonthlyToHourlyDemand:
"""
MonthlyToHourlyDemand class
"""
def __init__(self, building, conditioning_seasons):
self._hourly_heating = pd.DataFrame()
self._hourly_cooling = pd.DataFrame()
self._building = building
self._conditioning_seasons = conditioning_seasons
def hourly_heating(self, key):
"""
hourly distribution of the monthly heating of a building
:param key: string
:return: [hourly_heating]
"""
# todo: this method and the insel model have to be reviewed for more than one thermal zone
external_temp = self._building.external_temperature[cte.HOUR]
# todo: review index depending on how the schedules are defined, either 8760 or 24 hours
for usage in self._building.usages:
temp_set = float(usage.heating_setpoint)-3
temp_back = float(usage.heating_setback)-3
# todo: if these are data frames, then they should be called as (Occupancy should be in low case):
# usage.schedules.Occupancy
# self._conditioning_seasons.heating
occupancy = Occupant().get_complete_year_schedule(usage.schedules['Occupancy'])
heating_schedule = self._conditioning_seasons['heating']
hourly_heating = []
i = 0
j = 0
temp_grad_day = []
for month in range(1, 13):
temp_grad_month = 0
month_range = cal.monthrange(2015, month)[1]
for _ in range(1, month_range+1):
external_temp_med = 0
for hour in range(0, 24):
external_temp_med += external_temp[key][i]/24
for hour in range(0, 24):
if external_temp_med < temp_set and heating_schedule[month-1] == 1:
if occupancy[hour] > 0:
hdd = temp_set - external_temp[key][i]
if hdd < 0:
hdd = 0
temp_grad_day.append(hdd)
else:
hdd = temp_back - external_temp[key][i]
if hdd < 0:
hdd = 0
temp_grad_day.append(hdd)
else:
temp_grad_day.append(0)
temp_grad_month += temp_grad_day[i]
i += 1
for _ in range(1, month_range + 1):
for hour in range(0, 24):
monthly_demand = self._building.heating[cte.MONTH][month-1]
if monthly_demand == 'NaN':
monthly_demand = 0
if temp_grad_month == 0:
hourly_demand = 0
else:
hourly_demand = float(monthly_demand)*float(temp_grad_day[j])/float(temp_grad_month)
hourly_heating.append(hourly_demand)
j += 1
self._hourly_heating = pd.DataFrame(data=hourly_heating, columns=['monthly to hourly'])
return self._hourly_heating
def hourly_cooling(self, key):
"""
hourly distribution of the monthly cooling of a building
:param key: string
:return: [hourly_cooling]
"""
# todo: this method and the insel model have to be reviewed for more than one thermal zone
external_temp = self._building.external_temperature[cte.HOUR]
# todo: review index depending on how the schedules are defined, either 8760 or 24 hours
for usage in self._building.usages:
temp_set = float(usage.cooling_setpoint)
temp_back = 100
occupancy = Occupant().get_complete_year_schedule(usage.schedules['Occupancy'])
cooling_schedule = self._conditioning_seasons['cooling']
hourly_cooling = []
i = 0
j = 0
temp_grad_day = []
for month in range(1, 13):
temp_grad_month = 0
month_range = cal.monthrange(2015, month)[1]
for _ in range(1, month_range[1] + 1):
for hour in range(0, 24):
if external_temp[key][i] > temp_set and cooling_schedule[month - 1] == 1:
if occupancy[hour] > 0:
cdd = external_temp[key][i] - temp_set
if cdd < 0:
cdd = 0
temp_grad_day.append(cdd)
else:
cdd = external_temp[key][i] - temp_back
if cdd < 0:
cdd = 0
temp_grad_day.append(cdd)
else:
temp_grad_day.append(0)
temp_grad_month += temp_grad_day[i]
i += 1
for _ in range(1, month_range[1] + 1):
for hour in range(0, 24):
# monthly_demand = self._building.heating[cte.MONTH]['INSEL'][month-1]
monthly_demand = self._building.cooling[cte.MONTH][month - 1]
if monthly_demand == 'NaN':
monthly_demand = 0
if temp_grad_month == 0:
hourly_demand = 0
else:
hourly_demand = float(monthly_demand) * float(temp_grad_day[j]) / float(temp_grad_month)
hourly_cooling.append(hourly_demand)
j += 1
self._hourly_cooling = pd.DataFrame(data=hourly_cooling, columns=['monthly to hourly'])
return self._hourly_cooling

View File

@ -1,58 +0,0 @@
"""
OsmSubway module parses osm files and import the metro location into the city model structure
SPDX - License - Identifier: LGPL - 3.0 - or -later
Copyright © 2022 Concordia CERC group
Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
"""
import sys
import xmltodict
from pyproj import Transformer
from hub.city_model_structure.city import City
from hub.city_model_structure.subway_entrance import SubwayEntrance
class OsmSubway:
"""
Open street map subway
"""
def __init__(self, path):
self._city = None
self._subway_entrances = []
with open(path) as osm:
self._osm = xmltodict.parse(osm.read(), force_list='tag')
for node in self._osm['osm']['node']:
if 'tag' not in node:
continue
for tag in node['tag']:
if '@v' not in tag:
continue
if tag['@v'] == 'subway_entrance':
subway_entrance = SubwayEntrance(node['@id'], node['@lat'], node['@lon'])
self._subway_entrances.append(subway_entrance)
@property
def city(self) -> City:
"""
Get a city with subway entrances
"""
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857")
lower_corner = [sys.float_info.max, sys.float_info.max, 0]
upper_corner = [sys.float_info.min, sys.float_info.min, 0]
x = 0
y = 1
for subway_entrance in self._subway_entrances:
coordinate = transformer.transform(subway_entrance.longitude, subway_entrance.latitude)
if coordinate[x] >= upper_corner[x]:
upper_corner[x] = coordinate[x]
if coordinate[y] >= upper_corner[y]:
upper_corner[y] = coordinate[y]
if coordinate[x] < lower_corner[x]:
lower_corner[x] = coordinate[x]
if coordinate[y] < lower_corner[y]:
lower_corner[y] = coordinate[y]
city = City(lower_corner, upper_corner, 'unknown')
for subway_entrance in self._subway_entrances:
city.add_city_object(subway_entrance)
return city

View File

@ -9,7 +9,6 @@ import geopandas
from hub.city_model_structure.city import City
from hub.imports.geometry.citygml import CityGml
from hub.imports.geometry.obj import Obj
from hub.imports.geometry.osm_subway import OsmSubway
from hub.imports.geometry.rhino import Rhino
from hub.imports.geometry.gpandas import GPandas
from hub.imports.geometry.geojson import Geojson
@ -83,14 +82,6 @@ class GeometryFactory:
self._function_field,
self._function_to_hub).city
@property
def _osm_subway(self) -> City:
"""
Enrich the city by using OpenStreetMap information as data source
:return: City
"""
return OsmSubway(self._path).city
@property
def _rhino(self) -> City:
"""