diff --git a/city_model_structure/city.py b/city_model_structure/city.py index 4f661bd2..5002de02 100644 --- a/city_model_structure/city.py +++ b/city_model_structure/city.py @@ -1,6 +1,7 @@ from city_model_structure.city_object import CityObject from typing import List, Union import pyproj +from pyproj import Transformer import reverse_geocoder as rg @@ -25,8 +26,8 @@ class City: except pyproj.exceptions.CRSError: print('Invalid projection reference system, please check the input data. (e.g. in CityGML files: srs_name)') quit() - - coordinates = pyproj.transform(input_reference, gps, self.lower_corner[0], self.lower_corner[1]) + transformer = Transformer.from_crs(input_reference, gps) + coordinates = transformer.transform(self.lower_corner[0], self.lower_corner[1]) self._location = rg.search(coordinates) return self._location @@ -59,6 +60,7 @@ class City: return None def add_city_object(self, new_city_object): + print('add city object') if self._city_objects is None: self._city_objects = [] for city_object in self.city_objects: @@ -66,6 +68,7 @@ class City: for surface in city_object.surfaces: for surface2 in new_city_object.surfaces: surface.shared(surface2) + print('added city object') self._city_objects.append(new_city_object) @property diff --git a/city_model_structure/city_object.py b/city_model_structure/city_object.py index 0e9ca95a..47355115 100644 --- a/city_model_structure/city_object.py +++ b/city_model_structure/city_object.py @@ -14,7 +14,8 @@ from typing import Union, List class CityObject: - def __init__(self, name, lod, surfaces, terrains, year_of_construction, function, attic_heated=0, basement_heated=0): + def __init__(self, name, lod, surfaces, terrains, year_of_construction, function, lower_corner, attic_heated=0, + basement_heated=0): self._name = name self._lod = lod self._surfaces = surfaces @@ -24,6 +25,7 @@ class CityObject: self._terrains = terrains self._year_of_construction = year_of_construction self._function = function + self._lower_corner = lower_corner self._geometry = Geometry() self._average_storey_height = None self._storeys_above_ground = None @@ -42,6 +44,7 @@ class CityObject: t.bounded = [ThermalBoundary(s, [t]) for s in t.surfaces] surface_id = 0 for s in self._surfaces: + s.lower_corner = self._lower_corner s.parent(self, surface_id) surface_id += 1 diff --git a/city_model_structure/surface.py b/city_model_structure/surface.py index e7989e3d..b431caad 100644 --- a/city_model_structure/surface.py +++ b/city_model_structure/surface.py @@ -2,6 +2,7 @@ from __future__ import annotations import numpy as np import pyny3d.geoms as pn from helpers.geometry import Geometry +from typing import Union class Surface: @@ -14,8 +15,10 @@ class Surface: 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 @@ -25,6 +28,10 @@ class Surface: self._parent = None self._shapely = None self._projected_surface = None + self._lower_corner = 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) @@ -54,6 +61,54 @@ class Surface: 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): + if self._ground_points is None: + coordinates = '' + for point in self.points: + x = point[0] - self._lower_corner[0] + y = point[1] - self._lower_corner[1] + z = point[2] - self._lower_corner[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): if self._points_list is None: @@ -71,6 +126,16 @@ class Surface: self._polygon = None return self._polygon + @property + def ground_polygon(self): + 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: @@ -139,16 +204,19 @@ class Surface: self._type = 'Roof' return self._type - def add_shared(self, surface): - print(surface, self, 'are shared') - self._shared_surfaces.append((100, surface)) + 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': + if self.type is not 'Wall' or surface.type is not 'Wall': return if self._geometry.is_almost_same_surface(self, surface): - self._shared_surfaces.append((100, surface)) - surface.add_shared(self) + 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): @@ -174,15 +242,53 @@ class Surface: 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._lower_corner = (min_x, min_y, min_z) + surface._lower_corner = (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 = '' - for coordinate in self.shapely.exterior.coords: + 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' - self._projected_surface = Surface(coordinates) - return self._projected_surface + 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 diff --git a/geometry/geometry_feeders/city_gml.py b/geometry/geometry_feeders/city_gml.py index 57ee84f6..ad7c3104 100644 --- a/geometry/geometry_feeders/city_gml.py +++ b/geometry/geometry_feeders/city_gml.py @@ -98,7 +98,8 @@ class CityGml: year_of_construction = o['Building']['yearOfConstruction'] if 'function' in o['Building']: function = o['Building']['function'] - self._city.add_city_object(CityObject(name, lod, surfaces, terrains, year_of_construction, function)) + self._city.add_city_object(CityObject(name, lod, surfaces, terrains, year_of_construction, function, + self._lower_corner)) return self._city diff --git a/helpers/geometry.py b/helpers/geometry.py index ec6ed722..c99dabff 100644 --- a/helpers/geometry.py +++ b/helpers/geometry.py @@ -19,20 +19,25 @@ class Geometry: # s1 and s2 are at least almost parallel surfaces # calculate distance point to plane using all the vertex # select surface1 value for the point (X,Y,Z) where two of the values are 0 - minimum_distance = 200000 + minimum_distance = self._delta + 1 parametric = s2.polygon.get_parametric() n2 = s2.normal for point in s1.points: distance = abs( (point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3]) normal_module = math.sqrt(pow(n2[0], 2) + pow(n2[1], 2) + pow(n2[2], 2)) + if normal_module == 0: - print(normal_module) continue distance = distance / normal_module if distance < minimum_distance: minimum_distance = distance - return minimum_distance <= self._delta + if minimum_distance <= self._delta: + break + if minimum_distance > self._delta or s1.intersect(s2) is None: + return False + else: + return True @staticmethod def to_points_matrix(points, remove_last=False): diff --git a/physics/physics_feeders/us_physics_parameters.py b/physics/physics_feeders/us_physics_parameters.py index 79946968..34202f0e 100644 --- a/physics/physics_feeders/us_physics_parameters.py +++ b/physics/physics_feeders/us_physics_parameters.py @@ -5,6 +5,6 @@ from physics.physics_feeders.helpers.us_to_library_types import UsToLibraryTypes class UsPhysicsParameters(UsBasePhysicsParameters): def __init__(self, city): self._city = city - self._climate_zone = UsToLibraryTypes.city_to_climate_zone(city.city_name) + self._climate_zone = UsToLibraryTypes.city_to_climate_zone(city.name) super().__init__(self._climate_zone, self._city.city_objects, lambda function: function) diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 3e590c58..4633c6ea 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -1,6 +1,7 @@ from unittest import TestCase from pathlib import Path from geometry.geometry_factory import GeometryFactory +from city_model_structure.surface import Surface class TestGeometryFactory(TestCase): @@ -127,4 +128,21 @@ class TestGeometryFactory(TestCase): self.assertIsNone(thermal_opening.outside_reflectance, 'thermal_opening outside_reflectance was initialized') self.assertIsNone(thermal_opening.thickness_m, 'thermal_opening thickness_m was initialized') - self.assertRaises(Exception, lambda:thermal_opening.u_value, 'thermal_opening u_value was initialized') + self.assertRaises(Exception, lambda: thermal_opening.u_value, 'thermal_opening u_value was initialized') + + def test_rotation(self): + s = Surface('-4.330008674901071 -2.2499999999999982 0.0 0.12097547700892843 -2.2499999999999982 0.0 0.12097547700892843 -2.25 0.0 -4.330008674901074 -2.2499999999999996 0.0', remove_last=False) + print(s.inclination) + print(s.polygon) + surface3 = Surface('0 0 0 0 0 3 3 0 3 3 0 0 0 0 0') + surface4 = Surface('0 0 0 0 1 3 3 1 3 3 0 0 0 0 0') + + + surface4.intersect(surface3) + ''' + surface1 = Surface('301104.7614957978 56642.02970768605 9.0 301097.510030954 56646.7245319048 0.0 301091.984640329 56650.30216862355 9.0 301104.7614957978 56642.02970768605 9.0') + surface2 = Surface('301070.5715543915 56635.9657428423 0.0 301070.863546579 56617.6366412798 0.0 301067.7356168915 56631.5018756548 0.0 301070.5715543915 56635.9657428423 0.0') + surface1.intersect(surface2) + surface2.intersect(surface1) + ''' + diff --git a/tests_data/buildings.gml b/tests_data/buildings.gml index d00d2d3d..a86229e2 100644 --- a/tests_data/buildings.gml +++ b/tests_data/buildings.gml @@ -3,8 +3,8 @@ Gowanus 2050 Best Practice Scenario - 299606.4441129853 55348.37638737355 0 - 301879.9050504853 57594.05119206105 62.04879541695123 + 301067.2312223603 55182.50627018605 -15.371661394834565 + 301852.9778043915 57588.63517643605 75.67673757672314