diff --git a/src/Model/Geometry/PointXYZ.py b/src/Model/Geometry/PointXYZ.py
index 0db9ff3de26f7f5e7e83a662504257ecb4bd023c..9dd0a6e90ab13e055730b52f8177d512eb6e4367 100644
--- a/src/Model/Geometry/PointXYZ.py
+++ b/src/Model/Geometry/PointXYZ.py
@@ -204,3 +204,11 @@ class PointXYZ(Point, SQLSubModel):
             Euclidean distance between the two points
         """
         return dist((p1.x, p1.y, p1.z), (p2.x, p2.y, p2.z))
+
+    @staticmethod
+    def areatriangle3d(p1, p2, p3):
+        a = PointXYZ.distance(p1, p2)
+        b = PointXYZ.distance(p2, p3)
+        c = PointXYZ.distance(p3, p1)
+        s = (a + b + c) / 2
+        return (s*(s-a) * (s-b)*(s-c)) ** 0.5
diff --git a/src/Model/Geometry/ProfileXYZ.py b/src/Model/Geometry/ProfileXYZ.py
index f3ebe88fb220b018332ac02ed7a1033b1629be8f..0df81d3043d1fe0defd6f3c3e2c4e5986f82c5ed 100644
--- a/src/Model/Geometry/ProfileXYZ.py
+++ b/src/Model/Geometry/ProfileXYZ.py
@@ -530,3 +530,36 @@ class ProfileXYZ(Profile, SQLSubModel):
             pt_right = self.point(self.number_points - 1)
 
         return pt_left, pt_right
+
+    def purge(self, np_purge):
+        """
+        Remove points to keep at most np_purge points.
+        """
+
+        if (self.nb_points <= np_purge): return
+
+        nb_named = 2 # we consider the first and last point as named
+        area = [0.0]
+        for i in range(1, self.nb_points-1):
+            if self.point(i).point_is_named():
+                area.append(9999999.999)
+                nb_named += 1
+                print(self.point(i).name.strip())
+            else:
+                area.append(PointXYZ.areatriangle3d(self.point(i-1),self.point(i),self.point(i+1)))
+        area.append(0.0)
+        print(area)
+
+        while (self.nb_points > max(np_purge, nb_named)):
+            to_rm = np.argmin(area[1:self.nb_points-1])+1
+            print('to rm = ', to_rm)
+            self.delete_i([to_rm])
+            area.pop(to_rm)
+            for i in [to_rm-1, to_rm]:
+                if (i == 0): continue
+                if (i == self.nb_points - 1): continue
+                if self.point(i).point_is_named():
+                    area[i] = 9999999.999
+                else:
+                    area[i] = PointXYZ.areatriangle3d(self.point(i-1),self.point(i),self.point(i+1))
+                #self._study.river.edges()[0].reach.profile(0).purge(5)