city_retrofit/city_model_structure/surface.py

"""
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 import Geometry


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 = Geometry()
    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 = Geometry.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 = Geometry.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.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