diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py
index 011a3d11..a8ed3d42 100644
--- a/city_model_structure/attributes/polyhedron.py
+++ b/city_model_structure/attributes/polyhedron.py
@@ -65,6 +65,10 @@ class Polyhedron:
   def _point(coordinates):
     return coordinates[0], coordinates[1], coordinates[2]
 
+  def _ground_coordinates(self, surface):
+    for point in surface.points:
+      print(point[0] - self.min_x, point[1] - self.min_y, point[2] - self.min_z)
+
   def _triangulate(self, surface):
     triangles = []
     triangles_count = len(surface.points) - 2
@@ -78,7 +82,8 @@ class Polyhedron:
       rest_points = ' '.join(str(e) for e in [*points_list[0:point_index+3], *points_list[point_index+6:]])
       triangular_surface = Surface(triangle_points, remove_last=False)
       rest_surface = Surface(rest_points, remove_last=False)
-      if self._geometry.almost_same_area(area, (triangular_surface.area + rest_surface.area)):
+      t_area = triangular_surface.area
+      if t_area == 0 and self._geometry.almost_same_area(area, (t_area + rest_surface.area)):
         area = rest_surface.area
         triangles.append(triangular_surface)
         points_list = rest_surface.points_list
diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py
index a6fd3a80..e3bd1bfc 100644
--- a/city_model_structure/attributes/surface.py
+++ b/city_model_structure/attributes/surface.py
@@ -217,42 +217,46 @@ class Surface:
     Surface area in square meters
     :return: float
     """
+    # New method to calculate area
     if self._area is None:
+      if len(self.points) < 3:
+        print('Warning: the area of a line or point cannot be calculated. Area = 0')
+        return 0
       alpha = 0
       vec_1 = self.points[1] - self.points[0]
       for i in range(2, len(self.points)):
         vec_2 = self.points[i] - self.points[0]
         alpha += self.angle_between_vectors(vec_1, vec_2)
-      if len(self.points) < 3 or alpha == 0:
-        area = 0
+      if alpha == 0:
         print('Warning: the area of a line or point cannot be calculated. Area = 0')
-      else:
-        points_2d = np.array(self.rotate_surface_to_horizontal())
-        polygon_2d = pn.Polygon(points_2d)
-        area = 0
-        for i in range(0, len(polygon_2d.points)-1):
-          point = polygon_2d.points[i]
-          next_point = polygon_2d.points[i+1]
-          area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
-        next_point = polygon_2d.points[0]
-        point = polygon_2d.points[len(polygon_2d.points)-1]
+        return 0
+      points_2d = self.rotate_surface_to_horizontal(self)
+      polygon_2d = pn.Polygon(np.array(points_2d))
+      area = 0
+      for i in range(0, len(polygon_2d.points)-1):
+        point = polygon_2d.points[i]
+        next_point = polygon_2d.points[i+1]
         area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
+      next_point = polygon_2d.points[0]
+      point = polygon_2d.points[len(polygon_2d.points)-1]
+      area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
       self._area = abs(area)
     return self._area
 
-  def rotate_surface_to_horizontal(self):
+  @staticmethod
+  def rotate_surface_to_horizontal(surface):
     z_vector = [0, 0, 1]
-    normal_vector = self.normal
+    normal_vector = surface.normal
     points_2d = []
     x = normal_vector[0]
     y = normal_vector[1]
 
     if x == 0 and y == 0:
       # Already horizontal
-      for point in self.points:
+      for point in surface.points:
         points_2d.append([point[0], point[1], 0])
     else:
-      alpha = self.angle_between_vectors(normal_vector, z_vector)
+      alpha = surface.angle_between_vectors(normal_vector, z_vector)
       rotation_line = np.cross(normal_vector, z_vector)
       third_axis = np.cross(normal_vector, rotation_line)
       w_1 = rotation_line / np.linalg.norm(rotation_line)
@@ -268,7 +272,7 @@ class Surface:
       if rotation_base_matrix is None:
         print('Error processing rotation base matrix')
       else:
-        for point in self.points:
+        for point in surface.points:
           new_point = np.matmul(rotation_base_matrix, point)
           points_2d.append(new_point)
     return points_2d
diff --git a/helpers/geometry_helper.py b/helpers/geometry_helper.py
index 4d5bb477..1f57d1d8 100644
--- a/helpers/geometry_helper.py
+++ b/helpers/geometry_helper.py
@@ -41,6 +41,8 @@ class GeometryHelper:
     :param a2
     :return: Boolean
     """
+    if a1 == 0 or a2 == 0:
+      return False
     delta = math.fabs(a1 - a2)
     return delta <= self._area_delta