""" Surface module SPDX - License - Identifier: LGPL - 3.0 - or -later Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca """ from __future__ import annotations from typing import Union import numpy as np import pyny3d.geoms as pn from helpers.geometry_helper import GeometryHelper class Surface: """ Surface class """ def __init__(self, coordinates, surface_type=None, name=None, swr='0.2', remove_last=True, is_projected=False): self._coordinates = coordinates self._type = surface_type self._name = name self._swr = swr self._remove_last = remove_last self._is_projected = is_projected self._geometry_helper = GeometryHelper() self._polygon = None self._ground_polygon = None self._area = None self._points = None self._ground_points = None self._points_list = None self._normal = None self._azimuth = None self._inclination = None self._area_above_ground = None self._area_below_ground = None self._parent = None self._shapely = None self._projected_surface = None self._min_x = None self._min_y = None self._min_z = None self._shared_surfaces = [] self._global_irradiance_hour = np.zeros(8760) self._global_irradiance_month = np.zeros(12) self._ground_coordinates = (self.min_x, self.min_y, self.min_z) def parent(self, parent, surface_id): """ Assign a city object as surface parent and a surface id :param parent: CityObject :param surface_id: str :return: None """ self._parent = parent self._name = str(surface_id) @property def name(self): """ Surface name :return: str """ if self._name is None: raise Exception('surface has no name') return self._name @property def swr(self): """ Get surface short wave reflectance :return: float """ return self._swr @swr.setter def swr(self, value): """ Set surface short wave reflectance :param value: float :return: None """ self._swr = value @property def points(self) -> np.ndarray: """ Surface point matrix :return: np.ndarray """ if self._points is None: self._points = np.fromstring(self._coordinates, dtype=float, sep=' ') self._points = GeometryHelper.to_points_matrix(self._points, self._remove_last) return self._points def _min_coord(self, axis): if axis == 'x': axis = 0 elif axis == 'y': axis = 1 else: axis = 2 min_coordinate = '' for point in self.points: if min_coordinate == '': min_coordinate = point[axis] elif min_coordinate > point[axis]: min_coordinate = point[axis] return min_coordinate @property def min_x(self): """ Surface minimal x value :return: float """ if self._min_x is None: self._min_x = self._min_coord('x') return self._min_x @property def min_y(self): """ Surface minimal y value :return: float """ if self._min_y is None: self._min_y = self._min_coord('y') return self._min_y @property def min_z(self): """ Surface minimal z value :return: float """ if self._min_z is None: self._min_z = self._min_coord('z') return self._min_z @property def ground_points(self) -> np.ndarray: """ Surface grounded points matrix :return: np.ndarray """ if self._ground_points is None: coordinates = '' for point in self.points: x = point[0] - self._ground_coordinates[0] y = point[1] - self._ground_coordinates[1] z = point[2] - self._ground_coordinates[2] if coordinates != '': coordinates = coordinates + ' ' coordinates = coordinates + str(x) + ' ' + str(y) + ' ' + str(z) self._ground_points = np.fromstring(coordinates, dtype=float, sep=' ') self._ground_points = GeometryHelper.to_points_matrix(self._ground_points, False) return self._ground_points @property def points_list(self) -> np.ndarray: """ Surface point list :return: np.ndarray """ if self._points_list is None: s = self.points self._points_list = np.reshape(s, len(s) * 3) return self._points_list @property def polygon(self) -> Union[pn.Polygon, None]: """ Surface polygon :return: None or pyny3d.Polygon """ if self._polygon is None: try: self._polygon = pn.Polygon(self.points) except ValueError: # is not really a polygon but a line so just return none self._polygon = None return self._polygon @property def ground_polygon(self) -> Union[pn.Polygon, None]: """ Surface grounded polygon :return: None or pyny3d.Polygon """ if self._ground_polygon is None: try: self._ground_polygon = pn.Polygon(self.ground_points) except ValueError: # is not really a polygon but a line so just return none self._ground_polygon = None return self._ground_polygon @property def area(self): """ Surface area in square meters :return: float """ if self._area is None: self._area = self.polygon.get_area() return self._area def _is_almost_same_terrain(self, terrain_points, ground_points): equal = 0 for terrain_point in terrain_points: for ground_point in ground_points: if self._geometry.almost_equal(terrain_point, ground_point): equal += 1 return equal == len(terrain_points) @property def _is_terrain(self): for t_points in self._parent.terrains: if len(t_points) == len(self.points) and self._is_almost_same_terrain(t_points, self.points): return True return False @property def area_above_ground(self): """ Surface area above ground in square meters :return: float """ if self._area_above_ground is None: self._area_above_ground = self.area - self.area_below_ground return self._area_above_ground @property def area_below_ground(self): """ Surface area below ground in square meters :return: float """ if self._area_below_ground is None: self._area_below_ground = 0.0 if self._is_terrain: self._area_below_ground = self.area return self._area_below_ground @property def normal(self) -> np.ndarray: """ Surface normal vector :return: np.ndarray """ if self._normal is None: points = self.points cross_product = np.cross(points[1] - points[0], points[2] - points[0]) self._normal = cross_product / np.linalg.norm(cross_product) return self._normal @property def azimuth(self): """ Surface azimuth in radians :return: float """ if self._azimuth is None: normal = self.normal self._azimuth = np.arctan2(normal[1], normal[0]) return self._azimuth @property def inclination(self): """ Surface inclination in radians :return: float """ if self._inclination is None: self._inclination = np.arccos(self.normal[2]) return self._inclination @property def type(self): """ Surface type Ground, Wall or Roof :return: str """ if self._type is None: grad = np.rad2deg(self.inclination) if grad >= 170: self._type = 'Ground' elif 80 <= grad <= 100: self._type = 'Wall' else: self._type = 'Roof' return self._type def add_shared(self, surface, intersection_area): """ Add a given surface and shared area in percent to this surface. :param surface: :param intersection_area: :return: """ percent = intersection_area / self.area self._shared_surfaces.append((percent, surface)) def shared(self, surface): """ Check if given surface share some area with this surface :param surface: Surface :return: None """ if self.type != 'Wall' or surface.type != 'Wall': return if self._geometry_helper.is_almost_same_surface(self, surface): intersection_area = self.intersect(surface).area self.add_shared(surface, intersection_area) surface.add_shared(self, intersection_area) @property def global_irradiance_hour(self): """ Get surface global irradiance hour in Wh/m2 :return: float """ return self._global_irradiance_hour @global_irradiance_hour.setter def global_irradiance_hour(self, value): """ Set surface global irradiance per hour in Wh/m2 :param value: float :return: None """ self._global_irradiance_hour = value @property def global_irradiance_month(self): """ Get surface global irradiance per month in Wh/m2 :return: float """ return self._global_irradiance_month @global_irradiance_month.setter def global_irradiance_month(self, value): """ Set surface global irradiance per month in Wh/m2 :param value: float :return: None """ self._global_irradiance_month = value @property def shapely(self) -> Union[None, pn.Polygon]: """ Surface shapely (Z projection) :return: None or pyny3d.Polygon """ if self.polygon is None: return None if self._shapely is None: self._shapely = self.polygon.get_shapely() return self._shapely @staticmethod def _polygon_to_surface(polygon) -> Surface: coordinates = '' for coordinate in polygon.exterior.coords: if coordinates != '': coordinates = coordinates + ' ' coordinates = coordinates + str(coordinate[0]) + ' ' + str(coordinate[1]) + ' 0.0' return Surface(coordinates, remove_last=False) @property def projection(self) -> Surface: """ Projected surface (Z projection) :return: Surface """ if self._is_projected: return self if self._projected_surface is None: self._projected_surface = self._polygon_to_surface(self.shapely) return self._projected_surface def intersect(self, surface) -> Union[Surface, None]: """ Get the intersection surface, if any, between the given surface and this surface :param surface: Surface :return: None or Surface """ min_x = min(self.min_x, surface.min_x) min_y = min(self.min_y, surface.min_y) min_z = min(self.min_z, surface.min_z) self._ground_coordinates = (min_x, min_y, min_z) surface._ground_coordinates = (min_x, min_y, min_z) origin = (0, 0, 0) azimuth = self.azimuth - (np.pi / 2) while azimuth < 0: azimuth += (np.pi / 2) inclination = self.inclination - np.pi while inclination < 0: inclination += np.pi polygon1 = self.ground_polygon.rotate(azimuth, 'z', origin).rotate(inclination, 'x', origin) polygon2 = surface.ground_polygon.rotate(azimuth, 'z', origin).rotate(inclination, 'x', origin) try: coordinates = '' intersection = pn.Surface([polygon1]).intersect_with(polygon2) if len(intersection) == 0: return None for coordinate in pn.Surface([polygon1]).intersect_with(polygon2)[0]: if coordinates != '': coordinates = coordinates + ' ' coordinates = coordinates + str(coordinate[0]) + ' ' + str(coordinate[1]) + ' 0.0' if coordinates == '': return None intersect_surface = Surface(coordinates, remove_last=False) if intersect_surface.polygon is None: return None return Surface(coordinates, remove_last=False) except Exception as err: print('Error', err) return None