hub/city_model_structure/surface.py

294 lines
8.6 KiB
Python

from __future__ import annotations
import numpy as np
import pyny3d.geoms as pn
from helpers.geometry import Geometry
from typing import Union
class Surface:
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):
self._parent = parent
self._name = str(surface_id)
@property
def name(self):
if self._name is None:
raise Exception('surface has no name')
return self._name
@property
def swr(self):
return self._swr
@swr.setter
def swr(self, value):
self._swr = value
@property
def points(self) -> 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):
if self._min_x is None:
self._min_x = self._min_coord('x')
return self._min_x
@property
def min_y(self):
if self._min_y is None:
self._min_y = self._min_coord('y')
return self._min_y
@property
def min_z(self):
if self._min_z is None:
self._min_z = self._min_coord('z')
return self._min_z
@property
def ground_points(self) -> 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:
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) -> pn.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]:
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):
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 t in terrain_points:
for g in ground_points:
if self._geometry.almost_equal(t, g):
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):
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):
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:
if self._normal is None:
points = self.points
n = np.cross(points[1] - points[0], points[2] - points[0])
self._normal = n / np.linalg.norm(n)
return self._normal
@property
def azimuth(self):
if self._azimuth is None:
normal = self.normal
self._azimuth = np.arctan2(normal[1], normal[0])
return self._azimuth
@property
def inclination(self):
if self._inclination is None:
self._inclination = np.arccos(self.normal[2])
return self._inclination
@property
def type(self):
if self._type is None:
grad = np.rad2deg(self.inclination)
if 170 <= grad:
self._type = 'Ground'
elif 80 <= grad <= 100:
self._type = 'Wall'
else:
self._type = 'Roof'
return self._type
def add_shared(self, surface, intersection_area):
percent = intersection_area / self.area
self._shared_surfaces.append((percent, surface))
def shared(self, surface):
if self.type is not 'Wall' or surface.type is not 'Wall':
return
if self._geometry.is_almost_same_surface(self, surface):
intersection_area = self.intersect(surface).area
percent = intersection_area / self.area
self._shared_surfaces.append((percent, surface))
surface.add_shared(self, intersection_area)
@property
def global_irradiance_hour(self):
return self._global_irradiance_hour
@global_irradiance_hour.setter
def global_irradiance_hour(self, value):
self._global_irradiance_hour = value
@property
def global_irradiance_month(self):
return self._global_irradiance_month
@global_irradiance_month.setter
def global_irradiance_month(self, value):
self._global_irradiance_month = value
@property
def shapely(self):
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):
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:
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]:
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