diff --git a/src/Model/Geometry/ProfileXYZ.py b/src/Model/Geometry/ProfileXYZ.py
index 82ef682c..fe517fe6 100644
--- a/src/Model/Geometry/ProfileXYZ.py
+++ b/src/Model/Geometry/ProfileXYZ.py
@@ -749,8 +749,8 @@ class ProfileXYZ(Profile, SQLSubModel):
             if zz[i] >= z and zz[i+1] < z:
                 y = np.interp(
                     z,
-                    [zz[i], zz[i+1]],
-                    [station[i], station[i+1]]
+                    [zz[i+1], zz[i]],
+                    [station[i+1], station[i]]
                 )
                 line.append([y, z])
 
@@ -839,11 +839,13 @@ class ProfileXYZ(Profile, SQLSubModel):
 
         for i in range(self.number_points-1):
             if zz[i] > z and zz[i+1] <= z:
-                y = np.interp(
-                    z,
-                    [zz[i], zz[i+1]],
-                    [station[i], station[i+1]]
-                )
+                fact = (z - zz[i]) / (zz[i+1] - zz[i])
+                y = station[i] + fact * (station[i+1] - station[i])
+                #y = np.interp(
+                    #z,
+                    #[zz[i+1], zz[i]],
+                    #[station[i+1], station[i]]
+                #)
                 start.append(y)
 
         end = []
@@ -852,11 +854,13 @@ class ProfileXYZ(Profile, SQLSubModel):
 
         for i in reversed(range(self.number_points-1)):
             if zz[i] <= z and zz[i+1] > z:
-                y = np.interp(
-                    z,
-                    [zz[i], zz[i+1]],
-                    [station[i], station[i+1]]
-                )
+                fact = (z - zz[i]) / (zz[i+1] - zz[i])
+                y = station[i] + fact * (station[i+1] - station[i])
+                #y = np.interp(
+                    #z,
+                    #[zz[i], zz[i+1]],
+                    #[station[i], station[i+1]]
+                #)
                 end.append(y)
 
         if len(start) != len(end):
@@ -887,33 +891,39 @@ class ProfileXYZ(Profile, SQLSubModel):
 
         # Interpolate points at river left side
         if (i_left > 0):
-            x = np.interp(
-                z,
-                [self.point(i_left).z, self.point(i_left - 1).z],
-                [self.point(i_left).x, self.point(i_left - 1).x]
-            )
-            y = np.interp(
-                z,
-                [self.point(i_left).z, self.point(i_left - 1).z],
-                [self.point(i_left).y, self.point(i_left - 1).y]
-            )
-            pt_left = PointXYZ(x, y, z, name="wl_left")
+            fact = (z - self.point(i_left).z) / (self.point(i_left - 1).z - self.point(i_left).z)
+            x = self.point(i_left).x + fact * (self.point(i_left - 1).x - self.point(i_left).x)
+            y = self.point(i_left).y + fact * (self.point(i_left - 1).y - self.point(i_left).y)
+            #x = np.interp(
+                #z,
+                #[self.point(i_left - 1).z, self.point(i_left).z],
+                #[self.point(i_left - 1).x, self.point(i_left).x]
+            #)
+            #y = np.interp(
+                #z,
+                #[self.point(i_left).z, self.point(i_left - 1).z],
+                #[self.point(i_left).y, self.point(i_left - 1).y]
+            #)
+            pt_left = PointXYZ(x=x, y=y, z=z, name="wl_left")
         else:
             pt_left = self.point(0)
 
         # Interpolate points at river right side
         if (i_right < self.number_points - 1):
-            x = np.interp(
-                z,
-                [self.point(i_right).z, self.point(i_right + 1).z],
-                [self.point(i_right).x, self.point(i_right + 1).x]
-            )
-            y = np.interp(
-                z,
-                [self.point(i_right).z, self.point(i_right + 1).z],
-                [self.point(i_right).y, self.point(i_right + 1).y]
-            )
-            pt_right = PointXYZ(x, y, z, name="wl_right")
+            fact = (z - self.point(i_right).z) / (self.point(i_right + 1).z - self.point(i_right).z)
+            x = self.point(i_right).x + fact * (self.point(i_right + 1).x - self.point(i_right).x)
+            y = self.point(i_right).y + fact * (self.point(i_right + 1).y - self.point(i_right).y)
+            #x = np.interp(
+                #z,
+                #[self.point(i_right).z, self.point(i_right + 1).z],
+                #[self.point(i_right).x, self.point(i_right + 1).x]
+            #)
+            #y = np.interp(
+                #z,
+                #[self.point(i_right).z, self.point(i_right + 1).z],
+                #[self.point(i_right).y, self.point(i_right + 1).y]
+            #)
+            pt_right = PointXYZ(x=x, y=y, z=z, name="wl_right")
         else:
             pt_right = self.point(self.number_points - 1)