diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py
index 347b5c95..27970eeb 100644
--- a/city_model_structure/attributes/polyhedron.py
+++ b/city_model_structure/attributes/polyhedron.py
@@ -71,7 +71,6 @@ class Polyhedron:
     points_list = surface.points_list
     point_index = 0
     area = surface.area
-    print('AREA:', area)
     while len(triangles) < triangles_count:
       # select a triangle starting at point index
       triangle_points = ' '.join(str(e) for e in [*points_list[point_index:point_index + 9]])
@@ -79,8 +78,6 @@ 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)
-      print('triangular_surface:', triangular_surface.points)
-      print('triangular_surface_AREA:', triangular_surface.area)
       if self._geometry.almost_same_area(area, (triangular_surface.area + rest_surface.area)):
         area = rest_surface.area
         triangles.append(triangular_surface)
diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py
index 788dec4a..7e1ea463 100644
--- a/city_model_structure/attributes/surface.py
+++ b/city_model_structure/attributes/surface.py
@@ -219,33 +219,33 @@ class Surface:
     """
     # New method to calculate area
     if self._area is None:
-      print('NEW METHOD TO CALCULATE AREA')
-      print('original:')
-      print(self.points)
+#      print('NEW METHOD TO CALCULATE AREA')
+#      print('original:')
+#      print(self.points)
       if len(self.points) < 3:
         area = 0
       else:
         # 1. 3D -> 2D
         z_vector = [0, 0, 1]
         normal_vector = self.normal
-        print(normal_vector)
+#        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)
+#        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))
+#          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)
+#          print('alpha:', alpha)
           turning_line = np.cross(normal_vector, z_vector)
-          print('turning_line:', turning_line)
+#          print('turning_line:', turning_line)
           third_axis = np.cross(normal_vector, turning_line)
-          print('third_axis:', third_axis)
+#          print('third_axis:', third_axis)
           # orthonormal base
           w_1 = turning_line / np.linalg.norm(turning_line)
           w_2 = normal_vector
@@ -254,19 +254,19 @@ class Surface:
           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)
+#          print('turning_matrix:', turning_matrix)
           base_matrix = np.array([w_1, w_2, w_3])
-          print('base_matrix:', base_matrix)
+#          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)
+#          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)
+#              print('new_point:', new_point)
               points_2d.append(new_point)
   #            points_2d.append([new_point[0], new_point[1]])
 
@@ -275,8 +275,8 @@ class Surface:
             points_2d.append([point[0], point[1], 0])
 
         polygon_2d = pn.Polygon(np.array(points_2d))
-        print('2D:')
-        print(polygon_2d.points)
+#        print('2D:')
+#        print(polygon_2d.points)
         # 2. calculate area:
         area = 0
         for i in range(0, len(polygon_2d.points)-1):
@@ -287,7 +287,7 @@ class Surface:
         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)
-      print(self._area)
+#      print(self._area)
     return self._area
 
   def _is_almost_same_terrain(self, terrain_points, ground_points):