Surface and Polygon classes finished. Added test_surfaces to test_geometry_factory.py. Still not working version.

This commit is contained in:
Pilar 2021-03-03 16:33:33 -05:00
parent fa83e195c4
commit dbd32e4d20
7 changed files with 161 additions and 271 deletions

View File

@ -17,20 +17,14 @@ class Polygon:
def __init__(self, vertices):
self._vertices = vertices
self._remove_last = True
self._area = None
self._points = None
self._normal = None
@property
def points(self) -> np.ndarray:
"""
Surface point matrix
:return: np.ndarray
"""
if self._points is None:
self._points = np.fromstring(self._vertices, dtype=float, sep=' ')
self._points = gh.to_points_matrix(self._points, self._remove_last)
self._points = self._vertices
return self._points
@property

View File

@ -13,6 +13,7 @@ from city_model_structure.attributes.surface import Surface
class Polyhedron:
# todo: modify class to call polygons better than surfaces
"""
Polyhedron class
"""
@ -218,7 +219,7 @@ class Polyhedron:
triangle = Surface(points, remove_last=False)
error_sum = 0
for i in range(0, len(normal)):
error_sum += triangle.normal[i] - normal[i]
error_sum += triangle.perimeter_surface.normal[i] - normal[i]
if np.abs(error_sum) < accepted_error:
is_concave = True
return is_concave
@ -291,7 +292,7 @@ class Polyhedron:
:param points: [point]
:return: boolean
"""
area_ear = ear.area
area_ear = ear.perimeter_surface.area
for point in points:
area_points = 0
point_is_not_vertex = True
@ -306,7 +307,7 @@ class Polyhedron:
else:
new_points = ' '.join(str(e) for e in [*ear.points[i][:], *point[:], *ear.points[0][:]])
new_triangle = Surface(new_points, remove_last=False)
area_points += new_triangle.area
area_points += new_triangle.perimeter_surface.area
if abs(area_points - area_ear) < 1e-6:
# point_inside_ear = True
return False

View File

@ -2,58 +2,45 @@
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
contributors Pilar Monsalvete Álvarez de Uribarri pilar.monsalvete@concordia.ca
"""
from __future__ import annotations
from typing import Union
import sys
from __future__ import annotations
import numpy as np
import pyny3d.geoms as pn
import math
from helpers.geometry_helper import GeometryHelper as gh
from city_model_structure.attributes.polygon import Polygon
# todo: @Guille, review methods depending on pyny3d
# todo: remove pyny3d, seems to not be supported and has some issues
# todo substitute pyny3d.polygon with city_model_structure.attributes.polygon
class Surface:
"""
Surface class
"""
def __init__(self, coordinates, surface_type=None, name=None, swr='0.2', remove_last=True, is_projected=False):
def __init__(self, coordinates, holes_coordinates=None, surface_type=None, name=None, swr='0.2', remove_last=True):
self._coordinates = coordinates
self._holes_coordinates = holes_coordinates
self._type = surface_type
self._name = name
self._swr = swr
self._remove_last = remove_last
self._is_projected = is_projected
self._polygon = None
self._ground_polygon = None
self._points = None
self._ground_points = None
self._points_list = None
self._holes_points = None
self._holes_points_list = 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
self._perimeter = Polygon(self.perimeter_vertices)
# todo: rethink how to do this
self._holes = [Polygon(self.vertices)]
self._opaque = Polygon(self.opaque_vertices)
self._perimeter_surface = None
self._hole_surfaces = None
self._solid_surface = None
self._perimeter_vertices = None
def parent(self, parent, surface_id):
"""
@ -95,7 +82,7 @@ class Surface:
@property
def points(self) -> np.ndarray:
"""
Surface point matrix
Surface point coordinates list [x, y, z, x, y, z,...]
:return: np.ndarray
"""
if self._points is None:
@ -103,6 +90,45 @@ class Surface:
self._points = gh.to_points_matrix(self._points, self._remove_last)
return self._points
@property
def holes_points(self) -> [np.ndarray]:
"""
Surface point coordinates list [x, y, z, x, y, z,...]
:return: np.ndarray
"""
if self._holes_coordinates is not None:
self._holes_points = []
for hole_coordinates in self._holes_coordinates:
hole_points = np.fromstring(hole_coordinates, dtype=float, sep=' ')
hole_points = gh.to_points_matrix(hole_points, self._remove_last)
self._holes_points.append(hole_points)
return self._holes_points
@property
def points_list(self) -> np.ndarray:
"""
Surface point matrix [[x, y, z],[x, y, z],...]
: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 holes_points_list(self) -> np.ndarray:
"""
Surface point matrix [[x, y, z],[x, y, z],...]
:return: np.ndarray
"""
if self._holes_coordinates is not None:
self._holes_points_list = np.array([])
for hole_points in self.holes_points:
s = hole_points
hole_points_list = np.reshape(s, len(s) * 3)
np.add(self._holes_points_list, hole_points_list)
return self._holes_points_list
def _min_coord(self, axis):
if axis == 'x':
axis = 0
@ -148,80 +174,6 @@ class Surface:
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 = gh.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
@staticmethod
def _is_almost_same_terrain(terrain_points, ground_points):
equal = 0
for terrain_point in terrain_points:
for ground_point in ground_points:
if gh().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):
"""
@ -229,20 +181,17 @@ class Surface:
:return: float
"""
if self._area_above_ground is None:
self._area_above_ground = self._perimeter.area - self.area_below_ground
self._area_above_ground = self.perimeter_surface.area - self.area_below_ground
return self._area_above_ground
# todo: to be implemented
@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._perimeter.area
return self._area_below_ground
return 0.0
@property
def azimuth(self):
@ -251,7 +200,7 @@ class Surface:
:return: float
"""
if self._azimuth is None:
normal = self._perimeter.normal
normal = self.perimeter_surface.normal
self._azimuth = np.arctan2(normal[1], normal[0])
return self._azimuth
@ -262,7 +211,7 @@ class Surface:
:return: float
"""
if self._inclination is None:
self._inclination = np.arccos(self._perimeter.normal[2])
self._inclination = np.arccos(self.perimeter_surface.normal[2])
return self._inclination
@property
@ -288,24 +237,19 @@ class Surface:
:param intersection_area:
:return:
"""
percent = intersection_area / self._perimeter.area
percent = intersection_area / self.perimeter_surface.area
self._shared_surfaces.append((percent, surface))
# todo reimplement
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 gh().is_almost_same_surface(self, surface):
try:
intersection_area = self.intersect(surface)._perimeter.area
except ValueError:
intersection_area = 0
self.add_shared(surface, intersection_area)
surface.add_shared(self, intersection_area)
# intersection_area = 0
# surface.add_shared(self, intersection_area)
raise NotImplementedError
@property
def global_irradiance(self) -> dict:
@ -324,94 +268,76 @@ class Surface:
self._global_irradiance = value
@property
def shapely(self) -> Union[None, pn.Polygon]:
def perimeter_surface(self) -> Polygon:
"""
Surface shapely (Z projection)
:return: None or pyny3d.Polygon
total surface defined by the perimeter, merging solid and holes
:return: 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)
if self._perimeter_surface is None:
self._perimeter_surface = Polygon(self.perimeter_vertices)
return self._perimeter_surface
@property
def projection(self) -> Surface:
def solid_surface(self) -> Polygon:
"""
Projected surface (Z projection)
:return: Surface
solid surface
:return: Polygon
"""
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
if self._solid_surface is None:
self._solid_surface = Polygon(self.points)
return self._solid_surface
@property
def convex(self):
return pn.Polygon.is_convex(self.polygon.points)
def hole_surfaces(self) -> [Polygon]:
"""
hole surfaces, a list of hole polygons found in the surface
:return: None, [] or [Polygon]
None -> not known whether holes exist in reality or not due to low level of detail of input data
[] -> no holes in the surface
[Polygon] -> one or more holes in the surface
"""
if self._hole_surfaces is None:
if self.holes_points is None:
self._hole_surfaces = None
else:
self._hole_surfaces = []
for hole_vertices in self.holes_points:
self._hole_surfaces.append(Polygon(hole_vertices))
return self._hole_surfaces
@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
def perimeter_vertices(self) -> np.ndarray:
"""
vertices of the perimeter organized in the same order as in coordinates
:return:
"""
if self._perimeter_vertices is None:
if self.holes_points is None:
self._perimeter_vertices = self.points
else:
first_point = True
for point in self.points:
point_of_hole = False
for hole_points in self.holes_points:
for hole_point in hole_points:
if gh().almost_equal(0.0, point, hole_point):
point_of_hole = True
if not point_of_hole:
if first_point:
self._perimeter_vertices = np.array([point])
first_point = False
else:
self._perimeter_vertices = np.append(self._perimeter_vertices, [point], axis=0)
# remove duplicated points
pv = np.array([self._perimeter_vertices[0]])
for point in self._perimeter_vertices:
duplicated_point = False
for p in pv:
if gh().almost_equal(0.0, p, point):
duplicated_point = True
if not duplicated_point:
pv = np.append(pv, [point], axis=0)
self._perimeter_vertices = pv
if not self._remove_last:
self._perimeter_vertices = np.append(self._perimeter_vertices, [self._perimeter_vertices[0]], axis=0)
return self._perimeter_vertices

