From 9999c8a33188b6de071f9a3ab72b1d753b9d61e2 Mon Sep 17 00:00:00 2001
From: Guille <guillermo.gutierrezmorote@concordia.ca>
Date: Thu, 28 May 2020 15:22:46 -0400
Subject: [PATCH] shared wall implementation and code refactoring.

---
 city_model_structure/city.py    | 33 +++++++++++++++++++++++++++++----
 city_model_structure/surface.py |  5 +----
 helpers/geometry.py             | 30 +++++++++++++++++++++++++++---
 3 files changed, 57 insertions(+), 11 deletions(-)

diff --git a/city_model_structure/city.py b/city_model_structure/city.py
index 02b095f0..a41c771b 100644
--- a/city_model_structure/city.py
+++ b/city_model_structure/city.py
@@ -1,5 +1,7 @@
 from city_model_structure.city_object import CityObject
 from typing import List, Union
+import pyproj
+import reverse_geocoder as rg
 
 
 class City:
@@ -10,6 +12,33 @@ class City:
     self._upper_corner = upper_corner
     self._city_objects = city_objects
     self._srs_name = srs_name
+    # todo: right now extracted at city level, in the future should be extracted also at building level if exist
+    self._location = None
+
+  @property
+  def location(self):
+    if self._location is None:
+      gps = pyproj.CRS('EPSG:4326')  # LatLon with WGS84 datum used by GPS units and Google Earth
+      input_reference = None
+      try:
+        input_reference = pyproj.CRS(self.srs_name)  # Projected coordinate system from input data
+      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])
+      self._location = rg.search(coordinates)
+    return self._location
+
+  @property
+  def country_code(self):
+    return self.location[0]['cc']
+
+  @property
+  def name(self):
+    if self._name is None:
+      self._name = self.location[0]['name']
+    return self._name
 
   @property
   def city_objects(self) -> Union[List[CityObject], None]:
@@ -43,10 +72,6 @@ class City:
   def srs_name(self):
     return self._srs_name
 
-  @property
-  def name(self):
-    return self._name
-
   @name.setter
   def name(self, value):
     self._name = value
diff --git a/city_model_structure/surface.py b/city_model_structure/surface.py
index 10b17e76..5a43a649 100644
--- a/city_model_structure/surface.py
+++ b/city_model_structure/surface.py
@@ -77,9 +77,6 @@ class Surface:
       self._area = self.polygon.get_area()
     return self._area
 
-  def _is_shared(self, surface):
-    return False
-
   def _is_almost_same_terrain(self, terrain_points, ground_points):
     equal = 0
     for t in terrain_points:
@@ -145,7 +142,7 @@ class Surface:
   def shared(self, surface):
     if self.type is not 'Wall':
       return
-    if self._is_shared(surface):
+    if self._geometry.is_almost_same_surface(self, surface):
       self._shared_surfaces.append((100, surface))
       surface.shared(self)
 
diff --git a/helpers/geometry.py b/helpers/geometry.py
index 5474b531..bd77683e 100644
--- a/helpers/geometry.py
+++ b/helpers/geometry.py
@@ -7,13 +7,37 @@ class Geometry:
     self._delta = delta
 
   def almost_equal(self, v1, v2):
-    delta = math.sqrt(pow((v1[0]-v2[0]), 2) + pow((v1[1]-v2[1]), 2) + pow((v1[2]-v2[2]), 2))
+    delta = math.sqrt(pow((v1[0] - v2[0]), 2) + pow((v1[1] - v2[1]), 2) + pow((v1[2] - v2[2]), 2))
     return delta <= self._delta
 
+  def is_almost_same_surface(self, s1, s2):
+    # delta is grads an need to be converted into radians
+    delta = np.rad2deg(self._delta)
+    difference = (s1.inclination - s2.inclination) % math.pi
+    if abs(difference) > delta:
+      return False
+    # s1 and s2 are at least almost parallel surfaces
+    p1 = s1.polygon.get_parametric()
+    selected_coefficient = 0
+    for coefficient in p1:
+      if coefficient != 0:
+        break;
+      selected_coefficient += 1
+
+    # calculate distance point to plane
+    # select surface1 value for the point (X,Y,Z) where two of the values are 0
+
+    s1_coefficient = -p1[3] / p1[selected_coefficient]
+    p2 = s2.polygon.get_parametric()
+    n2 = s2.normal
+    parametric = abs(p2[2] * s1_coefficient + p2[3])
+    distance = parametric / math.sqrt(pow(n2[0], 2) + pow(n2[1], 2) + pow(n2[2], 2))
+    return distance <= self._delta
+
   @staticmethod
   def to_points_matrix(points, remove_last=False):
-    rows = points.size//3
+    rows = points.size // 3
     points = points.reshape(rows, 3)
     if remove_last:
-      points = np.delete(points, rows-1, 0)
+      points = np.delete(points, rows - 1, 0)
     return points