diff --git a/city_model_structure/attributes/polygon.py b/city_model_structure/attributes/polygon.py index fe8cb2da..08584db2 100644 --- a/city_model_structure/attributes/polygon.py +++ b/city_model_structure/attributes/polygon.py @@ -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 diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py index 0869c9c4..655361b0 100644 --- a/city_model_structure/attributes/polyhedron.py +++ b/city_model_structure/attributes/polyhedron.py @@ -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 diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py index de3377fb..ac3c5dd5 100644 --- a/city_model_structure/attributes/surface.py +++ b/city_model_structure/attributes/surface.py @@ -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 diff --git a/city_model_structure/building.py b/city_model_structure/building.py index 4607fadc..b134a677 100644 --- a/city_model_structure/building.py +++ b/city_model_structure/building.py @@ -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): diff --git a/helpers/geometry_helper.py b/helpers/geometry_helper.py index f29be90d..3f0e62f2 100644 --- a/helpers/geometry_helper.py +++ b/helpers/geometry_helper.py @@ -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): """ diff --git a/imports/geometry_factory.py b/imports/geometry_factory.py index 902e88f3..a9c55d85 100644 --- a/imports/geometry_factory.py +++ b/imports/geometry_factory.py @@ -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 diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 8d94bd46..3a31ccb8 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -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)