Merge pull request 'shared_surfaces_method' (#5) from shared_surfaces_method into main

Reviewed-on: https://nextgenerations-cities.encs.concordia.ca/gitea/CERC/hub/pulls/5

hub/city_model_structure/city_object.py

@@ -34,6 +34,7 @@ class CityObject:
     self._max_y = ConfigurationHelper().min_coordinate
     self._max_z = ConfigurationHelper().min_coordinate
     self._centroid = None
+    self._volume = None
     self._external_temperature = dict()
     self._global_horizontal = dict()
     self._diffuse = dict()
@@ -63,7 +64,13 @@ class CityObject:
     Get city object volume in cubic meters
     :return: float
     """
-    return self.simplified_polyhedron.volume
+    if self._volume is None:
+      self._volume = self.simplified_polyhedron.volume
+    return self._volume
+
+  @volume.setter
+  def volume(self, value):
+    self._volume = value
 
   @property
   def detailed_polyhedron(self) -> Polyhedron:

hub/exports/exports_factory.py

@@ -71,7 +71,7 @@ class ExportsFactory:
     Export the city geometry to obj with grounded coordinates
     :return: None
     """
-    return Obj(self._city, self._path).to_ground_points()
+    return Obj(self._city, self._path)
 
   @property
   def _sra(self):

hub/exports/formats/obj.py

@@ -11,24 +11,51 @@ from hub.exports.formats.triangular import Triangular
 from hub.imports.geometry_factory import GeometryFactory
 
 
-class Obj(Triangular):
+class Obj:
   """
   Export to obj format
   """
   def __init__(self, city, path):
-    super().__init__(city, path, 'obj')
+    self._city = city
+    self._path = path
+    self._export()
+
+  def _to_vertex(self, coordinate):
+    x = coordinate[0] - self._city.lower_corner[0]
+    y = coordinate[1] - self._city.lower_corner[1]
+    z = coordinate[2] - self._city.lower_corner[2]
+    return f'v {x} {y} {z}\n'
+
+  def _export(self):
+    if self._city.name is None:
+      self._city.name = 'unknown_city'
+    file_name = self._city.name + '.obj'
+    file_path = (Path(self._path).resolve() / file_name).resolve()
+    vertices = {}
+    with open(file_path, 'w') as obj:
+      obj.write("# cerc-hub export\n")
+      vertex_index = 0
+      faces = []
+      for building in self._city.buildings:
+        obj.write(f'# building {building.name}\n')
+        for surface in building.surfaces:
+          obj.write(f'# surface {surface.name}\n')
+          face = 'f '
+          for coordinate in surface.perimeter_polygon.coordinates:
+            vertex = self._to_vertex(coordinate)
+            if vertex not in vertices.keys():
+              vertex_index += 1
+              vertices[vertex] = vertex_index
+              current = vertex_index
+              obj.write(vertex)
+            else:
+              current = vertices[vertex]
+
+            face = f'{face} {current}'
+
+          faces.append(f'{face} {face.split(" ")[1]}\n')
+          obj.writelines(faces)
+          faces = []
+
+
 
-  def to_ground_points(self):
-    """
-    Move closer to the origin
-    """
-    file_name_in = self._city.name + '.' + self._triangular_format
-    file_name_out = self._city.name + '_ground.' + self._triangular_format
-    file_path_in = (Path(self._path).resolve() / file_name_in).resolve()
-    file_path_out = (Path(self._path).resolve() / file_name_out).resolve()
-    obj = GeometryFactory('obj', path=file_path_in)
-    scene = obj.scene
-    scene.rezero()
-    obj_file = trimesh.exchange.obj.export_obj(scene)
-    with open(file_path_out, 'w') as file:
-      file.write(obj_file)

hub/helpers/geometry_helper.py

@@ -16,6 +16,8 @@ from hub.city_model_structure.attributes.polygon import Polygon
 from hub.city_model_structure.attributes.polyhedron import Polyhedron
 from hub.helpers.location import Location
 
+from PIL import Image
+
 
 class MapPoint:
   def __init__(self, x, y):
@@ -63,7 +65,7 @@ class GeometryHelper:
     return MapPoint(((city.upper_corner[0] - coordinate[0]) * 0.5), ((city.upper_corner[1] - coordinate[1]) * 0.5))
 
   @staticmethod
-  def city_mapping(city, building_names=None):
+  def city_mapping(city, building_names=None, plot=False):
     """
 
     Returns a shared_information dictionary like
@@ -78,6 +80,8 @@ class GeometryHelper:
     y = int((city.upper_corner[1] - city.lower_corner[1]) * 0.5) + 1
     city_map = [['' for _ in range(y + 1)] for _ in range(x + 1)]
     map_info = [[{} for _ in range(y + 1)] for _ in range(x + 1)]
+    img = Image.new('RGB', (x + 1, y + 1), "black")  # create a new black image
+    city_image = img.load()  # create the pixel map
     for building_name in building_names:
       building = city.city_object(building_name)
       line = 0
@@ -103,13 +107,14 @@ class GeometryHelper:
                 'line_start': (coordinate[0], coordinate[1]),
                 'line_end': (next_coordinate[0], next_coordinate[1]),
               }
+              city_image[x, y] = (100, 0, 0)
             elif city_map[x][y] != building.name:
               neighbour = city.city_object(city_map[x][y])
               neighbour_info = map_info[x][y]
 
               # prepare the keys
-              neighbour_start_coordinate = f'{neighbour_info["line_start"][0]}_{neighbour_info["line_start"][1]}'
-              building_start_coordinate = f'{coordinate[0]}_{coordinate[1]}'
+              neighbour_start_coordinate = f'{GeometryHelper.coordinate_to_map_point(neighbour_info["line_start"], city)}'
+              building_start_coordinate = f'{GeometryHelper.coordinate_to_map_point(coordinate, city)}'
               neighbour_key = f'{neighbour.name}_{neighbour_start_coordinate}_{building_start_coordinate}'
               building_key = f'{building.name}_{building_start_coordinate}_{neighbour_start_coordinate}'
 
@@ -126,6 +131,10 @@ class GeometryHelper:
                   'line_end': (next_coordinate[0], next_coordinate[1]),
                   'neighbour_line_start': neighbour_info['line_start'],
                   'neighbour_line_end': neighbour_info['line_end'],
+                  'coordinate_start': f"{GeometryHelper.coordinate_to_map_point(coordinate, city)}",
+                  'coordinate_end': f"{GeometryHelper.coordinate_to_map_point(next_coordinate, city)}",
+                  'neighbour_start': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_start'], city)}",
+                  'neighbour_end': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_end'], city)}",
                   'shared_points': 1
                 }
 
@@ -142,6 +151,10 @@ class GeometryHelper:
                   'line_end': neighbour_info['line_end'],
                   'neighbour_line_start': (coordinate[0], coordinate[1]),
                   'neighbour_line_end': (next_coordinate[0], next_coordinate[1]),
+                  'neighbour_start': f"{GeometryHelper.coordinate_to_map_point(coordinate, city)}",
+                  'neighbour_end': f"{GeometryHelper.coordinate_to_map_point(next_coordinate, city)}",
+                  'coordinate_start': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_start'], city)}",
+                  'coordinate_end': f"{GeometryHelper.coordinate_to_map_point(neighbour_info['line_end'], city)}",
                   'shared_points': 1
                 }
 
@@ -154,6 +167,8 @@ class GeometryHelper:
               elif building not in neighbour.neighbours:
                 neighbour.neighbours.append(building)
           line += 1
+    if plot:
+      img.show()
     return lines_information
 
   @staticmethod

hub/imports/geometry/geojson.py

@@ -6,13 +6,15 @@ Project Coder Guillermo Gutierrez Guillermo.GutierrezMorote@concordia.ca
 """
 import json
 
+import numpy as np
 import trimesh.creation
 
 from pyproj import Transformer
 from shapely.geometry import Polygon as ShapelyPolygon
 
 import hub.helpers.constants as cte
-from hub.imports.geometry.helpers.geometry_helper import GeometryHelper
+from hub.helpers.geometry_helper import GeometryHelper
+from hub.imports.geometry.helpers.geometry_helper import GeometryHelper as igh
 from hub.city_model_structure.attributes.polygon import Polygon
 from hub.city_model_structure.building import Building
 from hub.city_model_structure.building_demand.surface import Surface
@@ -62,9 +64,11 @@ class Geojson:
     surfaces = []
     buildings = []
     for zone, surface_coordinates in enumerate(surfaces_coordinates):
-      points = GeometryHelper.points_from_string(GeometryHelper.remove_last_point_from_string(surface_coordinates))
+      points = igh.points_from_string(igh.remove_last_point_from_string(surface_coordinates))
+      # geojson provides the roofs, need to be transform into grounds
+      points = igh.invert_points(points)
       polygon = Polygon(points)
-      surfaces.append(Surface(polygon, polygon, surface_type=cte.GROUND))
+      surfaces.append(Surface(polygon, polygon))
       buildings.append(Building(f'{name}_zone_{zone}', surfaces, year_of_construction, function))
     return buildings
 
@@ -73,22 +77,41 @@ class Geojson:
     lod0_buildings = Geojson._create_buildings_lod0(name, year_of_construction, function, surface_coordinates)
     surfaces = []
     buildings = []
+
     for zone, lod0_building in enumerate(lod0_buildings):
-      for surface in lod0_building.surfaces:
-        shapely_polygon = ShapelyPolygon(surface.solid_polygon.coordinates)
-        if not shapely_polygon.is_valid:
-          print(surface.solid_polygon.area)
-          print('error?', name, surface_coordinates)
-          continue
-        mesh = trimesh.creation.extrude_polygon(shapely_polygon, height)
-        for face in mesh.faces:
-          points = []
-          for vertex_index in face:
-            points.append(mesh.vertices[vertex_index])
-          polygon = Polygon(points)
-          surface = Surface(polygon, polygon)
+      for surface in lod0_building.grounds:
+        volume = surface.solid_polygon.area * height
         surfaces.append(surface)
-        buildings.append(Building(f'{name}_zone_{zone}', surfaces, year_of_construction, function))
+        roof_coordinates = []
+        # adding a roof means invert the polygon coordinates and change the Z value
+        for coordinate in surface.solid_polygon.coordinates:
+          roof_coordinate = np.array([coordinate[0], coordinate[1], height])
+          # insert the roof rotated already
+          roof_coordinates.insert(0, roof_coordinate)
+        polygon = Polygon(roof_coordinates)
+        roof = Surface(polygon, polygon)
+        surfaces.append(roof)
+        # adding a wall means add the point coordinates and the next point coordinates with Z's height and 0
+        coordinates_length = len(roof.solid_polygon.coordinates)
+        for i, coordinate in enumerate(roof.solid_polygon.coordinates):
+          j = i + 1
+          if j == coordinates_length:
+            j = 0
+          next_coordinate = roof.solid_polygon.coordinates[j]
+          wall_coordinates = [
+            np.array([coordinate[0], coordinate[1], 0.0]),
+            np.array([next_coordinate[0], next_coordinate[1], 0.0]),
+            np.array([next_coordinate[0], next_coordinate[1], next_coordinate[2]]),
+            np.array([coordinate[0], coordinate[1], coordinate[2]])
+          ]
+          polygon = Polygon(wall_coordinates)
+          wall = Surface(polygon, polygon)
+          surfaces.append(wall)
+
+        building = Building(f'{name}_zone_{zone}', surfaces, year_of_construction, function)
+        building.volume = volume
+        buildings.append(building)
+
     return buildings
 
   def _get_polygons(self, polygons, coordinates):
@@ -106,6 +129,50 @@ class Geojson:
       polygons.append(transformed_coordinates.lstrip(' '))
       return polygons
 
+  @staticmethod
+  def _find_wall(line_1, line_2):
+    for i in range(0, 2):
+      point_1 = line_1[i]
+      point_2 = line_2[i]
+      distance = GeometryHelper.distance_between_points(point_1, point_2)
+      if distance > 1e-2:
+        return False
+    return True
+
+  def _store_shared_percentage_to_walls(self, city, city_mapped):
+    for building in city.buildings:
+      if building.name not in city_mapped.keys():
+        continue
+      building_mapped = city_mapped[building.name]
+      for wall in building.walls:
+        percentage = 0
+        ground_line = []
+        for point in wall.perimeter_polygon.coordinates:
+          if point[2] < 0.5:
+            ground_line.append(point)
+        # todo: erase when we have no triangulation
+        if len(ground_line) < 2:
+          continue
+        # todo: erase down to here
+        for entry in building_mapped:
+          if building_mapped[entry]['shared_points'] <= 5:
+            continue
+          line = [building_mapped[entry]['line_start'], building_mapped[entry]['line_end']]
+          neighbour_line = [building_mapped[entry]['neighbour_line_start'],
+                            building_mapped[entry]['neighbour_line_end']]
+          neighbour_height = city.city_object(building_mapped[entry]['neighbour_name']).max_height
+          if self._find_wall(line, ground_line):
+            line_shared = (GeometryHelper.distance_between_points(line[0], line[1]) +
+                           GeometryHelper.distance_between_points(neighbour_line[0], neighbour_line[1]) -
+                           GeometryHelper.distance_between_points(line[1], neighbour_line[0]) -
+                           GeometryHelper.distance_between_points(line[0], neighbour_line[1])) / 2
+            percentage_ground = line_shared / GeometryHelper.distance_between_points(line[0], line[1])
+            percentage_height = neighbour_height / building.max_height
+            if percentage_height > 1:
+              percentage_height = 1
+            percentage += percentage_ground * percentage_height
+        wall.percentage_shared = percentage
+
   @property
   def city(self) -> City:
     """
@@ -115,6 +182,7 @@ class Geojson:
       missing_functions = []
       buildings = []
       building_id = 0
+      lod = 1
       for feature in self._geojson['features']:
         extrusion_height = 0
         if self._extrusion_height_field is not None:
@@ -140,7 +208,6 @@ class Geojson:
           building_name = f'building_{building_id}'
           building_id += 1
         polygons = []
-        lod = 1
         for part, coordinates in enumerate(geometry['coordinates']):
           polygons = self._get_polygons(polygons, coordinates)
           for zone, polygon in enumerate(polygons):
@@ -163,6 +230,9 @@ class Geojson:
       for building in buildings:
         self._city.add_city_object(building)
       self._city.level_of_detail.geometry = lod
+      if lod == 1:
+        lines_information = GeometryHelper.city_mapping(self._city)
+        self._store_shared_percentage_to_walls(self._city, lines_information)
       if len(missing_functions) > 0:
         print(f'There are unknown functions {missing_functions}')
     return self._city

hub/imports/geometry/helpers/geometry_helper.py

@@ -45,3 +45,10 @@ class GeometryHelper:
     array = points.split(' ')
     res = " "
     return res.join(array[0:len(array) - 3])
+
+  @staticmethod
+  def invert_points(points):
+    res = []
+    for point in points:
+      res.insert(0,point)
+    return res

hub/requirements.txt

@@ -23,3 +23,4 @@ shapely
 geopandas
 triangle
 psycopg2-binary
+PIL

hub/unittests/test_geometry_factory.py

@@ -4,13 +4,11 @@ SPDX - License - Identifier: LGPL - 3.0 - or -later
 Copyright © 2022 Concordia CERC group
 Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
 """
-import datetime
+
 from pathlib import Path
 from unittest import TestCase
 from hub.helpers.geometry_helper import GeometryHelper
 
-from numpy import inf
-
 import hub.exports.exports_factory
 from hub.imports.construction_factory import ConstructionFactory
 from hub.imports.geometry_factory import GeometryFactory
@@ -118,7 +116,6 @@ class TestGeometryFactory(TestCase):
     city = self._get_city(file, 'rhino')
     self.assertIsNotNone(city, 'city is none')
     self.assertTrue(len(city.buildings) == 36)
-    i = 0
 
   def test_import_obj(self):
     """
@@ -142,7 +139,7 @@ class TestGeometryFactory(TestCase):
                           function_field='CODE_UTILI')
 
     hub.exports.exports_factory.ExportsFactory('obj', city, self._output_path).export()
-    self.assertEqual(207, len(city.buildings), 'wrong number of buildings')
+    self.assertEqual(195, len(city.buildings), 'wrong number of buildings')
     self._check_buildings(city)
 
   def test_map_neighbours(self):
@@ -150,13 +147,21 @@ class TestGeometryFactory(TestCase):
     Test neighbours map creation
     """
     file = 'neighbours.geojson'
+    city = self._get_city(file, 'geojson',
+                          year_of_construction_field='ANNEE_CONS',
+                          function_field='LIBELLE_UT')
+    info_lod0 = GeometryHelper.city_mapping(city, plot=False)
+
     city = self._get_city(file, 'geojson',
                           height_field='citygml_me',
                           year_of_construction_field='ANNEE_CONS',
                           function_field='LIBELLE_UT')
-    print(GeometryHelper.city_mapping(city))
+    info_lod1 = GeometryHelper.city_mapping(city, plot=False)
+    hub.exports.exports_factory.ExportsFactory('obj', city, self._output_path).export()
+    self.assertEqual(info_lod0, info_lod1)
     for building in city.buildings:
       self.assertEqual(2, len(building.neighbours))
+      print(building.volume)
 
     self.assertEqual('2_part_0_zone_0',city.city_object('1_part_0_zone_0').neighbours[0].name)
     self.assertEqual('3_part_0_zone_0',city.city_object('1_part_0_zone_0').neighbours[1].name)
@@ -164,4 +169,3 @@ class TestGeometryFactory(TestCase):
     self.assertEqual('3_part_0_zone_0',city.city_object('2_part_0_zone_0').neighbours[1].name)
     self.assertEqual('1_part_0_zone_0', city.city_object('3_part_0_zone_0').neighbours[0].name)
     self.assertEqual('2_part_0_zone_0', city.city_object('3_part_0_zone_0').neighbours[1].name)
-

hub/unittests/tests_data/neighbours.geojson

@@ -12,18 +12,18 @@
               -73.580414175680588,
               45.497641136608358
             ],