View File

@ -8,12 +8,6 @@ contributors: Pilar Monsalvete pilar_monsalvete@yahoo.es
from typing import List
import matplotlib.patches as patches
import numpy as np
from matplotlib import pylab
from shapely import ops
from shapely.geometry import MultiPolygon
from city_model_structure.attributes.surface import Surface
from city_model_structure.attributes.thermal_boundary import ThermalBoundary
from city_model_structure.attributes.thermal_zone import ThermalZone
@ -207,56 +201,14 @@ class Building(CityObject):
def _tuple_to_point(xy_tuple):
return [xy_tuple[0], xy_tuple[1], 0.0]
def _plot(self, polygon):
points = ()
for point_tuple in polygon.exterior.coords:
almost_equal = False
for point in points:
point_1 = Building._tuple_to_point(point)
point_2 = Building._tuple_to_point(point_tuple)
if self._geometry.almost_equal(point_1, point_2):
almost_equal = True
break
if not almost_equal:
points = points + (point_tuple,)
points = points + (points[0],)
pylab.scatter([point[0] for point in points], [point[1] for point in points])
pylab.gca().add_patch(patches.Polygon(points, closed=True, fill=True))
pylab.grid()
pylab.show()
# todo, to be implemented
@property
def foot_print(self) -> Surface:
"""
City object foot print surface
:return: Surface
"""
if self._foot_print is None:
shapelys = []
union = None
for surface in self.surfaces:
if surface.shapely.is_empty or not surface.shapely.is_valid:
continue
shapelys.append(surface.shapely)
union = ops.unary_union(shapelys)
shapelys = [union]
if isinstance(union, MultiPolygon):
Exception('foot print returns a multipolygon')
points_list = []
for point_tuple in union.exterior.coords:
# ToDo: should be Z 0.0 or min Z?
point = Building._tuple_to_point(point_tuple)
almost_equal = False
for existing_point in points_list:
if self._geometry.almost_equal(point, existing_point):
almost_equal = True
break
if not almost_equal:
points_list.append(point)
points_list = np.reshape(points_list, len(points_list) * 3)
points = np.array_str(points_list).replace('[', '').replace(']', '')
self._foot_print = Surface(points, remove_last=False, is_projected=True)
return self._foot_print
raise NotImplementedError
@property
def type(self):

