480 lines
14 KiB
Python
480 lines
14 KiB
Python
"""
|
|
Surface module
|
|
SPDX - License - Identifier: LGPL - 3.0 - or -later
|
|
Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
|
|
contributors: Pilar Monsalvete pilar_monsalvete@yahoo.es
|
|
"""
|
|
from __future__ import annotations
|
|
from typing import Union
|
|
|
|
import numpy as np
|
|
import pyny3d.geoms as pn
|
|
import math
|
|
|
|
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 = dict()
|
|
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
|
|
"""
|
|
# New method to calculate area
|
|
if self._area is None:
|
|
if len(self.points) < 3:
|
|
print('Warning: the area of a line or point cannot be calculated. Area = 0')
|
|
return 0
|
|
alpha = 0
|
|
vec_1 = self.points[1] - self.points[0]
|
|
for i in range(2, len(self.points)):
|
|
vec_2 = self.points[i] - self.points[0]
|
|
alpha += self.angle_between_vectors(vec_1, vec_2)
|
|
if alpha == 0:
|
|
print('Warning: the area of a line or point cannot be calculated. Area = 0')
|
|
return 0
|
|
horizontal_points = self.rotate_surface_to_horizontal(self)
|
|
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)
|
|
return self._area
|
|
|
|
@staticmethod
|
|
def rotate_surface_to_horizontal(surface):
|
|
z_vector = [0, 0, 1]
|
|
normal_vector = surface.normal
|
|
horizontal_points = []
|
|
x = normal_vector[0]
|
|
y = normal_vector[1]
|
|
|
|
if x == 0 and y == 0:
|
|
# Already horizontal
|
|
for point in surface.points:
|
|
horizontal_points.append([point[0], point[1], 0])
|
|
else:
|
|
alpha = surface.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, math.cos(alpha), -math.sin(alpha)],
|
|
[0, math.sin(alpha), math.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:
|
|
print('Error processing rotation base matrix')
|
|
else:
|
|
for point in surface.points:
|
|
new_point = np.matmul(rotation_base_matrix, point)
|
|
horizontal_points.append(new_point)
|
|
return horizontal_points
|
|
|
|
@staticmethod
|
|
def angle_between_vectors(vec_1, vec_2):
|
|
alpha = math.acos(np.dot(vec_1, vec_2) / np.linalg.norm(vec_1) / np.linalg.norm(vec_2))
|
|
if np.linalg.norm(vec_1) == 0 or np.linalg.norm(vec_2) == 0:
|
|
print('ERROR')
|
|
# todo: rise error
|
|
return alpha
|
|
|
|
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
|
|
|
|
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):
|
|
try:
|
|
intersection_area = self.intersect(surface).area
|
|
except ValueError:
|
|
intersection_area = 0
|
|
self.add_shared(surface, intersection_area)
|
|
surface.add_shared(self, intersection_area)
|
|
|
|
@property
|
|
def global_irradiance(self) -> dict:
|
|
"""
|
|
global irradiance on surface in Wh/m2
|
|
:return: dict{DataFrame(float)}
|
|
"""
|
|
return self._global_irradiance
|
|
|
|
@global_irradiance.setter
|
|
def global_irradiance(self, value):
|
|
"""
|
|
global irradiance on surface in Wh/m2
|
|
:param value: dict{DataFrame(float)}
|
|
"""
|
|
self._global_irradiance = 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:
|
|
shapely = self.shapely
|
|
if shapely is not None:
|
|
self._projected_surface = self._polygon_to_surface(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('Warning: intersecting surfaces', err)
|
|
return None
|
|
|
|
@property
|
|
def convex(self):
|
|
return pn.Polygon.is_convex(self.polygon.points)
|