-            [
-              -73.581414175680588,
-              45.497641136608358
-            ],
-            [
-              -73.581414175680588,
-              45.498641136608358
-            ],
             [
               -73.580414175680588,
               45.498641136608358
             ],
+            [
+              -73.581414175680588,
+              45.498641136608358
+            ],
+            [
+              -73.581414175680588,
+              45.497641136608358
+            ],
             [
               -73.580414175680588,
               45.497641136608358
@@ -204,19 +204,20 @@
             [
               -73.581414175680588,
               45.497641136608358
+            ]
+            ,
+            [
+              -73.581414175680588,
+              45.498441136608358
+            ],
+            [
+              -73.582214175680588,
+              45.498441136608358
             ],
             [
               -73.582214175680588,
               45.497641136608358
             ],
-            [
-              -73.582214175680588,
-              45.498441136608358
-            ],
-            [
-              -73.581414175680588,
-              45.498441136608358
-            ],
             [
               -73.581414175680588,
               45.497641136608358
@@ -399,31 +400,30 @@
               -73.581914175680588,
               45.498441136608358
             ],
-            [
-              -73.581914175680588,
-              45.499641136608358
-            ],
-            [
-              -73.580914175680588,
-              45.499641136608358
-            ],
-            [
-              -73.580914175680588,
-              45.498641136608358
-            ],
-            [
-              -73.581414175680588,
-              45.498641136608358
-            ],
             [
               -73.581414175680588,
               45.498441136608358
             ],
+            [
+              -73.581414175680588,
+              45.498641136608358
+            ],
+            [
+              -73.580914175680588,
+              45.498641136608358
+            ],
+            [
+              -73.580914175680588,
+              45.499641136608358
+            ],
+            [
+              -73.581914175680588,
+              45.499641136608358
+            ],
             [
               -73.581914175680588,
               45.498441136608358
             ]
-
           ]
         ]
       },