View File

@ -50,15 +50,16 @@ class GeometryHelper:
delta = math.fabs(a1 - a2)
return delta <= self._area_delta
def almost_equal(self, v1, v2):
def almost_equal(self, delta_max, v1, v2):
"""
Compare two points and decides if they are almost equal (quadratic error under delta)
Compare two points and decides if they are almost equal (distance under delta_max)
:param delta_max: maximum distance to be considered same point
:param v1: [x,y,z]
:param v2: [x,y,z]
:return: Boolean
"""
delta = self.distance_between_points(v1, v2)
return delta <= self._delta
return delta <= delta_max
def is_almost_same_surface(self, s1, s2):
"""

View File

@ -41,6 +41,10 @@ class GeometryFactory:
"""
return getattr(self, self._file_type, lambda: None)
@property
def _city_debug(self):
return CityGml(self._path).city
@property
def _osm_subway(self):
return OsmSubway(self._path).subway_entrances

View File

@ -6,9 +6,9 @@ Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@conc
import os
from pathlib import Path
from unittest import TestCase
import math
from city_model_structure.city import City
from imports.geometry_factory import GeometryFactory
from city_model_structure.attributes.surface import Surface
class TestGeometryFactory(TestCase):
@ -29,7 +29,7 @@ class TestGeometryFactory(TestCase):
def _get_citygml(self):
if self._city_gml is None:
file_path = (self._example_path / 'lod2_buildings.gml').resolve()
self._city_gml = GeometryFactory('citygml', file_path).city
self._city_gml = GeometryFactory('citygml', file_path)._city_debug
self.assertIsNotNone(self._city_gml, 'city is none')
return self._city_gml
@ -123,10 +123,10 @@ class TestGeometryFactory(TestCase):
:return: None
"""
city = self._get_citygml()
print(city)
for building in city.buildings:
for surface in building.surfaces:
self.assertIsNotNone(surface.name, 'surface name is none')
self.assertIsNotNone(surface.area, 'surface area is none')
self.assertIsNotNone(surface.type, 'surface type is none')
self.assertIsNotNone(surface.azimuth, 'surface azimuth is none')
self.assertIsNotNone(surface.inclination, 'surface inclination is none')
@ -134,18 +134,13 @@ class TestGeometryFactory(TestCase):
self.assertIsNotNone(surface.area_above_ground, 'surface area_above_ground is none')
self.assertIsNotNone(surface.points, 'surface points is none')
self.assertIsNotNone(surface.points_list, 'surface points_list is none')
self.assertIsNotNone(surface.polygon, 'surface polygon is none')
self.assertIsNotNone(surface.shapely, 'surface shapely is none')
self.assertIsNotNone(surface.global_irradiance, 'monthly irradiance is none')
self.assertIsNotNone(surface.normal, 'surface normal is none')
self.assertIsNotNone(surface.projection, 'surface projection is none')
self.assertIsNotNone(surface.swr, 'surface swr is none')
self.assertIsNotNone(surface.min_x, 'surface min_x is none')
self.assertIsNotNone(surface.min_y, 'surface min_y is none')
self.assertIsNotNone(surface.min_z, 'surface min_z is none')
self.assertIsNotNone(surface.ground_polygon, 'surface ground_polygon is none')
self.assertIsNotNone(surface.ground_points, 'surface ground_points is none')
self.assertIsNotNone(surface.intersect(surface), 'self intersection is none')
print('points', surface.points)
print('coordinates', surface._coordinates)
def test_citygml_thermal_zone(self):
"""
@ -243,10 +238,27 @@ class TestGeometryFactory(TestCase):
def test_divide_mesh_by_plane(self):
file_path = (self._example_path / 'FZK-Haus-LoD-all-KIT-IAI-KHH-B36-V1.gml').resolve()
# todo @Guille: this file has 5 lods (0, 1, 2, 3 and 4), all as one single city_object.
# Only lod1 is read ans saved
# Only lod1 is read and saved
city = GeometryFactory('citygml', file_path).city
for building in city.buildings:
print(building.name)
print(building.volume)
print(building.thermal_zones[0].floor_area)
print(len(building.storeys))
def test_surface(self):
coordinates = '0.0 0.0 0.0 0.0 4.0 0.0 4.0 4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 ' \
'1.0 1.0 0.0 2.0 1.0 0.0 2.0 2.0 0.0 1.0 2.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 ' \
'2.0 2.0 0.0 3.0 2.0 0.0 3.0 3.0 0.0 2.0 3.0 0.0 2.0 2.0 0.0 0.0 0.0 0.0'
holes_coordinates = ['1.0 1.0 0.0 2.0 1.0 0.0 2.0 2.0 0.0 1.0 2.0 0.0 1.0 1.0 0.0',
'2.0 2.0 0.0 3.0 2.0 0.0 3.0 3.0 0.0 2.0 3.0 0.0 2.0 2.0 0.0']
surface = Surface(coordinates, holes_coordinates=holes_coordinates, remove_last=True)
print(surface.points)
print(surface.holes_points)
print(surface.perimeter_vertices)
print(surface.hole_surfaces[1].area)
print(surface.perimeter_surface.area)
print(surface.solid_surface.area)
print(surface.hole_surfaces[1].normal)
print(surface.perimeter_surface.normal)
print(surface.solid_surface.normal)