add a basic st_to_m method

src/Model/Geometry/Profile.py

@@ -150,6 +150,19 @@ class Profile(object):
         self._status.modified()
 
 
+    def add_point(self, point:Point):
+        """append point at end of profile.
+
+        Args:
+            point: The point.
+
+        Returns:
+            Nothing.
+        """
+        self._points.append(point)
+        self._status.modified()
+
+
     def delete(self, indexes: int):
         """Delete points at index
 

src/Model/Geometry/ProfileXYZ.py

@@ -271,3 +271,28 @@ class ProfileXYZ(Profile):
                 constant = station[index_profile_z_min[0]]
 
             return list(map(lambda s: s - constant, station))
+
+    @timer
+    def length(self) -> list:
+        l = []
+        x = self.x()
+        y = self.y()
+        z = self.z()
+        for i in len(self.points)-1:
+            dx = x[i+1] - x[i]
+            dy = y[i+1] - y[i]
+            dz = z[i+1] - z[i]
+            l.append( np.sqrt(dx*dx + dy*dy + dz*dz) )
+        return l
+
+
+    @timer
+    def cumul_length(self) -> list:
+        lcum = []
+        lc  = 0.0
+        l = self.length()
+        for i in l:
+            lc = lc + i
+            lcum.append(lc)
+        return lcum
+

src/Model/Geometry/Reach.py

@@ -459,7 +459,7 @@ class Reach:
     # Interpolate
 
     @timer
-    def interpolate_step(self, lim: list[int], guideline1: str, guideline2: str, step: float, longmi: int, linear: bool)
+    def interpolate_step(self, lim: list, guideline1: str, guideline2: str, step: float, longmi: int, linear: bool):
         """interpolate with a given step between profile lim[0] and profile lim[1]
         Args:
             lim: index of profiles for interpolation
