shared wall implementation bug fixing

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()
-
+      transformer = Transformer.from_crs(input_reference, gps)
-      coordinates = pyproj.transform(input_reference, gps, self.lower_corner[0], self.lower_corner[1])
+      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

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
 

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):
+  def add_shared(self, surface, intersection_area):
-    print(surface, self, 'are shared')
+    percent = intersection_area / self.area
-    self._shared_surfaces.append((100, surface))
+    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))
+      intersection_area = self.intersect(surface).area
-      surface.add_shared(self)
+      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)
+      if coordinates == '':
-    return self._projected_surface
+        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

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
 
 

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):

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)
 

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)
+    '''
+

tests_data/buildings.gml

@@ -3,8 +3,8 @@
 	<name>Gowanus 2050 Best Practice Scenario</name>
 	<boundedBy>
 		<Envelope srsName="EPSG:32118" srsDimension="3" xmlns:brid="http://www.opengis.net/citygml/bridge/2.0" xmlns:tran="http://www.opengis.net/citygml/transportation/2.0" xmlns:frn="http://www.opengis.net/citygml/cityfurniture/2.0" xmlns:wtr="http://www.opengis.net/citygml/waterbody/2.0" xmlns:sch="http://www.ascc.net/xml/schematron" xmlns:veg="http://www.opengis.net/citygml/vegetation/2.0" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:tun="http://www.opengis.net/citygml/tunnel/2.0" xmlns:tex="http://www.opengis.net/citygml/texturedsurface/2.0" xmlns:gml="http://www.opengis.net/gml" xmlns:gen="http://www.opengis.net/citygml/generics/2.0" xmlns:dem="http://www.opengis.net/citygml/relief/2.0" xmlns:app="http://www.opengis.net/citygml/appearance/2.0" xmlns:luse="http://www.opengis.net/citygml/landuse/2.0" xmlns:xAL="urn:oasis:names:tc:ciq:xsdschema:xAL:2.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:smil20lang="http://www.w3.org/2001/SMIL20/Language" xmlns:pbase="http://www.opengis.net/citygml/profiles/base/2.0" xmlns:smil20="http://www.w3.org/2001/SMIL20/" xmlns:bldg="http://www.opengis.net/citygml/building/2.0" xmlns:core="http://www.opengis.net/citygml/2.0" xmlns:grp="http://www.opengis.net/citygml/cityobjectgroup/2.0">
-			<lowerCorner>299606.4441129853 55348.37638737355 0</lowerCorner>
+			<gml:lowerCorner>301067.2312223603 55182.50627018605 -15.371661394834565</gml:lowerCorner>
-			<upperCorner>301879.9050504853 57594.05119206105 62.04879541695123</upperCorner>
+            <gml:upperCorner>301852.9778043915 57588.63517643605 75.67673757672314</gml:upperCorner>
 		</Envelope>
 	</boundedBy>
 	<cityObjectMember>