-"""
-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_helper.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
-