@@ -476,9 +476,9 @@ class Reach:
         """
         if len(lim) != 2:
             print('TODO: CRITICAL ERROR!') # error for now, only one interval
-        if lim[0] < 0 .or. lim[0] >= lim[1] .or. lim[1] >= len(self._profiles):
+        if lim[0] < 0 or lim[0] >= lim[1] or lim[1] >= self.number_profiles:
             print('TODO: CRITICAL ERROR!')
-        if lim[1] < 0 .or. lim[0] > lim[1]:
+        if lim[1] < 0 or lim[0] > lim[1]:
             print('TODO: CRITICAL ERROR!')
         comp, inconmp = self.compute_guidelines()
         if not guideline1 in comp:
@@ -487,29 +487,29 @@ class Reach:
             print('TODO: CRITICAL ERROR!')
 
         nbsti = [] # number of profiles to add between self._profile[i] and self._profile[i+1]
-        for i in range(lim[0])
+        for i in range(lim[0]):
             nbsti.append(0)
         for i in range(lim[0],lim[1]-1):
-                points = self.profile(i).points
+            points = self.profile(i).points
-                point11 = next(filter(lambda p: p.name == guideline1, points))
+            point11 = next(filter(lambda p: p.name == guideline1, points))
-                point12 = next(filter(lambda p: p.name == guideline2, points))
+            point12 = next(filter(lambda p: p.name == guideline2, points))
-                points = self.profile(i+1).points
+            points = self.profile(i+1).points
-                point21 = next(filter(lambda p: p.name == guideline1, points))
+            point21 = next(filter(lambda p: p.name == guideline1, points))
-                point22 = next(filter(lambda p: p.name == guideline2, points))
+            point22 = next(filter(lambda p: p.name == guideline2, points))
-                d1 = point11.dist(point21)
+            d1 = point11.dist(point21)
-                d2 = point12.dist(point22) # warning: 3D distance
+            d2 = point12.dist(point22) # warning: 3D distance
 
-            if (longmi == 1) then
+            if (longmi == 1):
                 d = max(d1, d2)
-            else if (longmi == 2) then
+            elif (longmi == 2):
                 d = min(d1, d2)
-            else
+            else:
                 d = 0.5*(d1 + d2)
 
             nbsti[i] = nint(d / pas) - 1
-            if(nbsti[i].lt.0) nbsti[i] = 0
+            if(nbsti[i] < 0): nbsti[i] = 0
 
-        for i in range(lim[1]-1,len(self._profiles)-1)
+        for i in range(lim[1]-1,self.number_profiles-1):
             nbsti.append(0)
 
         return self.interpolate_nprof(nbsti, lineaire)
@@ -547,10 +547,10 @@ class Reach:
         sections_out = []
         kps = self.get_kp()
         for isect_in in range(len(nbsti)):
-            sections_out.append(self.section(isect_in))
+            sections_out.append(self.profiles(isect_in))
             dkp = kps[isect_in + 1] - kps[isect_in]
-            sect_in1 = self.section(isect_in)
+            sect_in1 = self.profiles(isect_in)
-            sect_in2 = self.section(isect_in+1)
+            sect_in2 = self.profiles(isect_in+1)
             for isect_loc in range(nbsti[isect_in]):
                 # RATIO entre les deux sections initiales
                 d=float(isect_loc)/float(nbsti(isect_in)+1)
@@ -565,11 +565,142 @@ class Reach:
                                               reach = self))
                 for i in range(len(x1)): # can do better...
                     # add the points to the created section
-                    sections_out[-1].add(PointXYZ(x = x1[i] + d*(x2[i] - x1[i]),
+                    sections_out[-1].add_point(PointXYZ(x = x1[i] + d*(x2[i] - x1[i]),
-                                                  y = y1[i] + d*(y2[i] - y1[i]),
+                                                        y = y1[i] + d*(y2[i] - y1[i]),
-                                                  z = z1[i] + d*(z2[i] - z1[i]),
+                                                        z = z1[i] + d*(z2[i] - z1[i]),
-                                                  name = tags[i]
+                                                        name = tags[i]))
-                                                 )
         # don't forget the last one:
-        sections_out.append(self.section(self.number_profiles()-1))
+        sections_out.append(self.profiles[self.number_profiles-1])
         return sections_out
+
+    @timer
+    def st_to_m(self):
+        """convert a ST bief to a M bief (same number of points on every section)
+
+        Args:
+            nbsti: number of profiles to add between self._profile[i] and self._profile[i+1]
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+        isectmax = 0
+        npoint_max = 0
+        new_profiles = []
+        np = []
+        nprof = len(self.number_profiles)
+        # find the profile with more points
+        for isect in range(nprof):
+            new_profiles.append(self.profile(isect))
+            points = new_profiles[isect].points
+            np.append(len(points))
+            if(npoint_max < len(points)):
+                npoint_max = len(points)
+                isectmax = isect
+
+        # let's go ! up from isectmax
+        for isect in range(isectmax+1,nprof):
+            # TODO case new_profiles[isect].length() = 0.0
+            # TODO case new_profiles[isect-1].length() = 0.0
+            l = new_profiles[isect-1].length() # length of each segment
+            lcum = new_profiles[isect-1].cumul_length() # lcum[-1] is length
+            alpha = [i/lcum[-1] for i in lcum]
+            l = new_profiles[isect].length()
+            lcum = new_profiles[isect].cumul_length()
+            beta = [i/lcum[-1] for i in lcum]
+            for i in range(np[isect], npoints_max):
+                beta.append(1.0)
+
+
+
+            for k in range(np[isect]): # boucle sur les points à rajouter
+                trouve = False
+                for j0 in range(np[isect]-1):
+                    if(trouve):break
+                    undeplus=True
+                    for j1 in range(j0, np[isect]-1):
+                        if(undeplus):
+                            if(beta[j1].le.alpha[j1]):
+                                undeplus=False
+                                break
+
+                    if (undeplus):
+                        if (j0 == np[isect]-2):
+                            trouve=True
+                        else:
+                            j1=j0+1
+                            if (beta[j0]-alpha[j0] > alpha[j1]-beta[j0]): # beta(j0) > (alpha(j1)+alpha(j2))/2
+                                trouve = True
+
+
+                    if (trouve):
+                    # on rajoute un point
+                        ratio=(alpha[j0]-beta[j0-1])/(beta[j0]-beta[j0-1])
+                        if (ratio < 0.0):
+                            # on double le point a gauche
+                            new_profiles[isect].insert_point(j0,new_profiles[isect].points[j0])
+                            beta[j0]=beta[j0-1]
+                        else:
+                            p = PointXYZ(x = (1.0-ratio)*new_profiles[isect].points[j0].x + ratio*new_profiles[isect].points[j0].x,
+                                         y = (1.0-ratio)*new_profiles[isect].points[j0].y + ratio*new_profiles[isect].points[j0].y,
+                                         z = (1.0-ratio)*new_profiles[isect].points[j0].z + ratio*new_profiles[isect].points[j0].z,
+                                        )
+                            new_profiles[isect].insert_point(j0,p)
+
+            n = len(new_profiles[isect].points)
+            if n <= npoint_max: # we add the missing points at the end of the profile
+                for k in range(n,npoints_max): new_profiles[isect].add_point(new_profiles[isect].points[-1])
+
+
+        # let's go ! down from isectmax
+        for isect in range(isectmax-1,-1,-1):
+            # TODO case new_profiles[isect].length() = 0.0
+            # TODO case new_profiles[isect-1].length() = 0.0
+            l = new_profiles[isect+1].length() # length of each segment
+            lcum = new_profiles[isect+1].cumul_length() # lcum[-1] is length
+            alpha = [i/lcum[-1] for i in lcum]
+            l = new_profiles[isect].length()
+            lcum = new_profiles[isect].cumul_length()
+            beta = [i/lcum[-1] for i in lcum]
+            for i in range(np[isect], npoints_max):
+                beta.append(1.0)
+
+
+
+            for k in range(np[isect]): # boucle sur les points à rajouter
+                trouve = False
+                for j0 in range(np[isect]-1):
+                    if(trouve):break
+                    undeplus=True
+                    for j1 in range(j0, np[isect]-1):
+                        if(undeplus):
+                            if(beta[j1].le.alpha[j1]):
+                                undeplus=False
+                                break
+
+                    if (undeplus):
+                        if (j0 == np[isect]-2):
+                            trouve=True
+                        else:
+                            j1=j0+1
+                            if (beta[j0]-alpha[j0] > alpha[j1]-beta[j0]): # beta(j0) > (alpha(j1)+alpha(j2))/2
+                                trouve = True
+
+
+                    if (trouve):
+                    # on rajoute un point
+                        ratio=(alpha[j0]-beta[j0-1])/(beta[j0]-beta[j0-1])
+                        if (ratio < 0.0):
+                            # on double le point a gauche
+                            new_profiles[isect].insert_point(j0,new_profiles[isect].points[j0])
+                            beta[j0]=beta[j0-1]
+                        else:
+                            p = PointXYZ(x = (1.0-ratio)*new_profiles[isect].points[j0].x + ratio*new_profiles[isect].points[j0].x,
+                                         y = (1.0-ratio)*new_profiles[isect].points[j0].y + ratio*new_profiles[isect].points[j0].y,
+                                         z = (1.0-ratio)*new_profiles[isect].points[j0].z + ratio*new_profiles[isect].points[j0].z,
+                                        )
+                            new_profiles[isect].insert_point(j0,p)
+
+            n = len(new_profiles[isect].points)
+            if n <= npoint_max: # we add the missing points at the end of the profile
+                for k in range(n,npoints_max): new_profiles[isect].add_point(new_profiles[isect].points[-1])
+

src/Scripts/plot_3DST.py

@@ -32,7 +32,7 @@ def set_axes_equal(ax):
 
     ax.set_xlim3d([x_middle - plot_radius, x_middle + plot_radius])
     ax.set_ylim3d([y_middle - plot_radius, y_middle + plot_radius])
-    ax.set_zlim3d([z_middle - plot_radius, z_middle + plot_radius])
+    ax.set_zlim3d([z_middle - 0.1*plot_radius, z_middle + 0.1*plot_radius])
 
 st_file = sys.argv[1]
 my_reach = Reach(None)