diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py
index 449737ab..44d9bfac 100644
--- a/city_model_structure/attributes/surface.py
+++ b/city_model_structure/attributes/surface.py
@@ -9,6 +9,7 @@ from typing import Union
 
 import numpy as np
 import pyny3d.geoms as pn
+import math
 
 from helpers.geometry_helper import GeometryHelper
 
@@ -198,7 +199,7 @@ class Surface:
     return self._ground_polygon
 
   @property
-  def area(self):
+  def area_geoms_class(self):
     """
     Surface area in square meters
     :return: float
@@ -210,6 +211,79 @@ class Surface:
         self._area = 0
     return self._area
 
+  @property
+  def area(self):
+    """
+    Surface area in square meters
+    :return: float
+    """
+    # New method to calculate area
+    if self._area is None:
+      print('NEW METHOD TO CALCULATE AREA')
+      print('original:')
+      print(self.points)
+
+      # 1. 3D -> 2D
+      z_vector = [0, 0, 1]
+      normal_vector = self.normal
+      print(normal_vector)
+      turning_base_matrix = None
+      points_2d = []
+      x = normal_vector[0]
+      y = normal_vector[1]
+      z = normal_vector[2]
+      print('x:', x)
+      print('y:', y)
+      print('z:', z)
+      if x != 0 or y != 0:
+        # cos(alpha) = n.z/|n|.|z|
+        print('dot_mult:', np.dot(normal_vector, z_vector))
+        alpha = math.acos(np.dot(normal_vector, z_vector) / np.linalg.norm(normal_vector) / np.linalg.norm(z_vector))
+        print('alpha:', alpha)
+        turning_line = np.cross(normal_vector, z_vector)
+        print('turning_line:', turning_line)
+        third_axis = np.cross(normal_vector, turning_line)
+        print('third_axis:', third_axis)
+        # orthonormal base
+        w_1 = turning_line / np.linalg.norm(turning_line)
+        w_2 = normal_vector
+        w_3 = third_axis / np.linalg.norm(third_axis)
+        # turning_base_matrix
+        turning_matrix = np.array([[1, 0, 0],
+                          [0, math.cos(alpha), -math.sin(alpha)],
+                          [0, math.sin(alpha), math.cos(alpha)]])
+        print('turning_matrix:', turning_matrix)
+        base_matrix = np.array([w_1, w_2, w_3])
+        print('base_matrix:', base_matrix)
+        turning_base_matrix = np.matmul(base_matrix.transpose(), turning_matrix.transpose())
+        turning_base_matrix = np.matmul(turning_base_matrix, base_matrix)
+        print('turning_base_matrix:', turning_base_matrix)
+
+        if turning_base_matrix is None:
+          print('ERROR')
+        else:
+          for point in self.points:
+            new_point = np.matmul(turning_base_matrix, point)
+            print('new_point:', new_point)
+            points_2d.append(new_point)
+#            points_2d.append([new_point[0], new_point[1]])
+
+      else:
+        for point in self.points:
+          points_2d.append([point[0], point[1], 0])
+
+      polygon_2d = pn.Polygon(np.array(points_2d))
+      print('2D:')
+      print(polygon_2d.points)
+      # 2. calculate area:
+      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])
+      self._area = area
+    return self._area
+
   def _is_almost_same_terrain(self, terrain_points, ground_points):
     equal = 0
     for terrain_point in terrain_points: