+"""
+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
+
+
+[docs]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)
+
+
[docs] 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
+
+
[docs] 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))
+
+
[docs] 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
+
+
[docs] 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
+