Merge branch 'master' of gitlab-ssh.irstea.fr:theophile.terraz/pamhyr

src/Model/Geometry/Profile.py

@@ -345,3 +345,9 @@ class Profile(object):
     # Abstract method, must be implemented for in non abstract class
     def get_station(self):
         raise NotImplementedMethodeError(self, self.get_station)
+
+    # Computation method
+
+    # Abstract method for width approximation
+    def width_approximation(self):
+        raise NotImplementedMethodeError(self, self.width_approximation)

src/Model/Geometry/ProfileXYZ.py

@@ -287,6 +287,27 @@ class ProfileXYZ(Profile, SQLSubModel):
         except IndexError:
             raise IndexError(f"Invalid point index: {index}")
 
+    def get_point_by_name(self, name: str) -> PointXYZ:
+        """Get point by name.
+
+        Args:
+            name: Point name.
+
+        Returns:
+            The point.
+        """
+        try:
+            n_name = name.lower().strip()
+            return next(
+                filter(
+                    lambda p: p.name.lower().strip() == n_name,
+                    self.points
+                )
+            )
+        except Exception as e:
+            logger.debug(f"{e}")
+            raise IndexError(f"Invalid point name: {name}")
+
     def add(self):
         """Add a new PointXYZ to profile.
 
@@ -415,3 +436,9 @@ class ProfileXYZ(Profile, SQLSubModel):
                 constant = station[index_profile_z_min[0]]
 
             return list(map(lambda s: s - constant, station))
+
+    def width_approximation(self):
+        rg = self.get_point_by_name("rg")
+        rd = self.get_point_by_name("rd")
+
+        return abs(rg.dist(rd))

src/Model/Geometry/Reach.py

@@ -578,3 +578,103 @@ class Reach(SQLSubModel):
 
                 file_st.write("    999.9990     999.9990     999.9990")
                 file_st.write("\n")
+
+    def get_incline(self):
+        first = self.profile(0)
+        last = self.profile(len(self) - 1)
+
+        z_first = first.z_min()
+        kp_first = first.kp
+
+        z_last = last.z_min()
+        kp_last = last.kp
+
+        return (
+            (z_last - z_first)
+            /
+            (kp_last - kp_first)
+        )
+
+    def get_incline_mean(self):
+        profiles = self.profiles
+        previous = profiles[0]
+
+        incline_acc = 0
+
+        for profile in profiles[1:]:
+            z_first = previous.z_min()
+            kp_first = previous.kp
+
+            z_last = profile.z_min()
+            kp_last = profile.kp
+
+            incline_acc += (
+                (z_last - z_first)
+                /
+                (kp_last - kp_first)
+            )
+
+            previous = profile
+
+        return (incline_acc / (len(profiles) - 1))
+
+    def get_incline_median(self):
+        profiles = self.profiles
+        previous = profiles[0]
+
+        incline_acc = []
+
+        for profile in profiles[1:]:
+            z_first = previous.z_min()
+            kp_first = previous.kp
+
+            z_last = profile.z_min()
+            kp_last = profile.kp
+
+            incline_acc += [
+                (z_last - z_first)
+                /
+                (kp_last - kp_first)
+            ]
+
+            previous = profile
+
+        incline_acc.sort()
+        return incline_acc[round(len(profiles)/2)]
+
+    def get_incline_median_mean(self):
+        profiles = self.profiles
+        previous = profiles[0]
+
+        incline_acc = []
+
+        for profile in profiles[1:]:
+            z_first = previous.z_min()
+            kp_first = previous.kp
+
+            z_last = profile.z_min()
+            kp_last = profile.kp
+
+            incline_acc += [
+                (z_last - z_first)
+                /
+                (kp_last - kp_first)
+            ]
+
+            previous = profile
+
+        incline_acc.sort()
+
+        marge = round(len(incline_acc) * 0.1)
+        incline_set = incline_acc[marge:-marge]
+
+        logger.debug(f"+{incline_acc}")
+        logger.debug(f"-{incline_set}")
+
+        return (
+            reduce(
+                lambda acc, x: acc + x,
+                incline_set,
+                0.0
+            ) / (len(incline_set))
+        )

src/Model/InitialConditions/InitialConditions.py

@@ -16,12 +16,16 @@
 
 # -*- coding: utf-8 -*-
 
+import logging
+
 from copy import copy, deepcopy
 from tools import trace, timer
 from functools import reduce
 
 from Model.Tools.PamhyrDB import SQLSubModel
 
+logger = logging.getLogger()
+
 
 class Data(SQLSubModel):
     def __init__(self, name: str = "",
@@ -358,37 +362,79 @@ class InitialConditions(SQLSubModel):
         profiles = self._reach.reach.profiles.copy()
         self._sort_by_z_and_kp(profiles)
 
-        prev = None
+        incline = self._reach.reach.get_incline_median_mean()
+        logger.debug(f"incline = {incline}")
+
+        previous_elevation = -99999.99
         for profile in profiles:
+            width = profile.width_approximation()
+            strickler = 25
+            discharge = (
+                ((width * 0.8)
+                 * strickler
+                 * (height ** (5/3))
+                 * (abs(incline) ** (0.5)))
+            )
+
+            elevation= max(
+                profile.z_min() + height,
+                previous_elevation
+            )
+
+            logger.debug(f"({profile.kp}):")
+            logger.debug(f"  width  = {width}")
+            logger.debug(f"  strickler = {strickler}")
+            logger.debug(f"  discharge = {discharge}")
+
             new = Data(reach=self._reach, status=self._status)
             new["kp"] = profile.kp
+            new["discharge"] = discharge
 
-            if prev is None:
-                new["elevation"] = profile.z_min() + height
-            else:
-                new["elevation"] = max(
-                    profile.z_min() + height,
-                    prev["elevation"]
-                )
+            new["elevation"] = elevation
 
             self._data.append(new)
-            prev = new
+            previous_elevation = elevation
+
 
-        is_reverse = False
-        if profiles[0].kp > profiles[-1].kp:
-            is_reverse = True
+    def generate_discharge(self, discharge: float):
+        self._data = []
 
-        self._data.sort(
-            reverse=not is_reverse,
-            key=lambda d: d['kp']
+        self._generate_height_estimation_from_discharge(
+            discharge
         )
 
-    def generate_discharge(self, discharge: float):
-        self._new = []
+    def _generate_height_estimation_from_discharge(self, discharge: float):
+        profiles = self._reach.reach.profiles.copy()
+        self._sort_by_z_and_kp(profiles)
 
-        for d in self._data:
-            n = d.copy()
-            n['discharge'] = discharge
-            self._new.append(n)
+        previous_elevation = -99999.99
+
+        incline = self._reach.reach.get_incline_median_mean()
+        logger.debug(f"incline = {incline}")
+
+        for profile in profiles:
+            width = profile.width_approximation()
+            strickler = 25
+            height = (
+                discharge
+                /
+                ((width * 0.8) * strickler * (abs(incline) ** (0.5)))
+            ) ** (0.6)
+
+            elevation= max(
+                profile.z_min() + height,
+                previous_elevation
+            )
+
+            logger.debug(f"({profile.kp}):")
+            logger.debug(f"  width  = {width}")
+            logger.debug(f"  strickler = {strickler}")
+            logger.debug(f"  hieght = {height}")
+
+            new = Data(reach=self._reach, status=self._status)
+            new["kp"] = profile.kp
+            new["discharge"] = discharge
+            new["elevation"] = elevation
 
-        self._data = self._new
+            previous_elevation = elevation
+            self._data.append(new)