"""
Main class controlling the wrapper.
This class is managing three main components:
1. Reading/writing inputs and output files
We used a generic formalation for the function's names based on the CATHY files.
Example:
-------
In order to update the soil file:
`update_soil(arguments)`
In order to update the atmbc file:
`update_atmbc(arguments)`
Note that the name of the arguments are similar to the name of the variable to update
In order to update the parm file with a new minimum time step DTMIN:
`update_parm(DTMIN=1e3)`
The update function is composed of three main actions:
- Set the defaut parameters (it is actually not reading the input file to set the defaut parameters)
- replace values by the new ones introduced via the function argument
- write the new file (it overwrites the existing file)
Note also that updating an input file may affect the CATHY.H control file.
When needed, the code takes care of updating values of the CATHY.H file retroactivelly.
Remenber to update all your prepo files before calling `run_preprocessor()`, and
all you input files before calling `run_processor()`
2. Compiling the fortran files
3. Running executable
Once all the files were updated, the `run_preprocessor()` or `run_processor()`
are taking care of recompiling the source files via bash cmd.
"""
import json
import os
import sys
import warnings
from pathlib import Path
import numpy as np
from git import Repo # In order to fetch the file directly from the repo
import tempfile
import pyCATHY.meshtools as mt
from pyCATHY.importers import cathy_inputs as in_CT
from pyCATHY.importers import cathy_outputs as out_CT
from pyCATHY.plotters import cathy_plots as plt_CT
from pyCATHY.cathy_utils import dictObs_2pd
from pyCATHY.config_manager import VersionConfigManager
warnings.filterwarnings("ignore", category=DeprecationWarning)
import glob
import pickle
import re
import shutil
import subprocess
import time # measure time of computation
from collections import OrderedDict
from os import listdir
from os.path import isfile, join
import pandas as pd
import pyvista as pv # read .vtk files
import rich.console
from rich import print
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf
from scipy.ndimage import distance_transform_edt
# -----------------------------------------------------------------------------
[docs]
def subprocess_run_multi(pathexe_list):
"""
Run multiple exe files in parallel.
Multiprocessor functions outside main CATHY object
# https://stackoverflow.com/questions/44144584/typeerror-cant-pickle-thread-lock-objects
"""
print(f"x= {pathexe_list}, PID = {os.getpid()}")
os.chdir(pathexe_list)
callexe = "./" + "cathy"
with warnings.catch_warnings():
warnings.simplefilter("ignore")
p = subprocess.run(
[callexe],
text=True,
# capture_output=True
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
)
return p
def make_console(verbose):
"""
Start up the :class:`rich.console.Console` instance we'll use.
Parameters
----------
verbose : bool
Whether or not to print status messages to stderr.
"""
return rich.console.Console(stderr=True, quiet=not verbose)
# -----------------------------------------------------------------------------
[docs]
class CATHY:
"""
Main CATHY object.
When instantiated it creates the tree project directories with 'prj_name' as root folder.
The src files are fetched from the online repository if not existing
(note that it is possible to call a specific version).
"""
[docs]
def __init__(
self,
dirName=None,
prj_name="my_cathy_prj",
notebook=False,
version="1.0.0",
verbose=True,
**kwargs,
):
"""
Create CATHY object.
..note:
All variables in CAPITAL LETTERS use the semantic from the CATHY legacy fortran codes,
while the others are created exclusively for the wrapper.
..note:
All file variable are self dictionnary objects; example: soil file is contained in self.soil
"""
self.t0 = time.time() # executation time estimate
self.console = make_console(verbose)
self.console.print(":checkered_flag: [b]Initiate CATHY object[/b]")
self.verbose = False
# flag for notebook execution
# ---------------------------------------------------------------------
self.notebook = notebook # flag if the script is run in a notebook
# version tag used to select version-specific behavior/config
# (see pyCATHY.config_manager.VersionConfigManager and
# version_config.py) -- e.g. "SCF_variable" enables the
# per-vegetation-type SCF soil-file format; "default"/anything
# else keeps the classic format unchanged.
# ---------------------------------------------------------------------
self.version = version
# create working and project dir
# ---------------------------------------------------------------------
if dirName is None:
self.workdir = os.path.join(os.getcwd())
else:
self.workdir = os.path.join(os.getcwd(), dirName)
# change working directory
if not os.path.exists(os.path.join(self.workdir)):
os.makedirs(self.workdir, exist_ok=True)
# os.chdir(self.workdir)
self.project_name = prj_name
if not os.path.exists(os.path.join(self.workdir, self.project_name)):
os.makedirs(os.path.join(self.workdir, self.project_name), exist_ok=True)
# convention for directory tree
# ---------------------------------------------------------------------
self.processor_name = "cathy"
self.input_dirname = "input"
self.output_dirname = "output"
# dict related to CATHY inputs files
# ---------------------------------------------------------------------
self.parm = {} # dict of parm input parameters
self.soil = {} # dict of soil input parameters
# self.zone = None # raster of zone input parameters
self.ic = {} # dict of ic input parameters
self.cathyH = {} # dict of cathyH C header parameters
# self.nudging = {'NUDN': 0} # Temporary
# THIS SHOULD BE REMOVED #
# dict related to Data assimilation
# ---------------------------------------------------------------------
self.DAFLAG = False # Flag to trigger data assimilation process
# (self.dict_obs is assimiliatiom times size)
self.dict_obs = OrderedDict() # dictionnary containing all the observation data
# (self.stacked_data_cov is data_size * assimiliatiom times size)
self.count_DA_cycle = None # counter of the Assimilation Cycle
# self.df_performance = pd.DataFrame() # pandas dataframe with performence evaluation metrics
# dict related to the mesh
# ---------------------------------------------------------------------
self.grid3d = {}
for key, value in kwargs.items():
# clear src files if required
# ---------------------------------------------------------------------
if key == "clear_src": # clear src files
if value == True:
if os.path.exists(
os.path.join(self.workdir, self.project_name, "src")
):
print("clear src files")
shutil.rmtree(
os.path.join(self.workdir, self.project_name, "src")
)
shutil.rmtree(
os.path.join(self.workdir, self.project_name, "input")
)
shutil.rmtree(
os.path.join(self.workdir, self.project_name, "prepro")
)
# fetch src files if not existing from Gitbucket repo
# ---------------------------------------------------------------------
if not os.path.exists(os.path.join(self.workdir, self.project_name, "src")):
self.console.print(":worried_face: [b]src files not found[/b]")
self.console.print("working directory is:" + str((self.workdir)))
if version == "1.0.0":
try:
# Repo.clone_from(
# "https://bitbucket.org/cathy1_0/cathy.git",
# os.path.join(self.workdir, self.project_name, "tmp_src"),
# branch="master",
#)
Repo.clone_from(
"https://github.com/CATHY-Org/CATHY_src_org.git",
os.path.join(self.workdir, self.project_name, "tmp_src"),
branch="main",
)
self.console.print(":inbox_tray: [b]Fetch cathy src files[/b]")
shutil.move(
os.path.join(self.workdir, self.project_name, "tmp_src/src"),
os.path.join(self.workdir, self.project_name, "src"),
)
pathsrc = os.path.join(
self.workdir, self.project_name, "tmp_src/runs/weilletal/"
)
onlyfiles = [
f for f in listdir(pathsrc) if isfile(join(pathsrc, f))
]
for (
file
) in (
onlyfiles
):
shutil.move(
os.path.join(pathsrc, file),
os.path.join(self.workdir, self.project_name, file),
)
except:
print("no internet connection to fetch the files")
sys.exit()
pass
if version == "G. Manoli":
print("fetch cathy G. Manoli src files")
path_manoli = "/home/ben/Documents/CATHY/CathyGitbucket/Test_Gabriele/1_Gabriele_Piante_NON_modificato/CATHY_RWU_ABL_1D/"
shutil.copytree(
path_manoli, os.path.join(self.workdir, self.project_name, "src")
)
if version == "SCF_variable":
print("fetch cathy SCF_variable src files")
dest = os.path.join(self.workdir, self.project_name)
with tempfile.TemporaryDirectory() as tmp:
Repo.clone_from(
"https://github.com/CATHY-Org/CATHY_SCF_LAI.git",
tmp,
branch="main",
depth=1,
)
shutil.copytree(
tmp,
dest,
dirs_exist_ok=True,
ignore=shutil.ignore_patterns(".git", ".gitattributes", ".gitignore"),
)
#path_SCF_variable = "/home/z0272571a@campus.csic.es/Nextcloud/BenCSIC/Codes/CATHY-org/CATHY_src_SCF_variable_dev"
#shutil.copytree(
# path_SCF_variable,
# os.path.join(self.workdir, self.project_name),
# dirs_exist_ok=True,
#)
if not os.path.exists(os.path.join(self.workdir, self.project_name, "prepro")):
self.console.print(
":inbox_tray: [b]Fetch cathy prepro src files[/b]"
)
shutil.move(
os.path.join(
self.workdir,
self.project_name,
"tmp_src/runs/weilletal/prepro",
),
os.path.join(self.workdir, self.project_name, "prepro"),
)
# os.remove(os.path.join(
# self.workdir, self.project_name, "prepro") + '/dem' )
if not os.path.exists(os.path.join(self.workdir, self.project_name, "input")):
self.console.print(":inbox_tray: [b]Fetch cathy input files[/b]")
shutil.move(
os.path.join(
self.workdir,
self.project_name,
"tmp_src/runs/weilletal/input",
),
os.path.join(self.workdir, self.project_name, "input"),
)
if not os.path.exists(os.path.join(self.workdir, self.project_name, "output")):
self.create_output(output_dirname="output")
shutil.rmtree(
os.path.join(self.workdir, self.project_name, "tmp_src")
)
os.remove(os.path.join(self.workdir, self.project_name, "readme.txt"))
# clear output files if required
# ---------------------------------------------------------------------
for key, value in kwargs.items():
if key == "clear_outputs":
if value == True:
if not os.path.exists(
os.path.join(self.workdir, self.project_name, "output")
):
self.create_output(output_dirname="output")
else:
shutil.rmtree(
os.path.join(self.workdir, self.project_name, "output")
)
shutil.rmtree(
os.path.join(self.workdir, self.project_name, "vtk")
)
self.create_output(output_dirname="output")
pass
# --------------------------------------------------------------------------- #
# replace console cmd by a subprocess call
# --------------------------------------------------------------------------- #
def run_preprocessor(self, KeepOutlet=True, verbose=False, **kwargs):
"""
Run cppp.exe
1. updates cathy parameters prepo
2. recompile the source files prepo
3. run the preprocessor using bash cmd
Running the executable file has allowed to generate a complete set of
files describing physiographic features of the drainage system, as shown in Table 2.
"""
self.console.print(":cooking: [b]gfortran compilation[/b]")
for loopi in range(2): # run it twice (to avoid the first error)
os.chdir(os.path.join(self.workdir, self.project_name, "prepro/src/"))
# clean all files compiled
for file in glob.glob("*.o"):
os.remove(file)
# if self.notebook==False:
# bashCommand = "gfortran -O -o pycppp mpar.f90 mbbio.f90 wbb_sr.f90 csort.f90 qsort.f90 depit.f90 cca.f90 smean.f90 dsf.f90 facet.f90 hg.f90 mrbb_sr.f90 bb2shp_sr.f90 shape.f90 dbase.f90 streamer.f90 cppp.f90"
# bashCommand = "gfortran -O -mcmodel=large -o pycppp mpar.f90 mbbio.f90 wbb_sr.f90 csort.f90 qsort.f90 depit.f90 cca.f90 smean.f90 dsf.f90 facet.f90 hg.f90 mrbb_sr.f90 bb2shp_sr.f90 shape.f90 dbase.f90 streamer.f90 cppp.f90"
bashCommand = "gfortran -O -mcmodel=large -fPIC -o pycppp mpar.f90 mbbio.f90 wbb_sr.f90 csort.f90 qsort.f90 depit.f90 cca.f90 smean.f90 dsf.f90 facet.f90 hg.f90 mrbb_sr.f90 bb2shp_sr.f90 shape.f90 dbase.f90 streamer.f90 cppp.f90"
try:
p = os.system(bashCommand + "> /dev/null 2>&1")
# run it twice (to avoid the first error)
p = os.system(bashCommand + "> /dev/null 2>&1")
except:
print("bash cmd not recognized")
os.chdir(self.workdir)
# move to the directory where the source FORTRAN files are contained (cathy_main.f)
# shutil.move(
# os.path.join(self.workdir, self.project_name, "prepro/src/pycppp"),
# os.path.join(self.workdir, self.project_name, "prepro/pycppp"),
# )
src = Path(self.workdir) / self.project_name / "prepro/src/pycppp"
dst = Path(self.workdir) / self.project_name / "prepro/pycppp"
if dst.exists():
dst.unlink() # Delete the existing destination file
src.rename(dst)
self.console.print(":athletic_shoe: [b]Run preprocessor[/b]")
os.chdir(os.path.join(self.workdir, self.project_name, "prepro"))
bashcmd = "./pycppp"
my_data = "2\n0\n1\n" # user input data
with warnings.catch_warnings():
warnings.simplefilter("ignore")
p = subprocess.run(
[bashcmd],
text=True,
input=my_data,
# capture_output=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
)
if verbose:
print(p.stdout)
print(p.stderr)
if "DEM resolution is not a multiple of the rivulet spacing!" in p.stdout:
print("DEM resolution is not a multiple of the rivulet spacing!")
sys.exit()
if "catchment with more than one outlet cell!" in p.stdout:
print("catchment with more than one outlet cell!")
sys.exit()
if KeepOutlet == False:
print("remove outlet")
idoutlet = np.where(self.DEM == min(np.unique(self.DEM)))
self.DEM[idoutlet[0], idoutlet[1]] = max(np.unique(self.DEM))
if hasattr(self, 'DEM') is False:
DEM_mat, DEM_header = in_CT.read_dem(
os.path.join(self.workdir, self.project_name, "prepro/dem"),
os.path.join(self.workdir, self.project_name, "prepro/dtm_13.val"),
)
self.DEM = DEM_mat
# if hasattr(self, 'hapin') is False:
# self.update_cathyH(
# ROWMAX=self.hapin['N'],
# COLMAX=self.hapin['M']
# )
# if np.shape(self.DEM)[1]==self.hapin['N']:
# self.console.rule(''':warning: !transposing DEM dimension!
# This should be avoided in a future version
# :warning:''',
# style="yellow")
# self.update_prepo_inputs(N=self.hapin['M'],
# M=self.hapin['N'])
# self.update_zone()
os.chdir(self.workdir)
pass
# def recompileSrc(self, verbose=False):
# """
# Other option is self.run_processor(runProcess=False)
# """
# ti = time.time() # executation time estimate
# self.console.print(
# ":hammer_and_wrench: [b] Recompile src files[/b] ["
# + str(int(abs(ti - self.t0)))
# + "s]"
# )
# # gfortran_executable = 'gfortran.exe'
# # # Execute the where command to find the path to gfortran
# # where_command = ['where', gfortran_executable]
# # result = subprocess.run(where_command, stdout=subprocess.PIPE, shell=True, check=True)
# # # Get the path to gfortran from the output of the where command
# # gfortran_path = result.stdout.decode('utf-8').strip().split('\r\n')[0]
# # clean all files previously compiled
# for file in glob.glob("*.o"):
# os.remove(file)
# # list all the fortran files to compile and compile
# for file in glob.glob("*.f"):
# # filepath = Path(self.workdir) / self.project_name / 'src' / str(file)
# filepath = str(file)
# bashCommand = "gfortran -c " + str(filepath)
# # gfortran -c *.f
# # gfortran *.o -L\MinGW\lib -llapack -lblas -o cathy
# # print(bashCommand)
# with warnings.catch_warnings():
# warnings.simplefilter("ignore")
# process = subprocess.Popen(
# bashCommand.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE
# )
# # if verbose==True:
# # output, error = process.communicate()
# output, error = process.communicate()
# # list all the fortran compiled files to compile and run
# files = ""
# for file in glob.glob("*.o"):
# # files += " " + os.path.join(self.workdir, self.project_name, 'prepro/src/' + str(file))
# files += " " + str(file)
# # bashCommand = "gfortran" + files + " -llapack -lblas -o " + self.processor_name
# bashCommand = "gfortran -mcmodel=large -fPIC" + files + " -llapack -lblas -o " + self.processor_name
# ti = time.time() # executation time estimate
# self.console.print(
# ":cooking: [b]gfortran compilation[/b] ["
# + str(int(abs(ti - self.t0)))
# + "s]"
# )
# with warnings.catch_warnings():
# warnings.simplefilter("ignore")
# process = subprocess.Popen(
# bashCommand.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE
# )
# output, error = process.communicate()
# print(error)
# try:
# shutil.move(
# os.path.join(
# self.workdir, self.project_name, "src", self.processor_name
# ),
# os.path.join(self.workdir, self.project_name, self.processor_name),
# )
# except:
# self.console.print(":pensive_face: [b]Cannot find the new processsor[/b]")
# pass
def recompileSrc(self, verbose=False):
"""
Other option is self.run_processor(runProcess=False)
"""
ti = time.time()
self.console.print(
":hammer_and_wrench: [b] Recompile src files[/b] ["
+ str(int(abs(ti - self.t0)))
+ "s]"
)
# clean all files previously compiled
for file in glob.glob("*.o"):
os.remove(file)
# Step 1: compile each .f file to .o WITH large memory model flags
compile_errors = []
for file in glob.glob("*.f"):
filepath = str(file)
bashCommand = f"gfortran -O -mcmodel=large -fPIC -c {filepath}"
with warnings.catch_warnings():
warnings.simplefilter("ignore")
process = subprocess.Popen(
bashCommand.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE
)
output, error = process.communicate()
if process.returncode != 0:
compile_errors.append((filepath, error.decode()))
if verbose:
self.console.print(f":x: [red]Error compiling {filepath}[/red]")
print(error.decode())
if compile_errors:
self.console.print(
f":x: [red]{len(compile_errors)} file(s) failed to compile[/red]"
)
for fname, err in compile_errors:
self.console.print(f" - {fname}: {err}")
return
# Step 2: link all .o files WITH the same flags
files = " ".join(glob.glob("*.o"))
bashCommand = (
f"gfortran -mcmodel=large -fPIC {files} "
f"-llapack -lblas -o {self.processor_name}"
)
ti = time.time()
self.console.print(
":cooking: [b]gfortran compilation[/b] ["
+ str(int(abs(ti - self.t0)))
+ "s]"
)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
process = subprocess.Popen(
bashCommand.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE
)
output, error = process.communicate()
if process.returncode != 0:
self.console.print(":x: [red]Linker error:[/red]")
print(error.decode())
return
else:
self.console.print(":white_check_mark: [green]Compilation successful![/green]")
if verbose:
print(output.decode())
try:
shutil.move(
os.path.join(self.workdir, self.project_name, "src", self.processor_name),
os.path.join(self.workdir, self.project_name, self.processor_name),
)
except Exception as e:
self.console.print(
f":pensive_face: [b]Cannot find the new processor[/b]: {e}"
)
def run_processor(self, recompile=True, runProcess=True, verbose=False, **kwargs):
"""
Run cathy.exe
1. updates cathy parameters and cathyH based on **kwargs
2. recompile the source files (set False for notebook) and create the executable
3. run the processor using bash cmd if runProcess is True
Parameters
----------
recompile : Bool, optional
recomplile CATHY src files. The default is True.
runProcess : Bool, optional
exectute compile CATHY exe. The default is True.
verbose : Bool, optional
Output CATHY log. The default is False.
"""
# set parallel flag for mutirun (used for DA ensemble)
# --------------------------------------------------------------------
parallel = False
if "parallel" in kwargs:
parallel = kwargs["parallel"]
# update parm and cathyH
# --------------------------------------------------------------------
# VERY VERY IMPORTANT NEVER COMMENT !
self.check_DEM_versus_inputs() # to uncomment
# if len(kwargs)>0:
self.update_parm(**kwargs) # to uncomment
self.update_cathyH(**kwargs) ### to uncomment
# self.cathyH
if recompile:
# recompile
# --------------------------------------------------------------------
os.chdir(os.path.join(self.workdir, self.project_name, "src"))
self.recompileSrc()
if "path_CATHY_folder" in kwargs:
os.chdir(kwargs["path_CATHY_folder"])
else:
# back to root dir
os.chdir(os.path.join(self.workdir, self.project_name))
# some checkings
# --------------------------------------------------------------------
try:
os.path.exists(os.path.join("inputs"))
except OSError:
print("input folder missing")
sys.exit()
# check if output folder
if not os.path.exists("output"):
os.makedirs("output", exist_ok=True)
# check if vtk folder
if not os.path.exists("vtk"):
os.makedirs("vtk", exist_ok=True)
# run the processor
# --------------------------------------------------------------------
if runProcess:
# t0 = time.time() # executation time estimate
self.console.print(":athletic_shoe: [b]Run processor[/b]")
callexe = "./" + self.processor_name
# case of simple simulation
# ----------------------------------------------------------------
with warnings.catch_warnings():
warnings.simplefilter("ignore")
p = subprocess.run(
[callexe],
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
)
# process = subprocess.Popen(
# [callexe], stdout=subprocess.PIPE, stderr=subprocess.PIPE
# )
# # if verbose==True:
# # output, error = process.communicate()
# output, error = process.communicate()
# Redirecting output to /dev/null or NUL
# output_redirection = "> /dev/null 2>&1" if os.name != 'nt' else "> NUL 2>&1"
# # Execute the command without displaying the shell
# os.system(f"{callexe} {output_redirection}")
# result = os.popen(f"{callexe} {output_redirection}").read()
# p = subprocess.run(
# [callexe],
# stdout=subprocess.DEVNULL,
# stderr=subprocess.DEVNULL,
# )
# p = subprocess.Popen([callexe],
# )
# p = subprocess.Popen([callexe],
# stdout=subprocess.PIPE,
# stderr=subprocess.PIPE
# )
# warnings.simplefilter("ignore")
# p = subprocess.Popen([callexe], stdout=subprocess.PIPE, stderr=subprocess.PIPE)
# p = subprocess.run(
# [callexe],
# text=True,
# capture_output=True,
# check=True
# )
# print('verbose='+ str(verbose))
if verbose:
# stdout, stderr = p.communicate() # This will block until the process finishes
# print(result)
print(p.stdout)
# print(p.stderr)
os.chdir(os.path.join(self.workdir))
try:
self.grid3d = out_CT.read_grid3d(
os.path.join(self.workdir, self.project_name, 'output', 'grid3d')
)
except Exception as exc:
# Previously this was a bare `except: self.run_processor(IPRT1=3)`
# with no retry limit and no visibility into the actual
# failure. If output/grid3d is missing/malformed for a
# reason that doesn't self-heal, every retry reproduces
# the same failure and calls run_processor again, forever:
# an infinite "Update parm file / Recompile / Run
# processor" loop with the real error silently swallowed
# every single time (bare `except:` even catches
# KeyboardInterrupt/SystemExit). Retry exactly once (for a
# genuine one-off race), then surface the real exception
# instead of recursing further.
if getattr(self, "_grid3d_read_retried", False):
self._grid3d_read_retried = False
raise RuntimeError(
"Failed to read output/grid3d after retrying "
"run_processor(IPRT1=3) once. Not retrying again "
"to avoid an infinite recompile/run loop - see "
"the CATHY stdout/stderr above, and confirm "
"output/grid3d is actually produced by this run "
"configuration."
) from exc
self._grid3d_read_retried = True
self.run_processor(IPRT1=3)
else:
self._grid3d_read_retried = False
# if not "nnod3" in self.grid3d.keys():
# self.run_processor(IPRT1=3)
# computation time
# ----------------------------------------------------------------
# t1 = time.time()
# self.total_computation_time = t1 - t0
return
def create_output(self, output_dirname="output"):
"""
Create output directories
Parameters
----------
output_dirname : str, optional
Name of the dir. The default is "output".
Returns
-------
None.
"""
# create project_name/output folders
if not os.path.exists(os.path.join(self.project_name, output_dirname)):
os.mkdir(os.path.join(self.workdir, self.project_name, output_dirname))
if not os.path.exists(os.path.join(self.project_name, "vtk")):
os.mkdir(os.path.join(self.workdir, self.project_name, "vtk"))
# compute the processor (make clean make) using the Makefile
pass
def display_time_run(self):
return self.total
#%% Checkings
def check_DEM_versus_inputs(self):
"""
Check shape consistency between attributes and DEM
"""
# try:
if hasattr(self, "veg_map"):
if np.shape(self.veg_map) != np.shape(self.DEM):
self.console.print(
f"Inconsistent shapes between vegetation map {np.shape(self.veg_map) } and DEM {np.shape(self.DEM)} - need to update veg first"
)
# Previously this always fell back to a fabricated
# single-class (all-ones) raster here — which, on a fresh
# session where self.MAXVEG hadn't been set yet, silently
# baked MAXVEG=1 into CATHY.H *and* overwrote a real,
# multi-class root_map on disk with garbage. Prefer
# reloading the real root_map that's actually on disk (its
# own header carries the true shape) before ever
# fabricating placeholder data.
root_map_path = os.path.join(
self.workdir, self.project_name, self.input_dirname, "root_map"
)
if os.path.exists(root_map_path):
self.console.print(
"Reloading existing root_map from disk instead of "
"fabricating a placeholder vegetation map."
)
self.update_veg_map()
else:
self.console.print(
"[yellow]No root_map found on disk either — "
"falling back to a single placeholder vegetation "
"class. MAXVEG will be 1 until a real vegetation "
"map is supplied via update_veg_map().[/yellow]"
)
self.update_veg_map(indice_veg=np.ones(np.shape(self.DEM)))
#%%
# --------------------------------------------------------------------------- #
# update files
# --------------------------------------------------------------------------- #
def update_cathyH(self, verbose=False, **kwargs):
"""
Update CathyH file input based on **kwargs parameters.
Parameters ---------- **kwargs : all variable parameters included into cathyH. For instance
ROWMAX # maximum NROW, with NROW = number of rows in the DEM.
Returns ------- type Write a modified cathyH file based on new parameters parsed
"""
indent = ""
for kk, value in kwargs.items():
if kk == "indent":
indent = value
if verbose:
self.console.print(
indent + ":arrows_counterclockwise: [b]Update cathyH files[/b]"
)
CATHYH_file = open(
os.path.join(self.workdir, self.project_name, "src", "CATHY.H"), "r"
)
Lines0_109 = CATHYH_file.readlines()
CATHYH_file.close()
# check if hapin and parm exist
# create them if not existing
if hasattr(self, "hapin") is False:
self.update_prepo_inputs()
if hasattr(self, "dem_parameters") is False:
self.update_dem_parameters()
if hasattr(self, "MAXVEG") is False:
# update_cathyH() is reached from early in the pipeline
# (update_prepo_inputs -> update_zone -> update_parm), well
# before the project's real vegetation map is built — that
# happens later, via an explicit
# update_veg_map(veg_map, maxveg=...) call in the driving
# script. Previously this branch called self.update_veg_map()
# here, which tried to read/validate root_map (and mask it
# against prepro/dem) using whatever happens to be sitting in
# input/root_map at this early stage — for a freshly created
# project that's just the default template file, unrelated to
# this project's real DEM/vegetation data, and could be the
# wrong shape entirely. That produced everything from a wrong
# MAXVEG=1 (silently) to an IndexError/ValueError (once that
# path was hardened to check shapes) — always at the wrong
# moment, because the real answer genuinely isn't known yet.
# Leave a provisional placeholder here instead, and do no file
# I/O or DEM validation; self.veg_map is deliberately left
# unset so the later, authoritative update_veg_map(...) call
# is still recognized as the first *real* vegetation setup and
# correctly overwrites both root_map and CATHY.H once the true
# vegetation data actually exists.
self.MAXVEG = 1
if "NR" not in self.parm:
self.console.rule(
":warning: warning messages above :warning:", style="yellow"
)
self.console.print(
r"""
The parm dictionnary is empty
Falling back to defaults to update CATHYH
This can have consequences !!
""",
style="yellow",
)
self.console.rule("", style="yellow")
self.update_parm()
DEMRES = 1
ROWMAX = self.hapin["M"]
COLMAX = self.hapin["N"]
if len(self.cathyH) == 0:
self.cathyH = {
# "ROWMAX": self.hapin["N"], # maximum NROW, with NROW = number of rows in the DEM
# "COLMAX": self.hapin["M"], # maximum NCOL, with NCOL = number of columns in the DEM
"ROWMAX": ROWMAX, # maximum NROW, with NROW = number of rows in the DEM
"COLMAX": COLMAX, # maximum NCOL, with NCOL = number of columns in the DEM
# 'COARSE': ,
"MAXCEL": int(self.hapin["N"]) * int(self.hapin["M"]),
"MAXRES": 1,
"DEMRES": DEMRES,
# "NODMAX": (int(self.hapin["N"]) / DEMRES + 1) * (int(self.hapin["M"]) / DEMRES + 1),
# "NODMAX": 1e6,
"NTRMAX": 2* (int(self.hapin["N"]) * int(self.hapin["M"]))/ (DEMRES * DEMRES),
"NP2MAX": 1,
"MAXSTR": self.dem_parameters["nstr"],
"NPMAX": 1,
"NQMAX": 1,
"NSFMAX": 1,
"NNSFMX": 1,
"MAXNUDN": 1, # self.nudging["NUDN"],
"MAXNUDT": 1,
"MAXNUDC": 1,
# 'MAXNENS': ,
"MAXZON": self.dem_parameters["nzone"], # maximum NZONE, with NZONE = number of material types in the porous medium
"MAXTRM": self.hapin["N"]*self.hapin["M"]*self.dem_parameters["nstr"]*30,
"MAXIT": 30,
"NRMAX": self.parm["NR"]+1,
# maximum NPRT (ref. parm file), with NPRT = number of time values for detailed output
"MAXPRT": self.parm["NPRT"],
# maximum NUMVP (ref. parm input file), NUMVP = number of surface nodes for vertical profile output
"MAXVP": int(self.parm["NUMVP"]),
"N1MAX": self.dem_parameters[
"n1"
], # maximum N1 (it is good to have N1 ≤ 20), N1 = maximum number of element connections to a node
"MAXBOT": 1,
"INTBOT": 1,
"NTAB": 100,
"MAXQOUT": 1,
"MAXVTKPRT": self.parm["NPRT"],
# Is related to data assimilation (DA) (to do not considered if you do not have DA)
"MAXENNOUT": 52560,
"MAXVEG": self.MAXVEG, #
# "MAXVEG": 1, #
}
for kk, value in kwargs.items():
if kk in self.cathyH.keys():
self.cathyH[kk] = value
# self.cathyH
# self.cathyH['MAXCEL'] = int(self.cathyH["ROWMAX"]) * int(self.cathyH["COLMAX"])
# self.cathyH['NODMAX'] = int((int(self.cathyH["ROWMAX"]) / DEMRES + 1) * (int(self.cathyH["COLMAX"]) / DEMRES + 1))
# self.cathyH['NODMAX'] = int(300)
# NODMAX = (ROWMAX+1) * (COLMAX+1) — number of nodes (corners), not cells
self.cathyH['NODMAX'] = int((self.cathyH["ROWMAX"] + 1) * (self.cathyH["COLMAX"] + 1))
# self.cathyH['NODMAX'] = 1e6
# MAXCEL = ROWMAX * COLMAX — number of cells
self.cathyH['MAXCEL'] = int(self.cathyH["ROWMAX"] * self.cathyH["COLMAX"])
# self.cathyH['MAXCEL'] = 1e6
# NTRMAX = 2 * MAXCEL — 2 triangles per cell
self.cathyH['NTRMAX'] = int(2 * self.cathyH['MAXCEL'])
# self.cathyH['NTRMAX'] = int((2*self.cathyH['MAXCEL'])/ (DEMRES * DEMRES))
self.cathyH['MAXTRM'] = self.cathyH["ROWMAX"]*self.cathyH["COLMAX"]*self.dem_parameters["nstr"]*30
self.cathyH['NFACEMAX'] = self.cathyH['NODMAX']*5
# self.cathyH['MAXTRM'] = 1599000
# MAXTRM=1599000
# cathyH_laC = {
# # "ROWMAX": 247, # maximum NROW, with NROW = number of rows in the DEM
# # "COLMAX": 221, # maximum NCOL, with NCOL = number of columns in the DEM
# # # 'COARSE': ,
# # "MAXRES": 1,
# # "DEMRES": 1,
# "NODMAX": 28440,
# "NTRMAX": 2*27848,
# "MAXTRM": 2364075,
# "NRMAX": 1,
# }
# for k in cathyH_laC:
# self.cathyH[k] = cathyH_laC[k]
# ---------------------------------------------------------------------
# Hard safety net: verify self.cathyH["MAXVEG"] against what's
# actually on disk in root_map (the file datin.f itself reads to
# compute NVEG = MAX(VEG_TYPE) inside the valid DEM domain) before
# ever writing CATHY.H. Several code paths feed into MAXVEG (an
# explicit maxveg=, an auto-derive on a fresh session, a stale
# self.cathyH dict from a much earlier call) and any one of them
# getting it wrong previously meant the *first* symptom was a
# cryptic Fortran "Bad integer for item 1 in list input" crash deep
# into a run. Checking here instead means a mismatch is caught
# immediately, in Python, with the actual numbers involved.
#
# Only meaningful once a real vegetation map has actually been
# established for this project (self.veg_map set, via a genuine
# update_veg_map(...) call) — before that, root_map on disk may
# still be the project template (wrong shape/content, unrelated to
# this project's real DEM), and MAXVEG is a known, temporary
# placeholder that the later, authoritative call will overwrite.
# Checking against the template at that stage would raise on a
# perfectly normal pipeline order, not catch a real bug.
if hasattr(self, "veg_map"):
root_map_path = os.path.join(
self.workdir, self.project_name, self.input_dirname, "root_map"
)
if os.path.exists(root_map_path):
veg_on_disk, _ = in_CT.read_root_map(root_map_path)
veg_on_disk = np.asarray(veg_on_disk)
if hasattr(self, "DEM"):
valid_mask = self.DEM != -9999
if valid_mask.shape == veg_on_disk.shape:
veg_on_disk = veg_on_disk[valid_mask]
actual_nveg = int(np.max(veg_on_disk)) if veg_on_disk.size else 0
if actual_nveg > 0 and actual_nveg != int(self.cathyH["MAXVEG"]):
raise ValueError(
f"MAXVEG mismatch: CATHY.H is about to be written with "
f"MAXVEG={self.cathyH['MAXVEG']}, but root_map on disk "
f"({root_map_path}) requires NVEG={actual_nveg} (max "
f"vegetation class id inside the valid DEM domain). "
f"Writing CATHY.H now would desync it from the SOIL "
f"table and crash the Fortran run at the IVGHU read. "
f"Call update_veg_map(..., maxveg={actual_nveg}) (or "
f"pass MAXVEG={actual_nveg} explicitly here) before "
f"regenerating CATHY.H."
)
with open(
os.path.join(self.workdir,
self.project_name,
"src",
"CATHY.H"
), "w+") as CATHYH_file:
for i, l in enumerate(Lines0_109):
if i < 109:
CATHYH_file.write(l)
CATHYH_file.write(
" PARAMETER (ROWMAX={},COLMAX={},DEMRES={})\n".format(
self.cathyH["ROWMAX"], self.cathyH["COLMAX"], self.cathyH["DEMRES"]
)
)
CATHYH_file.write(" PARAMETER (MAXCEL=ROWMAX*COLMAX,MAXRES=1)\n")
# CATHYH_file.write(
# " PARAMETER (NODMAX=(ROWMAX/DEMRES+1)*(COLMAX/DEMRES+1))\n"
# )
CATHYH_file.write(
" PARAMETER (NODMAX={})\n".format(self.cathyH["NODMAX"])
)
# CATHYH_file.write(
# " PARAMETER (NTRMAX=2*MAXCEL/(DEMRES*DEMRES))\n".format()
# )
CATHYH_file.write(
" PARAMETER (NTRMAX={})\n".format(self.cathyH["NTRMAX"])
)
CATHYH_file.write(
" PARAMETER (NP2MAX=1,MAXSTR={})\n".format(self.cathyH["MAXSTR"])
)
CATHYH_file.write(
# " PARAMETER (NFACEMAX={})\n".format(self.cathyH['NFACEMAX']))
" PARAMETER (NFACEMAX={})\n".format(self.cathyH['NFACEMAX']))
CATHYH_file.write(
" PARAMETER (NMAX=NODMAX*(MAXSTR + 1),NTEMAX=3*NTRMAX*MAXSTR)\n".format()
)
CATHYH_file.write(
" PARAMETER (NPMAX={},NPMAX_TRA=1,NQMAX={},NSFMAX={})\n".format(
self.cathyH["NPMAX"], self.cathyH["NQMAX"], self.cathyH["NSFMAX"]
)
)
CATHYH_file.write(
" PARAMETER (NNSFMX={},MAXDIR=NODMAX+NPMAX+NSFMAX*NNSFMX)\n".format(
self.cathyH["NNSFMX"]
)
)
CATHYH_file.write(
" PARAMETER (MAXNUDN={},MAXNUDT={},MAXNUDC={})\n".format(
self.cathyH["MAXNUDN"],
self.cathyH["MAXNUDT"],
self.cathyH["MAXNUDC"],
)
)
CATHYH_file.write(
" PARAMETER (MAXZON={},MAXTRM={},MAXIT={},MAXVEG={})\n".format(
self.cathyH["MAXZON"],
self.cathyH["MAXTRM"],
self.cathyH["MAXIT"],
self.cathyH["MAXVEG"],
)
)
CATHYH_file.write(
" PARAMETER (NRMAX={},MAXPRT={},MAXVP={})\n".format(
self.cathyH["NRMAX"], self.cathyH["MAXPRT"], self.cathyH["MAXVP"]
)
)
CATHYH_file.write(
" PARAMETER (N1MAX={},NTPMAX=N1MAX*NMAX)\n".format(
self.cathyH["N1MAX"]
)
)
CATHYH_file.write(
" PARAMETER (MAXBOT={},INTBOT={},MAXQOUT={})\n".format(
self.cathyH["MAXBOT"], self.cathyH["INTBOT"], self.cathyH["MAXQOUT"]
)
)
CATHYH_file.write(
"cxcx PARAMETER (NIAUXMAX=NFACEMAX + MAXTRM + 1)\n".format()
)
CATHYH_file.write(
"cxcx PARAMETER (NRAUXMAX=5*NFACEMAX + MAXTRM,NQMAX_TRA=NODMAX)\n".format()
)
CATHYH_file.write(" PARAMETER (NIAUXMAX=NMAX + MAXTRM + 1)\n".format())
CATHYH_file.write(
" PARAMETER (NRAUXMAX=5*NMAX + MAXTRM,NQMAX_TRA=NODMAX)\n".format()
)
CATHYH_file.write(
" PARAMETER (MAXVTKPRT={})\n".format(self.cathyH["MAXVTKPRT"])
)
# CATHYH_file.write(" PARAMETER (MAXVTKPRT=9)\n".format())
CATHYH_file.write(" PARAMETER (MAXFCONTONODE=100,MAXLKP=3)\n")
CATHYH_file.close()
pass
def update_hapin(self, Lines, hapin, tmp_lnb, tmp_param_value):
L = []
tmp_param_value_new = []
count = 0
# Loop over hap.in file lines
# --------------------------------------------------------------------
for i, line in enumerate(Lines):
xnew = line.strip().split("=")
if len(xnew) > 1: # skip lines with no =
# loop over hapin dictionnary
# -----------------------------------
for key, value in self.hapin.items():
if count < len(hapin):
if key == hapin[tmp_lnb[count]]:
if value != tmp_param_value[count]:
# print(key)
# print(value)
if isinstance(value, list):
value_str = "/t".join(value)
xnew[1] = value_str
else:
xnew[1] = value
line = xnew[0] + "= " + str(xnew[1]) + "\n"
tmp_param_value_new.append(value)
count += 1 # count line nb
L.append(line)
# # write the new hap.in file
hap_file = open(
os.path.join(self.workdir, self.project_name, "prepro/hap.in"), "w+"
)
hap_file.writelines(L)
hap_file.close()
def update_prepo_inputs(self, DEM=None, verbose=False, show=False, **kwargs):
"""
Update default prepro inputs i.e. hap.in and dtm_13.val files based on kwargs
Parameters
----------
DEM : TYPE, optional
DESCRIPTION. The default is None.
verbose : TYPE, optional
DESCRIPTION. The default is False.
show : TYPE, optional
DESCRIPTION. The default is False.
**kwargs : TYPE
DESCRIPTION.
"""
# print("update hap.in")
self.console.print(":arrows_counterclockwise: [b]Update hap.in file[/b]")
structural_parameter = [
"delta_x",
"delta_y",
"N",
"M",
"N_celle",
"xllcorner",
"yllcorner",
]
terrain_parameter = [
"pt",
"imethod",
"lambda",
"CC_threshold",
"ndcf",
"nchc",
"A_threshold",
"ASk_threshold",
"kas",
"DN_threshold",
"local_slope_t",
"p_outflow_vo",
"bcc",
"cqm",
"cqg",
]
# 'RIVULET NETWORK PARAMETERS (HYDRAULIC GEOMETRY OF THE SINGLE RIVULET)'
rivulet_parameter = [
"dr",
"As_rf",
"(Qsf_rf,w_rf)",
"(Wsf_rf,b1_rf,b2_rf)",
"(kSsf_rf,y1_rf,y2_rf)",
"Qsi_rf",
]
channel_parameter = [
"As_cf",
"(Qsf_cf,w_cf)",
"(Wsf_cf,b1_cf,b2_cf)",
"(kSsf_cf,y1_cf,y2_cf)",
"Qsi_cf",
]
# idhapin = np.ones(22) + [1,1,2,3,3,1] #+ [1,2,3,3,1]
hapin = (
structural_parameter
+ terrain_parameter
+ rivulet_parameter
+ channel_parameter
)
hap_file = open(
os.path.join(self.workdir, self.project_name, "prepro/hap.in"), "r"
)
Lines = hap_file.readlines()
hap_file.close()
# read current values in hap.in
tmp_param_value = []
tmp_lnb = [] # save line nb where values are existing
count = 0
for line in Lines:
x = line.strip().split("=")
if len(x) > 1: # skip lines with no =
xs = " ".join(x[1].split())
xsr = xs.replace(" ", ",")
l = xsr.split(",")
if isinstance(l, list):
l = "/t".join(l)
tmp_param_value.append(l)
tmp_lnb.append(count)
count += 1
self.hapin = {}
for i in range(len(tmp_param_value)):
if re.search("/t", tmp_param_value[i]) is None:
if tmp_param_value[i].isdigit():
tmp_param_value[i] = int(tmp_param_value[i])
else:
tmp_param_value[i] = float(tmp_param_value[i])
else:
list_tmp = tmp_param_value[i].split("/t")
tmp_param_value[i] = [float(x) for x in list_tmp]
self.hapin[hapin[i]] = tmp_param_value[i]
# if self.hapin["dr"] != self.hapin["delta_x"]:
# if self.hapin["dr"] % self.hapin["delta_x"] != 0:
# print("adapt rivulet param to dem resolution")
# self.hapin["dr"] = self.hapin["delta_x"]
# sys.exit()
for key, value in kwargs.items():
self.hapin[key] = value
self.update_hapin(Lines, hapin, tmp_lnb, tmp_param_value)
self.update_dem_parameters(**kwargs)
# dtm_13.val
# If we start with a DEM file ("dtm_13.val") for an already delineated
# catchment (i.e., a "catchment" DEM file instead of a "full" DEM file), then
# only the last step in the above procedure (i.e., run "cppp" just once) is
# needed (make sure that "Boundary channel construction" is set to 0 in
# "hap.in").
if DEM is not None:
self.DEM = DEM
self.hapin["N"] = np.shape(DEM)[1]
self.hapin["M"] = np.shape(DEM)[0]
self.update_hapin(Lines,
hapin,
tmp_lnb,
tmp_param_value
)
# check presence of the outlet
if len(np.unique(DEM)) == 1:
print("Error: outlet not defined")
DEM[0, 0] = 0
self.console.print(":arrows_counterclockwise: [b]Update dtm_13 file[/b]")
self.update_dem(DEM)
self.update_zone()
with open(
os.path.join(self.workdir, self.project_name, "prepro/dtm_13.val"), "w+"
) as f:
# use exponential notation
np.savetxt(f, DEM, fmt="%1.4e")
# np.shape(DEM)
self.update_cathyH(
ROWMAX=self.hapin["M"],
COLMAX=self.hapin["N"],
MAXSTR=self.dem_parameters['nstr'],
DEMRES=min([self.hapin["delta_x"], self.hapin["delta_y"]])
)
self.update_dem_parameters(
**kwargs)
if show == True:
plt_CT.show_dem(DEM, self.hapin)
pass
def update_dem_parameters(self, **kwargs):
"""
Update DEM parameters - if no args reset to default values
ivert
=0 each layer will be parallel to the surface, including the base of the 3‐d grid.
ZRATIO is applied to each vertical cross section.
=1 base of the 3‐d grid will be flat.
ZRATIO is applied to each vertical cross section.
=2 base of the 3‐d grid will be flat, as will the NSTR‐1 horizontal cross sections above it.
ZRATIO is applied only to the vertical cross section having the lowest elevation.
=3 for each cell of the dem a single depth value is read in file input IIN60 (basement).
ZRATIO is applied to each vertical cross section.
=4 the first NSTR‐1 layers from the surface will be parallel to the surface and the base of the 3‐d grid will be flat.
ZRATIO is applied only to the vertical cross section having the lowest elevation.
isp
=0 for flat surface layer (only one Z value is read in, and is replicated to all surface nodes);
otherwise surface layer is not flat (Z values read in for each surface node);
(for ISP=0, IVERT=0, 1, and 2 yield the same 3‐d mesh, given the same values of BASE and ZRATIO).
base
Value which defines the thickness or base of the 3‐d mesh.
For IVERT=0, BASE is subtracted from each surface elevation value,
so that each vertical cross section will be of thickness BASE,
and the base of the 3‐d mesh will be parallel to the surface.
For IVERT=1 or 2, BASE is subtracted from the lowest surface elevation value,
say ZMIN, so that each vertical cross section will be of thickness (Z ‐ ZMIN) + BASE,
where Z is the surface elevation for that cross section.
The base of the 3‐d mesh will thus be flat.
"""
self.console.print(
":arrows_counterclockwise: [b]Update dem_parameters file [/b]"
)
if hasattr(self, "dem_parameters") == False:
dem_parameters_dict = in_CT.read_dem_parameters(
os.path.join(
self.workdir,
self.project_name,
"input",
"dem_parameters",
)
)
self.dem_parameters = dem_parameters_dict
if 'zratio(i).i=1.nstr' in self.dem_parameters:
self.dem_parameters['zratio(i),i=1,nstr'] = self.dem_parameters.pop('zratio(i).i=1.nstr')
# create dictionnary from kwargs
for keykwargs, value in kwargs.items():
if keykwargs == "zratio":
key = "zratio(i),i=1,nstr"
value = list(value)
if (abs(float(sum(value))-float(1)))>1e-3:
self.console.rule(
":warning: warning messages above :warning:", style="yellow"
)
self.console.print(
r"The sum of all the layers is not equal to 1 but to {}".format(
np.round(sum(value), 2)
),
style="yellow",
)
self.console.rule("", style="yellow")
else:
key = keykwargs
try:
self.dem_parameters[key]
self.dem_parameters[key] = value
except:
pass
dem_parameters_tmp = {}
for key, value in self.dem_parameters.items():
if isinstance(value, list):
strlst = "\t".join(str(e) for e in value)
self.dem_parameters[key] = strlst
# if (key +'_list') in self.dem_parameters.keys() is False:
# dem_parameters_tmp[key +'_list'] = value
if 'zratio' in kwargs:
self.dem_parameters['nstr'] = len(kwargs['zratio'])
# self.dem_parameters = self.dem_parameters | dem_parameters_tmp
# write file
header_fmt = [1, 1, 1, 1, 3, 3, 1]
with open(
os.path.join(
self.workdir, self.project_name, self.input_dirname, "dem_parameters"
),
"w+",
) as dem_parametersfile:
# Write values first
# -----------------------
counth = 0
for h in header_fmt:
if h == 3:
dem_parametersfile.write(
str(list(self.dem_parameters.values())[counth])
+ "\t"
+ str(list(self.dem_parameters.values())[counth + 1])
+ "\t"
+ str(list(self.dem_parameters.values())[counth + 2])
+ "\t"
+ "\n"
)
counth += 3
if h == 1:
# print(str(list(self.dem_parameters.values())[counth]))
dem_parametersfile.write(
str(list(self.dem_parameters.values())[counth]) + "\t" + "\n"
)
counth += 1
# Write keys
# -----------------------
counth = 0
for h in header_fmt:
if h == 3:
dem_parametersfile.write(
str(list(self.dem_parameters.keys())[counth])
+ "\t"
+ str(list(self.dem_parameters.keys())[counth + 1])
+ "\t"
+ str(list(self.dem_parameters.keys())[counth + 2])
+ "\t"
+ "\n"
)
counth += 3
if h == 1:
dem_parametersfile.write(
str(list(self.dem_parameters.keys())[counth]) + "\t" + "\n"
)
counth += 1
dem_parametersfile.close()
pass
def update_lakes_map(self, lakes_map=[]):
"""
Update lakes_map file
Parameters
----------
dem : np.array([]), optional
2d numpy array describing the zone for x and y direction.
The array dim must be equal to DEM dim.
The default is []."""
with open(
os.path.join(self.workdir, self.project_name, "prepro/lakes_map"), "w+"
) as lakesfile:
lakesfile.write("north: 0" + "\n")
lakesfile.write("south: " + str(self.hapin["yllcorner"]) + "\n")
lakesfile.write("east: 0" + "\n")
lakesfile.write("west: " + str(self.hapin["xllcorner"]) + "\n")
lakesfile.write("rows: " + str(self.hapin["N"]) + "\n")
lakesfile.write("cols: " + str(self.hapin["M"]) + "\n")
np.savetxt(lakesfile, lakes_map, fmt="%i")
# np.shape(DEM)
lakesfile.close()
pass
def update_dem(self, DEM=[]):
"""
Update zone file
Parameters
----------
dem : np.array([]), optional
2d numpy array describing the zone for x and y direction.
The array dim must be equal to DEM dim.
The default is []."""
with open(
os.path.join(self.workdir, self.project_name, "prepro/dem"), "w+"
) as demfile:
demfile.write("north: 0" + "\n")
demfile.write("south: " + str(self.hapin["yllcorner"]) + "\n")
demfile.write("east: 0" + "\n")
demfile.write("west: " + str(self.hapin["xllcorner"]) + "\n")
demfile.write("rows: " + str(self.hapin["N"]) + "\n")
demfile.write("cols: " + str(self.hapin["M"]) + "\n")
np.savetxt(demfile, DEM, fmt="%1.4e")
# np.shape(DEM)
demfile.close()
with open(
os.path.join(self.workdir, self.project_name, "prepro/dtm_13.val"), "w+"
) as f:
np.savetxt(f, DEM, fmt="%1.4e")
self.DEM = DEM
pass
def update_zone(self, zone=[]):
"""
Update zone file
Parameters
----------
zone : np.array([]), optional
2d numpy array describing the zone for x and y direction.
The array dim must be equal to DEM dim.
The smallest numbering is 1 (not 0!)
If zone is empty, reset to homogeneous zone = 1
The default is [].
Notes
-----
.. note::
- Create object zone variable
- Update dem_parameters file (number of zones).
- Update parm file.
- Update CATHYH file.
- Update vtk mesh file.
"""
self.console.print(":arrows_counterclockwise: [b]update zone file [/b]")
with open(
os.path.join(self.workdir, self.project_name, "prepro/zone"), "w+"
) as zonefile:
zonefile.write("north: 0" + "\n")
zonefile.write("south: " + str(self.hapin["yllcorner"]) + "\n")
zonefile.write("east: 0" + "\n")
zonefile.write("west: " + str(self.hapin["xllcorner"]) + "\n")
zonefile.write("rows: " + str(self.hapin["M"]) + "\n")
zonefile.write("cols: " + str(self.hapin["N"]) + "\n")
if len(zone) == 0:
zone = np.c_[np.ones([self.hapin["M"], self.hapin["N"]])]
np.savetxt(zonefile, zone, fmt="%i")
else:
# if np.shape(zone)== :
np.savetxt(zonefile, zone, fmt="%i")
zonefile.close()
self.zone = zone
# update number of zone in the dem parameter file
self.update_dem_parameters(nzone=len(np.unique(zone)))
# self.dem_parameters
self.update_parm()
self.update_cathyH(MAXZON=len(np.unique(zone)))
# %% INPUT FILES
def update_parm(self, **kwargs):
"""
Update parm file
.. note:
- create object zone parm
- Updated CATHYH file (NPRT and NUMVP).
"""
self.console.print(":arrows_counterclockwise: [b]Update parm file [/b]")
warnings_parm = []
if len(self.parm) == 0:
#%%
dict_parm = in_CT.read_parm(os.path.join(self.workdir,
self.project_name,
"input",
"parm"
)
)
#%%
self.parm = dict_parm
if '(TIMPRT(I)I=1NPRT)' in self.parm:
self.parm['(TIMPRT(I),I=1,NPRT)'] = self.parm.pop('(TIMPRT(I)I=1NPRT)')
# create self.parm dictionnary from kwargs
# --------------------------------------------------------------------
for kk, value in kwargs.items():
if kk in self.parm.keys():
if self.verbose == True:
print(f"key: {kk} | value: {value}")
# times of interest TIMPRTi
# ----------------------------------------------------------------
if kk == "TIMPRTi":
key = "(TIMPRT(I),I=1,NPRT)"
TIMPRTinew = [int(lv) for lv in list(value)]
self.parm[key] = TIMPRTinew
# points of interest NODVP
# ----------------------------------------------------------------
elif kk == "NODVP":
key = "(NODVP(I),I=1,NUMVP)"
self.parm[key] = value
# check if consistency between node of interest and
# number of nodes of interest
# ------------------------------------------------------------
if len(value) != self.parm["NUMVP"]:
warnings_parm.append(
"Adjusting NUMVP with respect to NODVP requested" + "\n"
)
self.parm["NUMVP"] = len(value)
# points of interest NR
# ----------------------------------------------------------------
elif kk == "ID_NR":
key = "CONTR(I),I=1,NR"
self.parm[key] = value
# check if consistency between node of interest and
# number of nodes of interest
# ------------------------------------------------------------
if len(value) != self.parm["NR"]:
warnings_parm.append(
"Adjusting NR with respect to CONTR requested" + "\n"
)
self.parm["NR"] = len(value)
# points of interest NR
# ----------------------------------------------------------------
elif kk == "ID_QOUT":
key = "(ID_QOUT(I),I=1,NUM_QOUT)"
self.parm[key] = value
# check if consistency between node of interest and
# number of nodes of interest
# ------------------------------------------------------------
if len(value) != self.parm["NUM_QOUT"]:
warnings_parm.append(
"Adjusting NUM_QOUT with respect to ID_QOUT requested" + "\n"
)
self.parm["NUM_QOUT"] = len(value)
# other type of kwargs
# ----------------------------------------------------------------
else:
if kk in self.parm.keys():
self.parm[kk] = value
# check if consistency between times of interest and TMAX
# ------------------------------------------------------------
if "(TIMPRT(I),I=1,NPRT)" in kk:
if self.parm["TMAX"] < max(self.parm["(TIMPRT(I),I=1,NPRT)"]):
self.parm["TMAX"] = int(max(self.parm["(TIMPRT(I),I=1,NPRT)"]))
# check if consistency between times of interest and
# number of times of interest
# ------------------------------------------------------------
if len(self.parm["(TIMPRT(I),I=1,NPRT)"]) != self.parm["NPRT"]:
warnings_parm.append(
"Adjusting NPRT with respect to time of interests requested" + "\n"
)
self.parm["NPRT"] = len(self.parm["(TIMPRT(I),I=1,NPRT)"])
# check if consistency between DELTAT, DTMIN, and DTMAX
# ------------------------------------------------------------
if min(np.diff(self.parm["(TIMPRT(I),I=1,NPRT)"])) <= 0:
raise ValueError("Vtk time steps should be monotonically increasing")
if self.parm["DTMIN"] > min(np.diff(self.parm["(TIMPRT(I),I=1,NPRT)"])
# ) or self.parm["DELTAT"] > 1e-1 * min(
# np.diff(self.parm["(TIMPRT(I),I=1,NPRT)"])
):
warnings_parm.append(
"Adjusting DTMIN with respect to time of interests requested" + "\n"
)
self.parm["DTMIN"] = min(np.diff(self.parm["(TIMPRT(I),I=1,NPRT)"])) / 1e2
if self.parm["DELTAT"] < self.parm["DTMIN"]:
warnings_parm.append(
"adjusting 2*DTMIN == DELTAT: " + str(self.parm["DTMIN"]) + "\n"
)
self.parm["DTMIN"] = self.parm["DELTAT"] / 2
if self.parm["DELTAT"] >= self.parm["DTMAX"]:
warnings_parm.append("Adjusting DTMAX == 2*DELTAT" + "\n")
self.parm["DTMAX"] = 1e2 * self.parm["DELTAT"]
if self.parm["TMAX"] < max(self.parm["(TIMPRT(I),I=1,NPRT)"]):
warnings_parm.append(
"Adjusting TMAX with respect to time of interests requested" + "\n"
)
self.parm["TMAX"] = int(max(self.parm["(TIMPRT(I),I=1,NPRT)"]))
if len(warnings_parm) > 0:
self.console.rule(
":warning: warning messages above :warning:", style="yellow"
)
for ww in warnings_parm:
self.console.print(warnings_parm, style="yellow")
self.console.rule("", style="yellow")
# transform array args to list
# --------------------------------------------------------------------
for kk, value in self.parm.items():
if "numpy.array" in str(type(value)):
value = list(value)
self.parm[kk] = value
# update CATHYH
# --------------------------------------------------------------------
self.update_cathyH(
MAXPRT=self.parm["NPRT"],
MAXVP=self.parm["NUMVP"],
indent=" :arrow_right_hook:",
**kwargs,
)
# write parm file
# --------------------------------------------------------------------
if "filename" in kwargs:
file2write = kwargs["filename"]
else:
file2write = os.path.join(
self.workdir, self.project_name, self.input_dirname, "parm"
)
# backup file during DA scheme cycle
# --------------------------------------------------------------------
backup = False
if "backup" in kwargs:
backup = kwargs["backup"]
if backup == True:
if self.count_DA_cycle is not None:
dst_dir = file2write + str(self.count_DA_cycle)
shutil.copy(file2write, dst_dir)
self._write_parm_file(file2write)
pass
def _write_parm_file(self, file2write):
"""
Overwrite existing parm file
Returns
-------
New overwritten file.
"""
if int(self.parm['NR'])==0:
header_fmt_parm = [3, 3, 2, 4, 4, 3, 3, 2, 4, 3, 3, 4, 4, 4, 2, 1, 1, 1]
try:
del self.parm["CONTR(I),I=1,NR"]
except:
pass
else:
header_fmt_parm = [3, 3, 2, 4, 4, 3, 3, 2, 4, 3, 3, 4, 4, 4, 2, 1, 1, 1, 1]
counth = 0
# with open(
# os.path.join(self.workdir, self.project_name, self.input_dirname, "parm"),
# "w+",
# ) as parmfile:
with open(file2write, "w+") as parmfile:
# Write line by line according to header format
# ----------------------------------------------------------------
for i, h in enumerate(header_fmt_parm):
left = right = []
# left = values
# ------------------------------------------------------------
left = str(list(self.parm.values())[counth : counth + h])
left = left.strip("[]").replace(",", "")
left = left.strip("[]").replace("[", "")
# right = keys
# ------------------------------------------------------------
right = str(list(self.parm.keys())[counth : counth + h])
right = right.strip("[]").replace("',", "")
right = right.replace("'", "")
# add left + right
# ------------------------------------------------------------
line = left + "\t" + right + "\n"
counth += h
parmfile.write(str(line))
parmfile.close()
pass
def update_ic(self, INDP=2, IPOND=0, WTPOSITION=0, verbose=False, **kwargs):
"""
The initial conditions file contains the pressure heads distribution for the study area (INDP)
For example, to simulate a uniform water table depth or 0.5 m or 1.0 m from the ground surface,
INDP=3 and WTHEIGHT=4.5 are selected
Parameters
----------
INDP : int, optional
Flag for pressure head initial conditions (all nodes). The default is 2.
- =0 for input of uniform initial conditions (one value read in)
- =1 for input of non-uniform IC's (one value read in for each node)
- =2 for calculation of fully saturated vertical hydrostatic equilibrium IC's
(calculated in subroutine ICVHE). In the case of IPOND>0, the fully saturated
hydrostatic IC is calculated (in subroutine ICVHEPOND) starting from the ponding head
values at the surface nodes, rather than surface pressure heads of 0.
- =3 for calculation of partially saturated vertical hydrostatic equilibrium IC's
(calculated in subroutine ICVHWT) with the water table height (relative to the base
of the 3‐d grid) given by parameter WTHEIGHT
- =4 for calculation of partially saturated vertical hydrostatic equilibrium IC's
(calculated in subroutine ICVDWT) with the water table depth (relative to the surface
of the 3‐d grid) given by parameter WTHEIGHT
IPOND : int, optional
Flag for ponding head initial conditions (surface nodes). The default is 0.
- =0 no input of ponding head initial conditions; otherwise (IPOND = 1 or 2) ponding
head initial conditions are read into PONDNOD, and, where PONDNOD > 0, these values
are used to update the surface node values in PtimeP read in according to the
previous INDP flag
- =1 uniform ponding head initial conditions (one value read in)
- =2 non-uniform ponding head initial conditions (one value read in for each node)
WTPOSITION : float, optional
For the case INDP=3, specifies the initial water table height relative to
the base of the 3‐d grid. The default is 0.
kwargs :
- pressure_head_ini (int): uniform value of pressure head to assign (only if INDP=0)
Returns
-------
new ic file written/overwritten.
"""
if "shellprint_update" in kwargs:
if kwargs["shellprint_update"] is False:
pass
else:
self.console.print(":arrows_counterclockwise: [b]Update ic[/b]")
# check value of WTPOSITION
# --------------------------------------------------------------------
# For the case INDP=3, specifies the initial water table
# height relative to the base of the 3‐d grid
# if WTPOSITION>0:
# print('WTPOSITION must be negative - water table height
# relative to the base of the 3‐d grid')
# sys.exit()
# set default parameters
# --------------------------------------------------------------------
self.ic = {
"INDP": INDP,
"IPOND": IPOND,
# For the case INDP=3, specifies the initial water table
# height relative to the base of the 3‐d grid
"WTPOSITION": WTPOSITION,
}
# write ic file
# --------------------------------------------------------------------
if "filename" in kwargs:
filename = kwargs["filename"]
else:
filename = os.path.join(
self.workdir, self.project_name, self.input_dirname, "ic"
)
backup = False
if "backup" in kwargs:
backup = kwargs["backup"]
if backup == True:
if self.count_DA_cycle is not None:
dst_dir = filename + str(self.count_DA_cycle - 1)
shutil.copy(filename, dst_dir)
with open(filename, "w+") as icfile:
if INDP == 0:
icfile.write(str(INDP) + "\t" + str(IPOND) + "\t" + "INDP \n")
if "pressure_head_ini" in kwargs:
np.savetxt(icfile, [kwargs["pressure_head_ini"]], fmt="%1.3e")
# try:
self.map_prop2mesh({"ic": kwargs["pressure_head_ini"]})
# except:
# print('Impossible to map pressure head to mesh')
else:
raise ValueError("Missing initial pressure value")
elif INDP == 1:
pressure_head_ini = kwargs["pressure_head_ini"]
icfile.write(
str(INDP)
+ "\t"
+ str(IPOND)
+ "\t"
+ "INDP"
+ "\t"
+ "IPOND"
+ "\n"
)
np.savetxt(icfile, pressure_head_ini, fmt="%1.3e")
self.map_prop2mesh({"ic": kwargs["pressure_head_ini"]})
elif INDP==2:
icfile.write(
str(INDP)
+ "\t"
+ str(IPOND)
+ "\t"
+ "INDP"
+ "\t"
+ "IPOND"
+ "\n"
)
elif INDP in [3, 4]:
icfile.write(
str(INDP)
+ "\t"
+ str(IPOND)
+ "\t"
+ "INDP"
+ "\t"
+ "IPOND"
+ "\n"
)
icfile.write(str(WTPOSITION) + "\t" + "WTPOSITION" + "\n")
icfile.close()
pass
def update_atmbc(
self,
HSPATM=1,
IETO=0,
time=None,
VALUE=[None, None],
netValue=[],
show=False,
verbose=False,
**kwargs,
):
"""
Atmospheric forcing term (atmbc - IIN6)
..note:
1 1 HSPATM,IETO
0.0000000e+00 time
5.5e-06 VALUE
12.000000e+03 time
0.00 VALUE
18.000000e+03 time
0.00 VALUE
The values are those of a 200-min rainfall event at a uniform
intensity of 3.3·10-4 m/min, followed by 100 min of drainage.
..note:
In case of simultaneous precipitation and evaporation, we impose at
the surface the net flux, i.e., precipitation minus evaporation.
Parameters
----------
HSPATM : int, optional
- =0 for spatially variable atmospheric boundary condition inputs;
blank or =9999 if unit IIN6 input is to be ignored; otherwise atmospheric BC's are
homogeneous in space.
IETO : int, optional
- =0 for linear interpolation of the atmospheric boundary condition inputs between different different ATMTIM
- otherwise the inputs are assigned as a piecewise constant function (ietograph).
The default is 0.
time : list
ATMBC Times in seconds. The default is None.
VALUE : list
List of array. The default is [Precipitation, EvapoTranspiration].
show : bool, optional
Plot atmbc. The default is False.
verbose : bool, optional
Display. The default is False.
**kwargs
Returns
-------
..note:
- Update parm file (NPRT).
- Update CATHYH file (MAXPRT).
"""
if "shellprint_update" in kwargs:
if kwargs["shellprint_update"] is False:
pass
else:
self.console.print(":arrows_counterclockwise: [b]Update atmbc[/b]")
# type(netValue)
if len(netValue) > 0:
v_atmbc = netValue
else:
# atmbc are homoegenous
# -----------------------------------------------------------------
if HSPATM == 1:
if len(VALUE) == 2: # take the difference between Precipitation and EvapoTranspiration
v_atmbc = VALUE[0] - abs(VALUE[1])
else: # Assume it is already the net
v_atmbc = VALUE
print('Assuming it is already the net')
else:
print('spatial v_atmbc must be provided as net between Precipitation and EvapoTranspiration')
# set default parameters
# --------------------------------------------------------------------
self.atmbc = {"HSPATM": HSPATM, "IETO": IETO, "time": time, "VALUE": v_atmbc}
# len(self.atmbc['time'])
# overwrite existing input atmbc file
# --------------------------------------------------------------------
if "filename" in kwargs:
filename = kwargs["filename"]
else:
filename = os.path.join(
self.workdir, self.project_name, self.input_dirname, "atmbc"
)
backup = True
if "backup" in kwargs:
backup = kwargs["backup"]
if backup == True:
if self.count_DA_cycle is not None:
dst_dir = filename + str(self.count_DA_cycle - 1)
shutil.copy(filename, dst_dir)
with open(filename, "w+") as atmbcfile:
atmbcfile.write(
str(HSPATM) + "\t" + str(IETO) + "\t" + "HSPATM" + "\t" + "IETO" + "\n"
)
# atmbc are homogeneous
# -----------------------------------------------------------------
if HSPATM == 1:
for t, v in zip(time, v_atmbc):
if verbose:
print(t, v)
# atmbcfile.write("{:.0f}".format(t) + "\t" + "time" + "\n")
atmbcfile.write("%i\t%s\n" % (t, "time"))
# atmbcfile.close()
if isinstance(v, float) | isinstance(v, int):
atmbcfile.write("{:.3e}".format(v) + "\t" + "VALUE" + "\n")
else:
atmbcfile.write("{:.3e}".format(float(v)) + "\t" + "VALUE" + "\n")
# atmbcfile.write(str(v) + "\t" + "VALUE" + "\n")
# atmbc are heteregeneous
# -----------------------------------------------------------------
else:
for t, v in zip(time, v_atmbc):
atmbcfile.write("{:.0f}".format(t) + "\t" + "time" + "\n")
np.savetxt(atmbcfile, v, fmt="%.3e")
atmbcfile.close()
self.update_parm(TIMPRTi=[time[0],time[-1]], NPRT=2, TMAX=int(max(time)))
# self.update_parm(TIMPRTi=time, NPRT=len(time), TMAX=max(time))
# don't need to update if sequential DA as the cathy.exe is already created
# ---------------------------------------------------------------------
if "omit_cathyH" not in kwargs:
self.update_cathyH(MAXPRT=len(time))
if show:
if HSPATM ==0:
print('impossible to plot for non homogeneous atmbc')
# # sys.exit()
else:
# x_units = "sec"
# for key, value in kwargs.items():
# if key == "x_units":
# x_units = value
if len(netValue) > 0:
plt_CT.show_atmbc(time, v_atmbc, **kwargs)
else:
plt_CT.show_atmbc(time, VALUE, **kwargs)
pass
def _record_bc_in_mesh_bounds_df(self, BCtypName, time, nodes, values, no_flow):
"""
Keep self.mesh_bound_cond_df (consumed by show_bc()/plot_mesh_bounds)
in sync with what update_nansfdirbc/update_nansfneubc/update_sfbc
actually wrote to the BC file on disk.
Previously update_nansfdirbc/update_nansfneubc/update_sfbc only
wrote the raw BC file and never touched mesh_bound_cond_df, so
show_bc() either crashed with AttributeError (df never created)
or, if the df had been built by hand elsewhere via
create_mesh_bounds_df(), crashed with KeyError (that method never
actually populates a "nansfdirbc"/"nansfneubc"/"sfbc" column -
only the never-called-here assign_mesh_bc_df() did, and only for
its no_flow branch).
This builds mesh_bound_cond_df on first use and stores, per
(node, time), the value actually imposed - NaN where no BC is
imposed on that node at that time. plot_mesh_bounds() treats NaN
as "not applied" and colours applied nodes by their real value.
"""
if no_flow or nodes is None:
# Nothing explicit imposed at any node/time for this BC type -
# leave the column as NaN ("not applied") rather than the old
# placeholder behaviour of flagging every mesh node (the
# noflow_bound mask in create_mesh_bounds_df is unconditionally
# True for all nodes, not just boundary ones, so reusing it
# here would mislabel the whole interior mesh as boundary BC).
# NOTE: this check now runs BEFORE the grid3d/mesh_bound_cond_df
# lookup below, since a pure no-flow call needs no grid info at
# all - previously it ran after, so a no_flow=True call (the
# common case, e.g. baseline runs) still tried to build the df
# and crashed if grid3d wasn't available yet.
if not hasattr(self, "mesh_bound_cond_df"):
return
if BCtypName not in self.mesh_bound_cond_df.columns:
self.mesh_bound_cond_df[BCtypName] = np.nan
return
if not hasattr(self, "mesh_bound_cond_df"):
# self.grid3d is set to {} in __init__, so it always exists as
# an attribute - hasattr() alone can't tell "not populated yet"
# from "populated". Check its content instead (as done
# elsewhere in this file, e.g. read_inputs()'s
# "if len(self.grid3d) == 0" check).
if len(self.grid3d) == 0 or "mesh3d_nodes" not in self.grid3d:
try:
self.grid3d = self.read_outputs("grid3d")
except Exception as e:
# update_nansfdirbc/neubc/sfbc are normally called
# BEFORE run_processor(), while output/grid3d doesn't
# exist on disk yet. This bookkeeping only feeds
# show_bc()/plot_mesh_bounds() after a run - it is not
# needed to write the BC file itself, so skip it
# quietly instead of aborting the whole run. It will
# be rebuilt lazily (successfully) once grid3d is
# available, e.g. the first time show_bc() is called.
print(
f"[_record_bc_in_mesh_bounds_df] grid3d not available yet "
f"({e!r}); skipping mesh_bound_cond_df bookkeeping for "
f"'{BCtypName}' - this does not affect the BC file written "
f"to disk."
)
return
if not self.grid3d or "mesh3d_nodes" not in self.grid3d:
print(
f"[_record_bc_in_mesh_bounds_df] grid3d has no 'mesh3d_nodes' "
f"yet; skipping mesh_bound_cond_df bookkeeping for "
f"'{BCtypName}' - this does not affect the BC file written "
f"to disk."
)
return
self.create_mesh_bounds_df(
BCtypName, self.grid3d["mesh3d_nodes"], times=list(time)
)
if BCtypName not in self.mesh_bound_cond_df.columns:
self.mesh_bound_cond_df[BCtypName] = np.nan
for tt in time:
tt_nodes = nodes[tt] if isinstance(nodes, dict) else nodes
tt_nodes = np.atleast_1d(np.asarray(tt_nodes))
if isinstance(values, dict):
tt_values = np.atleast_1d(np.asarray(values[tt]))
elif np.isscalar(values):
tt_values = np.full(tt_nodes.shape, values, dtype=float)
else:
tt_values = np.atleast_1d(np.asarray(values))
value_by_node = dict(zip(tt_nodes.tolist(), tt_values.tolist()))
row_mask = (
(self.mesh_bound_cond_df["time"] == tt)
& self.mesh_bound_cond_df["id_node"].isin(tt_nodes)
)
self.mesh_bound_cond_df.loc[row_mask, BCtypName] = (
self.mesh_bound_cond_df.loc[row_mask, "id_node"].map(value_by_node)
)
def _resolve_bc_time(self, time, caller_name):
"""
Shared time-argument default for update_nansfdirbc,
update_nansfneubc, and update_sfbc: if `time` isn't given,
fall back to the times already established via update_atmbc
(self.atmbc["time"]), so these boundary-condition files line
up with the atmospheric BC time stepping without repeating the
time list at every call site.
Passing `time=[]` explicitly still means "write nothing" (an
empty file of BC blocks) and is left untouched - only `None`
triggers the atmbc fallback.
"""
if time is not None:
return time
if not hasattr(self, "atmbc") or self.atmbc.get("time") is None:
raise ValueError(
f"{caller_name}: no `time` provided and no atmbc time "
"steps are set yet to default to - call update_atmbc "
"first, or pass `time` explicitly."
)
return self.atmbc["time"]
def update_nansfdirbc(
self,
time=None,
nodes=None,
pressure_head=None,
no_flow=False,
**kwargs,
):
"""
Write the nansfdirbc file (non-atmospheric, non-seepage-face
Dirichlet boundary conditions).
This mirrors exactly what the Fortran side (BCONE/BCNXT ->
RDNDBC/READBC) expects to read, one block per time step:
TIME
NODIN2 NODINF
[node id list] <- only written if NODIN2*NODINF-derived NBC>0
[value list]
We always use the NODIN2=0 ("fixed node list") branch of
RDNDBC/READBC, i.e. NODIN2 is always written as 0 and NODINF is
the number of explicit nodes for that time step. This is the
most general form: it does not rely on CATHY auto-replicating
surface nodes down the mesh (NODIN2>0 branch), so any arbitrary
3-D node can be prescribed directly.
NOTE: Fortran node numbering is 1-based. The mesh dataframes
built elsewhere in this wrapper (id_node) are 0-based, so node
indices passed in here are assumed 0-based and are converted
(+1) when written to file.
Parameters
----------
time : list, optional
Simulation times (s) at which a BC block is written. One
block is written per entry, in order. If not given (None),
defaults to the times already set via update_atmbc
(self.atmbc["time"]), so the Dirichlet/Neumann/seepage-face
BC files line up with the atmospheric BC time stepping by
default without having to repeat the time list at every
call site.
nodes : array-like or dict, optional
- array-like of 0-based node indices: same nodes used for
every time step.
- dict {t: array-like}: time-varying node set, keyed by the
matching entries of `time`.
Required unless no_flow=True.
pressure_head : float, array-like or dict, optional
- scalar: same pressure head applied to every listed node,
every time step.
- array-like matching `nodes` (when `nodes` is array-like):
one value per node, same for every time step.
- dict {t: array-like matching nodes[t]}: time-varying
values.
Required unless no_flow=True.
no_flow : bool, optional
If True, write NODIN2=NODINF=0 for every time step (no
Dirichlet nodes imposed at all -> CATHY falls back to its
default no-flow/natural boundary). `nodes`/`pressure_head`
are ignored in this case.
"""
time = self._resolve_bc_time(time, "update_nansfdirbc")
dirbcfile = open(os.path.join(self.workdir, self.project_name,
"input", "nansfdirbc"), "w+")
max_nodinf = 0
for tt in time:
dirbcfile.write("{:.0f}".format(tt) + "\t" + "time" + "\n")
if no_flow:
dirbcfile.write("0" + "\t" + "0" + "\n")
continue
# --- resolve nodes/values for this time step ---
tt_nodes = nodes[tt] if isinstance(nodes, dict) else nodes
tt_nodes = np.atleast_1d(np.asarray(tt_nodes))
if isinstance(pressure_head, dict):
tt_values = np.atleast_1d(np.asarray(pressure_head[tt]))
elif np.isscalar(pressure_head):
tt_values = np.full(tt_nodes.shape, pressure_head, dtype=float)
else:
tt_values = np.atleast_1d(np.asarray(pressure_head))
if len(tt_values) != len(tt_nodes):
raise ValueError(
"nansfdirbc: pressure_head and nodes must have the "
f"same length at time {tt} "
f"({len(tt_values)} != {len(tt_nodes)})"
)
nodin2 = 0
nodinf = len(tt_nodes)
max_nodinf = max(max_nodinf, nodinf)
dirbcfile.write(f"{nodin2}\t{nodinf}\n")
if nodinf > 0:
# Fortran node numbering is 1-based
fortran_nodes = tt_nodes + 1
dirbcfile.write(
" ".join(str(int(n)) for n in fortran_nodes) + "\n"
)
dirbcfile.write(
" ".join(f"{v:.6e}" for v in tt_values) + "\n"
)
dirbcfile.close()
# CATHY.H's NPMAX bounds NDIRC (the fixed-node count we use via
# the NODIN2=0 branch) - keep it in sync so the compiled solver
# doesn't overflow the BCNOD/BCINP arrays.
if not no_flow and max_nodinf > 0:
self.update_cathyH(NPMAX=max(max_nodinf, 1), verbose=False)
self._record_bc_in_mesh_bounds_df(
"nansfdirbc", time, nodes, pressure_head, no_flow
)
pass
def update_nansfneubc(
self,
time=None,
nodes=None,
flux=None,
no_flow=False,
**kwargs,
):
"""
Write the nansfneubc file (non-atmospheric, non-seepage-face
Neumann boundary conditions). Same file grammar and same
NODIN2=0 explicit-node-list convention as update_nansfdirbc
(RDNDBC's ZERO/NQ header maps directly onto NODIN2/NODINF).
Parameters
----------
time : list, optional
Simulation times (s), one BC block written per entry. If
not given (None), defaults to self.atmbc["time"] (see
update_nansfdirbc).
nodes : array-like or dict, optional
0-based node indices with an imposed flux; same layout
rules as in update_nansfdirbc. Required unless no_flow=True.
flux : float, array-like or dict, optional
Imposed flux value(s); same layout rules as pressure_head
in update_nansfdirbc. Required unless no_flow=True.
no_flow : bool, optional
If True, write ZERO=NQ=0 for every time step (no Neumann
nodes imposed).
"""
time = self._resolve_bc_time(time, "update_nansfneubc")
neubcfile = open(os.path.join(self.workdir, self.project_name,
"input", "nansfneubc"), "w+")
max_nq = 0
for tt in time:
neubcfile.write("{:.0f}".format(tt) + "\t" + "time" + "\n")
if no_flow:
neubcfile.write("0" + "\t" + "0" + "\n")
continue
tt_nodes = nodes[tt] if isinstance(nodes, dict) else nodes
tt_nodes = np.atleast_1d(np.asarray(tt_nodes))
if isinstance(flux, dict):
tt_values = np.atleast_1d(np.asarray(flux[tt]))
elif np.isscalar(flux):
tt_values = np.full(tt_nodes.shape, flux, dtype=float)
else:
tt_values = np.atleast_1d(np.asarray(flux))
if len(tt_values) != len(tt_nodes):
raise ValueError(
"nansfneubc: flux and nodes must have the same "
f"length at time {tt} "
f"({len(tt_values)} != {len(tt_nodes)})"
)
zero = 0
nq = len(tt_nodes)
max_nq = max(max_nq, nq)
neubcfile.write(f"{zero}\t{nq}\n")
if nq > 0:
fortran_nodes = tt_nodes + 1
neubcfile.write(
" ".join(str(int(n)) for n in fortran_nodes) + "\n"
)
neubcfile.write(
" ".join(f"{v:.6e}" for v in tt_values) + "\n"
)
neubcfile.close()
# CATHY.H's NQMAX bounds NQ - keep it in sync.
if not no_flow and max_nq > 0:
self.update_cathyH(NQMAX=max(max_nq, 1), verbose=False)
self._record_bc_in_mesh_bounds_df(
"nansfneubc", time, nodes, flux, no_flow
)
pass
def update_sfbc(
self,
time=None,
no_flow=True,
**kwargs,
):
"""
Write the sfbc file (seepage face boundary conditions).
Only the no-flow / no-seepage-face case (NSF=0 at every time
step) is implemented here, matching the previous behaviour.
Writing non-trivial seepage faces requires knowing the exact
record layout expected by the Fortran seepage-face reader
(e.g. SFONE/SFNODE, analogous to RDNDBC/READBC for
nansfdirbc/nansfneubc) - that source isn't available yet, so
implementing it here would risk guessing a wrong file layout
and silently corrupting a run. Share that source and this can
be extended the same way update_nansfdirbc/update_nansfneubc
were.
Parameters
----------
time : list, optional
Simulation times (s), one block written per entry. If not
given (None), defaults to self.atmbc["time"] (see
update_nansfdirbc).
"""
time = self._resolve_bc_time(time, "update_sfbc")
if not no_flow:
raise NotImplementedError(
"update_sfbc: non-trivial seepage faces are not "
"implemented - the Fortran seepage-face reader "
"(SFONE/SFNODE or equivalent) hasn't been provided, "
"so the exact file layout for NSF>0 blocks is unknown. "
"Only no_flow=True (NSF=0) is currently supported."
)
sfbcfile = open(os.path.join(self.workdir, self.project_name,
"input", "sfbc"), "w+")
for tt in time:
sfbcfile.write("{:.0f}".format(tt) + "\n")
sfbcfile.write("0" + "\n")
sfbcfile.close()
self._record_bc_in_mesh_bounds_df(
"sfbc", time, nodes=None, values=None, no_flow=True
)
pass
def update_soil(
self,
IVGHU=[],
FP=[],
FP_map=[],
SPP=[],
SPP_map=[],
zone3d=[],
show=False,
**kwargs,
):
"""
Soil parameters (soil - IIN4). The porous media properties.
..note::
The first thing that must be decides is the type of relationship to describe the hydraulic
characteristics of the unsaturated soil (i.e. retention curves). This can be done through
the choice of the parameter **IVGHU** amongst the several options.
Parameters
----------
IVGHU : int, optional
= -1 table look up for moisture curves
= 0 for van Genuchten moisture curves
= 1 for extended van Genuchten moisture curves
= 2 for moisture curves from Huyakorn et al (WRR 20(8) 1984, WRR 22(13) 1986) with Kr=Se**n conductivity relationship
= 3 for moisture curves from Huyakorn et al (WRR 20(8) 1984, WRR 22(13) 1986)with conductivity relationship from Table 3 of 1984 paper (log_10 Kr(Se) curve)
= 4 for Brooks‐Corey moisture curves.
The default is [].
FP : list, optional
Feddes Parameters. The default is [].
[PCANA PCREF PCWLT ZROOT PZ OMGC]
- 'PCANA': float, anaerobiosis point --> h2 (≈ -0.5m)
- 'PCREF': float, field capacity --> h3 (≈ -4m)
- 'PCWLT': float, wilting point --> h4 (≈ -150m)
- 'ZROOT': float, root depth
- 'PZ': float, pz is an empirical parameter
- 'OMGC': float, 0<OMGC<1 Compensatory mechanisms for root water uptake
.. note:
Calculate actual transpiration Ta (m d−1).
A multiplicative reduction factor is defined by four pressure heads (0 > h1 > h2 > h3 > h4)
delimiting five phases of uptake
For details, see http://dx.doi.org/10.1002/2015WR017139
BETA = (1-DEPTH/ZROOT(VEG_TYPE(I)))*DEXP(-1.0d0*PZ(VEG_TYPE(I))*DEPTH/ZROOT(VEG_TYPE(I)))
BTRANI(K) = MAX(ZERO,BETA*DZ*GX)
BTRAN(I) = BTRAN(I) + BETA*DZ
OMG(I) = OMG(I) + GX*BETA*DZ
QTRANIE(K)=ETP(I)*BTRANI(K)/BTRAN(I)/MAX(OMG(I),OMGC(VEG_TYPE(I)))
FP : Dict, optional
Dictionnary containing the SPP properies per zones (root map, indice of vegetation)
SPP : np.array, optional
Array of size (nb of layers*nb of zones)*8 (8: Physical properties)
SPP_map : Dict or Pd.DataFrame, optional
Dictionnary containing the SPP properies per zones
Example for 2 zones:
{
'PERMX': [0.000188, 9.4e-05],
'PERMY': [0.000188, 0.000188],
'PERMZ': [0.000188, 0.000188],
'ELSTOR': [1e-05, 1e-05],
'POROS': [0.55, 0.55],
'VGNCELL': [1.46, 1.46],
'VGRMCCELL': [0.15, 0.15],
'VGPSATCELL': [0.03125, 0.03125]
}
Soil Physical Properties. The default is [].
- 'PERMX' (NSTR, NZONE): saturated hydraulic conductivity - xx
- 'PERMY' (NSTR, NZONE): saturated hydraulic conductivity - yy
- 'PERMZ' (NSTR, NZONE): saturated hydraulic conductivity - zz
- 'ELSTOR' (NSTR, NZONE): specific storage
- 'POROS' (NSTR, NZONE): porosity (moisture content at saturation) = \thetaS
retention curves parameters VGN, VGRMC, and VGPSAT
- 'VGNCELL' (NSTR, NZONE): van Genuchten curve exponent = n
- 'VGRMCCELL' (NSTR, NZONE): residual moisture content = \thetaR
- 'VGPSATCELL' (NSTR, NZONE): van Genuchten curve exponent -->
VGPSAT == -1/alpha (with alpha expressed in [L-1]);
Example for 1 zone:
PERMX PERMY PERMZ ... VGNCELL VGRMCCELL VGPSATCELL
str zone ...
0 0 0.000188 0.000188 0.000188 ... 1.46 0.15 0.03125
1 0 0.000188 0.000188 0.000188 ... 1.46 0.15 0.03125
2 0 0.000188 0.000188 0.000188 ... 1.46 0.15 0.03125
3 0 0.000188 0.000188 0.000188 ... 1.46 0.15 0.03125
4 0 0.000188 0.000188 0.000188 ... 1.46 0.15 0.03125
..note::
At the file level, this is an example of 3 layers and 2 zones; The inner reading cycle is by zone, and the outer one by layers, i.e.:
PERMX_z1_str1 PERMY_z1_str1 PERMZ_z1_str1 ELSTOR_z1_str1 POROS_z1_str1 VGNCELL_z1_str1 VGRMCCELL__z1_str1 VGPSATCELL__z1_str1
PERMX_z2_str1 PERMY_z2_str1 PERMZ_z2_str1 ELSTOR_z2_str1 POROS_z2_str1 VGNCELL_z2_str1 VGRMCCELL__z2_str1 VGPSATCELL__z2_str1
PERMX_z1_str2 PERMY_z1_str2 PERMZ_z1_str2 ELSTOR_z1_str2 POROS_z1_str2 VGNCELL_z1_str2 VGRMCCELL__z1_str2 VGPSATCELL__z1_str2
PERMX_z2_str2 PERMY_z2_str2 PERMZ_z2_str2 ELSTOR_z2_str2 POROS_z2_str2 VGNCELL_z2_str2 VGRMCCELL__z2_str2 VGPSATCELL__z2_str2
PERMX_z1_str3 PERMY_z1_str3 PERMZ_z1_str3 ELSTOR_z1_str3 POROS_z1_str3 VGNCELL_z1_str3 VGRMCCELL__z1_str3 VGPSATCELL__z1_str3
PERMX_z2_str3 PERMY_z2_str3 PERMZ_z2_str3 ELSTOR_z2_str3 POROS_z2_str3 VGNCELL_z2_str3 VGRMCCELL__z2_str3 VGPSATCELL__z2_str3
Make sure all 8 parameters for each layer/zone are on the same line.
Also, make sure NZONE = 2 in dem_parameters and MAXZON in CATHY.H is larger than or equal to NZONE.
Returns
-------
- update parm file
- update CATHY.H file
- update mesh vtk file
"""
if "shellprint_update" in kwargs:
if kwargs["shellprint_update"] is False:
pass
else:
self.console.print(":arrows_counterclockwise: [b]Update soil[/b]")
# set default parameters if SPP and/or FP args are not existing yet
# --------------------------------------------------------------------
if len(self.soil) == 0:
self.set_SOIL_defaults()
try:
df = self.read_inputs('soil',
MAXVEG=self.MAXVEG
)
if len(FP_map)==0:
# FP_map = df[1].reset_index(drop=True).to_dict(orient='list')
FP_map = df[1] #.reset_index(drop=True).to_dict(orient='list')
except:
pass
if len(SPP_map) == 0:
SPP_map = self.set_SOIL_defaults(SPP_map_default=True)
if len(FP_map) == 0:
FP_map = self.set_SOIL_defaults(FP_map_default=True)
# check size of the heteregeneity
# -----------------------------------
if len(zone3d) == 0:
print("homogeneous soil")
# check size of soil properties map versus nb of zones/ nb of layers
# --------------------------------------------------------------------
if isinstance(SPP_map["PERMX"], float):
if self.dem_parameters["nzone"] != 1:
raise ValueError("Wrong number of zones")
else:
self.zone3d=zone3d
self.dem_parameters['nzone'] = len(np.unique(SPP_map.index.get_level_values(0)))
self.update_dem_parameters()
self.update_cathyH(MAXZON=self.dem_parameters['nzone'])
# read function arguments kwargs and udpate soil and parm files
# --------------------------------------------------------------------
for kk, value in kwargs.items():
self.soil[kk] = value
self.parm[kk] = value
# loop over Feddes parameters
# --------------------------------------------------------------------
for fp in FP_map: # loop over fedded parameterssoil_het_dim
self.soil[fp] = FP_map[fp]
# check consistency between parameters
# --------------------------------------------------------------------
if isinstance(SPP_map["PERMX"], float):
for k in SPP_map:
SPP_map[k] = [SPP_map[k]]
if isinstance(SPP_map, dict):
SPP_map_dict = SPP_map
if hasattr(self, 'zone3d'):
num_rows, num_cols = np.shape(self.zone3d)[0],(np.shape(self.zone3d)[1] * np.shape(self.zone3d)[2])
df_SPP_map = self.init_soil_SPP_map_df(num_cols, num_rows)
for i, layersi_zones in enumerate(self.zone3d):
layersi_zones = np.hstack(layersi_zones)
unique_zones = np.unique(layersi_zones)
layersi_zones_prop = np.zeros(np.shape(layersi_zones))
for c in SPP_map_dict.keys():
for zi in unique_zones:
layersi_zones_prop[layersi_zones == zi] = SPP_map_dict[c][int(zi-1)]
df_SPP_map.loc[(slice(None), i + 1), c] = layersi_zones_prop
else:
nzones, nstr = self.dem_parameters["nzone"], self.dem_parameters["nstr"]
df_SPP_map = self.init_soil_SPP_map_df(nzones, nstr)
for key, values in SPP_map.items():
if len(values) == 1:
values_all_layers = values*nstr
self.console.rule(
":warning: Assuming that soil is homogeneous with depth :warning:", style="yellow"
)
df_SPP_map[key] = values_all_layers
SPP_map = df_SPP_map
if (SPP_map["VGRMCCELL"] >= SPP_map["POROS"]).any():
self.console.rule(
":warning: residual water content is > porosity :warning:", style="yellow"
)
# create prepro inputs if not existing (containing info about the DEM)
# --------------------------------------------------------------------
if hasattr(self, "dem_parameters") is False:
self.update_prepo_inputs()
# Soil Physical Properties strat by strat
# --------------------------------------------------------------------
self.soil_SPP = {}
SPP_map = SPP_map.sort_values(by=['layer', 'zone'])
self.soil_SPP["SPP_map"] = SPP_map # mapping with respect to zones
if len(SPP) > 0:
self.soil_SPP["SPP"] = SPP # matrice with respect to zones
else:
SoilPhysProp = SPP_map
self.soil_SPP["SPP"] = SoilPhysProp # matrice with respect to zones/layers
# Vegetation properties (PCANA,PCREF,PCWLT,ZROOT,PZ,OMGC)
# --------------------------------------------------------------------
FeddesParam, SCFParam = self._prepare_SOIL_vegetation_tb(FP_map)
self.soil_FP = {}
self.soil_FP["FP"] = FeddesParam
self.soil_FP["FP_map"] = FP_map # mapping with respect to zones
self.soil_FP["SCF"] = SCFParam # SCF-VEG: per-vegetation-type SCF
map_veg = np.ones(np.shape(self.veg_map))
for i, value in enumerate(FP_map['PCANA']):
map_veg[self.veg_map == i + 1] = i + 1
# Pass maxveg=self.MAXVEG explicitly: FP_map was just validated
# against self.MAXVEG in _prepare_SOIL_vegetation_tb() above, so
# MAXVEG is already correct at this point. Without this, the
# update_veg_map() call below would re-derive MAXVEG from whichever
# classes happen to appear in `map_veg` this month (which can be a
# strict subset of the full class set) and silently overwrite the
# just-validated value — corrupting CATHY.H right before the solve.
self.update_veg_map(map_veg, maxveg=self.MAXVEG)
if show:
update_map_veg = self.map_prop_veg(FP_map)
fig, ax = plt_CT.dem_plot_2d_top(update_map_veg, label="all")
fig.savefig(
os.path.join(self.workdir, self.project_name, "map_veg.png"), dpi=400
)
# write soil file
# --------------------------------------------------------------------
self._write_SOIL_file(self.soil_SPP["SPP"],
FeddesParam,
SCFParam,
**kwargs
)
# map SPP to the mesh
# --------------------------------------------------------------------
dem, dem_header = self.read_inputs('dem')
if len(zone3d) > 0:
saveMeshPath = os.path.join(self.workdir, self.project_name,
'vtk',
self.project_name + '.vtk'
)
if 'saveMeshPath' in kwargs:
saveMeshPath = kwargs.pop('saveMeshPath')
mt.add_markers_zone3d_2_mesh(
zone3d,
dem,
self.mesh_pv_attributes,
self.dem_parameters,
self.hapin,
self.grid3d,
to_nodes=True,
show=show,
saveMeshPath=saveMeshPath,
)
pass
def _get_soil_SPP_columnsNames(self):
# Define the column names
columns = ['PERMX', 'PERMY', 'PERMZ','ELSTOR',
'POROS',
'VGNCELL', 'VGRMCCELL', 'VGPSATCELL'
]
return columns
def _get_soil_FP_columnsNames(self):
# Define the column names
# SCF-VEG: SCF is added as a per-vegetation-type column only for
# versions that opt into the new soil-file format (see
# version_config.py's CATHY_SOIL_FORMAT); "default"/"withIrr" keep
# the classic 6-column Feddes table unchanged. When present, SCF
# is NOT part of the 6-column Feddes table itself -- it's written
# as its own separate block (see _prepare_SOIL_vegetation_tb /
# _write_SOIL_file), so older code reading FeddesParam's 6
# columns is unaffected either way.
columns = ['PCANA', 'PCREF', 'PCWLT',
'ZROOT',
'PZ','OMGC'
]
cfg = VersionConfigManager("cathy", "soil_format", self.version)
if cfg.get("scf_per_veg", False):
columns = columns + ['SCF']
return columns
def init_soil_FP_map_df(self, nveg):
columns = self._get_soil_FP_columnsNames()
FP_map = pd.DataFrame(index=np.arange(1,nveg+1),
columns=columns
)
FP_map.index.name = 'Veg nb'
return FP_map
def init_soil_SPP_map_df(self, nzones, nstr):
columns = self._get_soil_SPP_columnsNames()
# Generate the lists of integers for zones and strings
zones = [ni + 1 for ni in range(nzones)]
strings = [ni + 1 for ni in range(nstr)]
# Create multi-level index
multi_index = pd.MultiIndex.from_product([zones, strings],
names=['zone', 'layer']
)
# Create an empty DataFrame with multi-level index and specified columns
SPP_map = pd.DataFrame(index=multi_index, columns=columns)
return SPP_map
def set_SOIL_defaults(self,
FP_map_default=False,
SPP_map_default=False,
nveg = None
):
cfg = VersionConfigManager("cathy", "soil_format", self.version)
scf_per_veg = cfg.get("scf_per_veg", False)
self.soil = {
"PMIN": -5.0,
"IPEAT": 0,
}
# SCF-VEG: for classic versions (scf_per_veg=False), SCF stays
# exactly where it always was -- a single global scalar written
# on the same line as IPEAT -- so "default"/"withIrr" input decks
# and file format are completely unchanged. Only "SCF_variable"
# omits it here, since it becomes a per-vegetation-type value
# carried in FP_map/FeddesParam instead (see the FP_map_default
# block below and _prepare_SOIL_vegetation_tb), matching
# SCF(VEG_TYPE(I)) on the Fortran side.
if not scf_per_veg:
self.soil["SCF"] = 1.0 # 1 Feddes approach (assume all soil is covered by the vegetation), or 0 if PMIN approach
self.soil.update({
"CBETA0": 0.4,
"CANG": 0.225,
# Feddes parameters default values
"PCANA": [0.0],
"PCREF": [-4.0],
"PCWLT": [-150],
"ZROOT": [0.1],
"PZ": [1.0],
"OMGC": [1.0],
"IVGHU": 0, # The first thing that must be decided!
"HUALFA": 0.02,
"HUBETA": 2,
"HUGAMA": 2,
"HUPSIA": 0,
"HUSWR": 0.333,
"HUN": 1,
"HUA": -5,
"HUB": 1,
"BCBETA": 1.2,
"BCRMC": 0,
"BCPSAT": -0.345,
})
if FP_map_default:
values = dict(PCANA=0.0,
PCREF=-4.0,
PCWLT=-150,
ZROOT=1,
PZ=1,
OMGC=1,
)
# SCF-VEG: only add a per-vegetation-type SCF default for
# versions using the new format. Default value (1.0) matches
# the old global self.soil["SCF"] default (Feddes approach,
# full canopy cover assumed), so behavior is unchanged unless
# the caller sets FP_map["SCF"] explicitly per vegetation type.
if scf_per_veg:
values["SCF"] = 1.0
if nveg is None:
nveg = len(np.unique(self.veg_map))
FP_map = self.init_soil_FP_map_df(nveg)
FP_map.update(pd.DataFrame([values]*len(FP_map), index=FP_map.index))
return FP_map
# # set Soil Physical Properties defaults parameters
# # --------------------------------------------------------------------
if SPP_map_default:
vals = dict(PERMX=1.88e-4, PERMY=1.88e-4, PERMZ=1.88e-4,
ELSTOR=1e-5,
POROS=0.55,
VGNCELL=1.46, VGRMCCELL=0.15, VGPSATCELL=0.03125
)
SPP_map = self.init_soil_SPP_map_df(self.dem_parameters["nzone"],
self.dem_parameters["nstr"]
)
# SPP_map.loc[:, vals.keys()] = pd.DataFrame([vals] * len(SPP_map)
# )
SPP_map.update(pd.DataFrame([vals]*len(SPP_map), index=SPP_map.index))
return SPP_map
pass
def _prepare_SPP_tb(self, SPP, zone3d):
"""
prepare SOIL Physical Properties table
Parameters
----------
SPP : TYPE
DESCRIPTION.
Returns
-------
np.array describing the SoilPhysProp with rows corresponding to the layer.
"""
# check number of zones
if self.dem_parameters["nzone"] > 1 or len(zone3d) > 0:
if len(SPP["PERMX"]) <= 1:
for i, spp in enumerate(SPP):
SPP[spp] = SPP[spp] * np.ones(self.dem_parameters["nzone"])
else:
pass
# check size of the heteregeneity of SPP
# ----------------------------------------
# 1d --> uniform
# 2d --> lateral variations due to zones defined in surface
# 3d --> lateral + vertical variations due to zones and strates
SoilPhysProp = []
if len(zone3d) == 0:
# loop over strates
# -----------------------------------------------------------------
for istr in range(self.dem_parameters["nstr"]):
# loop over zones (defined in the zone file)
# --------------------------------------------------------------
LayeriZonei = np.zeros([self.dem_parameters["nzone"], 8])
for izone in range(self.dem_parameters["nzone"]):
for i, spp in enumerate(SPP):
LayeriZonei[izone, i] = SPP[spp][izone]
SoilPhysProp.append(LayeriZonei)
else:
zones3d_def = np.arange(1, self.hapin["M"] * self.hapin["N"] + 1)
zones3d_def = np.reshape(zones3d_def, [self.hapin["M"], self.hapin["N"]])
self.update_zone(zones3d_def)
# loop over strates
# -----------------------------------------------------------------
for istr in range(self.dem_parameters["nstr"]):
# loop over zones (defined in the zone file)
# --------------------------------------------------------------
LayeriZonei = np.zeros([self.dem_parameters["nzone"], 8])
for izone in range(self.hapin["N"] * self.hapin["M"]):
for i, spp in enumerate(SPP):
flag_zone = int(np.ravel(zone3d[istr])[izone])
LayeriZonei[izone, i] = SPP[spp][flag_zone - 1]
SoilPhysProp.append(LayeriZonei)
SoilPhysProp = np.vstack(SoilPhysProp)
# np.shape(SoilPhysProp)
# case if there is only one zone in the mesh
# -------------------------------------------------------------
else:
if len(SPP["PERMX"]) <= 1:
izoneSoil = []
for spp in SPP:
izoneSoil.append(SPP[spp])
izoneSoil = np.hstack(izoneSoil)
SoilPhysProp = np.tile(izoneSoil, (self.dem_parameters["nstr"], 1))
# case where it is heterogeneous in depth (z)
# ------------------------------------------------------
else:
izoneSoil_per_layer = []
for stri in range(self.dem_parameters["nstr"]):
izoneSoil = []
for spp in SPP:
izoneSoil.append(SPP[spp][stri])
izoneSoil_per_layer.append(np.hstack(izoneSoil))
SoilPhysProp = np.vstack(izoneSoil_per_layer)
return SoilPhysProp
def _prepare_SOIL_vegetation_tb(self, FP_map):
"""
_prepare_SOIL_vegetation_tb
Parameters
----------
FP_map : dict
dict containing Feddes parameters.
- 'PCANA': anaerobiosis point
- 'PCREF': field capacity
- 'PCWLT': wilting point
- 'ZROOT': float root depth
- 'PZ': float, can be negative, defines the root profile
- 'OMGC': float, compensatory mechanism
- 'SCF': float, 0<SCF<1, canopy/soil cover fraction, per
vegetation type (SCF-VEG). Written to a SEPARATE block in
the soil file from the other six -- see _write_SOIL_file --
not as a 7th column of the Feddes table, so the Feddes
table's on-disk shape is unchanged for any code still
expecting 6 columns there.
For details, see http://dx.doi.org/10.1002/2015WR017139
Returns
-------
FeddesParam: numpy array, shape (MAXVEG, 6)
table describing the six Feddes parameters for a given DEM
SCFParam: numpy array, shape (MAXVEG, 1), or None
per-vegetation-type canopy/soil cover fraction (SCF-VEG).
None for versions that don't use the per-vegetation-type SCF
format (see version_config.py's CATHY_SOIL_FORMAT) -- in that
case SCF stays the classic single global scalar and this
function's behavior is otherwise unchanged from before
SCF-VEG was introduced.
"""
cfg = VersionConfigManager("cathy", "soil_format", self.version)
scf_per_veg = cfg.get("scf_per_veg", False)
# Vegetation properties (PCANA,PCREF,PCWLT,ZROOT,PZ,OMGC) [+ SCF]
# --------------------------------------------------------------------
# Check if root_map file exist and is updated
# -------------------------------------------
if hasattr(self, "veg_map") is False:
if len(FP_map[list(FP_map)[0]]) == 1:
self.update_veg_map()
else:
raise ValueError(
"Found multiple values of Feddes zones"
+ "but vegetation map is not defined"
)
# Check vegetation heterogeneity dimension
# ----------------------------------------
if self.cathyH["MAXVEG"] != len(FP_map[list(FP_map)[0]]):
raise ValueError(
"Wrong number of vegetations: PCANA size is "
+ str(len(FP_map[list(FP_map)[0]]))
+ " while MAXVEG is "
+ str(self.cathyH["MAXVEG"])
)
# SCF-VEG: the six Feddes columns are named explicitly here
# (rather than iterating "for sfp in FP_map") so FeddesParam keeps
# its fixed (MAXVEG, 6) shape even when scf_per_veg=True and
# FP_map carries an extra 'SCF' column -- iterating every FP_map
# column into a single hstack, as the original code did, would
# then silently try to pack 7 values into a 6-wide row.
feddes_cols = ['PCANA', 'PCREF', 'PCWLT', 'ZROOT', 'PZ', 'OMGC']
SCFParam = None
# check number of vegetation
# --------------------------------------------------------------------
if self.cathyH["MAXVEG"] > 1:
FeddesParam = np.zeros([self.cathyH["MAXVEG"], 6])
if scf_per_veg:
SCFParam = np.zeros([self.cathyH["MAXVEG"], 1])
# FP_map is indexed 1..MAXVEG (see init_soil_FP_map_df: "Veg nb"
# index starts at 1, matching the 1..MAXVEG vegetation class
# labels used everywhere else, e.g. in the root_map raster).
# Loop over those same 1-based labels for .loc, and only drop
# to a 0-based offset for the FeddesParam array itself.
for iveg in range(1, self.cathyH["MAXVEG"] + 1): # loop over veg zones within a strate
izoneVeg_tmp = [FP_map[sfp].loc[iveg] for sfp in feddes_cols]
izoneVeg_tmp = np.hstack(izoneVeg_tmp)
FeddesParam[iveg - 1, :] = izoneVeg_tmp
if scf_per_veg:
SCFParam[iveg - 1, 0] = FP_map['SCF'].loc[iveg]
# case where unique vegetation type
else:
FeddesParam = np.c_[
self.soil["PCANA"],
self.soil["PCREF"],
self.soil["PCWLT"],
self.soil["ZROOT"],
self.soil["PZ"],
self.soil["OMGC"],
]
if scf_per_veg:
# SCF-VEG: SCF is no longer a fixed-position key in
# self.soil (see set_SOIL_defaults), so fall back to
# FP_map directly, then to the old global default (1.0,
# full canopy cover) if neither is available, to preserve
# behavior for any caller that never sets SCF at all.
if 'SCF' in FP_map:
SCFParam = np.c_[FP_map["SCF"]]
else:
SCFParam = np.c_[[self.soil.get("SCF", 1.0)]]
return FeddesParam, SCFParam
def _write_SOIL_file(self, SoilPhysProp, FeddesParam, SCFParam, **kwargs):
"""
_write_SOIL_file
Parameters
----------
SoilPhysProp : Soil physical properties
Numpy array of Soil physical properties.
FeddesParam : Feddes Parameters
Dictionnatry of Feddes Parameters
SCFParam : SCF-VEG
numpy array, shape (MAXVEG, 1), or None. When provided
(scf_per_veg version), written as its own block right after
the Feddes table (matching the Fortran side's
READ(IIN4,*) (SCF(I),I=1,NVEG), read immediately after the
6-column Feddes table and before IVGHU). When None (classic
versions), no such block is written and the file format is
byte-for-byte the original one (SCF stays on the IPEAT line).
"""
# backup file during DA scheme cycle
# --------------------------------------------------------------------
backup = True
if "backup" in kwargs:
backup = kwargs["backup"]
# number of side header for each row
# SCF-VEG: classic format (SCFParam is None) keeps the original
# [1, 2, 2, 6, 1, 5, 1, 2, 3] -- IPEAT,SCF on one line, no extra
# block, byte-for-byte the file format from before SCF-VEG. The
# scf_per_veg format instead reads IPEAT alone (h=1) and inserts
# a new index 4 (h=0) right after the Feddes table (index 3) to
# write SCFParam -- h=0 because it doesn't consume any self.soil
# dict entries, it's special-cased below exactly like the Feddes
# table is.
if SCFParam is None:
header_fmt_soil = [1, 2, 2, 6, 1, 5, 1, 2, 3]
else:
header_fmt_soil = [1, 1, 2, 6, 0, 1, 5, 1, 2, 3]
# open soil file
# --------------------------------------------------------------------
if self.DAFLAG:
soil_filepath = os.path.join(os.getcwd(), self.input_dirname, "soil")
else:
soil_filepath = os.path.join(
self.workdir, self.project_name, self.input_dirname, "soil"
)
if "filename" in kwargs:
soil_filepath = os.path.join(kwargs["filename"]) #, "soil")
if backup:
if self.count_DA_cycle is not None:
dst_dir = soil_filepath + str(self.count_DA_cycle)
shutil.copy(soil_filepath, dst_dir)
# print('backup soil ' + dst_dir)
with open(os.path.join(soil_filepath), "w+") as soilfile:
counth = 0 # count header index
# Write line by line according to header format
# ----------------------------------------------------------------
for i, h in enumerate(header_fmt_soil):
# left = values
# ------------------------------------------------------------
left = right = []
left = str(list(self.soil.values())[counth : counth + h])
left = left.strip("[]").replace(",", "")
# right = keys
# ------------------------------------------------------------
right = str(list(self.soil.keys())[counth : counth + h])
right = right.strip("[]").replace(",", "")
right = right.replace("'", "")
# Line 4: write Feddes Parameters Table
# ------------------------------------------------------------
if i == 3: # Feddes parameters
np.savetxt(soilfile, FeddesParam, fmt="%1.3e")
counth += h
# SCF-VEG: write the per-vegetation-type SCF block
# immediately after the Feddes table, as its own separate
# lines -- NOT merged into the Feddes table's columns --
# so the Feddes table's own shape/format is unchanged for
# anything else that reads it. Guarded on SCFParam is not
# None: in the classic 9-slot header_fmt_soil, index 4 is
# the ordinary IVGHU line, not this block, and must fall
# through to the generic "else" branch below.
elif i == 4 and SCFParam is not None: # SCF, per vegetation type
np.savetxt(soilfile, SCFParam, fmt="%1.3e")
counth += h
else:
line = left + "\t" + right + "\n"
counth += h
soilfile.write(str(line))
# End of soil file: write Soil Properties Table
# ------------------------------------------------------------
np.savetxt(soilfile, SoilPhysProp, fmt="%1.3e")
soilfile.write(
"PERMX PERMY PERMZ ELSTOR POROS,VGNCELL,VGRMCCELL,VGPSATCELL" + "\n"
)
soilfile.close()
def _get_DEM_for_shape(self, shape):
"""
Return self.DEM, reloading it from disk if it's missing OR its shape
doesn't match `shape`.
self.DEM is cached on first use (`if not hasattr(self, "DEM")`) in
several places (update_veg_map, _check_outside_DEM). That caching is
unsafe across a batch run that reuses the same `simu` object (or
carries over cached attributes) for multiple projects/grids: a DEM
loaded for an earlier, differently-shaped project silently survives
into the next project's calls. Boolean-masking a differently-shaped
raster with it then either raises IndexError (shape mismatch caught)
or, worse, silently produces a wrong-sized mask if the shapes happen
to be broadcastable — corrupting MAXVEG/MAXZON without any error at
all. Always validating against the raster actually being processed,
and reloading when it disagrees, closes both failure modes.
"""
if not hasattr(self, "DEM") or np.shape(self.DEM) != tuple(shape):
DEM_mat, DEM_header = in_CT.read_dem(
os.path.join(self.workdir, self.project_name, "prepro/dem"),
os.path.join(self.workdir, self.project_name, "prepro/dtm_13.val"),
)
self.DEM = DEM_mat
if np.shape(self.DEM) != tuple(shape):
raise ValueError(
f"DEM on disk for project '{self.project_name}' has "
f"shape {np.shape(self.DEM)}, but the raster being "
f"checked against it has shape {tuple(shape)}. These "
f"must match — check that prepro/dem for this project "
f"is actually the one you expect (a stale/wrong "
f"project directory would produce exactly this)."
)
return self.DEM
def update_veg_map(self, indice_veg=None, show=False, maxveg=None, **kwargs):
"""
Update the vegetation raster (root_map / IIN3) describing which
vegetation class each cell belongs to, and keep MAXVEG in sync.
Parameters
----------
indice_veg : raster or int, optional
New vegetation distribution. An int broadcasts that single class
over the whole DEM grid. If None, the existing root_map file is
re-read from disk. Default is None.
show : bool, optional
Plot the vegetation raster. Default is False.
maxveg : int, optional
Fix the number of vegetation zones (MAXVEG) explicitly instead of
deriving it from this raster's classes. Useful for a time series
of rasters that share a fixed set of semantic classes (e.g. 1..5)
but don't all spatially contain every class.
When given, this method also guarantees that every class in
1..maxveg is actually present in root_map inside the valid DEM
footprint, injecting a single placeholder pixel for any class
that's spatially absent this raster. This matters because
CATHY's grid-generation step (and update_cathyH's consistency
check) derive NVEG from the classes actually found in root_map
inside the valid DEM domain, not from self.MAXVEG — without
this, a raster missing a class would leave NVEG < maxveg,
desyncing root_map from the SOIL table (built for the fixed
maxveg) and crashing the Fortran run at the IVGHU read.
If None and this method already ran once for this project
(self.veg_map is set), the resulting self.MAXVEG is reused
rather than recomputed — MAXVEG is "sticky" once a *real*
vegetation map has been established, so no call site can
silently shrink or grow it mid-run. Auto-derivation only
happens on the first genuine call for a project. Note this is
deliberately NOT keyed off hasattr(self, "MAXVEG") alone:
update_cathyH() parks a provisional self.MAXVEG = 1 before any
real vegetation map exists (and leaves self.veg_map unset while
doing so) — treating that placeholder as sticky would lock in
the wrong MAXVEG the first time a real map is supplied. To
force a recompute after a genuine call, pass `maxveg`
explicitly or `del self.MAXVEG` first.
Returns
-------
indice_veg : raster
The vegetation distribution raster that was written.
"""
if maxveg is None and hasattr(self, "veg_map"):
maxveg = self.MAXVEG
if indice_veg is None:
indice_veg, _ = in_CT.read_root_map(
os.path.join(
self.workdir, self.project_name, self.input_dirname, "root_map"
)
)
if not hasattr(self, "hapin"):
self.update_prepo_inputs()
indice_veg = self._veg_raster_from_input(indice_veg)
if maxveg is not None:
indice_veg = self._inject_missing_veg_classes(indice_veg, maxveg)
self.veg_map = indice_veg
self._write_root_map(indice_veg)
self._check_veg_classes(indice_veg, maxveg)
self.MAXVEG = self._resolve_maxveg(indice_veg, maxveg)
self.update_cathyH(MAXVEG=self.MAXVEG)
if show:
ax = plt_CT.show_indice_veg(self.veg_map, **kwargs)
return indice_veg, ax
return indice_veg
def _veg_raster_from_input(self, indice_veg):
"""
Normalize `indice_veg` into a full (M, N) raster, broadcasting a
single int class over the whole grid if needed.
Always returns a fresh copy: `indice_veg` may be the exact array
object the caller is still holding onto (e.g. a script that reuses
its `veg_map` variable right after calling `update_veg_map` for
plotting or history), so mutating it here would leak changes back
to the caller.
"""
if isinstance(indice_veg, int):
return (
np.c_[np.ones([int(self.hapin["M"]), int(self.hapin["N"])])]
* indice_veg
)
return np.array(indice_veg, copy=True)
def _inject_missing_veg_classes(self, indice_veg, maxveg, margin=3, patch_size=2):
"""
Guarantee every class in 1..maxveg is physically present in
`indice_veg`, inside the valid DEM footprint (DEM != -9999),
injecting a small placeholder block for any class that's
spatially absent.
root_map (IIN3) is what CATHY's grid-generation step reads to
derive NVEG for the "grid" file (IIN2), and that derivation counts
the classes it actually finds in root_map inside the valid DEM
domain — it does NOT know about self.MAXVEG. Without this, a
raster missing one class (e.g. only 1 of 2 present this month)
would produce NVEG=1 in the grid file while the SOIL file — sized
from the fixed maxveg=2 — still has 2 vegetation rows. datin.f's
sequential reads then desynchronize: the vegetation-parameter loop
(bounded by NVEG=1) stops one SOIL line short, and the very next
READ (IVGHU) consumes that leftover SOIL row instead, crashing
with "Bad integer for item 1 in list input" (the row's first token
is a real number, not an integer).
Two properties matter for *where* and *how much* we inject, and an
earlier version of this method got both wrong:
- A single pixel is not safe. CATHY's DEM-mode preprocessing
(RAST_INPUT_DEM / TRIANGOLI) turns the raw raster into an FE
mesh; a lone placeholder pixel is exactly the kind of
salt-and-pepper noise that boundary handling / cell-to-node
resampling in that step can dilute or drop before NVEG is ever
computed from the mesh. We now inject a `patch_size` x
`patch_size` block instead.
- Placement must not be "the first valid cell found". Scanning
`np.where(valid_mask)` in row-major order returns cells starting
at the very edge of the valid DEM footprint (wherever validity
first appears scanning from the array's origin) — i.e. cells
directly adjacent to the novalue boundary, which is the single
riskiest place to put a placeholder, since edge nodes are most
exposed to whatever boundary/novalue handling the mesh-generation
step applies. We now require placeholders to sit `margin` pixels
away from any novalue cell (via a distance transform on
valid_mask), so they land deep in the interior instead.
Mutates and returns `indice_veg` in place — callers must pass an
array they own (e.g. the fresh copy from `_veg_raster_from_input`),
never the caller's original array.
"""
DEM_here = self._get_DEM_for_shape(np.shape(indice_veg))
valid_mask = DEM_here != -9999
present = set(np.unique(indice_veg[valid_mask]).astype(int).tolist())
missing = [c for c in range(1, int(maxveg) + 1) if c not in present]
if not missing:
return indice_veg
interior_mask = distance_transform_edt(valid_mask) > margin
interior_rows, interior_cols = np.where(interior_mask)
if len(interior_rows) == 0:
# Fall back to the raw valid mask if the domain is too small/
# thin for any cell to be `margin` pixels from a novalue cell —
# still better than raising outright, but flag it clearly.
print(
f"WARNING: no DEM cell is >{margin} pixels from a novalue "
f"cell — falling back to the valid-mask boundary to inject "
f"missing vegetation class(es) {missing}. Placeholders may "
f"still be lost by mesh-generation boundary handling; "
f"consider lowering `margin`."
)
interior_rows, interior_cols = np.where(valid_mask)
if len(interior_rows) == 0:
print(
f"WARNING: no valid DEM cells at all — cannot inject "
f"missing vegetation class(es) {missing}. root_map/SOIL "
f"NVEG mismatch is still possible for this raster."
)
return indice_veg
# Spread placeholder blocks evenly across the interior candidates
# instead of clustering them all at the first cells found, so
# distinct classes don't collide in the same block.
order = np.argsort(interior_rows * indice_veg.shape[1] + interior_cols)
interior_rows, interior_cols = interior_rows[order], interior_cols[order]
step = max(1, len(interior_rows) // (len(missing) + 1))
n_rows, n_cols = indice_veg.shape
for k, cls in enumerate(missing, start=1):
idx = min(k * step, len(interior_rows) - 1)
r, c = interior_rows[idx], interior_cols[idx]
r0, r1 = r, min(r + patch_size, n_rows)
c0, c1 = c, min(c + patch_size, n_cols)
indice_veg[r0:r1, c0:c1] = cls
print(
f" → injected {patch_size}x{patch_size} placeholder block(s) "
f"for vegetation class(es) {missing} (absent this raster), "
f"{margin}+ pixels inside the DEM interior, so root_map keeps "
f"NVEG == MAXVEG == {maxveg}"
)
return indice_veg
def _write_root_map(self, indice_veg):
"""Write the vegetation raster to the root_map input file (IIN3)."""
root_map_path = os.path.join(
self.workdir, self.project_name, self.input_dirname, "root_map"
)
with open(root_map_path, "w+") as rootmapfile:
rootmapfile.write("north: 0\n")
rootmapfile.write(f"south: {self.hapin['yllcorner']}\n")
rootmapfile.write("east: 0\n")
rootmapfile.write(f"west: {self.hapin['xllcorner']}\n")
rootmapfile.write(f"rows: {self.hapin['M']}\n")
rootmapfile.write(f"cols: {self.hapin['N']}\n")
np.savetxt(rootmapfile, indice_veg, fmt="%i")
def _check_veg_classes(self, indice_veg, maxveg):
"""
Sanity-check vegetation class labels (warns via print, does not
raise, matching the previous behaviour).
- If `maxveg` is fixed, this raster may contain any subset of
{1..maxveg} — it need not include every class, nor be consecutive
from 1 (e.g. only classes {2, 4} present this month is fine).
Only the label range is checked.
- Otherwise, classes must start at 1 and be consecutive.
"""
unique_veg = np.sort(np.unique(indice_veg))
if maxveg is not None:
if unique_veg.min() < 1 or unique_veg.max() > maxveg:
print(
f"Vegetation classes must be within 1..{maxveg} (maxveg).\n"
f"Found: {unique_veg.tolist()}"
)
else:
expected = np.arange(1, len(unique_veg) + 1)
if not np.array_equal(unique_veg, expected):
print(
f"Vegetation classes must start at 1 and be consecutive.\n"
f"Found: {unique_veg.tolist()}\n"
f"Expected: {expected.tolist()}"
)
def _resolve_maxveg(self, indice_veg, maxveg):
"""
Determine MAXVEG the same way CATHY's Fortran side does: datin.f
sets NVEG = MAX(VEG_TYPE) over cells inside the valid DEM domain
(DEM != -9999) — the *maximum class id*, not a count of distinct
classes present.
If `maxveg` is given, it's treated as authoritative for the whole
run and returned as-is (this is what keeps MAXVEG, and every SOIL
vegetation table built against it, consistent across a time series
where a class can be temporarily absent from a given raster).
Otherwise it's derived from this raster via the same max-over-valid-
domain rule datin.f uses, so it can't desync from what CATHY itself
will compute.
"""
if maxveg is not None:
return maxveg
DEM_here = self._get_DEM_for_shape(np.shape(indice_veg))
valid_mask = DEM_here != -9999
veg_valid = indice_veg[valid_mask]
if veg_valid.size == 0:
raise ValueError(
"update_veg_map: no valid (non -9999) DEM cells to "
"derive MAXVEG from."
)
return int(np.max(veg_valid))
def _check_outside_DEM(self, raster2check):
DEM_here = self._get_DEM_for_shape(np.shape(raster2check))
valid_mask = DEM_here != -9999
veg_all = np.unique(raster2check)
veg_valid = np.unique(raster2check[valid_mask])
# Number of classes present in the full raster but absent from the
# valid (non -9999) DEM area — i.e. classes that only occur outside
# the domain. Previously this was `int(len(veg_valid) < len(veg_all))`,
# a boolean capped at 1 regardless of how many classes were actually
# excluded, which silently under-counted whenever more than one
# class fell entirely outside the DEM mask.
exclude_out_ind = len(veg_all) - len(veg_valid)
# # Remap vegetation IDs to consecutive values starting at 1
# veg_mapping = {old: new for new, old in enumerate(sorted(veg_valid), start=1)}
# raster2check_new = raster2check.copy()
# for old_id, new_id in veg_mapping.items():
# raster2check_new[raster2check == old_id] = new_id
# return exclude_out_ind, veg_valid, veg_mapping, raster2check_new
return exclude_out_ind, veg_valid
#%% Add inputs and outputs attributes to the mesh
def init_boundary_conditions(self, BCtypName, time, **kwargs):
"""
.. note:
The boundary conditions are defined in the nansfdirbc (Dirichlet),
nansfneubc (Neumann), and sfbc (seepage face) files.
We have two types of boundary conditions (BC):
- Neumann BC (or specifed flux)
- Dirichlet BC (or pressure).
.. note:
- Pioggia: condizioni di Neumann. Quando non ci può più essere
infiltrazione metto Dirichlet.
- Evaporazione: si indica un limite di pressione minimo ( Pmin ) al di
sotto del quale si ha uno switch da Neumann a Dirichlet
(in quanto al di sotto di questo valore non si ha più evapotraspirazione).
.. note:
The boundary condition for any given surface node can switch between a
Dirichlet condition and a Neumann condition depending on the saturation
(or pressure) state of that node.
.. note:
A Neumann (or specified flux) boundary condition corresponds to
atmosphere-controlled infiltration or exfiltration, with the flux equal
to the rainfall or potential evaporation rate given by the atmospheric input data.
When the surface node reaches a threshold level of saturation or moisture deficit,
the boundary condition is switched to a Dirichlet (specified head) condition,
and the infiltration or exfiltration process becomes soil limited [1].
Returns
-------
None.
"""
# self.update_nansfdirbc()
# self.update_nansfneubc()
# self.update_sfbc()
try:
self.console.print(
":orange_square: [b] init boundary condition dataframe [/b]"
)
# self.create_mesh_vtk()
if not "nnod3" in self.grid3d.keys():
self.run_processor(IPRT1=3)
if hasattr(self, "mesh_bound_cond_df") is False:
self.create_mesh_bounds_df(
BCtypName,
self.grid3d["mesh3d_nodes"],
time,
**kwargs
)
except:
raise ValueError("cannot init boundary conditions dataframe")
pass
pass
def check_for_inconsistent_BC(self):
pass
def set_BC_laterals(self, time, BC_type="", val=0):
"""
Set all sides expect surface one
"""
nnodes = len(
self.mesh_bound_cond_df[self.mesh_bound_cond_df["time (s)"] == 0]["id_node"]
)
for tt in time:
BC_bool_name = []
BC_bool_val = []
for id_node in range(nnodes):
if self.mesh_bound_cond_df["bound"].loc[int(id_node)] == True:
BC_bool_name.append(BC_type)
BC_bool_val.append(0)
else:
BC_bool_name.append(None)
BC_bool_val.append(val)
self.update_mesh_boundary_cond(
time=tt, BC_bool_name=BC_bool_name, BC_bool_val=BC_bool_val
)
self.check_for_inconsistent_BC()
# self.update_mesh_vtk(BC_bool_name,BC_bool_val)
pass
def get_outer_nodes(self, x, y, z):
x_min, x_max = np.min(x), np.max(x)
y_min, y_max = np.min(y), np.max(y)
z_min, z_max = np.min(z), np.max(z)
outer_mask = np.logical_or.reduce((x == x_min, x == x_max, y == y_min, y == y_max, z == z_min, z == z_max))
outer_nodes = np.column_stack((x[outer_mask], y[outer_mask], z[outer_mask]))
return outer_nodes, outer_mask
def create_mesh_bounds_df(self, BCtypName, grid3d, times, **kwargs):
"""
Create a dataframe with flag for different boundary condtions assigned to each nodes
Parameters
----------
grid3d : TYPE
DESCRIPTION.
Returns
-------
None.
"""
mesh_bound_np = np.c_[np.arange(0,len(grid3d)),grid3d]
self.mesh_bound_cond_df = pd.DataFrame(mesh_bound_np)
self.mesh_bound_cond_df = self.mesh_bound_cond_df.rename(
columns={
0: "id_node",
1: "x",
2: "y",
3: "z",
}
)
self.mesh_bound_cond_df["id_node"] = self.mesh_bound_cond_df["id_node"].astype(int)
self.mesh_bound_cond_df["time"] = 0
top_id = np.arange(0, self.grid3d['nnod'], 1)
top_mask = np.zeros(len(self.grid3d['mesh3d_nodes']), dtype=bool)
top_mask[top_id.astype(int)] = True
bot_id = np.arange(self.grid3d['nnod3']-self.grid3d['nnod'],
self.grid3d['nnod3'],
1
)
bot_mask = np.zeros(len(self.grid3d['mesh3d_nodes']), dtype=bool)
bot_mask[bot_id.astype(int)] = True
x_threshold = self.grid3d['mesh3d_nodes'][:,0].min()
xmin_side_mask = self.mesh_bound_cond_df['x'] <= x_threshold
x_threshold = self.grid3d['mesh3d_nodes'][:,0].max()
xmax_side_mask = self.mesh_bound_cond_df['x'] >= x_threshold
y_threshold = self.grid3d['mesh3d_nodes'][:,1].min()
ymin_side_mask = self.mesh_bound_cond_df['y'] <= y_threshold
y_threshold = self.grid3d['mesh3d_nodes'][:,1].max()
ymax_side_mask = self.mesh_bound_cond_df['y'] >= y_threshold
all_sides_mask = np.c_[xmin_side_mask,
xmax_side_mask,
ymin_side_mask,
ymax_side_mask,
]
all_sides_mask = np.any(all_sides_mask, axis=1)
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
# # Scatter plot in 3D
# sc = ax.scatter(self.grid3d['mesh3d_nodes'][:, 0],
# self.grid3d['mesh3d_nodes'][:, 1],
# self.grid3d['mesh3d_nodes'][:, 2],
# # c=self.grid3d['mesh3d_nodes'][:, 2], # Coloring by z-axis values
# c='k', # Coloring by z-axis values
# # cmap='viridis'
# ) # Optional: Change color map
# # Scatter plot in 3D
# sc = ax.scatter(self.grid3d['mesh3d_nodes'][all_sides_mask, 0],
# self.grid3d['mesh3d_nodes'][all_sides_mask, 1],
# self.grid3d['mesh3d_nodes'][all_sides_mask, 2],
# # c=self.grid3d['mesh3d_nodes'][left_side_mask, 2], # Coloring by z-axis values
# c='r', # Coloring by z-axis values
# # cmap='viridis'
# ) # Optional: Change color map
# Step 4: Fill the dataframe with flag for outer nodes
# -------------------------------------------------------------------
# print('Fill the dataframe with flag for outer nodes')
self.mesh_bound_cond_df["noflow_bound"] = np.ones(len(self.grid3d['mesh3d_nodes']))*True
self.mesh_bound_cond_df["top_bound"] = top_mask
self.mesh_bound_cond_df["bot_bound"] = bot_mask
self.mesh_bound_cond_df["xmin_bound"] = xmin_side_mask
self.mesh_bound_cond_df["ymin_bound"] = ymin_side_mask
self.mesh_bound_cond_df["xmax_bound"] = xmax_side_mask
self.mesh_bound_cond_df["ymax_bound"] = ymax_side_mask
self.mesh_bound_cond_df["all_sides"] = all_sides_mask
# Step 4: Replicate df for all given times
# -------------------------------------------------------------------
if len(times) > 1:
mesh_bound_cond_df_withtimes = pd.concat(
[self.mesh_bound_cond_df.assign(time=ti) for ti in times]
).reset_index(drop=True)
self.mesh_bound_cond_df = mesh_bound_cond_df_withtimes
# ------------------------------------------------------------------
# "Open outlet" Neumann BC helpers
#
# Thin, self-contained geometry helpers used by an "open outlet on
# one lateral side + optional free-bottom-drainage" Neumann BC
# (see e.g. Agramon_withLAI_withETp.py's configure_boundary_conditions()).
# All five assume create_mesh_bounds_df() has already been called
# (mesh_bound_cond_df must exist) except check_neumann_flux_vs_ks,
# which only needs self.soil_SPP / self.soil.
# ------------------------------------------------------------------
def _first_time_bounds_df(self):
"""mesh_bound_cond_df filtered to its first time block only -
the boundary flags (top/bot/xmin/.../all_sides) don't vary
with time, so any one block has the full node-flag picture."""
if not hasattr(self, "mesh_bound_cond_df"):
raise RuntimeError(
"mesh_bound_cond_df not built yet - call "
"create_mesh_bounds_df() first."
)
return self.mesh_bound_cond_df[
self.mesh_bound_cond_df["time"]
== self.mesh_bound_cond_df["time"].iloc[0]
]
def get_outlet_side(self, outlet_side=None):
"""
Resolve which lateral mesh side ("xmin", "xmax", "ymin",
"ymax") acts as the catchment outlet for an open-outlet
Neumann BC.
If `outlet_side` is given (any of the 4 labels, case-
insensitive), it's validated and returned as-is - no auto-
detection performed.
Otherwise the side is auto-detected as the lateral boundary
whose top-surface (DEM) nodes have the lowest mean elevation,
i.e. the topographically lowest edge of the domain - a
reasonable proxy for "where water would naturally leave" when
the true pour point isn't an interior node reachable by
re-shaping a whole mesh side.
Returns
-------
str : one of "xmin", "xmax", "ymin", "ymax"
"""
valid = {"xmin", "xmax", "ymin", "ymax"}
if outlet_side is not None:
side = str(outlet_side).lower()
if side not in valid:
raise ValueError(
f"outlet_side must be one of {sorted(valid)} "
f"(got '{outlet_side}')"
)
return side
df0 = self._first_time_bounds_df()
surface = df0[df0["top_bound"]]
means = {}
for side, col in (("xmin", "xmin_bound"), ("xmax", "xmax_bound"),
("ymin", "ymin_bound"), ("ymax", "ymax_bound")):
side_nodes = surface[surface[col]]
if len(side_nodes) == 0:
continue
means[side] = float(side_nodes["z"].mean())
if not means:
raise RuntimeError("get_outlet_side: no lateral boundary nodes found.")
resolved = min(means, key=means.get)
others = ", ".join(f"{s}={z:.3f}" for s, z in means.items() if s != resolved)
print(f" Outlet side auto-detected: '{resolved}' "
f"(mean top elevation {means[resolved]:.3f} m"
+ (f"; other sides: {others}" if others else "") + ")")
return resolved
def get_outlet_node(self):
"""
Diagnostic surface node id (0-based) at the catchment's
lowest-elevation top-surface point. Used only for
reporting/plotting (e.g. tracking pressure head "at the
outlet" through a spin-up) - independent of whichever lateral
side get_outlet_side() resolves to, since the pour point
itself need not sit exactly on that side.
Returns
-------
int : 0-based node id
"""
df0 = self._first_time_bounds_df()
surface = df0[df0["top_bound"]]
if len(surface) == 0:
raise RuntimeError("get_outlet_node: no top-surface nodes found.")
outlet_row = surface.loc[surface["z"].idxmin()]
return int(outlet_row["id_node"])
def compute_nodal_face_areas(self, side):
"""
Per-node vertical face area [m^2] for every 3-D mesh node on
the given lateral side ("xmin"/"xmax"/"ymin"/"ymax"), used to
convert a flux density [m/s] into a volumetric rate [m^3/s]
per node for the Neumann BC (nansfneubc expects m^3/s, not
m/s - see update_nansfneubc's docstring).
Geometry assumption: the mesh is a structured DEM-extrusion -
`mesh3d_nodes` is ordered by node-layer (create_mesh_bounds_df
flags the first `nnod` rows as "top_bound", the last `nnod`
as "bot_bound"), with `nnod3 / nnod` node-layers sharing the
same (x, y) per surface column, i.e. node id `i` and node id
`i + nnod` are vertically stacked in the same column. For
every boundary column, each node is given half the vertical
distance to its neighbour above plus half the distance to its
neighbour below (the full distance at the very top/bottom of
the column - a standard 1-D control-volume split), multiplied
by the along-side horizontal spacing (`delta_y` for xmin/xmax
faces, `delta_x` for ymin/ymax faces, from `self.hapin`).
Simplification: every boundary column - including the two
corner columns shared with an adjacent side - gets the full
along-side spacing rather than half of it at the domain's far
ends. This over-counts the side's total area by up to one
node-spacing at each end; fine for a tunable Neumann outflow,
not meant to be an exact quadrature.
Returns
-------
dict {0-based node id: area [m^2]}
"""
side = side.lower()
side_col = {
"xmin": "xmin_bound", "xmax": "xmax_bound",
"ymin": "ymin_bound", "ymax": "ymax_bound",
}[side]
along_spacing = (self.hapin["delta_y"] if side in ("xmin", "xmax")
else self.hapin["delta_x"])
df0 = self._first_time_bounds_df()
side_ids = df0.loc[df0[side_col], "id_node"].to_numpy(dtype=int)
if len(side_ids) == 0:
raise RuntimeError(f"compute_nodal_face_areas: no nodes flagged "
f"'{side_col}'.")
nnod = int(self.grid3d["nnod"])
nnod3 = int(self.grid3d["nnod3"])
n_layers = nnod3 // nnod
z_all = self.grid3d["mesh3d_nodes"][:, 2]
side_id_set = set(side_ids.tolist())
areas = {}
# nnod/nnod3 can come back as numpy floats from some grid3d
# readers, which would silently upcast side_ids % nnod (and
# everything derived from it) to float - force int here so
# these stay valid array indices.
columns = np.unique(side_ids % nnod).astype(np.int64)
for col in columns:
layer_node_ids = (col + np.arange(n_layers) * nnod).astype(np.int64)
z = z_all[layer_node_ids]
dz = np.abs(np.diff(z))
if len(dz) == 0:
thickness = np.array([0.0])
else:
thickness = np.zeros(n_layers)
thickness[0] = dz[0] / 2
thickness[-1] = dz[-1] / 2
if n_layers > 2:
thickness[1:-1] = (dz[:-1] + dz[1:]) / 2
for node_id, th in zip(layer_node_ids, thickness):
if int(node_id) in side_id_set:
areas[int(node_id)] = float(th * along_spacing)
return areas
def compute_bottom_face_area(self, res_x, res_y):
"""
Per-node horizontal face area [m^2] for a bottom-layer node,
used to convert a free-drainage flux density [m/s] into a
volumetric rate [m^3/s] per bottom node.
Assumes one bottom node per DEM raster cell (the standard
CATHY DEM-extrusion mesh), so each bottom node's tributary
area is simply the raster cell area `res_x * res_y`. Edge/
corner bottom nodes are given the same area as interior ones
(same simplification as compute_nodal_face_areas) rather than
a half- or quarter-cell - fine for a tunable Neumann outflow.
Parameters
----------
res_x, res_y : float
DEM pixel size in m (x and y directions).
Returns
-------
float : area [m^2], same for every bottom node.
"""
return float(res_x) * float(res_y)
def check_neumann_flux_vs_ks(self, flux, label="flux"):
"""
Sanity check only - never blocks or raises. Warns (prints) if
|flux| [m/s] exceeds the saturated hydraulic conductivity
(PERMX) of the soil, which is a red flag that the prescribed
Neumann flux is physically unrealistic (more water leaving or
entering than the soil could ever conduct). Looks first in
`self.soil_SPP["SPP_map"]` (per-zone/-layer soil property
table, set by update_soil), falling back to the scalar
`self.soil["PERMX"]` if that isn't available yet; if neither
is found, prints a note and skips the check rather than
guessing.
"""
ks = None
try:
spp = self.soil_SPP["SPP_map"]
if "PERMX" in spp.columns:
ks = float(spp["PERMX"].astype(float).max())
except Exception:
pass
if ks is None:
try:
permx = self.soil.get("PERMX")
if permx is not None:
ks = float(np.max(np.atleast_1d(permx)))
except Exception:
pass
if ks is None:
print(f" NOTE: could not determine Ks (PERMX) to sanity-check "
f"{label}={flux:.2e} m/s against - skipping check.")
return
if abs(flux) > ks:
print(f" WARNING: |{label}| = {abs(flux):.2e} m/s exceeds soil "
f"Ks (PERMX max = {ks:.2e} m/s) - this Neumann flux may "
f"be physically unrealistic. Not blocking the run.")
def assign_mesh_bc_df(self, BCtypName, times=0, **kwargs):
# step 5 add flag for each type of BC
# ---------------------------------------
# specified pressure
if "nansfdirbc" in BCtypName:
if "no_flow" in kwargs:
self.mesh_bound_cond_df[BCtypName] = -99
mask = self.mesh_bound_cond_df["noflow_bound"] == True
self.mesh_bound_cond_df.loc[mask, BCtypName] = 0
# specified flux
if "nansfneubc" in BCtypName:
if "no_flow" in kwargs:
self.mesh_bound_cond_df[BCtypName] = -99
mask = self.mesh_bound_cond_df["noflow_bound"] == True
self.mesh_bound_cond_df.loc[mask, BCtypName] = 0
# specified seepage
if "sfbc" in BCtypName:
if "no_flow" in kwargs:
self.mesh_bound_cond_df[BCtypName] = -99
mask = self.mesh_bound_cond_df["noflow_bound"] == True
self.mesh_bound_cond_df.loc[mask, BCtypName] = 0
# print(
# "SKip time dependence init boundary condition dataframe - consequences (?)"
# )
pass
def update_mesh_boundary_cond(
self,
time,
BC_name=[],
BC_val=[],
nodesId=[],
):
"""
update_mesh_bounds
Parameters
----------
bound_type : str
Neumann or Dirichlet.
bound_bool : TYPE
Boolean for bound cond.
"""
self.console.print(":sponge: [b]update boundary condition dataframe[/b]")
if 'BC_type' not in self.mesh_bound_cond_df.columns:
self.mesh_bound_cond_df.loc[
self.mesh_bound_cond_df["time"] == time,
BC_name
] = np.nan
self.mesh_bound_cond_df.loc[
(self.mesh_bound_cond_df["time"] == time) &
(self.mesh_bound_cond_df["id_node"].isin(nodesId)),
BC_name
] = BC_val
pass
def create_mesh_vtkris3d_vtk9(self):
"""
Create mesh for vtk format version 9
# vtk DataFile Version 9.0
3D Unstructured Grid of Linear Triangles
ASCII
DATASET STRUCTURED_GRID
FIELD FieldData 1
TIME 1 1 double
0.00000
POINTS 7056 float
TETRA (5,NT) - element connectivities in 3-d mesh (TETRA(5,I)
# C indicates material type for 3-d element I)
"""
if not "nnod3" in self.grid3d.keys():
self.run_processor(IPRT1=3)
with open(
os.path.join(self.workdir, self.project_name, "vtk/mesh_tmp.vtk"), "w+"
) as vtkmesh:
vtkmesh.write("# vtk DataFile Version 2.0\n")
vtkmesh.write("3D Unstructured Grid of Linear Triangles\n")
vtkmesh.write("ASCII\n")
vtkmesh.write("DATASET UNSTRUCTURED_GRID\n")
vtkmesh.write("FIELD FieldData 1\n")
vtkmesh.write("TIME 1 1 double\n")
vtkmesh.write(" 0.00000\n")
vtkmesh.write(
"POINTS " + "{:3.0f}".format(self.grid3d["nnod3"]) + " float\n"
)
np.savetxt(vtkmesh, self.grid3d["mesh3d_nodes"], fmt="%1.6e")
len(self.grid3d["mesh3d_nodes"])
ntetra = len(self.grid3d["mesh_tetra"])
numb = ntetra * 5
mesh_tretra_m = self.grid3d["mesh_tetra"].T[0:4] - 1
# np.shape(mesh_tretra_m)
tetra_mesh = np.vstack([4 * np.ones(ntetra), mesh_tretra_m]).T
vtkmesh.write(
"CELLS " + "{:d}".format(ntetra) + "\t" + "{:d}".format(numb) + "\n"
)
np.savetxt(vtkmesh, tetra_mesh, fmt="%d", delimiter="\t")
vtkmesh.write("CELL_TYPES " + "{:d}".format(ntetra) + "\n")
np.savetxt(vtkmesh, 10 * np.ones(ntetra).T, fmt="%d")
vtkmesh.close()
pass
def create_mesh_vtkris3d_vtk2(self, verbose):
"""
Create mesh for vtk format version 2
# vtk DataFile Version 2.0
3D Unstructured Grid of Linear Triangles
ASCII
DATASET STRUCTURED_GRID
FIELD FieldData 1
TIME 1 1 double
0.00000
POINTS 7056 float
TETRA (5,NT) - element connectivities in 3-d mesh (TETRA(5,I)
# C indicates material type for 3-d element I)
"""
with open(
os.path.join(self.workdir, self.project_name, "vtk/mesh_tmp.vtk"), "w+"
) as vtkmesh:
vtkmesh.write("# vtk DataFile Version 2.0\n")
vtkmesh.write("3D Unstructured Grid of Linear Triangles\n")
vtkmesh.write("ASCII\n")
vtkmesh.write("DATASET UNSTRUCTURED_GRID\n")
vtkmesh.write("FIELD FieldData 1\n")
vtkmesh.write("TIME 1 1 double\n")
vtkmesh.write(" 0.00000\n")
vtkmesh.write(
"POINTS " + "{:3.0f}".format(self.grid3d["nnod3"]) + " float\n"
)
np.savetxt(vtkmesh, self.grid3d["mesh3d_nodes"], fmt="%1.6e")
len(self.grid3d["mesh3d_nodes"])
ntetra = len(self.grid3d["mesh_tetra"])
numb = ntetra * 5
mesh_tretra_m = self.grid3d["mesh_tetra"].T[0:4] - 1
# np.shape(mesh_tretra_m)
tetra_mesh = np.vstack([4 * np.ones(ntetra), mesh_tretra_m]).T
vtkmesh.write(
"CELLS " + "{:d}".format(ntetra) + "\t" + "{:d}".format(numb) + "\n"
)
np.savetxt(vtkmesh, tetra_mesh, fmt="%d", delimiter="\t")
vtkmesh.write("CELL_TYPES " + "{:d}".format(ntetra) + "\n")
np.savetxt(vtkmesh, 10 * np.ones(ntetra).T, fmt="%d")
vtkmesh.close()
pass
def create_mesh_vtk(self, verbose=False):
"""
Create custum mesh
THIS SHOULD BE MOVED TO MESHTOOLS
"""
# if not 'nnod3' in self.grid3d.keys():
self.run_preprocessor(verbose=verbose)
self.run_processor(IPRT1=3, verbose=verbose)
self.create_mesh_vtkris3d_vtk2(verbose)
# self.create_mesh_vtkris3d_vtk9()
self.mesh_pv_attributes = pv.read(
os.path.join(self.workdir, self.project_name, "vtk/mesh_tmp.vtk")
)
self.mesh_pv_attributes.save(
os.path.join(
self.workdir, self.project_name, "vtk/", self.project_name + ".vtk"
),
binary=False,
)
# THIS IS A TEMPORARY IMPLEMENTATION OF THE PYVISTA MESH
# https://docs.pyvista.org/examples/00-load/create-structured-surface.html
pass
def update_mesh_vtk(self, prop="", prop_value=[], replaceVTK=True, **kwargs):
"""
https://docs.pyvista.org/api/core/_autosummary/pyvista.ExplicitStructuredGrid.add_field_data.html#pyvista.ExplicitStructuredGrid.add_field_data
Parameters
----------
prop : TYPE, optional
DESCRIPTION. The default is ''.
prop_value : TYPE, optional
DESCRIPTION. The default is [].
Returns
-------
None.
"""
self.mesh_pv_attributes.add_field_data(prop_value, prop)
if replaceVTK:
self.mesh_pv_attributes.save(
os.path.join(
self.workdir, self.project_name, "vtk/", self.project_name + ".vtk"
),
binary=False,
)
pass
def map_prop_2mesh_markers(self,
prop_name,
prop_map,
zones_markers_3d=None,
to_nodes=False,
# save_mesh=False,
**kwargs):
"""
Map a physical property to the CATHY mesh nodes/cells.
The mapping length should be equal to the mesh node markers
unique value length.
If no markers (i.e zones_markers_3d = None) are defined
for the mesh nodes then each layers
is associated with a unique marker
Parameters
----------
prop_name : str
property name i.e. ic, POROS, ... .
prop_map : list or pd.dataframe
Values of the property.
to_nodes : bool, optional
Map to the mesh nodes. The default is False.
Returns
-------
pv.Mesh
Updated pyvista mesh with new property.
"""
if hasattr(self, "mesh_pv_attributes") == False:
self.create_mesh_vtk()
saveMeshPath = None
if 'saveMeshPath' in kwargs:
saveMeshPath = kwargs.pop('saveMeshPath')
# mt.add_attribute_2mesh(data, mesh, kwargs)
if hasattr(self, 'zone3d'):
zones_markers_3d = self.zone3d
if zones_markers_3d is None:
zones_markers_3d = []
for l in range(self.dem_parameters['nstr']):
zones_markers_3d.append(np.ones([self.hapin["M"], self.hapin["N"]])*l)
np.unique(zones_markers_3d)
mt.add_markers_zone3d_2_mesh(
zones_markers_3d,
self.DEM,
self.mesh_pv_attributes,
self.dem_parameters,
self.hapin,
self.grid3d,
to_nodes=False,
show=False,
saveMeshPath = saveMeshPath
)
# Get unique markers in the mesh
unique_markers = np.unique(self.mesh_pv_attributes["cell_markers_zone3d"])
# Check that the number of properties matches the number of unique markers
# print('Skip Error Temporary!')
if len(prop_map) != len(unique_markers):
raise ValueError(
f"Mismatch between property map (len={len(prop_map)}) "
f"and unique mesh cell markers (len={len(unique_markers)}). "
f"Markers: {unique_markers}"
)
np.unique(self.mesh_pv_attributes["node_markers_zone3d"])
# np.unique(prop_mesh_nodes)
if to_nodes:
prop_mesh_nodes = np.zeros(len(self.mesh_pv_attributes["node_markers_zone3d"]))
for m in range(len(prop_map)):
print(m, prop_map[m])
prop_mesh_nodes[
self.mesh_pv_attributes["node_markers_zone3d"] == m
] = prop_map[m]
self.mesh_pv_attributes[f'{prop_name}_nodes'] = prop_mesh_nodes
if saveMeshPath is not None:
self.mesh_pv_attributes.save(saveMeshPath,
binary=False,
)
return list(prop_mesh_nodes)
else:
# Initialize property values for each mesh cell
prop_mesh_cells = np.zeros(len(self.mesh_pv_attributes["cell_markers_zone3d"]))
# Assign property values to mesh cells based on marker indices
for m in range(len(prop_map)):
prop_mesh_cells[self.mesh_pv_attributes["cell_markers_zone3d"] == m] = prop_map[m]
if np.any(prop_mesh_cells == 0):
raise ValueError("prop_mesh_cells contains 0 values. Check your property mapping!")
# Store cell property values in the mesh attributes
self.mesh_pv_attributes[f'{prop_name}_cells'] = prop_mesh_cells
# Set active scalars for proper visualization and conversion
self.mesh_pv_attributes.set_active_scalars(f'{prop_name}_cells')
# Convert cell data to point data (interpolation to nodes)
# prop_mesh_nodes = self.mesh_pv_attributes.cell_data_to_point_data(pass_cell_data=True)
prop_mesh_nodes = self.mesh_pv_attributes.cell_data_to_point_data()
# Assign interpolated property to mesh node attributes
self.mesh_pv_attributes[f'{prop_name}_nodes'] = prop_mesh_nodes[f'{prop_name}_cells']
if saveMeshPath is not None:
self.mesh_pv_attributes.save(saveMeshPath,
binary=False,
)
# self.mesh_pv_attributes.save("mesh_output.vtk")
self.mesh_pv_attributes['POROS_cells']
return list(prop_mesh_cells), list(self.mesh_pv_attributes[f'{prop_name}_nodes'])
def map_prop2mesh(self, dict_props):
"""
Add a given physical property to the CATHY mesh
"""
if hasattr(self, "mesh_pv_attributes") == False:
self.create_mesh_vtk()
# Guard against a stale cached mesh: if the CATHY grid has been
# (re)built since mesh_pv_attributes was created (different DEM/run),
# the node count silently drifts out of sync and any per-node
# property assignment below would be wrong. Rebuild rather than
# broadcast against the wrong length.
if "nnod3" in self.grid3d and len(self.mesh_pv_attributes.points) != self.grid3d["nnod3"]:
warnings.warn(
"mesh_pv_attributes ({} points) is out of sync with the "
"current grid3d (nnod3={}); rebuilding mesh_pv_attributes.".format(
len(self.mesh_pv_attributes.points), self.grid3d["nnod3"]
)
)
self.create_mesh_vtk()
for dp in dict_props.keys():
# FIX: `~isinstance(x, list)` is a bitwise-not on a bool, which
# is always truthy in Python (~True == -2, ~False == -1) — so
# this branch used to run unconditionally, even for per-node
# arrays (e.g. psi_ini), and broadcast them as if they were a
# single scalar value. Use an actual scalar check instead.
if np.isscalar(dict_props[dp]):
print(
"Single value detected for "
+ str(dp)
+ " ==> assumming it homogeneous"
)
self.update_mesh_vtk(
prop=dp,
prop_value=np.ones(len(self.mesh_pv_attributes.points))
* dict_props[dp],
)
else:
dict_props[dp] = np.asarray(dict_props[dp])
if len(dict_props[dp]) != len(self.mesh_pv_attributes.points):
raise ValueError(
"map_prop2mesh: '{}' has length {} but the mesh has {} "
"points — check that the property was computed on the "
"same grid as mesh_pv_attributes.".format(
dp, len(dict_props[dp]), len(self.mesh_pv_attributes.points)
)
)
self.update_mesh_vtk(prop=dp, prop_value=dict_props[dp])
pass
def map_prop_veg(self, dict_props):
if hasattr(self, "veg_map") == False:
warnings.warn("no known existing vegetation map.")
pass
else:
map_veg_dict = {}
update_map_veg = {}
for d in dict_props.keys():
map_veg = np.zeros(np.shape(self.veg_map))
for i, value in enumerate(dict_props[d]):
map_veg[self.veg_map == i+1] = value
update_map_veg[d] = map_veg
return update_map_veg
def map_prop2zone(self, dict_props, prop):
if hasattr(self, "zone") == False:
warnings.warn("no known existing zones.")
pass
else:
prop_zones = np.zeros(np.shape(self.zone))
for z in range(len(np.unique(self.zone))):
prop_zones[self.zone == z + 1] = dict_props[prop][z]
return prop_zones
# ------------------------------------------------------------------------
#%% Plot call
# ------------------------------------------------------------------------
def show(self, prop="hgsfdet", **kwargs):
"""
Call and parse to cathy.plotter from the main CATHY class
Parameters
----------
prop : str
property to plot.
**kwargs : kwargs
Returns
-------
None.
"""
try:
df = self.read_outputs(filename=prop)
except:
pass
if prop == "hgsfdet":
plt_CT.show_hgsfdet(df, **kwargs)
elif prop == "hgraph":
plt_CT.show_hgraph(df, **kwargs)
elif prop == "cumflowvol":
plt_CT.show_COCumflowvol(df, **kwargs)
elif prop == "dtcoupling":
plt_CT.show_dtcoupling(df, **kwargs)
elif prop == "wtdepth":
plt_CT.show_wtdepth(df, **kwargs)
elif prop == "spatialET":
df_fort777 = out_CT.read_fort777(os.path.join(self.workdir,
self.project_name,
'fort.777'),
)
cmap = plt_CT.show_spatialET(df_fort777, **kwargs)
return cmap
elif prop == "WTD": # water table depth
xyz = self.read_outputs('xyz')
df_psi = self.read_outputs('psi')
grid3d = self.read_outputs('grid3d')
nstr = self.dem_parameters['nstr']+1
nnod = int(grid3d['nnod'])
NPRT = np.shape(df_psi.values)[0]
XYZsurface= xyz[['x','y','z']].iloc[0:nnod].to_numpy()
WT, FLAG = self.infer_WTD_from_psi(df_psi.values,nnod,nstr,NPRT,xyz,XYZsurface)
cmap = plt_CT.plot_WTD(XYZsurface,WT,**kwargs)
return cmap
else:
print("no proxy to plot")
# elif filename == 'psi':
# plt_CT.show_psi(path)
# elif filename == 'sw':
# plt_CT.show_sw(path)
pass
def show_bc(self, BCtypName=None, time=0, ax=None, **kwargs):
"""Show bc"""
# mesh_bound_cond_df is normally populated as a side effect of
# update_nansfdirbc/update_nansfneubc/update_sfbc (see
# _record_bc_in_mesh_bounds_df). Fall back to building an empty
# skeleton here so show_bc() never raises AttributeError, e.g.
# when called on a project that was only ever loaded from disk.
if not hasattr(self, "mesh_bound_cond_df"):
if not hasattr(self, "grid3d"):
self.grid3d = self.read_outputs("grid3d")
self.create_mesh_bounds_df(
"nansfdirbc", self.grid3d["mesh3d_nodes"], times=[time]
)
if BCtypName is None:
fig = plt.figure()
# Plot nansfdirbc
# ------------------------------------
ax = fig.add_subplot(1, 3, 1, projection='3d')
cmap = plt_CT.plot_mesh_bounds('nansfdirbc',
self.mesh_bound_cond_df,
time,
ax
)
# Plot nansfneubc
# ------------------------------------
ax = fig.add_subplot(1, 3, 2, projection='3d')
cmap = plt_CT.plot_mesh_bounds('nansfneubc',
self.mesh_bound_cond_df,
time,
ax
)
# Plot sfbc
# ------------------------------------
ax = fig.add_subplot(1, 3, 3, projection='3d')
cmap = plt_CT.plot_mesh_bounds('sfbc',
self.mesh_bound_cond_df,
time,
ax
)
plt.tight_layout()
else:
cmap = plt_CT.plot_mesh_bounds(BCtypName,
self.mesh_bound_cond_df,
time,
ax
)
pass
def show_input(self, prop="hgsfdet", ax=None, **kwargs):
"""
Call and parse to cathy.plotter from the main CATHY class
Parameters
----------
prop : str
property to plot.
**kwargs : kwargs
Returns
-------
None.
"""
df = self.read_inputs(filename=prop, **kwargs)
if prop == "atmbc":
plt_CT.show_atmbc(df["time"], df["value"], ax=ax, **kwargs)
elif prop == "root_map":
plt_CT.show_indice_veg(df[0], ax=ax, **kwargs)
elif prop == "dem":
if hasattr(self, "hapin") is False:
self.update_prepo_inputs()
hapin = self.hapin
plt_CT.show_dem(df[0], hapin, ax=ax, **kwargs)
elif prop == "zone":
plt_CT.show_zone(df[0], ax=ax)
elif prop == "soil":
SPP_colname = self._get_soil_SPP_columnsNames()
FP_colname = self._get_soil_FP_columnsNames()
yprop = "PERMX"
if "yprop" in kwargs:
yprop = kwargs["yprop"]
if yprop in SPP_colname:
# in 2 dimensions
# -------------
zone_mat = in_CT.read_zone(
os.path.join(self.workdir, self.project_name, "prepro/zone")
)
layer_nb = 0
if "layer_nb" in kwargs:
layer_nb = kwargs["layer_nb"]
soil_map_prop = zone_mat[0]
exclude_zone, _ = self._check_outside_DEM(zone_mat[0])
NZONES = len(np.unique(zone_mat[0])) - exclude_zone
if NZONES-1>1:
for z in range(NZONES):
soil_map_prop[zone_mat[0] == z+1] = df[0][yprop].xs(
(z+1,layer_nb)
)
else:
soil_map_prop[zone_mat[0] == 1] = df[0][yprop].xs((1,
layer_nb)
)
cmap = plt_CT.show_soil(soil_map_prop, ax=ax,
**kwargs)
return cmap
elif yprop in FP_colname:
FP_map_prop = np.copy(self.veg_map)
for vegi in range(len(np.unique(FP_map_prop))):
FP_map_prop[FP_map_prop == vegi+1] = df[1][yprop][vegi+1]
cmap = plt_CT.show_soil(FP_map_prop, ax=ax,
**kwargs)
return cmap
elif ("dtm_" in prop) | ("lakes_map" in prop):
raster_mat, header_raster = in_CT.read_raster(
os.path.join(self.workdir, self.project_name, "prepro/" + prop)
)
if hasattr(self, "hapin") is False:
self.update_prepo_inputs()
hapin = self.hapin
plt_CT.show_raster(
raster_mat, header_raster, prop=prop, hapin=hapin, ax=ax, **kwargs
)
else:
print("no proxy to plot")
pass
#%% Read outputs/inputs
# ------------------------------------------------------------------------
def read_outputs(self, filename, **kwargs):
"""
Read CATHY format output file
Parameters
----------
filename : str
name of the output file to read.
Returns
-------
A dataframe or a dict describing file data/entries.
"""
path = os.path.join(self.workdir, self.project_name, "output", filename)
if "path" in kwargs:
path = kwargs["path"] + "/" + filename
if filename == "vp":
df = out_CT.read_vp(path)
return df
elif filename == "hgraph":
df = out_CT.read_hgraph(path)
return df
elif filename == "dtcoupling":
df = out_CT.read_dtcoupling(path)
return df
elif filename == "hgsfdet":
df = out_CT.read_hgsfdet(path)
return df
elif filename == "psi":
df = out_CT.read_psi(path)
return df
elif filename == "sw":
df = out_CT.read_sw(path)
return df
elif filename == "mbeconv":
df = out_CT.read_mbeconv(path)
return df
elif filename == "cumflowvol":
df = out_CT.read_cumflowvol(path)
return df
elif filename == "wtdepth":
df = out_CT.read_wtdepth(path)
elif filename == "xyz":
df = out_CT.read_xyz(path)
return df
elif filename == "recharge":
df = out_CT.read_recharge(path)
return df
elif filename == "grid3d":
df = out_CT.read_grid3d(path)
return df
elif filename == "ET":
path = os.path.join(self.workdir,
self.project_name,
'fort.777'
)
df = out_CT.read_fort777(path)
return df
else:
print("no file specified")
pass
def read_inputs(self, filename, **kwargs):
"""
Read CATHY format input file
Parameters
----------
filename : str
name of the input CATHY file to read.
Returns
-------
A dataframe or a dict describing file data/entries.
"""
if filename == "atmbc":
if len(self.grid3d) == 0:
# self.run_processor(IPRT1=3, DAFLAG=0)
self.grid3d = out_CT.read_grid3d(os.path.join(self.workdir,
self.project_name,
'output', 'grid3d')
)
df, HSPATM, IETO, = in_CT.read_atmbc(
os.path.join(self.workdir, self.project_name, "input", filename),
grid=self.grid3d
)
self.atmbc = {"HSPATM": HSPATM,
"IETO": IETO,
"time": df['time'].unique(),
"VALUE": df['value'],
"atmbc_df": df
}
return df
elif filename == "dem":
df = in_CT.read_dem(
os.path.join(self.workdir, self.project_name, "prepro/dem"),
os.path.join(self.workdir, self.project_name, "prepro/dtm_13.val"),
)
return df
elif filename == "root_map":
df = in_CT.read_root_map(
os.path.join(self.workdir, self.project_name, "input", filename)
)
return df
elif filename == "soil":
dem_parm = in_CT.read_dem_parameters(
os.path.join(self.workdir, self.project_name, "input", "dem_parameters")
)
if 'MAXVEG' in kwargs:
MAXVEG = kwargs['MAXVEG']
else:
MAXVEG = self.MAXVEG
# SCF-VEG: must match whatever the file was written with --
# driven by the same version-gated config used in
# _write_SOIL_file, so classic versions read the classic
# format and "SCF_variable" reads its own extra block.
cfg = VersionConfigManager("cathy", "soil_format", self.version)
scf_per_veg = cfg.get("scf_per_veg", False)
df = in_CT.read_soil(
os.path.join(self.workdir, self.project_name, "input", filename),
dem_parm,
MAXVEG=MAXVEG,
scf_per_veg=scf_per_veg,
)
# df[0]
return df
elif filename == "zone":
df = in_CT.read_zone(
os.path.join(self.workdir, self.project_name, "prepro/zone")
)
return df
elif ("dtm_" in filename) | ("lakes_map" in filename):
raster_mat, header_raster = in_CT.read_raster(
os.path.join(self.workdir, self.project_name, "prepro/" + filename)
)
return raster_mat, header_raster
else:
print("unknown file requested")
pass
# -------------------------------------------------------------------#
# %% utils
# -------------------------------------------------------------------#
def infer_WTD_from_psi(self,psi,nnod,nstr,NPRT,xyz,XYZsurface):
# Vertical profiles in rows and layers in columns,
# define topography Z
Z = np.zeros([nnod,nstr])
for l in range(nstr):
for nn in range(nnod):
Z[nn,l] = xyz['z'].iloc[l*nnod+nn]
Z = np.fliplr(Z)
FLAG = np.zeros([nnod,NPRT]);
WT=[]
for nprti in range(NPRT): # loop over NPRT
vpPSI = np.zeros([nnod,nstr])
for l in range(nstr):
for nn in range(nnod):
vpPSI[nn,l] = psi[nprti][l*nnod+nn]
vpPSI = np.fliplr(vpPSI)
# % Z0 contains for every vertical profile the height of the watertable,
Z0 = XYZsurface[:,2].copy()
for ni in range(nnod): ## loop over all mesh nodes
for nstri in range(nstr-1): ## loop over mesh layers
# print(nstri)
if vpPSI[ni, nstri] > 0 and vpPSI[ni, nstri+1] < 0 and FLAG[ni, nprti] == 0:
# Watertable, interpolate linearly (Z=rc*VPPSI+z0)
rc = (Z[ni, nstri] - Z[ni, nstri+1]) / (vpPSI[ni, nstri] - vpPSI[ni, nstri+1])
Z0[ni] = Z[ni, nstri] - rc * vpPSI[ni, nstri]
FLAG[ni, nprti] = 1
elif vpPSI[ni, nstri] > 0 and vpPSI[ni, nstri+1] < 0 and FLAG[ni, nprti] == 1:
# Second watertable
FLAG[ni, nprti] = 2
elif nstri == nstr - 2 and vpPSI[ni, nstri+1] >= 0 and FLAG[ni, nprti] == 0:
# Watertable not encountered and nodes still saturated
FLAG[ni, nprti] = 3
elif nstri == nstr -2 and FLAG[ni, nprti] == 0:
# Watertable not encountered, and not fully saturated profile. Set watertable to lowest node
Z0[ni] = Z[ni, nstri - nstr + 2]
FLAG[ni, nprti] = 4
Z0 = XYZsurface[:,2]-Z0;
WT.append(Z0)
WT = np.vstack(WT)
return WT, FLAG
def find_nearest_node(self, node_coords, grid3d=[]):
"""
Find nearest mesh node
Parameters
----------
node_coords : list
List of coordinates [[x,y,z],[x2,y2,z2]].
Returns
-------
closest_idx : list
Node indice in the grid3d.
closest : list
Node coordinate in the grid3d.
"""
if np.array(node_coords).ndim <= 1:
node_coords = [node_coords]
if len(grid3d) == 0:
grid3d = out_CT.read_grid3d(os.path.join(self.workdir,
self.project_name,
'output', 'grid3d')
)
closest_idx = []
closest = []
for i, nc in enumerate(node_coords):
# euclidean distance
if len(nc)==3:
d = (
(grid3d["mesh3d_nodes"][:, 0] - nc[0]) ** 2
+ (grid3d["mesh3d_nodes"][:, 1] - nc[1]) ** 2
# + (abs(grid3d["mesh3d_nodes"][:, 2]) - abs(nc[2])) ** 2
+ (grid3d["mesh3d_nodes"][:, 2] - nc[2]) ** 2
) ** 0.5
else:
d = (
(grid3d["mesh3d_nodes"][:, 0] - nc[0]) ** 2
+ (grid3d["mesh3d_nodes"][:, 1] - nc[1]) ** 2
) ** 0.5
closest_idx.append(np.argmin(d))
closest.append(grid3d["mesh3d_nodes"][closest_idx[i],:])
threshold = 5e-1
if d[np.argmin(d)] > threshold:
self.console.print(
":warning: [b]No node close to the required points![/b]"
)
return closest_idx, closest
def rich_display(self, title="Star Wars Movies", **kwargs):
"""
Describe the variable state and fate during the simulation with a rich table
Returns
-------
None.
"""
self.console.print(eval("self." + str(title)))
pass
def backup_simu(self):
"""
Save a copy of the simulation for reuse within python
Returns
-------
project_filename.pkl
"""
with open(
os.path.join(self.workdir, self.project_name, self.project_name + ".pkl"),
"wb",
) as f:
pickle.dump(CATHY, f)
f.close()
pass
def backup_results_DA(self, meta_DA=[]):
file_path = os.path.join(self.workdir, self.project_name, f"{self.project_name}")
try:
self.df_obs = dictObs_2pd(self.dict_obs) # Convert observations to DataFrame
# Save collumns that are not serialized are they are too large
if hasattr(self.df_obs, 'elecs'):
meta_DA['mesh_nodes_modif'] = self.df_obs['mesh_nodes_modif']
meta_DA['data_cov'] = self.df_obs['data_cov'][0]
meta_DA['elecs'] = self.df_obs['elecs'][0]
# Drop collumns that are not serialized are they are too large
self.df_obs = self.df_obs.drop(columns=["mesh_nodes_modif",'elecs','data_cov'], errors="ignore")
self.df_obs['data'] = self.df_obs['data'].apply(lambda x: x.to_json() if isinstance(x, pd.DataFrame) else x)
# attributes = ["dict_parm_pert", "df_DA", "df_obs", "df_performance", "ET_DA_xr", "df_Archie"]
attributes = ["dict_parm_pert", "df_DA", "df_obs", "ET_DA_xr", "df_Archie"]
for attr in attributes:
if hasattr(self, attr):
data = getattr(self, attr)
save_path = f"{file_path}_{attr}"
if isinstance(data, pd.DataFrame):
try:
data.to_parquet(f"{save_path}.parquet",
engine="pyarrow")
except:
with open(f"{save_path}.pkl", "wb") as f:
pickle.dump(data, f)
else:
with open(f"{save_path}.pkl", "wb") as f:
pickle.dump(data, f)
except Exception as e:
print(f"Error saving results: {e}")
def serialize(self,obj):
return json.dumps(obj.tolist() if isinstance(obj, (list, np.ndarray)) else obj)
def deserialize(self,obj):
try:
# Check if the object is a JSON string and convert back to list or np.ndarray
return json.loads(obj) if isinstance(obj, str) else obj
except:
return obj
def load_backup(self, filename="",**kwargs):
if not filename:
filename = os.path.join(self.workdir, self.project_name, f"{self.project_name}")
# List of all names (attributes) to load
all_names = [
"meta_DA",
"dict_parm_pert",
"df_DA",
"df_obs",
"df_performance",
"ET_DA_xr",
"df_Archie",
]
dict_backup = {}
for name in all_names:
if name in kwargs.keys():
if kwargs[name] == False:
continue
else:
file_base_path = f"{filename}_{name}"
# Try to find the file with either .pkl or .parquet extensions
file_path = None
for ext in [".pkl", ".parquet"]:
if os.path.exists(file_base_path + ext):
file_path = file_base_path + ext
break
if file_path is None:
print(f"File not found for {name}: {file_base_path} (with .pkl or .parquet extension)")
continue
else:
print(f"Reading File {name}: {file_base_path}")
# Load based on file extension
if file_path.endswith(".pkl"):
# Load from pickle file
try:
with open(file_path, "rb") as f:
obj = pickle.load(f)
# Deserialize DataFrame columns if needed
if isinstance(obj, pd.DataFrame):
if "data" in obj.columns:
obj["data"] = obj["data"].apply(lambda x: pd.read_json(x) if isinstance(x, str) else x)
dict_backup[name] = obj
except Exception as e:
print(f"Error loading pickle file for {name}: {e}")
elif file_path.endswith(".parquet"):
# Load from parquet file
try:
df = pd.read_parquet(file_path, engine="pyarrow")
# Deserialize specific columns if necessary
if "data_cov" in df.columns:
df["data_cov"] = df["data_cov"].apply(self.deserialize)
if "data" in df.columns:
df["data"] = df["data"].apply(lambda x: pd.read_json(x) if isinstance(x, str) else x)
dict_backup[name] = df
except Exception as e:
print(f"Error loading parquet file for {name}: {e}")
else:
print(f"Unsupported file type for {name}: {file_path}")
continue
return dict_backup
def DEPRECATED_load_pickle_backup(self,filename=""):
if len(filename) == 0:
filename = os.path.join(
self.workdir,
self.project_name,
self.project_name + "_df.pkl"
)
# filename = os.path.join(
# self.workdir, self.project_name, str(idsimu) + "_df.pkl"
# )
backup_list = []
all_names = [
"meta_DA",
"dict_parm_pert",
"df_DA",
"dict_obs",
"df_performance",
"ET_DA_xr",
"df_Archie",
]
names = []
# i = 0
# pickle_off = open(filename, "rb")
# df = pd.read_pickle(pickle_off)
with open(filename, "rb") as f:
for name in all_names:
try:
# Attempt to load the next object
obj = pickle.load(f)
backup_list.append(obj)
names.append(name)
except:
pass
# break
# Create a dictionary mapping names to loaded objects
dict_backup = {name: backup_list[i] for i, name in enumerate(names)}
return dict_backup
# def backup_results_DA(self, meta_DA=[]):
# """
# Save minimal dataframes of the simulation for result visualization within Python.
# Saves the specified attributes to a pickle file for later use.
# """
# file_path = os.path.join(
# self.workdir, self.project_name, f"{self.project_name}_df.pkl"
# )
# try:
# with open(file_path, "wb") as f:
# # Save the metadata
# pickle.dump(meta_DA, f)
# # Save attributes in a predefined order
# attributes = [
# "dict_parm_pert",
# "df_DA",
# "dict_obs",
# "df_performance",
# "ET_DA_xr",
# "df_Archie",
# ]
# for attr in attributes:
# if hasattr(self, attr):
# pickle.dump(getattr(self, attr), f)
# except Exception as e:
# print(f"An error occurred while saving results: {e}")
# raise
# def backup_results_DA(self, meta_DA=[]):
# """
# Save minimal dataframes of the simulation for result visualization within Python.
# Saves each variable in a separate pickle file.
# Returns
# -------
# None
# """
# save_dir = os.path.join(self.workdir, self.project_name)
# os.makedirs(save_dir, exist_ok=True) # Ensure the directory exists
# # Dictionary of attributes to save
# objects_to_save = {
# "meta_DA.pkl": meta_DA,
# "dict_parm_pert.pkl": self.dict_parm_pert,
# "df_DA.pkl": self.df_DA,
# "dict_obs.pkl": self.dict_obs,
# }
# # Add optional attributes dynamically
# # if hasattr(self, "df_performance"):
# # objects_to_save["df_performance.pkl"] = self.df_performance
# if hasattr(self, "ET_DA_xr"):
# objects_to_save["ET_DA_xr.pkl"] = self.ET_DA_xr
# if hasattr(self, "df_Archie"):
# objects_to_save["df_Archie.pkl"] = self.df_Archie
# # Save each object to its respective file
# for filename, obj in objects_to_save.items():
# file_path = os.path.join(save_dir, filename)
# try:
# with open(file_path, "wb") as f:
# pickle.dump(obj, f)
# print(f"Saved {filename} to {file_path}")
# except Exception as e:
# print(f"Error saving {filename}: {e}")
# def backup_results_DA(self, meta_DA=[]):
# """
# Save minimal dataframes of the simulation for result visualisation within python
# Returns
# -------
# project_filename_df.pkl
# """
# with open(
# os.path.join(
# self.workdir, self.project_name, self.project_name + "_df.pkl"
# ),
# "wb",
# ) as f:
# pickle.dump(meta_DA, f)
# pickle.dump(self.dict_parm_pert, f)
# pickle.dump(self.df_DA, f)
# pickle.dump(self.dict_obs, f)
# if hasattr(self, "df_performance"):
# pickle.dump(self.df_performance, f)
# if hasattr(self, "ET_DA_xr"):
# pickle.dump(self.ET_DA_xr, f)
# if hasattr(self, "df_Archie"):
# pickle.dump(self.df_Archie, f)
# f.close()
# def load_pickle_backup(self, filename=""):
# """
# Load a pickle backup file created by `backup_results_DA`.
# Parameters
# ----------
# filename : str, optional
# Path to the pickle file. Defaults to the project's standard backup location.
# Returns
# -------
# dict_backup : dict
# Dictionary mapping attribute names to loaded objects.
# """
# if not filename:
# filename = os.path.join(
# self.workdir, self.project_name, f"{self.project_name}_df.pkl"
# )
# all_names = [
# "meta_DA",
# "dict_parm_pert",
# "df_DA",
# "dict_obs",
# "df_performance",
# "ET_DA_xr",
# "df_Archie",
# ]
# try:
# with open(filename, "rb") as f:
# # Attempt to load objects in the same order they were saved
# backup_list = []
# for name in all_names:
# try:
# backup_list.append(pickle.load(f))
# except EOFError:
# # Reached the end of the file
# break
# # Create a dictionary of loaded objects
# return {name: obj for name, obj in zip(all_names, backup_list)}
# except (FileNotFoundError, pickle.UnpicklingError) as e:
# print(f"Error loading backup: {e}")
# return {}
def load_parquet_backup(self, filename=""):
if len(filename) == 0:
filename = os.path.join(
self.workdir,
self.project_name,
self.project_name + ".parquet"
)
backup_list = []
all_names = [
"meta_DA",
"dict_parm_pert",
"df_DA",
"df_obs",
"df_performance",
"ET_DA_xr",
"df_Archie",
]
names = []
# Open the Parquet file for each attribute (assuming each attribute is saved as a separate Parquet file)
for name in all_names:
try:
# Read the corresponding Parquet file for each attribute
file_path = filename.replace(".parquet", f"_{name}.parquet")
obj = pd.read_parquet(file_path)
backup_list.append(obj)
names.append(name)
except Exception as e:
print(f"Failed to load {name} from {file_path}: {e}")
pass
# Create a dictionary mapping names to loaded objects
dict_backup = {name: backup_list[i] for i, name in enumerate(names)}
return dict_backup