From 596cd24cc33977197181b11017b3674896b35603 Mon Sep 17 00:00:00 2001
From: Pierre-Antoine Rouby <pierre-antoine.rouby@inrae.fr>
Date: Mon, 21 Aug 2023 16:18:28 +0200
Subject: [PATCH] Sediment: Add GRA file reading.
---
src/Model/Results/River/River.py | 3 +
src/Solver/Mage.py | 162 ++++++++++++++++++++++++++++++-
2 files changed, 163 insertions(+), 2 deletions(-)
diff --git a/src/Model/Results/River/River.py b/src/Model/Results/River/River.py
index 5ae1e90f..9f72f731 100644
--- a/src/Model/Results/River/River.py
+++ b/src/Model/Results/River/River.py
@@ -98,6 +98,9 @@ class River(object):
def reach(self, id):
return self._reachs[id]
+ def has_reach(self, id):
+ return 0 <= id < len(self._reachs)
+
def add(self, reach_id):
reachs = self._study.river.enable_edges()
diff --git a/src/Solver/Mage.py b/src/Solver/Mage.py
index 5e7909d7..9bc9331e 100644
--- a/src/Solver/Mage.py
+++ b/src/Solver/Mage.py
@@ -101,7 +101,6 @@ class Mage(AbstractSolver):
l = super(Mage, self).cmd_args(study)
l.append("-r")
- l.append("-fp=1")
return l
@@ -475,7 +474,6 @@ class Mage8(Mage):
f.write(f"{name} {value}\n")
-
return files
@timer
@@ -555,6 +553,7 @@ class Mage8(Mage):
# RESULTS #
###########
+ @timer
def read_bin(self, study, repertory, results, qlog = None):
logger.info(f"read_bin: Start ...")
@@ -675,3 +674,162 @@ class Mage8(Mage):
logger.debug(reachs[0].profiles[0]._data)
results.set("timestamps", ts)
logger.info(f"read_bin: ... end with {len(ts)} timestamp read")
+
+ @timer
+ def read_gra(self, study, repertory, results, qlog = None):
+ # if not study.river.has_sediment():
+ # return
+
+ logger.info(f"read_gra: Start ...")
+
+ with mage_file_open(os.path.join(repertory, f"0.GRA"), "r") as f:
+ newline = lambda: np.fromfile(f, dtype=np.int32, count=1)
+ endline = lambda: np.fromfile(f, dtype=np.int32, count=1)
+
+ read_int = lambda size: np.fromfile(f, dtype=np.int32, count=size)
+ read_float = lambda size: np.fromfile(f, dtype=np.float32, count=size)
+ read_float64 = lambda size: np.fromfile(f, dtype=np.float64, count=size)
+
+ # Meta data (1st line)
+ newline()
+
+ data = read_int(3)
+
+ nb_reach = data[0]
+ nb_profile = data[1]
+ mage_version = data[2]
+
+ logger.debug(f"read_gra: nb_reach = {nb_reach}")
+ logger.debug(f"read_gra: nb_profile = {nb_profile}")
+ logger.debug(f"read_gra: mage_version = {mage_version}")
+
+ if mage_version <= 80:
+ msg = (
+ "Read GRA files: " +
+ f"Possible incompatible mage version '{mage_version}', " +
+ "please check your solver configuration..."
+ )
+ logger.warning(msg)
+
+ if qlog is not None:
+ qlog.put("[WARNING] " + msg)
+
+ results.set("gra_solver_version", f"Mage8 ({mage_version})")
+ results.set("gra_nb_reach", f"{nb_reach}")
+ results.set("gra_nb_profile", f"{nb_profile}")
+
+ endline()
+
+ # Reach information (2nd line)
+ newline()
+
+ reachs = []
+ iprofiles = {}
+ reach_offset = {}
+
+ data = read_int(2*nb_reach)
+
+ for i in range(nb_reach):
+ # Get results reach to reach list
+ r = results.river.reach(i)
+ reachs.append(r)
+
+ # ID of first and last reach profiles
+ i1 = data[2*i] - 1
+ i2 = data[2*i+1] - 1
+
+ # Add profile id correspondance to reach
+ key = (i1, i2)
+ iprofiles[key] = r
+
+ # Profile ID offset
+ reach_offset[r] = i1
+
+ logger.debug(f"read_gra: iprofiles = {iprofiles}")
+
+ endline()
+
+ # X (3rd line)
+ newline()
+ _ = read_float(nb_profile)
+ endline()
+
+ # npts (4th line)
+ newline()
+ _ = read_int(nb_profile)
+ endline()
+
+ # Data
+ ip_to_r = lambda i: iprofiles[
+ next(
+ filter(
+ lambda k: k[0] <= i <= k[1],
+ iprofiles
+ )
+ )
+ ]
+ ip_to_ri = lambda r, i: i - reach_offset[r]
+
+ ts = set()
+ ind = 0
+ end = False
+
+ newline()
+ while not end:
+ n = read_int(1)[0]
+ timestamp = read_float64(1)[0]
+
+ logger.debug(f"read_gra: timestamp = {timestamp} sec")
+ ts.add(timestamp)
+
+ endline()
+
+ for i in range(n):
+ newline()
+ nsl = read_int(1)[0]
+ endline()
+
+ # Get current profile id
+ reach = ip_to_r(i)
+ ri = ip_to_ri(reach, i)
+
+ if nsl > 1:
+ logger.warning(
+ "read_gra: " +
+ "Multiple sediment layers for one profile " +
+ "is not implemented yet..."
+ )
+
+ for j in range(nsl):
+ newline()
+ nl = read_int(1)[0]
+ endline()
+
+ sl = []
+
+ for k in range(nl):
+ newline()
+ data = read_float64(3)
+ endline()
+
+ h = data[0]
+ d50 = data[1]
+ sigma = data[2]
+
+ sl.append((h, d50, sigma))
+
+ reach.set(ri, timestamp, "sl", sl)
+
+ ind += 1
+ end = newline().size <= 0
+
+ logger.debug(reachs[0].profiles[0]._data)
+ results.set("timestamps", ts)
+ logger.info(f"read_gra: ... end with {len(ts)} timestamp read")
+
+ @timer
+ def results(self, study, repertory, qlog = None):
+ results = super(Mage8, self).results(study, repertory, qlog)
+ self.read_gra(study, repertory, results, qlog)
+
+ return results
--
GitLab