city_retrofit/city_model_structure/attributes/surface.py

"""
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 sys

import numpy as np
import pyny3d.geoms as pn
import math

from helpers.geometry_helper import GeometryHelper

# todo: remove pyny3d, seems to not be supported and has some issues


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)
    self._is_planar = None

  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:
        sys.stderr.write('Warning: the area of a line or point cannot be calculated 1. Area = 0\n')
        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 += GeometryHelper.angle_between_vectors(vec_1, vec_2)
      if alpha == 0:
        sys.stderr.write('Warning: the area of a line or point cannot be calculated 2. Area = 0\n')
        return 0
      horizontal_points = self.rotate_surface_to_horizontal
      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

  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
      # todo: IF THE FIRST ONE IS 0, START WITH THE NEXT
      point_origin = points[len(points)-2]
      vector_1 = points[len(points)-1] - point_origin
      vector_2 = points[0] - point_origin
      vector_3 = points[1] - point_origin
      cross_product = np.cross(vector_1, vector_2)
      if np.linalg.norm(cross_product) != 0:
        cross_product = cross_product / np.linalg.norm(cross_product)
        alpha = GeometryHelper.angle_between_vectors(vector_1, vector_2)
      else:
        # todo modify here
        cross_product = [0, 0, 0]
        alpha = 0
      if len(points) == 3:
        return cross_product
      alpha += self._angle(vector_2, vector_3, cross_product)
      for i in range(0, len(points)-4):
        vector_1 = points[i+1] - point_origin
        vector_2 = points[i+2] - point_origin
        alpha += self._angle(vector_1, vector_2, cross_product)
      vector_1 = points[len(points) - 1] - point_origin
      vector_2 = points[0] - point_origin
      if alpha < 0:
        cross_product = np.cross(vector_2, vector_1)
      else:
        cross_product = np.cross(vector_1, vector_2)
      self._normal = cross_product / np.linalg.norm(cross_product)
    return self._normal

  @staticmethod
  def _angle(vector_1, vector_2, cross_product):
    accepted_normal_difference = 0.01
    cross_product_next = np.cross(vector_1, vector_2)
    if np.linalg.norm(cross_product_next) != 0:
      cross_product_next = cross_product_next / np.linalg.norm(cross_product_next)
      alpha = GeometryHelper.angle_between_vectors(vector_1, vector_2)
    else:
      cross_product_next = [0, 0, 0]
      alpha = 0
    delta_normals = 0
    for j in range(0, 3):
      delta_normals += cross_product[j] - cross_product_next[j]
    if np.abs(delta_normals) < accepted_normal_difference:
      return alpha
    else:
      return -alpha

  @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:
      sys.stderr.write('Warning: intersecting surfaces ' + str(err))
      return None

  @property
  def convex(self):
    return pn.Polygon.is_convex(self.polygon.points)

  @property
  def is_planar(self) -> bool:
    if self._is_planar is None:
      self._is_planar = True
      vectors = []
      for i in range(1,len(self.points)):
        vectors.append(self.points[i] - self.points[0])
      for i in range(2, len(vectors)):
        product = np.dot(np.cross(vectors[0], vectors[1]), vectors[i])
        if math.fabs(product) > 1e-4:
          self._is_planar = False
          break
    return self._is_planar

  @property
  def rotate_surface_to_horizontal(self):
    z_vector = [0, 0, 1]
    normal_vector = self.normal
    horizontal_points = []
    x = normal_vector[0]
    y = normal_vector[1]

    if x == 0 and y == 0:
      # Already horizontal
      for point in self.points:
        horizontal_points.append([point[0], point[1], 0])
    else:
      alpha = GeometryHelper.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:
        sys.stderr.write('Warning: rotation base matrix returned None\n')
      else:
        for point in self.points:
          new_point = np.matmul(rotation_base_matrix, point)
          horizontal_points.append(new_point)
    return horizontal_points