From 654ac7f9d61a68c399415c8dd5eb8abc4737c51e Mon Sep 17 00:00:00 2001 From: Paul Duda Date: Fri, 14 Mar 2025 18:12:43 -0400 Subject: [PATCH 1/5] experimenting with in-memory version of HSP2, reading directly from UCI/WDM into dictionary of pandas dataframes and running off that dictionary, using no HDF5 file --- HSP2_DriverInMem.py | 17 + HSP2_DriverInMemUCIonly.py | 36 + src/hsp2/hsp2/HSP2mainInMem.py | 846 +++++++++++++++++ src/hsp2/hsp2/HSP2utilitiesInMem.py | 67 ++ src/hsp2/hsp2tools/UciDictToObject.py | 74 ++ src/hsp2/hsp2tools/readUCIinMem.py | 1233 +++++++++++++++++++++++++ src/hsp2/hsp2tools/readWDMinMem.py | 786 ++++++++++++++++ 7 files changed, 3059 insertions(+) create mode 100644 HSP2_DriverInMem.py create mode 100644 HSP2_DriverInMemUCIonly.py create mode 100644 src/hsp2/hsp2/HSP2mainInMem.py create mode 100644 src/hsp2/hsp2/HSP2utilitiesInMem.py create mode 100644 src/hsp2/hsp2tools/UciDictToObject.py create mode 100644 src/hsp2/hsp2tools/readUCIinMem.py create mode 100644 src/hsp2/hsp2tools/readWDMinMem.py diff --git a/HSP2_DriverInMem.py b/HSP2_DriverInMem.py new file mode 100644 index 00000000..a8609717 --- /dev/null +++ b/HSP2_DriverInMem.py @@ -0,0 +1,17 @@ +# sample script to illustrate running HSP2 with inMem version + +uciName = "C://dev//HSPsquared//tests//test10//HSP2results//test10.uci" +wdmName = "C://dev//HSPsquared//tests//test10//HSP2results//test10.wdm" + +from readers.utilities.readUCIinMem import readUCI +from readers.utilities.readWDMinMem import readWDM + +uciDict = readUCI(uciName,'') +wdmDict = readWDM(wdmName,'') +# combine the 2 dicts +uciDict.update(wdmDict) + +from utilities.HSP2mainInMem import main +main(uciDict, saveall=True, jupyterlab=False) + +x = uciDict diff --git a/HSP2_DriverInMemUCIonly.py b/HSP2_DriverInMemUCIonly.py new file mode 100644 index 00000000..8b166a1b --- /dev/null +++ b/HSP2_DriverInMemUCIonly.py @@ -0,0 +1,36 @@ +# sample script to illustrate running HSP2 with inMem version + +# reads UCI and WDMs into dictionary of pandas dataframes, +# executes HSP2 off the in-memory dictionary of pandas dataframes, +# writes output to same dictionary of pandas dataframes + +uciName = "C://dev//HSPsquared//tests//test10//HSP2results//test10.uci" + +from hsp2.hsp2tools.readUCIinMem import readUCI +from hsp2.hsp2tools.readWDMinMem import readWDM +import os + +uciDict = readUCI(uciName,'') +uci_path = os.path.dirname(uciName) + +input_file_missing = False +filesDf = uciDict['/CONTROL/FILES'] +for index, row in filesDf.iterrows(): + if row['TYPE'][0:3] == 'WDM': + wdmName = row['NAME'] + wdmName = os.path.join(uci_path, wdmName) + + if os.path.isfile(wdmName): + wdmDict = readWDM(wdmName,'') + # combine the 2 dicts + uciDict.update(wdmDict) + else: + input_file_missing = True + +if input_file_missing: + print("Cant run " + uciName + ", wdm file missing") +else: + from hsp2.hsp2.HSP2mainInMem import main + main(uciDict, saveall=True, jupyterlab=False) + +x = uciDict diff --git a/src/hsp2/hsp2/HSP2mainInMem.py b/src/hsp2/hsp2/HSP2mainInMem.py new file mode 100644 index 00000000..d4dc91ac --- /dev/null +++ b/src/hsp2/hsp2/HSP2mainInMem.py @@ -0,0 +1,846 @@ +"""Copyright (c) 2020 by RESPEC, INC. +Author: Robert Heaphy, Ph.D. +License: LGPL2 + + # This table defines the expansion to INFLOW, ROFLOW, OFLOW for RCHRES networks + d = [ + ['IVOL', 'ROVOL', 'OVOL', 'HYDRFG', 'HYDR'], + ['ICON', 'ROCON', 'OCON', 'CONSFG', 'CONS'], + ['IHEAT', 'ROHEAT', 'OHEAT', 'HTFG', 'HTRCH'], + ['ISED', 'ROSED', 'OSED', 'SEDFG', 'SEDTRN'], + ['IDQAL', 'RODQAL', 'ODQAL', 'GQALFG', 'GQUAL'], + ['ISQAL', 'ROSQAL', 'OSQAL', 'GQALFG', 'GQUAL'], + ['OXIF', 'OXCF1', 'OXCF2', 'OXFG', 'OXRX'], + ['NUIF1', 'NUCF1', 'NUCF1', 'NUTFG', 'NUTRX'], + ['NUIF2', 'NUCF2', 'NUCF9', 'NUTFG', 'NUTRX'], + ['PKIF', 'PKCF1', 'PKCH2', 'PLKFG', 'PLANK'], + ['PHIF', 'PHCF1', 'PHCF2', 'PHFG', 'PHCARB']] + df = pd.DataFrame(d, columns=['INFLOW', 'ROFLOW', 'OFLOW', 'Flag', 'Name']) + df.to_hdf(h2name, '/FLOWEXPANSION', format='t', data_columns=True) +""" + +import os +from datetime import datetime as dt +from re import S +from typing import Union + +from numpy import float32, float64 +from pandas import DataFrame, date_range +from pandas.tseries.offsets import Minute + +from hsp2.hsp2.configuration import activities, expand_masslinks, noop +from hsp2.hsp2.om import om_init_state, state_load_dynamics_om, state_om_model_run_prep +from hsp2.hsp2.SPECL import specl_load_state +from hsp2.hsp2.state import ( + init_state_dicts, + state_context_hsp2, + state_init_hsp2, + state_load_dynamics_hsp2, + state_siminfo_hsp2, +) +from hsp2.hsp2.utilities import ( + expand_timeseries_names, + get_gener_timeseries, + get_timeseries, + save_timeseries, + versions, +) +from hsp2.hsp2io.hdf import HDF5 +from hsp2.hsp2io.io import Category, IOManager, SupportsReadTS + +# special in-memory version of these +from .HSP2utilitiesInMem import get_timeseries_InMem, save_timeseries_InMem +from hsp2.hsp2tools.UciDictToObject import convert_uci_InMem + + +def main( + io_manager: Union[str, IOManager], saveall: bool = False, jupyterlab: bool = True +) -> None: + """ + Run main HSP2 program. + Parameters + ---------- + io_manager + An instance of IOManager class. + saveall: bool, default=False + Saves all calculated data ignoring SAVE tables. + jupyterlab: bool, default=True + Flag for specific output behavior for jupyter lab. + + Return + ------------ + None + + """ + if isinstance(io_manager, str): + hdf5_instance = HDF5(io_manager) + io_manager = IOManager(hdf5_instance) + elif isinstance(io_manager, dict): + hdfname = '' + else: + hdfname = io_manager._input.file_path + if not os.path.exists(hdfname): + raise FileNotFoundError(f"{hdfname} HDF5 File Not Found") + + msg = messages() + msg(1, f"Processing started for file {hdfname}; saveall={saveall}") + + # read user control, parameters, states, and flags uci and map to local variables + if isinstance(io_manager, dict): + uci_obj = convert_uci_InMem(io_manager) + else: + uci_obj = io_manager.read_uci() + opseq = uci_obj.opseq + ddlinks = uci_obj.ddlinks + ddmasslinks = uci_obj.ddmasslinks + ddext_sources = uci_obj.ddext_sources + ddgener = uci_obj.ddgener + uci = uci_obj.uci + siminfo = uci_obj.siminfo + ftables = uci_obj.ftables + specactions = uci_obj.specactions + monthdata = uci_obj.monthdata + + start, stop = siminfo["start"], siminfo["stop"] + + copy_instances = {} + gener_instances = {} + + ####################################################################################### + # initialize STATE dicts + ####################################################################################### + # Set up Things in state that will be used in all modular activities like SPECL + state = init_state_dicts() + state_siminfo_hsp2(uci_obj, siminfo) + # Add support for dynamic functions to operate on STATE + if not isinstance(io_manager, dict): + # - Load any dynamic components if present, and store variables on objects + state_load_dynamics_hsp2(state, io_manager, siminfo) + else: + state['state_step_hydr'] = 'disabled' + state['hsp2_local_py'] = False + state["state_step_om"] = "disabled" + # Iterate through all segments and add crucial paths to state + # before loading dynamic components that may reference them + state_init_hsp2(state, opseq, activities) + # - finally stash specactions in state, not domain (segment) dependent so do it once + state["specactions"] = specactions # stash the specaction dict in state + om_init_state(state) # set up operational model specific state entries + specl_load_state(state, io_manager, siminfo) # traditional special actions + if not isinstance(io_manager, dict): + state_load_dynamics_om( + state, io_manager, siminfo + ) # operational model for custom python + # finalize all dynamically loaded components and prepare to run the model + state_om_model_run_prep(state, io_manager, siminfo) + ####################################################################################### + + # main processing loop + msg(1, f"Simulation Start: {start}, Stop: {stop}") + tscat = {} + for _, operation, segment, delt in opseq.itertuples(): + msg(2, f"{operation} {segment} DELT(minutes): {delt}") + siminfo["delt"] = delt + siminfo["tindex"] = date_range(start, stop, freq=Minute(delt))[1:] + siminfo["steps"] = len(siminfo["tindex"]) + + if operation == "COPY": + copy_instances[segment] = activities[operation]( + io_manager, siminfo, ddext_sources[(operation, segment)] + ) + elif operation == "GENER": + try: + if not isinstance(io_manager, dict): + ts = get_timeseries( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + else: + ts = get_timeseries_InMem( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + ts = get_gener_timeseries( + ts, gener_instances, ddlinks[segment], ddmasslinks + ) + get_flows( + io_manager, + ts, + {}, + uci, + segment, + ddlinks, + ddmasslinks, + siminfo["steps"], + msg, + ) + gener_instances[segment] = activities[operation]( + segment, + siminfo, + copy_instances, + gener_instances, + ddlinks, + ddmasslinks, + ts, + ddgener, + ) + except NotImplementedError as e: + print(f"GENER '{segment}' may not function correctly. '{e}'") + else: + # now conditionally execute all activity modules for the op, segment + if not isinstance(io_manager, dict): + ts = get_timeseries( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + else: + ts = get_timeseries_InMem( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + ts = get_gener_timeseries( + ts, gener_instances, ddlinks[segment], ddmasslinks + ) + flags = uci[(operation, "GENERAL", segment)]["ACTIVITY"] + if operation == "RCHRES": + # Add nutrient adsorption flags: + if flags["NUTRX"] == 1: + flags["TAMFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "NH3FG" + ] + flags["ADNHFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "ADNHFG" + ] + flags["PO4FG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "PO4FG" + ] + flags["ADPOFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "ADPOFG" + ] + + get_flows( + io_manager, + ts, + flags, + uci, + segment, + ddlinks, + ddmasslinks, + siminfo["steps"], + msg, + ) + + for activity, function in activities[operation].items(): + if function == noop: # or not flags[activity]: + continue + + if (activity in flags) and (not flags[activity]): + continue + + if ( + (activity == "RQUAL") + and (not flags["OXRX"]) + and (not flags["NUTRX"]) + and (not flags["PLANK"]) + and (not flags["PHCARB"]) + ): + continue + + msg(3, f"{activity}") + # Set context for dynamic executables and special actions + state_context_hsp2(state, operation, segment, activity) + + ui = uci[(operation, activity, segment)] # ui is a dictionary + if operation == "PERLND" and activity == "SEDMNT": + # special exception here to make CSNOFG available + ui["PARAMETERS"]["CSNOFG"] = uci[(operation, "PWATER", segment)][ + "PARAMETERS" + ]["CSNOFG"] + if operation == "PERLND" and activity == "PSTEMP": + # special exception here to make AIRTFG available + ui["PARAMETERS"]["AIRTFG"] = flags["ATEMP"] + if operation == "PERLND" and activity == "PWTGAS": + # special exception here to make CSNOFG available + ui["PARAMETERS"]["CSNOFG"] = uci[(operation, "PWATER", segment)][ + "PARAMETERS" + ]["CSNOFG"] + if operation == "RCHRES": + if not "PARAMETERS" in ui: + ui["PARAMETERS"] = {} + ui["PARAMETERS"]["NEXITS"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["NEXITS"] + if activity == "ADCALC": + ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] + ui["PARAMETERS"]["KS"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["KS"] + ui["PARAMETERS"]["VOL"] = uci[(operation, "HYDR", segment)][ + "STATES" + ]["VOL"] + ui["PARAMETERS"]["ROS"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["ROS"] + nexits = uci[(operation, "HYDR", segment)]["PARAMETERS"][ + "NEXITS" + ] + for index in range(nexits): + ui["PARAMETERS"]["OS" + str(index + 1)] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["OS" + str(index + 1)] + if activity == "HTRCH": + ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] + if activity == "CONS": + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + if activity == "SEDTRN": + ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] + ui["PARAMETERS"]["HTFG"] = flags["HTRCH"] + ui["PARAMETERS"]["AUX3FG"] = 0 + if flags["HYDR"]: + ui["PARAMETERS"]["LEN"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["LEN"] + ui["PARAMETERS"]["DELTH"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DELTH"] + ui["PARAMETERS"]["DB50"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DB50"] + ui["PARAMETERS"]["AUX3FG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["AUX3FG"] + if activity == "GQUAL": + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + ui["PARAMETERS"]["HTFG"] = flags["HTRCH"] + ui["PARAMETERS"]["SEDFG"] = flags["SEDTRN"] + # ui['PARAMETERS']['REAMFG'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['REAMFG'] + ui["PARAMETERS"]["HYDRFG"] = flags["HYDR"] + if flags["HYDR"]: + ui["PARAMETERS"]["LKFG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LKFG"] + ui["PARAMETERS"]["AUX1FG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["AUX1FG"] + ui["PARAMETERS"]["AUX2FG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["AUX2FG"] + ui["PARAMETERS"]["LEN"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["LEN"] + ui["PARAMETERS"]["DELTH"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DELTH"] + if flags["OXRX"]: + ui["PARAMETERS"]["LKFG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LKFG"] + ui["PARAMETERS"]["CFOREA"] = uci[ + (operation, "OXRX", segment) + ]["PARAMETERS"]["CFOREA"] + if flags["SEDTRN"]: + ui["PARAMETERS"]["SSED1"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED1"] + ui["PARAMETERS"]["SSED2"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED2"] + ui["PARAMETERS"]["SSED3"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED3"] + if flags["HTRCH"]: + ui["PARAMETERS"]["CFSAEX"] = uci[ + (operation, "HTRCH", segment) + ]["PARAMETERS"]["CFSAEX"] + elif flags["PLANK"]: + if ( + "CFSAEX" + in uci[(operation, "PLANK", segment)]["PARAMETERS"] + ): + ui["PARAMETERS"]["CFSAEX"] = uci[ + (operation, "PLANK", segment) + ]["PARAMETERS"]["CFSAEX"] + + if activity == "RQUAL": + # RQUAL inputs: + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + if flags["HYDR"]: + ui["PARAMETERS"]["LKFG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LKFG"] + + ui["FLAGS"]["HTFG"] = flags["HTRCH"] + ui["FLAGS"]["SEDFG"] = flags["SEDTRN"] + ui["FLAGS"]["GQFG"] = flags["GQUAL"] + ui["FLAGS"]["OXFG"] = flags["OXFG"] + ui["FLAGS"]["NUTFG"] = flags["NUTRX"] + ui["FLAGS"]["PLKFG"] = flags["PLANK"] + ui["FLAGS"]["PHFG"] = flags["PHCARB"] + if flags["CONS"]: + if "PARAMETERS" in uci[(operation, "CONS", segment)]: + if ( + "NCONS" + in uci[(operation, "CONS", segment)]["PARAMETERS"] + ): + ui["PARAMETERS"]["NCONS"] = uci[ + (operation, "CONS", segment) + ]["PARAMETERS"]["NCONS"] + + # OXRX module inputs: + ui_oxrx = uci[(operation, "OXRX", segment)] + + if flags["HYDR"]: + ui_oxrx["PARAMETERS"]["LEN"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LEN"] + ui_oxrx["PARAMETERS"]["DELTH"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DELTH"] + + if flags["HTRCH"]: + ui_oxrx["PARAMETERS"]["ELEV"] = uci[ + (operation, "HTRCH", segment) + ]["PARAMETERS"]["ELEV"] + + if flags["SEDTRN"]: + ui["PARAMETERS"]["SSED1"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED1"] + ui["PARAMETERS"]["SSED2"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED2"] + ui["PARAMETERS"]["SSED3"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED3"] + + # PLANK module inputs: + if flags["HTRCH"]: + ui["PARAMETERS"]["CFSAEX"] = uci[ + (operation, "HTRCH", segment) + ]["PARAMETERS"]["CFSAEX"] + + # NUTRX, PLANK, PHCARB module inputs: + ui_nutrx = uci[(operation, "NUTRX", segment)] + ui_plank = uci[(operation, "PLANK", segment)] + ui_phcarb = uci[(operation, "PHCARB", segment)] + + ############ calls activity function like snow() ############## + if operation not in ["COPY", "GENER"]: + if activity == "HYDR": + errors, errmessages = function( + io_manager, siminfo, ui, ts, ftables, state + ) + elif activity == "SEDTRN": + errors, errmessages = function( + io_manager, siminfo, ui, ts, state + ) + elif activity != "RQUAL": + errors, errmessages = function(io_manager, siminfo, ui, ts) + else: + errors, errmessages = function( + io_manager, + siminfo, + ui, + ui_oxrx, + ui_nutrx, + ui_plank, + ui_phcarb, + ts, + monthdata, + ) + ############################################################### + + for errorcnt, errormsg in zip(errors, errmessages): + if errorcnt > 0: + msg(4, f"Error count {errorcnt}: {errormsg}") + + # default to hourly output + outstep = 2 + outstep_oxrx = 2 + outstep_nutrx = 2 + outstep_plank = 2 + outstep_phcarb = 2 + if "BINOUT" in uci[(operation, "GENERAL", segment)]: + if activity in uci[(operation, "GENERAL", segment)]["BINOUT"]: + outstep = uci[(operation, "GENERAL", segment)]["BINOUT"][ + activity + ] + elif activity == "RQUAL": + outstep_oxrx = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "OXRX" + ] + outstep_nutrx = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "NUTRX" + ] + outstep_plank = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "PLANK" + ] + outstep_phcarb = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "PHCARB" + ] + + if "SAVE" in ui: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui["SAVE"], + siminfo, + saveall, + operation, + segment, + activity, + jupyterlab, + outstep, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui["SAVE"], + siminfo, + saveall, + operation, + segment, + activity, + jupyterlab, + outstep, + ) + + if activity == "RQUAL": + if "SAVE" in ui_oxrx: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_oxrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "OXRX", + jupyterlab, + outstep_oxrx, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_oxrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "OXRX", + jupyterlab, + outstep_oxrx, + ) + if "SAVE" in ui_nutrx and flags["NUTRX"] == 1: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_nutrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "NUTRX", + jupyterlab, + outstep_nutrx, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_nutrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "NUTRX", + jupyterlab, + outstep_nutrx, + ) + if "SAVE" in ui_plank and flags["PLANK"] == 1: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_plank["SAVE"], + siminfo, + saveall, + operation, + segment, + "PLANK", + jupyterlab, + outstep_plank, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_plank["SAVE"], + siminfo, + saveall, + operation, + segment, + "PLANK", + jupyterlab, + outstep_plank, + ) + if "SAVE" in ui_phcarb and flags["PHCARB"] == 1: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_phcarb["SAVE"], + siminfo, + saveall, + operation, + segment, + "PHCARB", + jupyterlab, + outstep_phcarb, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_phcarb["SAVE"], + siminfo, + saveall, + operation, + segment, + "PHCARB", + jupyterlab, + outstep_phcarb, + ) + + msglist = msg(1, "Done", final=True) + + df = DataFrame(msglist, columns=["logfile"]) + if isinstance(io_manager, dict): + io_manager["logfile"] = df + else: + io_manager.write_log(df) + if jupyterlab: + df = versions(["jupyterlab", "notebook"]) + io_manager.write_versioning(df) + print("\n\n", df) + return + + +def messages(): + """Closure routine; msg() prints messages to screen and run log""" + start = dt.now() + mlist = [] + + def msg(indent, message, final=False): + now = dt.now() + m = str(now)[:22] + " " * indent + message + if final: + mn, sc = divmod((now - start).seconds, 60) + ms = (now - start).microseconds // 100_000 + m = "; ".join((m, f"Run time is about {mn:02}:{sc:02}.{ms} (mm:ss)")) + print(m) + mlist.append(m) + return mlist + + return msg + + +def get_flows( + io_manager: SupportsReadTS, + ts, + flags, + uci, + segment, + ddlinks, + ddmasslinks, + steps, + msg, +): + # get inflows to this operation + for x in ddlinks[segment]: + if x.SVOL != "GENER": # gener already handled in get_gener_timeseries + recs = [] + if x.MLNO == "": # Data from NETWORK part of Links table + rec = {} + rec["MFACTOR"] = x.MFACTOR + rec["SGRPN"] = x.SGRPN + rec["SMEMN"] = x.SMEMN + rec["SMEMSB1"] = x.SMEMSB1 + rec["SMEMSB2"] = x.SMEMSB2 + rec["TMEMN"] = x.TMEMN + rec["TMEMSB1"] = x.TMEMSB1 + rec["TMEMSB2"] = x.TMEMSB2 + rec["SVOL"] = x.SVOL + recs.append(rec) + else: # Data from SCHEMATIC part of Links table + mldata = ddmasslinks[x.MLNO] + for dat in mldata: + if dat.SMEMN != "": + rec = {} + rec["MFACTOR"] = dat.MFACTOR + rec["SGRPN"] = dat.SGRPN + rec["SMEMN"] = dat.SMEMN + rec["SMEMSB1"] = dat.SMEMSB1 + rec["SMEMSB2"] = dat.SMEMSB2 + rec["TMEMN"] = dat.TMEMN + rec["TMEMSB1"] = dat.TMEMSB1 + rec["TMEMSB2"] = dat.TMEMSB2 + rec["SVOL"] = dat.SVOL + recs.append(rec) + else: + # this is the kind that needs to be expanded + if dat.SGRPN == "ROFLOW" or dat.SGRPN == "OFLOW": + recs = expand_masslinks(flags, uci, dat, recs) + + for rec in recs: + mfactor = rec["MFACTOR"] + sgrpn = rec["SGRPN"] + smemn = rec["SMEMN"] + smemsb1 = rec["SMEMSB1"] + smemsb2 = rec["SMEMSB2"] + tmemn = rec["TMEMN"] + tmemsb1 = rec["TMEMSB1"] + tmemsb2 = rec["TMEMSB2"] + + if x.AFACTR != "": + afactr = x.AFACTR + factor = afactr * mfactor + else: + factor = mfactor + + # KLUDGE until remaining HSP2 modules are available. + if tmemn not in { + "IVOL", + "ICON", + "IHEAT", + "ISED", + "ISED1", + "ISED2", + "ISED3", + "IDQAL", + "ISQAL1", + "ISQAL2", + "ISQAL3", + "OXIF", + "NUIF1", + "NUIF2", + "PKIF", + "PHIF", + "ONE", + "TWO", + }: + continue + if (sgrpn == "OFLOW" and smemn == "OVOL") or ( + sgrpn == "ROFLOW" and smemn == "ROVOL" + ): + sgrpn = "HYDR" + if (sgrpn == "OFLOW" and smemn == "OHEAT") or ( + sgrpn == "ROFLOW" and smemn == "ROHEAT" + ): + sgrpn = "HTRCH" + if (sgrpn == "OFLOW" and smemn == "OSED") or ( + sgrpn == "ROFLOW" and smemn == "ROSED" + ): + sgrpn = "SEDTRN" + if (sgrpn == "OFLOW" and smemn == "ODQAL") or ( + sgrpn == "ROFLOW" and smemn == "RODQAL" + ): + sgrpn = "GQUAL" + if (sgrpn == "OFLOW" and smemn == "OSQAL") or ( + sgrpn == "ROFLOW" and smemn == "ROSQAL" + ): + sgrpn = "GQUAL" + if (sgrpn == "OFLOW" and smemn == "OXCF2") or ( + sgrpn == "ROFLOW" and smemn == "OXCF1" + ): + sgrpn = "OXRX" + if ( + sgrpn == "OFLOW" + and (smemn == "NUCF9" or smemn == "OSNH4" or smemn == "OSPO4") + ) or (sgrpn == "ROFLOW" and (smemn == "NUCF1" or smemn == "NUFCF2")): + sgrpn = "NUTRX" + if (sgrpn == "OFLOW" and smemn == "PKCF2") or ( + sgrpn == "ROFLOW" and smemn == "PKCF1" + ): + sgrpn = "PLANK" + if (sgrpn == "OFLOW" and smemn == "PHCF2") or ( + sgrpn == "ROFLOW" and smemn == "PHCF1" + ): + sgrpn = "PHCARB" + + if tmemn == "ISED" or tmemn == "ISQAL": + tmemn = tmemn + tmemsb1 # need to add sand, silt, clay subscript + if (sgrpn == "HYDR" and smemn == "OVOL") or ( + sgrpn == "HTRCH" and smemn == "OHEAT" + ): + smemsb2 = "" + if sgrpn == "GQUAL" and smemsb2 == "": + smemsb2 = "1" + + smemn, tmemn = expand_timeseries_names( + sgrpn, smemn, smemsb1, smemsb2, tmemn, tmemsb1, tmemsb2 + ) + + path = f"RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}" + MFname = f"{x.SVOL}{x.SVOLNO}_MFACTOR" + AFname = f"{x.SVOL}{x.SVOLNO}_AFACTR" + data = f"{smemn}{smemsb1}{smemsb2}" + + if not isinstance(io_manager, dict): + data_frame = io_manager.read_ts( + Category.RESULTS, x.SVOL, x.SVOLNO, sgrpn + ) + else: + path = f'/RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' + data_frame = io_manager[path] + + try: + if data in data_frame.columns: + t = data_frame[data].astype(float64).to_numpy()[0:steps] + else: + t = data_frame[smemn].astype(float64).to_numpy()[0:steps] + + if MFname in ts and AFname in ts: + t *= ts[MFname][:steps] * ts[AFname][0:steps] + msg(4, f"MFACTOR modified by timeseries {MFname}") + msg(4, f"AFACTR modified by timeseries {AFname}") + elif MFname in ts: + t *= afactr * ts[MFname][0:steps] + msg(4, f"MFACTOR modified by timeseries {MFname}") + elif AFname in ts: + t *= mfactor * ts[AFname][0:steps] + msg(4, f"AFACTR modified by timeseries {AFname}") + else: + t *= factor + + # if poht to iheat, imprecision in hspf conversion factor requires a slight adjustment + if (smemn == "POHT" or smemn == "SOHT") and tmemn == "IHEAT": + t *= 0.998553 + if (smemn == "PODOXM" or smemn == "SODOXM") and tmemn == "OXIF1": + t *= 1.000565 + + # ??? ISSUE: can fetched data be at different frequency - don't know how to transform. + if tmemn in ts: + ts[tmemn] += t + else: + ts[tmemn] = t + + except KeyError: + print("ERROR in FLOWS, cant resolve ", path + " " + smemn) + + return diff --git a/src/hsp2/hsp2/HSP2utilitiesInMem.py b/src/hsp2/hsp2/HSP2utilitiesInMem.py new file mode 100644 index 00000000..fbf48b62 --- /dev/null +++ b/src/hsp2/hsp2/HSP2utilitiesInMem.py @@ -0,0 +1,67 @@ +# in memory versions of HSP2 utilities get_timeseries and save_timeseries + +from typing import Protocol, Dict, Any, List, Union, runtime_checkable +from numba import types +from numba.typed import Dict +from enum import Enum +from hsp2.hsp2.utilities import transform, clean_name +import numpy as np +import pandas as pd + +def get_timeseries_InMem(timeseries_inputs, ext_sourcesdd, siminfo): + ''' makes timeseries for the current timestep and trucated to the sim interval''' + # explicit creation of Numba dictionary with signatures + ts = Dict.empty(key_type=types.unicode_type, value_type=types.float64[:]) + for row in ext_sourcesdd: + path = f'/TIMESERIES/{row.SVOLNO}' + if path in timeseries_inputs: + data_frame = timeseries_inputs[path] + temp_data_frame = data_frame.copy() + + if row.MFACTOR != 1.0: + temp_data_frame *= row.MFACTOR + t = transform(temp_data_frame, row.TMEMN, row.TRAN, siminfo) + + tname = clean_name(row.TMEMN,row.TMEMSB) + if tname in ts: + ts[tname] += t + else: + ts[tname] = t + return ts + +def save_timeseries_InMem(timeseries, ts, savedict, siminfo, saveall, operation, segment, activity, compress=True, outstep=2): + df = pd.DataFrame(index=siminfo['tindex']) + if (operation == 'IMPLND' and activity == 'IQUAL') or (operation == 'PERLND' and activity == 'PQUAL'): + for y in savedict.keys(): + for z in set(ts.keys()): + if '/' + y in z: + zrep = z.replace('/','_') + zrep2 = zrep.replace(' ', '') + df[zrep2] = ts[z] + if '_' + y in z: + df[z] = ts[z] + elif (operation == 'RCHRES' and (activity == 'CONS' or activity == 'GQUAL')): + for y in savedict.keys(): + for z in set(ts.keys()): + if '_' + y in z: + df[z] = ts[z] + for y in (savedict.keys() & set(ts.keys())): + df[y] = ts[y] + else: + for y in (savedict.keys() & set(ts.keys())): + df[y] = ts[y] + df = df.astype(np.float32).sort_index(axis='columns') + + if saveall: + save_columns = df.columns + else: + save_columns = [key for key,value in savedict.items() if value or saveall] + + if not df.empty: + path = f'{operation}_{segment}/{activity}' + #if category: + path = '/RESULTS/' + path + timeseries[path] = df + else: + print(f'DataFrame Empty for {operation}|{activity}|{segment}') + return \ No newline at end of file diff --git a/src/hsp2/hsp2tools/UciDictToObject.py b/src/hsp2/hsp2tools/UciDictToObject.py new file mode 100644 index 00000000..b2fec5d8 --- /dev/null +++ b/src/hsp2/hsp2tools/UciDictToObject.py @@ -0,0 +1,74 @@ +''' +Convert a dict of UCI data frames into the uci object as used in HSP2 +''' + +from hsp2.hsp2.uci import UCI +import pandas as pd + +def convert_uci_InMem(uciDict): + """ + Convert a dict of UCI data frames into the uci object as used in HSP2 + Cousin of read_uci in HSP2 hdf.py, which converts HDF5 format UCI into uci object + + Parameters + ---------- + uciDict : dict of pandas dataframes as returned by readUci + + Returns + ------- + uci object as used by HSP2 + """ + + uci = UCI() + for path in uciDict.keys(): + op, module, *other = path[1:].split(sep="/", maxsplit=3) + s = "_".join(other) + if op == "CONTROL": + if module == "GLOBAL": + temp = uciDict[path].to_dict()["Info"] + uci.siminfo["start"] = pd.Timestamp(temp["Start"]) + uci.siminfo["stop"] = pd.Timestamp(temp["Stop"]) + uci.siminfo["units"] = 1 + if "Units" in temp: + if int(temp["Units"]): + uci.siminfo["units"] = int(temp["Units"]) + elif module == "LINKS": + for row in uciDict[path].fillna("").itertuples(): + if row.TVOLNO != "": + uci.ddlinks[f"{row.TVOLNO}"].append(row) + else: + uci.ddlinks[f"{row.TOPFST}"].append(row) + + elif module == "MASS_LINKS": + for row in uciDict[path].replace("na", "").itertuples(): + uci.ddmasslinks[row.MLNO].append(row) + elif module == "EXT_SOURCES": + for row in uciDict[path].replace("na", "").itertuples(): + uci.ddext_sources[(row.TVOL, row.TVOLNO)].append(row) + elif module == "OP_SEQUENCE": + uci.opseq = uciDict[path] + elif op in {"PERLND", "IMPLND", "RCHRES"}: + for id, vdict in uciDict[path].to_dict("index").items(): + uci.uci[(op, module, id)][s] = vdict + elif op == "GENER": + for row in uciDict[path].itertuples(): + if len(row) > 1: + if len(row.OPNID.split()) == 1: + start = int(row.OPNID) + stop = start + else: + start, stop = row.OPNID.split() + for i in range(int(start), int(stop) + 1): + uci.ddgener[module][f"G{i:03d}"] = row[2] + else: + pass # for cases not yet implemented + elif op == "FTABLES": + uci.ftables[module] = uciDict[path] + elif op == "SPEC_ACTIONS": + uci.specactions[module] = uciDict[path] + elif op == "MONTHDATA": + if not uci.monthdata: + uci.monthdata = {} + uci.monthdata[f"{op}/{module}"] = uciDict[path] + + return uci diff --git a/src/hsp2/hsp2tools/readUCIinMem.py b/src/hsp2/hsp2tools/readUCIinMem.py new file mode 100644 index 00000000..ef5ed438 --- /dev/null +++ b/src/hsp2/hsp2tools/readUCIinMem.py @@ -0,0 +1,1233 @@ +""" +Cousin version of readUCI.py, if no hdf5 file name passed in, +instead puts all tables into 'store' as dictionary +and returns store for access by calling routine. +""" + +import os.path +from collections import defaultdict + +import pandas as pd + +from hsp2 import hsp2tools + +pd.set_option("io.hdf.default_format", "table") + +Lapse = pd.Series( + [ + 0.0035, + 0.0035, + 0.0035, + 0.0035, + 0.0035, + 0.0035, + 0.0037, + 0.0040, + 0.0041, + 0.0043, + 0.0046, + 0.0047, + 0.0048, + 0.0049, + 0.0050, + 0.0050, + 0.0048, + 0.0046, + 0.0044, + 0.0042, + 0.0040, + 0.0038, + 0.0037, + 0.0036, + ] +) + +Seasons = pd.Series([0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0]).astype(bool) + +Svp = pd.Series( + [ + 1.005, + 1.005, + 1.005, + 1.005, + 1.005, + 1.005, + 1.005, + 1.005, + 1.005, + 1.005, + 1.01, + 1.01, + 1.015, + 1.02, + 1.03, + 1.04, + 1.06, + 1.08, + 1.1, + 1.29, + 1.66, + 2.13, + 2.74, + 3.49, + 4.40, + 5.55, + 6.87, + 8.36, + 10.1, + 12.2, + 14.6, + 17.5, + 20.9, + 24.8, + 29.3, + 34.6, + 40.7, + 47.7, + 55.7, + 64.9, + ] +) + + +def reader(filename): + # simple reader to return non blank, non comment and proper length lines + with open(filename) as uci_file: + for line in uci_file: + nline = line[:81].rstrip() # UCI max line length is 80 + if "***" in nline or not nline: + continue + yield f"{nline: <90}" # prevent short line problems + + +def getlines(f): + lines = [] + for line in f: + if line.startswith("END"): + break + lines.append(line) + return lines + + +convert = {"C": str, "I": int, "R": float} + + +def parseD(line, parse): + d = {} + for name, type_, start, end, default in parse: + field = line[start:end].strip() + d[name] = convert[type_](field) if field else convert[type_](default) + return d + + +def parseD2(line, parse, d): + icnt = 0 + for name, type_, start, end, default in parse: + field = line[start:end].strip() + icnt += 1 + # don't do anything with the first 7 values + if icnt > 8: + d[name] = convert[type_](field) if field else convert[type_](default) + return d + + +def parseD3(line, parse, d): + icnt = 0 + for name, type_, start, end, default in parse: + field = line[start:end].strip() + icnt += 1 + # don't do anything with the first 14 values + if icnt > 15: + d[name] = convert[type_](field) if field else convert[type_](default) + return d + + +def get_opnid(opnidstr, operation): + first, *last = opnidstr.split() + b = int(last[0]) if last else int(first) + a = int(first) + for x in range(a, b + 1): + yield f"{operation[0]}{x:03d}" + + +def fix_df(df, op, save, ddfaults, valid): + """fix NANs and excess ids, missing indicies, bad names""" + extra_ids = set(df.index) - valid + if extra_ids: + df = df.drop(index=extra_ids) # drop unnecessary ids + + missing_ids = valid - set(df.index) + for name1 in missing_ids: + # add missing ids with NaNs + df = pd.concat([df, pd.DataFrame(index=[name1])]) + + if df.isna().any().any(): # replace NaNs with defaults + for col in df.columns: + try: + df[col] = df[col].fillna(ddfaults[op, save, col]) + except KeyError: + pass + cols = [c.replace("(", "").replace(")", "") for c in df.columns] + df.columns = cols + df = df.apply(pd.to_numeric, errors="ignore") # todo: 'ignore' is deprecated. + return df + + +# Ignore these tables during processing, not used by HSP2 +skip = { + ("PERLND", "PRINT-INFO"), + ("PERLND", "BINARY-INFO"), + ("IMPLND", "PRINT-INFO"), + ("IMPLND", "BINARY-INFO"), + ("RCHRES", "PRINT-INFO"), + ("RCHRES", "BINARY-INFO"), +} + + +ops = {"PERLND", "IMPLND", "RCHRES", "COPY", "GENER"} +conlike = {"CONS": "NCONS", "PQUAL": "NQUAL", "IQUAL": "NQUAL", "GQUAL": "NQUAL"} + + +def readUCI(uciname, hdfname, overwrite=True): + """ + Read data from a UCI file and create a series of dataframes with the data, + if hdf5 file name passed, write to hdf5 file. + + Parameters + ---------- + uciname : str + The name of the UCI file to read. + hdfname : str + The name of the HDF file to store the data. + overwrite : bool, optional + Whether to overwrite existing data in the HDF file. Defaults to True. + + Returns + ------- + If HDF5 file, returns None. + If no HDF5 file name, returns dict of dataframes + """ + if len(hdfname) > 0: + if overwrite is True and os.path.exists(hdfname): + os.remove(hdfname) + + # create lookup dictionaries from 'ParseTable.csv' and 'rename.csv' + parse = defaultdict(list) + defaults = {} + cat = {} + path = {} + hsp_paths = {} + datapath = os.path.join(hsp2tools.__path__[0], "data", "ParseTable.csv") + for row in pd.read_csv(datapath).itertuples(): + parse[row.OP, row.TABLE].append( + (row.NAME, row.TYPE, row.START, row.STOP, row.DEFAULT) + ) + + defaults[row.OP, row.SAVE, row.NAME] = convert[row.TYPE](row.DEFAULT) + + cat[row.OP, row.TABLE] = row.CAT + path[row.OP, row.TABLE] = row.SAVE + + # store paths for checking defaults: + hsp_path = f"/{row.OP}/{row.SAVE}/{row.CAT}" + if not hsp_path in hsp_paths: + hsp_paths[hsp_path] = {} + + hsp_paths[hsp_path][row.NAME] = defaults[row.OP, row.SAVE, row.NAME] + + rename = {} + extendlen = {} + datapath = os.path.join(hsp2tools.__path__[0], "data", "rename.csv") + for row in pd.read_csv(datapath).itertuples(): + if row.LENGTH != 1: + extendlen[row.OPERATION, row.TABLE] = row.LENGTH + rename[row.OPERATION, row.TABLE] = row.RENAME + + net = None + sc = None + if len(hdfname) > 0: + store = pd.HDFStore(hdfname, mode='a') + else: + store = {} + + info = (store, parse, path, defaults, cat, rename, extendlen) + + f = reader(uciname) + for line in f: + if line.startswith("GLOBAL"): + global_(info, getlines(f)) + elif line.startswith("FILES"): + files(info, getlines(f)) + elif line.startswith("OPN"): + opn(info, getlines(f)) + elif line.startswith("NETWORK"): + net = network(info, getlines(f)) + elif line.startswith("SCHEMATIC"): + sc = schematic(info, getlines(f)) + elif line.startswith("MASS-LINK"): + masslink(info, getlines(f)) + elif line.startswith("FTABLES"): + ftables(info, getlines(f)) + elif line.startswith("EXT"): + ext(info, getlines(f)) + elif line.startswith("GENER"): + gener(info, getlines(f)) + elif line.startswith("PERLND"): + operation(info, getlines(f), "PERLND") + elif line.startswith("IMPLND"): + operation(info, getlines(f), "IMPLND") + elif line.startswith("RCHRES"): + operation(info, getlines(f), "RCHRES") + elif line.startswith("MONTH-DATA"): + monthdata(info, getlines(f)) + elif line.startswith("SPEC-ACTIONS"): + specactions(info, getlines(f)) + + colnames = ( + "AFACTR", + "MFACTOR", + "MLNO", + "SGRPN", + "SMEMN", + "SMEMSB", + "SVOL", + "SVOLNO", + "TGRPN", + "TMEMN", + "TMEMSB", + "TRAN", + "TVOL", + "TVOLNO", + "COMMENTS", + ) + if not ( (net is None) and (sc is None) ): + linkage = pd.concat((net, sc), ignore_index=True, sort=True) + for cname in colnames: + if cname not in linkage.columns: + linkage[cname] = "" + linkage = linkage.sort_values(by=["TVOLNO"]).replace("na","") + if not isinstance(store, dict): + linkage.to_hdf(store, key="/CONTROL/LINKS", data_columns=True) + else: + store["/CONTROL/LINKS"] = linkage + + if not isinstance(store, dict): + Lapse.to_hdf(store, key="TIMESERIES/LAPSE_Table") + Seasons.to_hdf(store, key="TIMESERIES/SEASONS_Table") + Svp.to_hdf(store, key="TIMESERIES/Saturated_Vapor_Pressure_Table") + else: + store["/TIMESERIES/LAPSE_Table"] = Lapse + store["/TIMESERIES/SEASONS_Table"] = Seasons + store["/TIMESERIES/Saturated_Vapor_Pressure_Table"] = Svp + + if not isinstance(store, dict): + keys = set(store.keys()) + else: + keys = store.keys() + # rename needed for restart. NOTE issue with line 157 in PERLND SNOW HSPF + # where PKSNOW = PKSNOW + PKICE at start - ONLY + path = "/PERLND/SNOW/STATES" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + df = (df.rename( + columns={"PKSNOW":"PACKF","PKICE":"PACKI","PKWATR":"PACKW"})) + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/IMPLND/SNOW/STATES" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + df = df.rename( + columns={"PKSNOW":"PACKF","PKICE":"PACKI","PKWATR":"PACKW"} + ) + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/PERLND/SNOW/FLAGS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "SNOPFG" not in df.columns: # didn't read SNOW-FLAGS table + df["SNOPFG"] = 0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/IMPLND/SNOW/FLAGS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "SNOPFG" not in df.columns: # didn't read SNOW-FLAGS table + df["SNOPFG"] = 0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + # Need to fixup missing data + path = "/IMPLND/IWATER/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "PETMIN" not in df.columns: # didn't read IWAT-PARM2 table + df["PETMIN"] = 0.35 + df["PETMAX"] = 40.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/IMPLND/IWTGAS/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "SDLFAC" not in df.columns: # didn't read LAT-FACTOR table + df["SDLFAC"] = 0.0 + df["SLIFAC"] = 0.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/IMPLND/IQUAL/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "SDLFAC" not in df.columns: # didn't read LAT-FACTOR table + df["SDLFAC"] = 0.0 + df["SLIFAC"] = 0.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/PERLND/PWTGAS/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "SDLFAC" not in df.columns: # didn't read LAT-FACTOR table + df["SDLFAC"] = 0.0 + df["SLIFAC"] = 0.0 + df["ILIFAC"] = 0.0 + df["ALIFAC"] = 0.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + if "SOTMP" not in df.columns: # didn't read PWT-TEMPS table + df["SOTMP"] = 60.0 + df["IOTMP"] = 60.0 + df["AOTMP"] = 60.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + if "SODOX" not in df.columns: # didn't read PWT-GASES table + df["SODOX"] = 0.0 + df["SOCO2"] = 0.0 + df["IODOX"] = 0.0 + df["IOCO2"] = 0.0 + df["AODOX"] = 0.0 + df["AOCO2"] = 0.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/PERLND/PWATER/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "FZG" not in df.columns: # didn't read PWAT-PARM5 table + df["FZG"] = 1.0 + df["FZGL"] = 0.1 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/PERLND/PQUAL/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "SDLFAC" not in df.columns: # didn't read LAT-FACTOR table + df["SDLFAC"] = 0.0 + df["SLIFAC"] = 0.0 + df["ILIFAC"] = 0.0 + df["ALIFAC"] = 0.0 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/RCHRES/GENERAL/INFO" + if path in keys: + if not isinstance(store, dict): + dfinfo = pd.read_hdf(store, path) + else: + dfinfo = store[path] + path = "/RCHRES/HYDR/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + df["NEXITS"] = dfinfo["NEXITS"] + df["LKFG"] = dfinfo["LKFG"] + if "IREXIT" not in df.columns: # didn't read HYDR-IRRIG table + df["IREXIT"] = 0 + df["IRMINV"] = 0.0 + df["FTBUCI"] = df["FTBUCI"].map(lambda x: f"FT{int(x):03d}") + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + del dfinfo["NEXITS"] + del dfinfo["LKFG"] + if not isinstance(store, dict): + dfinfo.to_hdf(store, key="RCHRES/GENERAL/INFO", data_columns=True) + else: + store["/RCHRES/GENERAL/INFO"] = dfinfo + + path = "/RCHRES/HTRCH/FLAGS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "BEDFLG" not in df.columns: # didn't read HT-BED-FLAGS table + df["BEDFLG"] = 0 + df["TGFLG"] = 2 + df["TSTOP"] = 55 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/RCHRES/HTRCH/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "ELEV" not in df.columns: # didn't read HEAT-PARM table + df["ELEV"] = 0.0 + df["ELDAT"] = 0.0 + df["CFSAEX"]= 1.0 + df["KATRAD"]= 9.37 + df["KCOND"] = 6.12 + df["KEVAP"] = 2.24 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/RCHRES/HTRCH/PARAMETERS" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + if "MUDDEP" not in df.columns: # didn't read HT-BED-PARM table + df["MUDDEP"]= 0.33 + df["TGRND"] = 59.0 + df["KMUD"] = 50.0 + df["KGRND"] = 1.4 + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + path = "/RCHRES/HTRCH/STATES" + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + #if 'TW' not in df.columns: # didn't read HEAT-INIT table + # df['TW'] = 60.0 + # df['AIRTMP']= 60.0 + + # apply defaults: + for path in hsp_paths: + if path in keys: + if not isinstance(store, dict): + df = pd.read_hdf(store, path) + else: + df = store[path] + dct_params = hsp_paths[path] + + new_columns = {} + for par_name in dct_params: + if par_name == "CFOREA": + ichk = 0 + + if par_name not in df.columns: # missing value in HDF5 path + def_val = dct_params[par_name] + if def_val != "None": + # df[par_name] = def_val + new_columns[par_name] = def_val + + new_columns = pd.DataFrame(new_columns, index=df.index) + df1 = pd.concat([df, new_columns], axis=1) + + if not isinstance(store, dict): + df1.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df1 + else: + if path[-6:] == "STATES": + # need to add states if it doesn't already exist to save initial state variables + # such as the case where entire IWAT-STATE1 table is being defaulted + if "df" in locals(): + for column in df.columns: # clear out existing data frame columns + df = df.drop([column], axis=1) + dct_params = hsp_paths[path] + for par_name in dct_params: + def_val = dct_params[par_name] + if def_val != "None": + df[par_name] = def_val + if not isinstance(store, dict): + df.to_hdf(store, key=path, data_columns=True) + else: + store[path] = df + + if not isinstance(store, dict): + return + else: + return store + + +def global_(info, lines): + store, parse, path, *_ = info + d = parseD(lines[1], parse["GLOBAL", "START"]) + start = str( + pd.Timestamp(f"{d['SYR']}-{d['SMO']}-{d['SDA']}") + + pd.Timedelta(int(d["SHR"]), "h") + + pd.Timedelta(int(d["SMI"]), "min") + )[0:16] + stop = str( + pd.Timestamp(f"{d['EYR']}-{d['EMO']}-{d['EDA']}") + + pd.Timedelta(int(d["EHR"]), "h") + + pd.Timedelta(int(d["EMI"]), "min") + )[0:16] + units = lines[3].strip()[56:60] + data = [lines[0].strip(), start, stop, units] + dfglobal = pd.DataFrame( + data, index=["Comment", "Start", "Stop", "Units"], columns=["Info"] + ) + if not isinstance(store, dict): + dfglobal.to_hdf(store, key="/CONTROL/GLOBAL", data_columns=True) + else: + store["/CONTROL/GLOBAL"] = dfglobal + + +def files(info, lines): + store, parse, path, *_ = info + lst = [] + for line in lines: + ftype = line[0:6].strip() + funit = line[7:12].strip() + fname = line[16:80].strip() + lst.append((ftype, funit, fname)) + dffiles = pd.DataFrame(lst, columns=["TYPE", "UNIT", "NAME"]) + if not isinstance(store, dict): + dffiles.to_hdf(store, key="/CONTROL/FILES", data_columns=True) + else: + store["/CONTROL/FILES"] = dffiles + + +def opn(info, lines): + store, parse, path, *_ = info + lst = [] + for line in lines: + tokens = line.split() + if tokens[0] == "INGRP" and tokens[1] == "INDELT": + s = tokens[2].split(":") + indelt = int(s[0]) if len(s) == 1 else 60 * int(s[0]) + int(s[1]) + elif tokens[0] in ops: + s = f"{tokens[0][0]}{int(tokens[1]):03d}" + lst.append((tokens[0], s, indelt)) + dfopn = pd.DataFrame(lst, columns=["OPERATION", "SEGMENT", "INDELT_minutes"]) + if not isinstance(store, dict): + dfopn.to_hdf(store, key="/CONTROL/OP_SEQUENCE", data_columns=True) + else: + store["/CONTROL/OP_SEQUENCE"] = dfopn + + +def network(info, lines): + store, parse, path, *_ = info + lst = [] + for line in lines: + d = parseD(line, parse["NETWORK", "na"]) + if d["SVOL"] in ops and d["TVOL"] in ops: + d["SVOLNO"] = f"{d['SVOL'][0]}{int(d['SVOLNO']):03d}" + if "TVOLNO" in d: + d["TVOLNO"] = f"{d['TVOL'][0]}{int(d['TVOLNO']):03d}" + elif "TOPFST" in d: + d["TOPFST"] = f"{d['TVOL'][0]}{int(d['TOPFST']):03d}" + lst.append(d) + return pd.DataFrame(lst, columns=d) if lst else pd.DataFrame() + + +def schematic(info, lines): + store, parse, path, *_ = info + lst = [] + for line in lines: + d = parseD(line, parse["SCHEMATIC", "na"]) + if d["SVOL"] in ops and d["TVOL"] in ops: + d["MLNO"] = f"ML{int(d['MLNO']):03d}" + d["SVOLNO"] = f"{d['SVOL'][0]}{int(d['SVOLNO']):03d}" + d["TVOLNO"] = f"{d['TVOL'][0]}{int(d['TVOLNO']):03d}" + lst.append(d) + return pd.DataFrame(lst, columns=d) if lst else pd.DataFrame() + + +def masslink(info, lines): + store, parse, path, *_ = info + lst = [] + for line in lines: + if line[2:11] == "MASS-LINK": + name = line[12:].rstrip() + elif line[2:5] != "END": + d = parseD(line, parse["MASS-LINK", "na"]) + d["MLNO"] = f"ML{int(name):03d}" + lst.append(d) + if lst: + dfmasslink = pd.DataFrame(lst, columns=d).replace("na", "") + del dfmasslink["TGRPN"] + dfmasslink["COMMENTS"] = "" + if not isinstance(store, dict): + dfmasslink.to_hdf(store, key="/CONTROL/MASS_LINKS", data_columns=True) + else: + store["/CONTROL/MASS_LINKS"] = dfmasslink + + +def ftables(info, llines): + store, parse, path, *_ = info + header = [ + "Depth", + "Area", + "Volume", + "Disch1", + "Disch2", + "Disch3", + "Disch4", + "Disch5", + ] + lines = iter(llines) + for line in lines: + if line[2:8] == "FTABLE": + unit = int(line[8:]) + name = f"FT{unit:03d}" + rows, cols = next(lines).split() + lst = [] + elif line[2:5] == "END": + dfftable = pd.DataFrame(lst, columns=header[0 : int(cols)]) + if not isinstance(store, dict): + dfftable.to_hdf(store, key=f"/FTABLES/{name}", data_columns=True) + else: + store[f"/FTABLES/{name}"] = dfftable + else: + lst.append(parseD(line, parse["FTABLES", "FTABLE"])) + + +def specactions(info, llines): + store, parse, path, *_ = info + lines = iter(llines) + # Notes: + # - Only "classic" special actions are handled here. + # - Other type of SA are recognized by the parser, but not stored in hdf5 + # - The CURLVL code is a place-holder. This has not been thought through. + # - Each type of actions "head_[action type]" should include an "CURLVL" + # column to match with conditional expression if applicable + # - The value of CURLVL matches with an expression + sa_actions = [] # referred to as "classic" in old HSPF code comments + head_actions = [ + "OPTYP", + "RANGE1", + "RANGE2", + "DC", + "DS", + "YR", + "MO", + "DA", + "HR", + "MN", + "D", + "T", + "VARI", + "S1", + "S2", + "AC", + "VALUE", + "TC", + "TS", + "NUM", + "CURLVL", + ] + sa_mult = [] + head_mult = [] + sa_uvquan = [] + head_uvquan = [] + sa_conditional = [] + head_conditional = [] + sa_distrb = [] + head_distrb = [] + sa_uvname = [] + head_uvname = [] + sa_if = [] + head_if = [] + in_if = False # are we in an if block? + curlvl = 0 + for line in lines: + if line[2:5] == "MULT": + sa_mult.append(line) + elif line[2:8] == "UVQUAN": + sa_uvquan.append(line) + elif line[2:13] == "CONDITIONAL": + sa_conditional.append(line) + elif line[2:8] == "DISTRB": + sa_distrb.append(line) + elif line[2:8] == "UVNAME": + sa_uvname.append(line) + # This CURLVL code is a place-holder. This has not been thought through. + elif line[0:2] == "IF": + sa_if.append(line) + curlvl = len(sa_if) + elif line[0:7] == "END IF": + sa_if.append(line) + curlvl = curlvl - 1 + else: + # ACTIONS block + if not '=' in line: # not complete yet + d = parseD(line, parse["SPEC-ACTIONS", "ACTIONS"]) + d["CURLVL"] = curlvl + sa_actions.append(d.copy()) + if sa_actions: + dfftable = pd.DataFrame(sa_actions, columns=head_actions).replace("na", "") + if not isinstance(store, dict): + dfftable.to_hdf(store, key=f"/SPEC_ACTIONS/ACTIONS", data_columns=True) + else: + store[f"/SPEC_ACTIONS/ACTIONS"] = dfftable + +def ext(info, lines): + store, parse, path, *_ = info + lst = [] + lst_cols = {} + for line in lines: + d = parseD(line, parse["EXT SOURCES", "na"]) + if d["TVOL"] in ops: + d["SVOLNO"] = f"TS{int(d['SVOLNO']):03d}" + # d["SVOL"] = "*" # keep the WDMn designation + if d["TGRPN"] == "EXTNL": + d["TGRPN"] = "" + toplst = int(d["TOPFST"]) if d["TOPLST"] == "na" else int(d["TOPLST"]) + for i in range(int(d["TOPFST"]), toplst + 1): + d["TVOLNO"] = f"{d['TVOL'][0]}{i:03d}" + lst.append(d.copy()) + lst_cols = d + + if lst: + dfext = pd.DataFrame(lst, columns=lst_cols).replace("na", "") + dfext["COMMENT"] = "" + del dfext["TOPFST"] + del dfext["TOPLST"] + dfext = dfext.sort_values(by=["TVOLNO"]) + if not isinstance(store, dict): + dfext.to_hdf(store, key="/CONTROL/EXT_SOURCES", data_columns=True) + else: + store["/CONTROL/EXT_SOURCES"] = dfext + + +Months = ( + "JAN", + "FEB", + "MAR", + "APR", + "MAY", + "JUN", + "JUL", + "AUG", + "SEP", + "OCT", + "NOV", + "DEC", +) + + +def operation(info, llines, op): + store, parse, dpath, ddfaults, dcat, rename, extendlen = info + counter = set() + history = defaultdict(list) + lines = iter(llines) + gcount = 0 + canCount = 0 + cancovCount = 0 + for line in lines: + tokens = line.split() + if len(tokens) == 1: + table = tokens[0] + rows = {} + if dcat[op, table] == "EXTENDED": + extended_line = 0 + for line in lines: + extended_line += 1 + if (op, table) not in parse or line[2:5] == "END": + break + if extended_line == 1: + d = parseD(line, parse[op, table]) + elif extended_line == 2: + d = parseD2(line, parse[op, table], d) + elif extended_line == 3: + d = parseD3(line, parse[op, table], d) + for opnid in get_opnid(d.pop("OPNID"), op): + rows[opnid] = d + else: + for line in lines: + if (op, table) not in parse or line[2:5] == "END": + break + d = parseD(line, parse[op, table]) + for opnid in get_opnid(d.pop("OPNID"), op): + rows[opnid] = d + df = pd.DataFrame.from_dict(rows, orient="index") + cat = dcat[op, table] + if cat.startswith("GQUAL"): + if table == "GQ-QALDATA": + gcount += 1 + cat = cat + str(gcount) + if table == "PWAT-CANOPY": + canId = int(d["CANOPY"]) + if canId > 1: + df = df.rename( + columns={ + "CANOPY": "CANOPY" + str(canId), + "CANNAM": "CANNAM" + str(canId), + "CEPSC": "CEPSC" + str(canId), + "CEPS": "CEPS" + str(canId), + } + ) + newtable = "" + if table == "MON-INTERCEP": # special cases for canopy tables + canCount += 1 + if canCount > 1: + newtable = table + str(canCount) + if table == "MON-COVER": + cancovCount += 1 + if cancovCount > 1: + newtable = table + str(cancovCount) + if len(newtable) == 0: + history[dpath[op, table], cat].append((table, df)) + else: + history[dpath[op, table], cat].append((newtable, df)) + if len(history["GENERAL", "INFO"]) > 0: + (_, df) = history["GENERAL", "INFO"][0] + valid = set(df.index) + for path, cat in history: + counter.add(path) + if cat == "SKIP": + continue + if cat in { + "PARAMETERS", + "STATES", + "FLAGS", + "ACTIVITY", + "INFO", + "BINOUT", + "CANOPY", + }: + df = pd.concat( + [temp[1] for temp in history[path, cat]], axis="columns", sort=False + ) + df = fix_df(df, op, path, ddfaults, valid) + if cat == "ACTIVITY" and op == "PERLND": + df = df.rename( + columns={ + "AIRTFG": "ATEMP", + "SNOWFG": "SNOW", + "PWATFG": "PWATER", + "SEDFG": "SEDMNT", + "PSTFG": "PSTEMP", + "PWGFG": "PWTGAS", + "PQALFG": "PQUAL", + "MSTLFG": "MSTLAY", + "PESTFG": "PEST", + "NITRFG": "NITR", + "PHOSFG": "PHOS", + "TRACFG": "TRACER", + } + ) + if cat == "ACTIVITY" and op == "IMPLND": + df = df.rename( + columns={ + "ATMPFG": "ATEMP", + "SNOWFG": "SNOW", + "IWATFG": "IWATER", + "SLDFG": "SOLIDS", + "IWGFG": "IWTGAS", + "IQALFG": "IQUAL", + } + ) + if cat == "ACTIVITY" and op == "RCHRES": + df = df.rename( + columns={ + "HYDRFG": "HYDR", + "ADFG": "ADCALC", + "CONSFG": "CONS", + "HTFG": "HTRCH", + "SEDFG": "SEDTRN", + "GQALFG": "GQUAL", + "OXFG": "OXRX", + "NUTFG": "NUTRX", + "PLKFG": "PLANK", + "PHFG": "PHCARB", + } + ) + if not isinstance(store, dict): + df.to_hdf(store, key=f"{op}/{path}/{cat}", data_columns=True) + else: + store[f"/{op}/{path}/{cat}"] = df + elif cat == "MONTHLYS": + for table, df in history[path, cat]: + df = fix_df(df, op, path, ddfaults, valid) + df.columns = Months + if table[:-1] == "MON-INTERCEP": # special cases for canopy tables + name = "CEPSC" + table[-1:] + elif table[:-1] == "MON-COVER": + name = "COVER" + table[-1:] + else: + name = rename[(op, table)] + if not isinstance(store, dict): + df.to_hdf( + store, key=f"{op}/{path}/MONTHLY/{name}", data_columns=True + ) + else: + store[f"/{op}/{path}/MONTHLY/{name}"] = df + elif cat == "EXTENDED": + temp = defaultdict(list) + for table, df in history[path, cat]: + temp[table].append(df) + for table, lst in temp.items(): + df = pd.concat(lst, axis="columns") + length = extendlen[op, table] + name = rename[op, table] + df.columns = [name + str(i) for i in range(len(df.columns))] + df = df[df.columns[0:length]] + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf( + store, key=f"{op}/{path}/EXTENDEDS/{name}", data_columns=True + ) + else: + store[f"/{op}/{path}/EXTENDEDS/{name}"] = df + elif cat == "SILTCLAY": + table, df = history[path, cat][0] + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf(store, key=f"{op}/{path}/SILT", data_columns=True) + else: + store[f"/{op}/{path}/SILT"] = df + table, df = history[path, cat][1] + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf(store, key=f"{op}/{path}/CLAY", data_columns=True) + else: + store[f"/{op}/{path}/CLAY"] = df + elif cat == "CONS": + count = 0 + for table, df in history[path, cat]: + if table == "NCONS": + temp_path = "/RCHRES/CONS/PARAMETERS" + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf(store, key=temp_path, data_columns=True) + else: + store[temp_path] = df + elif table == "CONS-DATA": + count += 1 + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf( + store, key=f"{op}/{path}/{cat}{count}", data_columns=True + ) + else: + store[f"/{op}/{path}/{cat}{count}"] = df + elif cat == "PQUAL" or cat == "IQUAL": + count = 0 + for table, df in history[path, cat]: + if table == "NQUALS": + if cat == "IQUAL": + temp_path = "/IMPLND/IQUAL/PARAMETERS" + else: + temp_path = "/PERLND/PQUAL/PARAMETERS" + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf(store, key=temp_path, data_columns=True) + else: + store[temp_path] = df + elif table.startswith("MON"): + name = rename[(op, table)] + df = fix_df(df, op, path, ddfaults, valid) + if df.size == 12: + df.columns = Months + if not isinstance(store, dict): + df.to_hdf( + store, + key=f"{op}/{path}/{cat}{count}/MONTHLY/{name}", + data_columns=True + ) + else: + store[f"/{op}/{path}/{cat}{count}/MONTHLY/{name}"] = df + else: + if table == "QUAL-PROPS": + count += 1 + tag = "FLAGS" + else: + tag = "PARAMETERS" + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf( + store, + key=f"{op}/{path}/{cat}{count}/{tag}", + data_columns=True + ) + else: + store[f"/{op}/{path}/{cat}{count}/{tag}"] = df + elif cat.startswith("GQUAL"): + count = 0 + for table, df in history[path, cat]: + if table.startswith("MON"): + name = rename[(op, table)] + df = fix_df(df, op, path, ddfaults, valid) + df.columns = Months + if not isinstance(store, dict): + df.to_hdf( + store, + key=f"{op}/{path}/{cat}/MONTHLY/{name}", + data_columns=True + ) + else: + store[f"/{op}/{path}/{cat}/MONTHLY/{name}"] = df + else: + if table == "GQ-QALDATA": + count += 1 + df = pd.concat( + [temp[1] for temp in history[path, cat]], axis="columns" + ) + df = fix_df(df, op, path, ddfaults, valid) + if not isinstance(store, dict): + df.to_hdf(store, key=f"{op}/{path}/{cat}", data_columns=True) + else: + store[f"/{op}/{path}/{cat}"] = df + else: + print("UCI TABLE is not understood (yet) by readUCI", op, cat) + + savetable = defaultdict(dict) + datapath = os.path.join(hsp2tools.__path__[0], "data", "SaveTable.csv") + for row in pd.read_csv(datapath).itertuples(): + savetable[row.OPERATION, row.ACTIVITY][row.NAME] = row.VALUE + for activity in counter - {"GENERAL"}: + df = pd.DataFrame(index=sorted(valid)) + for name, value in savetable[op, activity].items(): + df[name] = int(value) + if not isinstance(store, dict): + df.to_hdf(store, key=f"{op}/{activity}/SAVE", data_columns=True) + else: + if df.columns.size > 0: + store[f"/{op}/{activity}/SAVE"] = df + + +def copy(info, lines): + # placeholder - PRT - no sure I actually need to implement this method + pass + + +def gener(info, lines): + store, parse, path, *_ = info + lst = [] + sub_blocks = ["OPCODE", "PARM", "NTERMS", "COEFFS"] + current_block = "" + d = {} + for line in lines: + if line[2:5] == "END": + df = pd.DataFrame(lst, columns=d) + if not isinstance(store, dict): + df.to_hdf(store, key=f"GENER/{current_block}", data_columns=True) + else: + store[f"/GENER/{current_block}"] = df + lst.clear() + elif any(s in line for s in sub_blocks): + current_block = [s for s in sub_blocks if s in line][0] + else: + d = parseD(line, parse["GENER", current_block]) + lst.append(d) + + +def monthdata(info, llines): + store, parse, path, *_ = info + header = [ + "JAN", + "FEB", + "MAR", + "APR", + "MAY", + "JUN", + "JUL", + "AUG", + "SEP", + "OCT", + "NOV", + "DEC", + ] + lines = iter(llines) + for line in lines: + if line[2:12] == "MONTH-DATA": + unit = line[12:].strip() + name = "MONTHDATA" + unit + lst = [] + elif line[2:5] == "END": + dfftable = pd.DataFrame(lst, columns=header[0:12]) + if not isinstance(store, dict): + dfftable.to_hdf(store, f"/MONTHDATA/{name}", data_columns=True) + else: + store[f"/MONTHDATA/{name}"] = dfftable + else: + vals = [] + line = line.strip() + while line: + try: + vals.append(float(line[:6])) + except: + if '-' in line[:6]: + # assume this exception was caused by the E missing, as is allowed in HSPF + negpos = line[:6].index("-") + vals.append(float(line[:negpos] + 'E' + line[negpos:6])) + line = line[6:] + lst.append(vals) diff --git a/src/hsp2/hsp2tools/readWDMinMem.py b/src/hsp2/hsp2tools/readWDMinMem.py new file mode 100644 index 00000000..780b8d53 --- /dev/null +++ b/src/hsp2/hsp2tools/readWDMinMem.py @@ -0,0 +1,786 @@ +""" +Cousin version of readWDM.py, if no hdf5 file name passed in, +instead puts all tables into 'store' as dictionary +and returns store for access by calling routine. +""" + +import datetime +import warnings + +import numpy as np +import pandas as pd +from numba import jit, njit + +# look up attributes NAME, data type (Integer; Real; String) and data length by attribute number +attrinfo = { + 1: ("TSTYPE", "S", 4), + 2: ("STAID", "S", 16), + 11: ("DAREA", "R", 1), + 17: ("TCODE", "I", 1), + 27: ("TSBYR", "I", 1), + 28: ("TSBMO", "I", 1), + 29: ("TSBDY", "I", 1), + 30: ("TSBHR", "I", 1), + 32: ("TFILL", "R", 1), + 33: ("TSSTEP", "I", 1), + 34: ("TGROUP", "I", 1), + 45: ("STNAM", "S", 48), + 83: ("COMPFG", "I", 1), + 84: ("TSFORM", "I", 1), + 85: ("VBTIME", "I", 1), + 444: ("DATMOD", "S", 12), + 443: ("DATCRE", "S", 12), + 22: ("DCODE", "I", 1), + 10: ("DESCRP", "S", 80), + 7: ("ELEV", "R", 1), + 8: ("LATDEG", "R", 1), + 9: ("LNGDEG", "R", 1), + 288: ("SCENARIO", "S", 8), + 289: ("CONSTITUENT", "S", 8), + 290: ("LOCATION", "S", 8), +} + +freq = { + 7: "100YS", + 6: "YS", + 5: "MS", + 4: "D", + 3: "h", + 2: "min", + 1: "S", +} # pandas date_range() frequency by TCODE, TGROUP + + +def readWDM(wdmfile, hdffile, compress_output=False): + iarray = np.fromfile(wdmfile, dtype=np.int32) + farray = np.fromfile(wdmfile, dtype=np.float32) + + date_epoch = np.datetime64(0, "Y") + dt_year = np.timedelta64(1, "Y") + dt_month = np.timedelta64(1, "M") + dt_day = np.timedelta64(1, "D") + dt_hour = np.timedelta64(1, "h") + dt_minute = np.timedelta64(1, "m") + dt_second = np.timedelta64(1, "s") + + if iarray[0] != -998: + raise ValueError( + "Provided file does not match WDM format. First int32 should be -998." + ) + nrecords = iarray[28] # first record is File Definition Record + ntimeseries = iarray[31] + + dsnlist = [] + for index in range(512, nrecords * 512, 512): + if ( + not ( + iarray[index] == 0 + and iarray[index + 1] == 0 + and iarray[index + 2] == 0 + and iarray[index + 3] == 0 + ) + and iarray[index + 5] == 1 + ): + dsnlist.append(index) + if len(dsnlist) != ntimeseries: + raise RuntimeError( + f"Wrong number of Time Series Records found expecting:{ntimeseries} found:{len(dsnlist)}" + ) + + if len(hdffile) > 0: + store = pd.HDFStore(hdffile) + else: + store = {} + + summary = [] + summaryindx = [] + + # check to see which extra attributes are on each dsn + columns_to_add = [] + search = ["STAID", "STNAM", "SCENARIO", "CONSTITUENT", "LOCATION"] + for att in search: + found_in_all = True + for index in dsnlist: + dattr = {} + psa = iarray[index + 9] + if psa > 0: + sacnt = iarray[index + psa - 1] + for i in range(psa + 1, psa + 1 + 2 * sacnt, 2): + id = iarray[index + i] + ptr = iarray[index + i + 1] - 1 + index + if id not in attrinfo: + continue + name, atype, length = attrinfo[id] + if atype == "I": + dattr[name] = iarray[ptr] + elif atype == "R": + dattr[name] = farray[ptr] + else: + dattr[name] = "".join( + [ + _inttostr(iarray[k]) + for k in range(ptr, ptr + length // 4) + ] + ).strip() + if att not in dattr: + found_in_all = False + if found_in_all: + columns_to_add.append(att) + + for index in dsnlist: + # get layout information for TimeSeries Dataset frame + dsn = iarray[index + 4] + psa = iarray[index + 9] + if psa > 0: + sacnt = iarray[index + psa - 1] + pdat = iarray[index + 10] + pdatv = iarray[index + 11] + frepos = iarray[index + pdat] + + print(f"{dsn} reading from wdm") + # get attributes + dattr = { + "TSBDY": 1, + "TSBHR": 1, + "TSBMO": 1, + "TSBYR": 1900, + "TFILL": -999.0, + } # preset defaults + for i in range(psa + 1, psa + 1 + 2 * sacnt, 2): + id = iarray[index + i] + ptr = iarray[index + i + 1] - 1 + index + if id not in attrinfo: + # print('PROGRAM ERROR: ATTRIBUTE INDEX not found', id, 'Attribute pointer', iarray[index + i+1]) + continue + + name, atype, length = attrinfo[id] + if atype == "I": + dattr[name] = iarray[ptr] + elif atype == "R": + dattr[name] = farray[ptr] + else: + dattr[name] = "".join( + [_inttostr(iarray[k]) for k in range(ptr, ptr + length // 4)] + ).strip() + + # Get timeseries timebase data + records = [] + offsets = [] + for i in range(pdat + 1, pdatv - 1): + a = iarray[index + i] + if a != 0: + record, offset = _splitposition(a) + records.append(record) + offsets.append(offset) + if len(records) == 0: + continue + + # calculate number of data points in each group, tindex is final index for storage + tgroup = dattr["TGROUP"] + tstep = dattr["TSSTEP"] + tcode = dattr["TCODE"] + + records = np.asarray(records) + offsets = np.asarray(offsets) + + dates, values, stop_datetime = _process_groups( + iarray, farray, records, offsets, tgroup + ) + stop_datetime = datetime.datetime(*_bits_to_date(stop_datetime)) + dates = np.array(dates) + dates_converted = _date_convert( + dates, + date_epoch, + dt_year, + dt_month, + dt_day, + dt_hour, + dt_minute, + dt_second, + ) + series = pd.Series(values, index=dates_converted) + try: + series.index.freq = str(tstep) + freq[tcode] + except ValueError: + series.index.freq = None + + dsname = f"/TIMESERIES/TS{dsn:03d}" + + if not isinstance(store, dict): + if compress_output: + series.to_hdf(store, key=dsname, complib="blosc", complevel=9) + else: + series.to_hdf(store, key=dsname, format="t", data_columns=True) + else: + store[dsname] = series + + data = [ + str(series.index[0]), + str(stop_datetime), + str(tstep) + freq[tcode], + len(series), + dattr["TSTYPE"], + dattr["TFILL"], + ] + columns = ["Start", "Stop", "Freq", "Length", "TSTYPE", "TFILL"] + for x in columns_to_add: + if x in dattr: + data.append(dattr[x]) + columns.append(x) + + summary.append(data) + summaryindx.append(dsname[11:]) + + if not isinstance(store, dict): + dfsummary = pd.DataFrame(summary, index=summaryindx, columns=columns) + store.put("/TIMESERIES/SUMMARY", dfsummary, format="t", data_columns=True) + return dfsummary + else: + return store + +@njit +def _splitdate(x): + year = np.int64(x >> 14) + month = np.int64(x >> 10 & 0xF) + day = np.int64(x >> 5 & 0x1F) + hour = np.int64(x & 0x1F) + return _correct_date(year, month, day, hour, 0, 0) + + +@njit +def _splitcontrol(x): + nval = x >> 16 + ltstep = x >> 10 & 0x3F + ltcode = x >> 7 & 0x7 + comp = x >> 5 & 0x3 + qual = x & 0x1F + return nval, ltstep, ltcode, comp, qual + + +@njit +def _splitposition(x): + return ((x >> 9) - 1, (x & 0x1FF) - 1) # args: record, offset + + +@njit +def _inttostr(i): + return ( + chr(i & 0xFF) + chr(i >> 8 & 0xFF) + chr(i >> 16 & 0xFF) + chr(i >> 24 & 0xFF) + ) + + +@njit +def _bits_to_date(x): + year = x >> 26 + month = x >> 22 & 0xF + day = x >> 17 & 0x1F + hour = x >> 12 & 0x1F + minute = x >> 6 & 0x3F + second = x & 0x3F + return year, month, day, hour, minute, second + + +@njit +def _date_to_bits(year, month, day, hour, minute, second): + x = year << 26 | month << 22 | day << 17 | hour << 12 | minute << 6 | second + return x + + +@njit +def _increment_date(date, timecode, timestep): + year, month, day, hour, minute, second = _bits_to_date(date) + + if timecode == 7: + year += 100 * timestep + elif timecode == 6: + year += timestep + elif timecode == 5: + month += timestep + elif timecode == 4: + day += timestep + elif timecode == 3: + hour += timestep + elif timecode == 2: + minute += timestep + elif timecode == 1: + second += timestep + + return _correct_date(year, month, day, hour, minute, second) + + +@njit +def _correct_date(year, month, day, hour, minute, second): + while second >= 60: + second -= 60 + minute += 1 + while minute >= 60: + minute -= 60 + hour += 1 + while hour >= 24: + hour -= 24 + day += 1 + while day > _days_in_month(year, month): + day -= _days_in_month(year, month) + month += 1 + while month > 12: + month -= 12 + year += 1 + return _date_to_bits(year, month, day, hour, minute, second) + + +@njit +def _days_in_month(year, month): + if month > 12: + month %= 12 + + if month in (1, 3, 5, 7, 8, 10, 12): + return 31 + elif month in (4, 6, 9, 11): + return 30 + elif month == 2: + if _is_leapyear(year): + return 29 + else: + return 28 + + +@njit +def _is_leapyear(year): + if year % 400 == 0: + return True + if year % 100 == 0: + return False + if year % 4 == 0: + return True + else: + return False + + +@njit +def _date_convert( + dates, date_epoch, dt_year, dt_month, dt_day, dt_hour, dt_minute, dt_second +): + converted_dates = [] + for x in dates: + year, month, day, hour, minute, second = _bits_to_date(x) + date = date_epoch + date += (year - 1970) * dt_year + date += (month - 1) * dt_month + date += (day - 1) * dt_day + date += hour * dt_hour + date += minute * dt_minute + date += second * dt_second + converted_dates.append(date) + return converted_dates + + +@njit +def _process_groups(iarray, farray, records, offsets, tgroup): + date_array = [0] # need initialize with a type for numba + value_array = [0.0] + + for i in range(0, len(records)): + record = records[i] + offset = offsets[i] + index = record * 512 + offset + pscfwr = iarray[record * 512 + 3] # should be 0 for last record in timeseries + current_date = _splitdate(iarray[index]) + group_enddate = _increment_date(current_date, tgroup, 1) + offset += 1 + index += 1 + + while current_date < group_enddate: + nval, ltstep, ltcode, comp, qual = _splitcontrol(iarray[index]) + # compressed - only has single value which applies to full range + if comp == 1: + for i in range(0, nval, 1): + date_array.append(current_date) + current_date = _increment_date(current_date, ltcode, ltstep) + value_array.append(farray[index + 1]) + index += 2 + offset += 2 + else: + for i in range(0, nval, 1): + date_array.append(current_date) + current_date = _increment_date(current_date, ltcode, ltstep) + value_array.append(farray[index + 1 + i]) + index += 1 + nval + offset += 1 + nval + + if offset >= 511: + offset = 4 + index = (pscfwr - 1) * 512 + offset + record = pscfwr + pscfwr = iarray[ + (record - 1) * 512 + 3 + ] # should be 0 for last record in timeseries + + date_array = date_array[1:] + value_array = value_array[1:] + + return date_array, value_array, group_enddate + + +def get_wdm_data_set(wdmfile, attributes): + """ + Get single time series data from a WDM file + based on a collection of attributes (name-value pairs) + """ + if attributes == None: + return None + + search_loc = attributes["location"] + search_cons = attributes["constituent"] + search_dsn = attributes["dsn"] + + iarray = np.fromfile(wdmfile, dtype=np.int32) + farray = np.fromfile(wdmfile, dtype=np.float32) + + if iarray[0] != -998: + print("Not a WDM file, magic number is not -990. Stopping!") + return None + nrecords = iarray[28] # first record is File Definition Record + ntimeseries = iarray[31] + + dsnlist = [] + for index in range(512, nrecords * 512, 512): + if ( + not ( + iarray[index] == 0 + and iarray[index + 1] == 0 + and iarray[index + 2] == 0 + and iarray[index + 3] == 0 + ) + and iarray[index + 5] == 1 + ): + dsnlist.append(index) + if len(dsnlist) != ntimeseries: + print("PROGRAM ERROR, wrong number of DSN records found") + + summary = [] + summaryindx = [] + + # check to see which extra attributes are on each dsn + columns_to_add = [] + search = ["STAID", "STNAM", "SCENARIO", "CONSTITUENT", "LOCATION"] + """ + for att in search: + found_in_all = True + for index in dsnlist: + dattr = {} + psa = iarray[index + 9] + if psa > 0: + sacnt = iarray[index + psa - 1] + for i in range(psa + 1, psa + 1 + 2 * sacnt, 2): + id = iarray[index + i] + ptr = iarray[index + i + 1] - 1 + index + if id not in attrinfo: + continue + name, atype, length = attrinfo[id] + if atype == 'I': + dattr[name] = iarray[ptr] + elif atype == 'R': + dattr[name] = farray[ptr] + else: + dattr[name] = ''.join([itostr(iarray[k]) for k in range(ptr, ptr + length // 4)]).strip() + if att not in dattr: + found_in_all = False + if found_in_all: + columns_to_add.append(att) + """ + + for index in dsnlist: + # get layout information for TimeSeries Dataset frame + dsn = iarray[index + 4] + psa = iarray[index + 9] + if psa > 0: + sacnt = iarray[index + psa - 1] + pdat = iarray[index + 10] + pdatv = iarray[index + 11] + frepos = iarray[index + pdat] + + print(f"{dsn} reading from wdm") + + # get attributes + dattr = { + "TSBDY": 1, + "TSBHR": 1, + "TSBMO": 1, + "TSBYR": 1900, + "TFILL": -999.0, + } # preset defaults + for i in range(psa + 1, psa + 1 + 2 * sacnt, 2): + id = iarray[index + i] + ptr = iarray[index + i + 1] - 1 + index + if id not in attrinfo: + # print('PROGRAM ERROR: ATTRIBUTE INDEX not found', id, 'Attribute pointer', iarray[index + i+1]) + continue + + name, atype, length = attrinfo[id] + if atype == "I": + dattr[name] = iarray[ptr] + elif atype == "R": + dattr[name] = farray[ptr] + else: + dattr[name] = "".join( + [itostr(iarray[k]) for k in range(ptr, ptr + length // 4)] + ).strip() + + if search_dsn > 0 and search_dsn == dsn: + pass + else: + # could do more attribute based filtering here such as constituent, location etc + if search_cons == dattr["TSTYPE"]: + pass + else: + continue + + # Get timeseries timebase data + records = [] + for i in range(pdat + 1, pdatv - 1): + a = iarray[index + i] + if a != 0: + records.append(splitposition(a)) + if len(records) == 0: + continue # WDM preallocation, but nothing saved here yet + + srec, soffset = records[0] + start = splitdate(iarray[srec * 512 + soffset]) + + sprec, spoffset = splitposition(frepos) + finalindex = sprec * 512 + spoffset + + # calculate number of data points in each group, tindex is final index for storage + tgroup = dattr["TGROUP"] + tstep = dattr["TSSTEP"] + tcode = dattr["TCODE"] + cindex = pd.date_range(start=start, periods=len(records) + 1, freq=freq[tgroup]) + tindex = pd.date_range( + start=start, end=cindex[-1], freq=str(tstep) + freq[tcode] + ) + counts = np.diff(np.searchsorted(tindex, cindex)) + + ## Get timeseries data + floats = np.zeros(sum(counts), dtype=np.float32) + findex = 0 + for (rec, offset), count in zip(records, counts): + findex = getfloats( + iarray, + farray, + floats, + findex, + rec, + offset, + count, + finalindex, + tcode, + tstep, + ) + + ts = pd.Series(floats[:findex], index=tindex[:findex]) + df = pd.DataFrame({"ts": ts}) + return df + + return None + + +######################## +### legacy functions ### +######################## + + +def todatetime(yr=1900, mo=1, dy=1, hr=0): + """takes yr,mo,dy,hr information then returns its datetime64""" + warnings.warn( + "use '_date_convert' instead; Removed for numba compatible datetime handling; reference commit 1b52a1736e45a497ccdf78cd6e2eab8c0b7a444f", + DeprecationWarning, + ) + if hr == 24: + return datetime.datetime(yr, mo, dy, 23) + pd.Timedelta(1, "h") + else: + return datetime.datetime(yr, mo, dy, hr) + + +def splitdate(x): + """splits WDM int32 DATWRD into year, month, day, hour -> then returns its datetime64""" + warnings.warn( + "use '_splitdate' instead; naming updated to indicate internal function", + DeprecationWarning, + ) + return todatetime( + x >> 14, x >> 10 & 0xF, x >> 5 & 0x1F, x & 0x1F + ) # args: year, month, day, hour + + +def splitcontrol(x): + """splits int32 into (qual, compcode, units, tstep, nvalues)""" + warnings.warn( + "use '_splitcontrol' instead; naming updated to indicate internal function", + DeprecationWarning, + ) + return (x & 0x1F, x >> 5 & 0x3, x >> 7 & 0x7, x >> 10 & 0x3F, x >> 16) + + +def splitposition(x): + """splits int32 into (record, offset), converting to Pyton zero based indexing""" + warnings.warn( + "use '_spiltposition' instead; naming updated to indicate internal function", + DeprecationWarning, + ) + return ((x >> 9) - 1, (x & 0x1FF) - 1) + + +def itostr(i): + warnings.warn( + "use '_inttostr' instead; naming updated to indicate internal function; method also handles integer argments so updated name from 'i' to 'int' for additonal clarity", + DeprecationWarning, + ) + return ( + chr(i & 0xFF) + chr(i >> 8 & 0xFF) + chr(i >> 16 & 0xFF) + chr(i >> 24 & 0xFF) + ) + + +def getfloats( + iarray, farray, floats, findex, rec, offset, count, finalindex, tcode, tstep +): + warnings.warn( + "discontinue use and replace with new 'process_groups' instead; Function replaced by incompatible group/block processing approach; reference commit c5b2a1cdd6a55eccc0db67d7840ec3eaf904dcec .", + DeprecationWarning, + ) + index = rec * 512 + offset + 1 + stop = (rec + 1) * 512 + cntr = 0 + while cntr < count and findex < len(floats): + if index == stop - 1: + print( + "Problem?", str(rec) + ) # perhaps not, block cannot start at word 512 of a record because not spot for values + if index >= stop - 1: + rec = ( + iarray[rec * 512 + 3] - 1 + ) # 3 is forward data pointer, -1 is python indexing + print("Process record ", str(rec)) + index = rec * 512 + 4 # 4 is index of start of new data + stop = (rec + 1) * 512 + + x = iarray[index] # block control word or maybe date word at start of group + nval = x >> 16 + ltstep = x >> 10 & 0x3F + ltcode = x >> 7 & 0x7 + comp = x >> 5 & 0x3 + qual = x & 0x1F + ldate = todatetime() # dummy + if ltstep != tstep or ltcode != tcode: + nval = adjustNval(ldate, ltstep, tstep, ltcode, tcode, comp, nval) + if nval == -1: # unable to convert block + try: + ldate = splitdate(x) + if isinstance(ldate, datetime.date): + print("Problem resolved - date found ", ldate) + nval = 1 + else: + print( + "BlockConversionFailure at ", + str(rec + 1), + " ", + str(index % 512), + ) + except: + print( + "BlockConversionFailure at ", + str(rec + 1), + " ", + str(index % 512), + ) + # try next word + comp = -1 + index += 1 + if comp == 0: + for k in range(nval): + if findex >= len(floats): + return findex + floats[findex] = farray[index + k] + findex += 1 + index += nval + elif comp > 0: + for k in range(nval): + if findex >= len(floats): + return findex + floats[findex] = farray[index] + findex += 1 + index += 1 + cntr += nval + return findex + + +def adjustNval(ldate, ltstep, tstep, ltcode, tcode, comp, nval): + warnings.warn( + "supporting function for deprecated 'get_floats' function;", DeprecationWarning + ) + lnval = nval + if comp != 1: + nval = -1 # only can adjust compressed data + else: + if tcode == 2: # minutes + if ltcode == 6: # from years + if leap_year(ldate.year): + ldays = 366 + else: + ldays = 365 + nval = ldays * lnval * 60 / ltstep + elif ltcode == 5: # from months + from calendar import monthrange + + ldateRange = monthrange(ldate.year, ldate.month) + print("month block ", ldateRange) + nval = ldateRange[2] * lnval * 60 / ltstep + elif ltcode == 4: # from days + nval = lnval * 60 / ltstep + elif ltcode == 3: # from hours + nval = lnval * 1440 / ltstep + else: # dont know how to convert + nval = -1 + elif tcode == 3: # hours + if ltcode == 6: # from years + nval = -1 + elif ltcode == 5: # from months + nval = -1 + elif ltcode == 4: # from days + nval = lnval * 24 / ltstep + else: # dont know how to convert + nval = -1 + else: + nval = -1 # dont know how to convert + + nval = int(nval) + if nval == -1: # conversion problem + print( + "Conversion problem (tcode ", + str(tcode), + ", ", + str(ltcode), + "), (tstep ", + str(tstep), + ",", + str(ltstep), + "), (comp ", + str(comp), + ")", + ) + else: + print( + "Conversion complete (tcode ", + str(tcode), + ", ", + str(ltcode), + "), (tstep ", + str(tstep), + ",", + str(ltstep), + "), (nval ", + str(nval) + ",", + str(lnval), + ")", + ) + + return nval From fff224141dd60f38ac877e5771f7589ca47d53b0 Mon Sep 17 00:00:00 2001 From: Paul Duda Date: Wed, 26 Mar 2025 10:19:10 -0400 Subject: [PATCH 2/5] HSP2mainInMem.py -- add try catch for inactive sections --- src/hsp2/hsp2/HSP2mainInMem.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/hsp2/hsp2/HSP2mainInMem.py b/src/hsp2/hsp2/HSP2mainInMem.py index d4dc91ac..5564d2b4 100644 --- a/src/hsp2/hsp2/HSP2mainInMem.py +++ b/src/hsp2/hsp2/HSP2mainInMem.py @@ -24,6 +24,7 @@ from re import S from typing import Union +import pandas as pd from numpy import float32, float64 from pandas import DataFrame, date_range from pandas.tseries.offsets import Minute @@ -807,7 +808,10 @@ def get_flows( ) else: path = f'/RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' - data_frame = io_manager[path] + try: + data_frame = io_manager[path] + except KeyError: + data_frame = pd.DataFrame() try: if data in data_frame.columns: From 234ea1a4dc65ba16a0512f5f2ec5ec34da2fde5d Mon Sep 17 00:00:00 2001 From: Burgholzer Date: Fri, 5 Sep 2025 08:40:25 -0400 Subject: [PATCH 3/5] replace main.py for better diff --- src/hsp2/hsp2/HSP2mainInMem.py | 850 ---------------------------- src/hsp2/hsp2/main.py | 980 +++++++++++++++++++++++---------- 2 files changed, 687 insertions(+), 1143 deletions(-) delete mode 100644 src/hsp2/hsp2/HSP2mainInMem.py diff --git a/src/hsp2/hsp2/HSP2mainInMem.py b/src/hsp2/hsp2/HSP2mainInMem.py deleted file mode 100644 index 5564d2b4..00000000 --- a/src/hsp2/hsp2/HSP2mainInMem.py +++ /dev/null @@ -1,850 +0,0 @@ -"""Copyright (c) 2020 by RESPEC, INC. -Author: Robert Heaphy, Ph.D. -License: LGPL2 - - # This table defines the expansion to INFLOW, ROFLOW, OFLOW for RCHRES networks - d = [ - ['IVOL', 'ROVOL', 'OVOL', 'HYDRFG', 'HYDR'], - ['ICON', 'ROCON', 'OCON', 'CONSFG', 'CONS'], - ['IHEAT', 'ROHEAT', 'OHEAT', 'HTFG', 'HTRCH'], - ['ISED', 'ROSED', 'OSED', 'SEDFG', 'SEDTRN'], - ['IDQAL', 'RODQAL', 'ODQAL', 'GQALFG', 'GQUAL'], - ['ISQAL', 'ROSQAL', 'OSQAL', 'GQALFG', 'GQUAL'], - ['OXIF', 'OXCF1', 'OXCF2', 'OXFG', 'OXRX'], - ['NUIF1', 'NUCF1', 'NUCF1', 'NUTFG', 'NUTRX'], - ['NUIF2', 'NUCF2', 'NUCF9', 'NUTFG', 'NUTRX'], - ['PKIF', 'PKCF1', 'PKCH2', 'PLKFG', 'PLANK'], - ['PHIF', 'PHCF1', 'PHCF2', 'PHFG', 'PHCARB']] - df = pd.DataFrame(d, columns=['INFLOW', 'ROFLOW', 'OFLOW', 'Flag', 'Name']) - df.to_hdf(h2name, '/FLOWEXPANSION', format='t', data_columns=True) -""" - -import os -from datetime import datetime as dt -from re import S -from typing import Union - -import pandas as pd -from numpy import float32, float64 -from pandas import DataFrame, date_range -from pandas.tseries.offsets import Minute - -from hsp2.hsp2.configuration import activities, expand_masslinks, noop -from hsp2.hsp2.om import om_init_state, state_load_dynamics_om, state_om_model_run_prep -from hsp2.hsp2.SPECL import specl_load_state -from hsp2.hsp2.state import ( - init_state_dicts, - state_context_hsp2, - state_init_hsp2, - state_load_dynamics_hsp2, - state_siminfo_hsp2, -) -from hsp2.hsp2.utilities import ( - expand_timeseries_names, - get_gener_timeseries, - get_timeseries, - save_timeseries, - versions, -) -from hsp2.hsp2io.hdf import HDF5 -from hsp2.hsp2io.io import Category, IOManager, SupportsReadTS - -# special in-memory version of these -from .HSP2utilitiesInMem import get_timeseries_InMem, save_timeseries_InMem -from hsp2.hsp2tools.UciDictToObject import convert_uci_InMem - - -def main( - io_manager: Union[str, IOManager], saveall: bool = False, jupyterlab: bool = True -) -> None: - """ - Run main HSP2 program. - Parameters - ---------- - io_manager - An instance of IOManager class. - saveall: bool, default=False - Saves all calculated data ignoring SAVE tables. - jupyterlab: bool, default=True - Flag for specific output behavior for jupyter lab. - - Return - ------------ - None - - """ - if isinstance(io_manager, str): - hdf5_instance = HDF5(io_manager) - io_manager = IOManager(hdf5_instance) - elif isinstance(io_manager, dict): - hdfname = '' - else: - hdfname = io_manager._input.file_path - if not os.path.exists(hdfname): - raise FileNotFoundError(f"{hdfname} HDF5 File Not Found") - - msg = messages() - msg(1, f"Processing started for file {hdfname}; saveall={saveall}") - - # read user control, parameters, states, and flags uci and map to local variables - if isinstance(io_manager, dict): - uci_obj = convert_uci_InMem(io_manager) - else: - uci_obj = io_manager.read_uci() - opseq = uci_obj.opseq - ddlinks = uci_obj.ddlinks - ddmasslinks = uci_obj.ddmasslinks - ddext_sources = uci_obj.ddext_sources - ddgener = uci_obj.ddgener - uci = uci_obj.uci - siminfo = uci_obj.siminfo - ftables = uci_obj.ftables - specactions = uci_obj.specactions - monthdata = uci_obj.monthdata - - start, stop = siminfo["start"], siminfo["stop"] - - copy_instances = {} - gener_instances = {} - - ####################################################################################### - # initialize STATE dicts - ####################################################################################### - # Set up Things in state that will be used in all modular activities like SPECL - state = init_state_dicts() - state_siminfo_hsp2(uci_obj, siminfo) - # Add support for dynamic functions to operate on STATE - if not isinstance(io_manager, dict): - # - Load any dynamic components if present, and store variables on objects - state_load_dynamics_hsp2(state, io_manager, siminfo) - else: - state['state_step_hydr'] = 'disabled' - state['hsp2_local_py'] = False - state["state_step_om"] = "disabled" - # Iterate through all segments and add crucial paths to state - # before loading dynamic components that may reference them - state_init_hsp2(state, opseq, activities) - # - finally stash specactions in state, not domain (segment) dependent so do it once - state["specactions"] = specactions # stash the specaction dict in state - om_init_state(state) # set up operational model specific state entries - specl_load_state(state, io_manager, siminfo) # traditional special actions - if not isinstance(io_manager, dict): - state_load_dynamics_om( - state, io_manager, siminfo - ) # operational model for custom python - # finalize all dynamically loaded components and prepare to run the model - state_om_model_run_prep(state, io_manager, siminfo) - ####################################################################################### - - # main processing loop - msg(1, f"Simulation Start: {start}, Stop: {stop}") - tscat = {} - for _, operation, segment, delt in opseq.itertuples(): - msg(2, f"{operation} {segment} DELT(minutes): {delt}") - siminfo["delt"] = delt - siminfo["tindex"] = date_range(start, stop, freq=Minute(delt))[1:] - siminfo["steps"] = len(siminfo["tindex"]) - - if operation == "COPY": - copy_instances[segment] = activities[operation]( - io_manager, siminfo, ddext_sources[(operation, segment)] - ) - elif operation == "GENER": - try: - if not isinstance(io_manager, dict): - ts = get_timeseries( - io_manager, ddext_sources[(operation, segment)], siminfo - ) - else: - ts = get_timeseries_InMem( - io_manager, ddext_sources[(operation, segment)], siminfo - ) - ts = get_gener_timeseries( - ts, gener_instances, ddlinks[segment], ddmasslinks - ) - get_flows( - io_manager, - ts, - {}, - uci, - segment, - ddlinks, - ddmasslinks, - siminfo["steps"], - msg, - ) - gener_instances[segment] = activities[operation]( - segment, - siminfo, - copy_instances, - gener_instances, - ddlinks, - ddmasslinks, - ts, - ddgener, - ) - except NotImplementedError as e: - print(f"GENER '{segment}' may not function correctly. '{e}'") - else: - # now conditionally execute all activity modules for the op, segment - if not isinstance(io_manager, dict): - ts = get_timeseries( - io_manager, ddext_sources[(operation, segment)], siminfo - ) - else: - ts = get_timeseries_InMem( - io_manager, ddext_sources[(operation, segment)], siminfo - ) - ts = get_gener_timeseries( - ts, gener_instances, ddlinks[segment], ddmasslinks - ) - flags = uci[(operation, "GENERAL", segment)]["ACTIVITY"] - if operation == "RCHRES": - # Add nutrient adsorption flags: - if flags["NUTRX"] == 1: - flags["TAMFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ - "NH3FG" - ] - flags["ADNHFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ - "ADNHFG" - ] - flags["PO4FG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ - "PO4FG" - ] - flags["ADPOFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ - "ADPOFG" - ] - - get_flows( - io_manager, - ts, - flags, - uci, - segment, - ddlinks, - ddmasslinks, - siminfo["steps"], - msg, - ) - - for activity, function in activities[operation].items(): - if function == noop: # or not flags[activity]: - continue - - if (activity in flags) and (not flags[activity]): - continue - - if ( - (activity == "RQUAL") - and (not flags["OXRX"]) - and (not flags["NUTRX"]) - and (not flags["PLANK"]) - and (not flags["PHCARB"]) - ): - continue - - msg(3, f"{activity}") - # Set context for dynamic executables and special actions - state_context_hsp2(state, operation, segment, activity) - - ui = uci[(operation, activity, segment)] # ui is a dictionary - if operation == "PERLND" and activity == "SEDMNT": - # special exception here to make CSNOFG available - ui["PARAMETERS"]["CSNOFG"] = uci[(operation, "PWATER", segment)][ - "PARAMETERS" - ]["CSNOFG"] - if operation == "PERLND" and activity == "PSTEMP": - # special exception here to make AIRTFG available - ui["PARAMETERS"]["AIRTFG"] = flags["ATEMP"] - if operation == "PERLND" and activity == "PWTGAS": - # special exception here to make CSNOFG available - ui["PARAMETERS"]["CSNOFG"] = uci[(operation, "PWATER", segment)][ - "PARAMETERS" - ]["CSNOFG"] - if operation == "RCHRES": - if not "PARAMETERS" in ui: - ui["PARAMETERS"] = {} - ui["PARAMETERS"]["NEXITS"] = uci[(operation, "HYDR", segment)][ - "PARAMETERS" - ]["NEXITS"] - if activity == "ADCALC": - ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] - ui["PARAMETERS"]["KS"] = uci[(operation, "HYDR", segment)][ - "PARAMETERS" - ]["KS"] - ui["PARAMETERS"]["VOL"] = uci[(operation, "HYDR", segment)][ - "STATES" - ]["VOL"] - ui["PARAMETERS"]["ROS"] = uci[(operation, "HYDR", segment)][ - "PARAMETERS" - ]["ROS"] - nexits = uci[(operation, "HYDR", segment)]["PARAMETERS"][ - "NEXITS" - ] - for index in range(nexits): - ui["PARAMETERS"]["OS" + str(index + 1)] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["OS" + str(index + 1)] - if activity == "HTRCH": - ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] - ui["advectData"] = uci[(operation, "ADCALC", segment)][ - "adcalcData" - ] - # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - if activity == "CONS": - ui["advectData"] = uci[(operation, "ADCALC", segment)][ - "adcalcData" - ] - if activity == "SEDTRN": - ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] - ui["advectData"] = uci[(operation, "ADCALC", segment)][ - "adcalcData" - ] - # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - ui["PARAMETERS"]["HTFG"] = flags["HTRCH"] - ui["PARAMETERS"]["AUX3FG"] = 0 - if flags["HYDR"]: - ui["PARAMETERS"]["LEN"] = uci[(operation, "HYDR", segment)][ - "PARAMETERS" - ]["LEN"] - ui["PARAMETERS"]["DELTH"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["DELTH"] - ui["PARAMETERS"]["DB50"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["DB50"] - ui["PARAMETERS"]["AUX3FG"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["AUX3FG"] - if activity == "GQUAL": - ui["advectData"] = uci[(operation, "ADCALC", segment)][ - "adcalcData" - ] - ui["PARAMETERS"]["HTFG"] = flags["HTRCH"] - ui["PARAMETERS"]["SEDFG"] = flags["SEDTRN"] - # ui['PARAMETERS']['REAMFG'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['REAMFG'] - ui["PARAMETERS"]["HYDRFG"] = flags["HYDR"] - if flags["HYDR"]: - ui["PARAMETERS"]["LKFG"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["LKFG"] - ui["PARAMETERS"]["AUX1FG"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["AUX1FG"] - ui["PARAMETERS"]["AUX2FG"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["AUX2FG"] - ui["PARAMETERS"]["LEN"] = uci[(operation, "HYDR", segment)][ - "PARAMETERS" - ]["LEN"] - ui["PARAMETERS"]["DELTH"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["DELTH"] - if flags["OXRX"]: - ui["PARAMETERS"]["LKFG"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["LKFG"] - ui["PARAMETERS"]["CFOREA"] = uci[ - (operation, "OXRX", segment) - ]["PARAMETERS"]["CFOREA"] - if flags["SEDTRN"]: - ui["PARAMETERS"]["SSED1"] = uci[ - (operation, "SEDTRN", segment) - ]["STATES"]["SSED1"] - ui["PARAMETERS"]["SSED2"] = uci[ - (operation, "SEDTRN", segment) - ]["STATES"]["SSED2"] - ui["PARAMETERS"]["SSED3"] = uci[ - (operation, "SEDTRN", segment) - ]["STATES"]["SSED3"] - if flags["HTRCH"]: - ui["PARAMETERS"]["CFSAEX"] = uci[ - (operation, "HTRCH", segment) - ]["PARAMETERS"]["CFSAEX"] - elif flags["PLANK"]: - if ( - "CFSAEX" - in uci[(operation, "PLANK", segment)]["PARAMETERS"] - ): - ui["PARAMETERS"]["CFSAEX"] = uci[ - (operation, "PLANK", segment) - ]["PARAMETERS"]["CFSAEX"] - - if activity == "RQUAL": - # RQUAL inputs: - ui["advectData"] = uci[(operation, "ADCALC", segment)][ - "adcalcData" - ] - if flags["HYDR"]: - ui["PARAMETERS"]["LKFG"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["LKFG"] - - ui["FLAGS"]["HTFG"] = flags["HTRCH"] - ui["FLAGS"]["SEDFG"] = flags["SEDTRN"] - ui["FLAGS"]["GQFG"] = flags["GQUAL"] - ui["FLAGS"]["OXFG"] = flags["OXFG"] - ui["FLAGS"]["NUTFG"] = flags["NUTRX"] - ui["FLAGS"]["PLKFG"] = flags["PLANK"] - ui["FLAGS"]["PHFG"] = flags["PHCARB"] - if flags["CONS"]: - if "PARAMETERS" in uci[(operation, "CONS", segment)]: - if ( - "NCONS" - in uci[(operation, "CONS", segment)]["PARAMETERS"] - ): - ui["PARAMETERS"]["NCONS"] = uci[ - (operation, "CONS", segment) - ]["PARAMETERS"]["NCONS"] - - # OXRX module inputs: - ui_oxrx = uci[(operation, "OXRX", segment)] - - if flags["HYDR"]: - ui_oxrx["PARAMETERS"]["LEN"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["LEN"] - ui_oxrx["PARAMETERS"]["DELTH"] = uci[ - (operation, "HYDR", segment) - ]["PARAMETERS"]["DELTH"] - - if flags["HTRCH"]: - ui_oxrx["PARAMETERS"]["ELEV"] = uci[ - (operation, "HTRCH", segment) - ]["PARAMETERS"]["ELEV"] - - if flags["SEDTRN"]: - ui["PARAMETERS"]["SSED1"] = uci[ - (operation, "SEDTRN", segment) - ]["STATES"]["SSED1"] - ui["PARAMETERS"]["SSED2"] = uci[ - (operation, "SEDTRN", segment) - ]["STATES"]["SSED2"] - ui["PARAMETERS"]["SSED3"] = uci[ - (operation, "SEDTRN", segment) - ]["STATES"]["SSED3"] - - # PLANK module inputs: - if flags["HTRCH"]: - ui["PARAMETERS"]["CFSAEX"] = uci[ - (operation, "HTRCH", segment) - ]["PARAMETERS"]["CFSAEX"] - - # NUTRX, PLANK, PHCARB module inputs: - ui_nutrx = uci[(operation, "NUTRX", segment)] - ui_plank = uci[(operation, "PLANK", segment)] - ui_phcarb = uci[(operation, "PHCARB", segment)] - - ############ calls activity function like snow() ############## - if operation not in ["COPY", "GENER"]: - if activity == "HYDR": - errors, errmessages = function( - io_manager, siminfo, ui, ts, ftables, state - ) - elif activity == "SEDTRN": - errors, errmessages = function( - io_manager, siminfo, ui, ts, state - ) - elif activity != "RQUAL": - errors, errmessages = function(io_manager, siminfo, ui, ts) - else: - errors, errmessages = function( - io_manager, - siminfo, - ui, - ui_oxrx, - ui_nutrx, - ui_plank, - ui_phcarb, - ts, - monthdata, - ) - ############################################################### - - for errorcnt, errormsg in zip(errors, errmessages): - if errorcnt > 0: - msg(4, f"Error count {errorcnt}: {errormsg}") - - # default to hourly output - outstep = 2 - outstep_oxrx = 2 - outstep_nutrx = 2 - outstep_plank = 2 - outstep_phcarb = 2 - if "BINOUT" in uci[(operation, "GENERAL", segment)]: - if activity in uci[(operation, "GENERAL", segment)]["BINOUT"]: - outstep = uci[(operation, "GENERAL", segment)]["BINOUT"][ - activity - ] - elif activity == "RQUAL": - outstep_oxrx = uci[(operation, "GENERAL", segment)]["BINOUT"][ - "OXRX" - ] - outstep_nutrx = uci[(operation, "GENERAL", segment)]["BINOUT"][ - "NUTRX" - ] - outstep_plank = uci[(operation, "GENERAL", segment)]["BINOUT"][ - "PLANK" - ] - outstep_phcarb = uci[(operation, "GENERAL", segment)]["BINOUT"][ - "PHCARB" - ] - - if "SAVE" in ui: - if not isinstance(io_manager, dict): - save_timeseries( - io_manager, - ts, - ui["SAVE"], - siminfo, - saveall, - operation, - segment, - activity, - jupyterlab, - outstep, - ) - else: - save_timeseries_InMem( - io_manager, - ts, - ui["SAVE"], - siminfo, - saveall, - operation, - segment, - activity, - jupyterlab, - outstep, - ) - - if activity == "RQUAL": - if "SAVE" in ui_oxrx: - if not isinstance(io_manager, dict): - save_timeseries( - io_manager, - ts, - ui_oxrx["SAVE"], - siminfo, - saveall, - operation, - segment, - "OXRX", - jupyterlab, - outstep_oxrx, - ) - else: - save_timeseries_InMem( - io_manager, - ts, - ui_oxrx["SAVE"], - siminfo, - saveall, - operation, - segment, - "OXRX", - jupyterlab, - outstep_oxrx, - ) - if "SAVE" in ui_nutrx and flags["NUTRX"] == 1: - if not isinstance(io_manager, dict): - save_timeseries( - io_manager, - ts, - ui_nutrx["SAVE"], - siminfo, - saveall, - operation, - segment, - "NUTRX", - jupyterlab, - outstep_nutrx, - ) - else: - save_timeseries_InMem( - io_manager, - ts, - ui_nutrx["SAVE"], - siminfo, - saveall, - operation, - segment, - "NUTRX", - jupyterlab, - outstep_nutrx, - ) - if "SAVE" in ui_plank and flags["PLANK"] == 1: - if not isinstance(io_manager, dict): - save_timeseries( - io_manager, - ts, - ui_plank["SAVE"], - siminfo, - saveall, - operation, - segment, - "PLANK", - jupyterlab, - outstep_plank, - ) - else: - save_timeseries_InMem( - io_manager, - ts, - ui_plank["SAVE"], - siminfo, - saveall, - operation, - segment, - "PLANK", - jupyterlab, - outstep_plank, - ) - if "SAVE" in ui_phcarb and flags["PHCARB"] == 1: - if not isinstance(io_manager, dict): - save_timeseries( - io_manager, - ts, - ui_phcarb["SAVE"], - siminfo, - saveall, - operation, - segment, - "PHCARB", - jupyterlab, - outstep_phcarb, - ) - else: - save_timeseries_InMem( - io_manager, - ts, - ui_phcarb["SAVE"], - siminfo, - saveall, - operation, - segment, - "PHCARB", - jupyterlab, - outstep_phcarb, - ) - - msglist = msg(1, "Done", final=True) - - df = DataFrame(msglist, columns=["logfile"]) - if isinstance(io_manager, dict): - io_manager["logfile"] = df - else: - io_manager.write_log(df) - if jupyterlab: - df = versions(["jupyterlab", "notebook"]) - io_manager.write_versioning(df) - print("\n\n", df) - return - - -def messages(): - """Closure routine; msg() prints messages to screen and run log""" - start = dt.now() - mlist = [] - - def msg(indent, message, final=False): - now = dt.now() - m = str(now)[:22] + " " * indent + message - if final: - mn, sc = divmod((now - start).seconds, 60) - ms = (now - start).microseconds // 100_000 - m = "; ".join((m, f"Run time is about {mn:02}:{sc:02}.{ms} (mm:ss)")) - print(m) - mlist.append(m) - return mlist - - return msg - - -def get_flows( - io_manager: SupportsReadTS, - ts, - flags, - uci, - segment, - ddlinks, - ddmasslinks, - steps, - msg, -): - # get inflows to this operation - for x in ddlinks[segment]: - if x.SVOL != "GENER": # gener already handled in get_gener_timeseries - recs = [] - if x.MLNO == "": # Data from NETWORK part of Links table - rec = {} - rec["MFACTOR"] = x.MFACTOR - rec["SGRPN"] = x.SGRPN - rec["SMEMN"] = x.SMEMN - rec["SMEMSB1"] = x.SMEMSB1 - rec["SMEMSB2"] = x.SMEMSB2 - rec["TMEMN"] = x.TMEMN - rec["TMEMSB1"] = x.TMEMSB1 - rec["TMEMSB2"] = x.TMEMSB2 - rec["SVOL"] = x.SVOL - recs.append(rec) - else: # Data from SCHEMATIC part of Links table - mldata = ddmasslinks[x.MLNO] - for dat in mldata: - if dat.SMEMN != "": - rec = {} - rec["MFACTOR"] = dat.MFACTOR - rec["SGRPN"] = dat.SGRPN - rec["SMEMN"] = dat.SMEMN - rec["SMEMSB1"] = dat.SMEMSB1 - rec["SMEMSB2"] = dat.SMEMSB2 - rec["TMEMN"] = dat.TMEMN - rec["TMEMSB1"] = dat.TMEMSB1 - rec["TMEMSB2"] = dat.TMEMSB2 - rec["SVOL"] = dat.SVOL - recs.append(rec) - else: - # this is the kind that needs to be expanded - if dat.SGRPN == "ROFLOW" or dat.SGRPN == "OFLOW": - recs = expand_masslinks(flags, uci, dat, recs) - - for rec in recs: - mfactor = rec["MFACTOR"] - sgrpn = rec["SGRPN"] - smemn = rec["SMEMN"] - smemsb1 = rec["SMEMSB1"] - smemsb2 = rec["SMEMSB2"] - tmemn = rec["TMEMN"] - tmemsb1 = rec["TMEMSB1"] - tmemsb2 = rec["TMEMSB2"] - - if x.AFACTR != "": - afactr = x.AFACTR - factor = afactr * mfactor - else: - factor = mfactor - - # KLUDGE until remaining HSP2 modules are available. - if tmemn not in { - "IVOL", - "ICON", - "IHEAT", - "ISED", - "ISED1", - "ISED2", - "ISED3", - "IDQAL", - "ISQAL1", - "ISQAL2", - "ISQAL3", - "OXIF", - "NUIF1", - "NUIF2", - "PKIF", - "PHIF", - "ONE", - "TWO", - }: - continue - if (sgrpn == "OFLOW" and smemn == "OVOL") or ( - sgrpn == "ROFLOW" and smemn == "ROVOL" - ): - sgrpn = "HYDR" - if (sgrpn == "OFLOW" and smemn == "OHEAT") or ( - sgrpn == "ROFLOW" and smemn == "ROHEAT" - ): - sgrpn = "HTRCH" - if (sgrpn == "OFLOW" and smemn == "OSED") or ( - sgrpn == "ROFLOW" and smemn == "ROSED" - ): - sgrpn = "SEDTRN" - if (sgrpn == "OFLOW" and smemn == "ODQAL") or ( - sgrpn == "ROFLOW" and smemn == "RODQAL" - ): - sgrpn = "GQUAL" - if (sgrpn == "OFLOW" and smemn == "OSQAL") or ( - sgrpn == "ROFLOW" and smemn == "ROSQAL" - ): - sgrpn = "GQUAL" - if (sgrpn == "OFLOW" and smemn == "OXCF2") or ( - sgrpn == "ROFLOW" and smemn == "OXCF1" - ): - sgrpn = "OXRX" - if ( - sgrpn == "OFLOW" - and (smemn == "NUCF9" or smemn == "OSNH4" or smemn == "OSPO4") - ) or (sgrpn == "ROFLOW" and (smemn == "NUCF1" or smemn == "NUFCF2")): - sgrpn = "NUTRX" - if (sgrpn == "OFLOW" and smemn == "PKCF2") or ( - sgrpn == "ROFLOW" and smemn == "PKCF1" - ): - sgrpn = "PLANK" - if (sgrpn == "OFLOW" and smemn == "PHCF2") or ( - sgrpn == "ROFLOW" and smemn == "PHCF1" - ): - sgrpn = "PHCARB" - - if tmemn == "ISED" or tmemn == "ISQAL": - tmemn = tmemn + tmemsb1 # need to add sand, silt, clay subscript - if (sgrpn == "HYDR" and smemn == "OVOL") or ( - sgrpn == "HTRCH" and smemn == "OHEAT" - ): - smemsb2 = "" - if sgrpn == "GQUAL" and smemsb2 == "": - smemsb2 = "1" - - smemn, tmemn = expand_timeseries_names( - sgrpn, smemn, smemsb1, smemsb2, tmemn, tmemsb1, tmemsb2 - ) - - path = f"RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}" - MFname = f"{x.SVOL}{x.SVOLNO}_MFACTOR" - AFname = f"{x.SVOL}{x.SVOLNO}_AFACTR" - data = f"{smemn}{smemsb1}{smemsb2}" - - if not isinstance(io_manager, dict): - data_frame = io_manager.read_ts( - Category.RESULTS, x.SVOL, x.SVOLNO, sgrpn - ) - else: - path = f'/RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' - try: - data_frame = io_manager[path] - except KeyError: - data_frame = pd.DataFrame() - - try: - if data in data_frame.columns: - t = data_frame[data].astype(float64).to_numpy()[0:steps] - else: - t = data_frame[smemn].astype(float64).to_numpy()[0:steps] - - if MFname in ts and AFname in ts: - t *= ts[MFname][:steps] * ts[AFname][0:steps] - msg(4, f"MFACTOR modified by timeseries {MFname}") - msg(4, f"AFACTR modified by timeseries {AFname}") - elif MFname in ts: - t *= afactr * ts[MFname][0:steps] - msg(4, f"MFACTOR modified by timeseries {MFname}") - elif AFname in ts: - t *= mfactor * ts[AFname][0:steps] - msg(4, f"AFACTR modified by timeseries {AFname}") - else: - t *= factor - - # if poht to iheat, imprecision in hspf conversion factor requires a slight adjustment - if (smemn == "POHT" or smemn == "SOHT") and tmemn == "IHEAT": - t *= 0.998553 - if (smemn == "PODOXM" or smemn == "SODOXM") and tmemn == "OXIF1": - t *= 1.000565 - - # ??? ISSUE: can fetched data be at different frequency - don't know how to transform. - if tmemn in ts: - ts[tmemn] += t - else: - ts[tmemn] = t - - except KeyError: - print("ERROR in FLOWS, cant resolve ", path + " " + smemn) - - return diff --git a/src/hsp2/hsp2/main.py b/src/hsp2/hsp2/main.py index 41d24e49..5564d2b4 100644 --- a/src/hsp2/hsp2/main.py +++ b/src/hsp2/hsp2/main.py @@ -1,25 +1,62 @@ -''' Copyright (c) 2020 by RESPEC, INC. +"""Copyright (c) 2020 by RESPEC, INC. Author: Robert Heaphy, Ph.D. License: LGPL2 -''' + # This table defines the expansion to INFLOW, ROFLOW, OFLOW for RCHRES networks + d = [ + ['IVOL', 'ROVOL', 'OVOL', 'HYDRFG', 'HYDR'], + ['ICON', 'ROCON', 'OCON', 'CONSFG', 'CONS'], + ['IHEAT', 'ROHEAT', 'OHEAT', 'HTFG', 'HTRCH'], + ['ISED', 'ROSED', 'OSED', 'SEDFG', 'SEDTRN'], + ['IDQAL', 'RODQAL', 'ODQAL', 'GQALFG', 'GQUAL'], + ['ISQAL', 'ROSQAL', 'OSQAL', 'GQALFG', 'GQUAL'], + ['OXIF', 'OXCF1', 'OXCF2', 'OXFG', 'OXRX'], + ['NUIF1', 'NUCF1', 'NUCF1', 'NUTFG', 'NUTRX'], + ['NUIF2', 'NUCF2', 'NUCF9', 'NUTFG', 'NUTRX'], + ['PKIF', 'PKCF1', 'PKCH2', 'PLKFG', 'PLANK'], + ['PHIF', 'PHCF1', 'PHCF2', 'PHFG', 'PHCARB']] + df = pd.DataFrame(d, columns=['INFLOW', 'ROFLOW', 'OFLOW', 'Flag', 'Name']) + df.to_hdf(h2name, '/FLOWEXPANSION', format='t', data_columns=True) +""" + +import os +from datetime import datetime as dt from re import S -from numpy import float64, float32 +from typing import Union + +import pandas as pd +from numpy import float32, float64 from pandas import DataFrame, date_range from pandas.tseries.offsets import Minute -from datetime import datetime as dt -from typing import Union -import os -from hsp2.hsp2io.hdf import HDF5 -from hsp2.hsp2.utilities import versions, get_timeseries, expand_timeseries_names, save_timeseries, get_gener_timeseries -from hsp2.hsp2.configuration import activities, noop, expand_masslinks -from hsp2.hsp2.state import init_state_dicts, state_siminfo_hsp2, state_load_dynamics_hsp2, state_init_hsp2, state_context_hsp2 -from hsp2.hsp2.om import om_init_state, state_om_model_run_prep, state_load_dynamics_om + +from hsp2.hsp2.configuration import activities, expand_masslinks, noop +from hsp2.hsp2.om import om_init_state, state_load_dynamics_om, state_om_model_run_prep from hsp2.hsp2.SPECL import specl_load_state +from hsp2.hsp2.state import ( + init_state_dicts, + state_context_hsp2, + state_init_hsp2, + state_load_dynamics_hsp2, + state_siminfo_hsp2, +) +from hsp2.hsp2.utilities import ( + expand_timeseries_names, + get_gener_timeseries, + get_timeseries, + save_timeseries, + versions, +) +from hsp2.hsp2io.hdf import HDF5 +from hsp2.hsp2io.io import Category, IOManager, SupportsReadTS -from hsp2.hsp2io.io import IOManager, SupportsReadTS, Category +# special in-memory version of these +from .HSP2utilitiesInMem import get_timeseries_InMem, save_timeseries_InMem +from hsp2.hsp2tools.UciDictToObject import convert_uci_InMem -def main(io_manager:Union[str, IOManager], saveall:bool=False, jupyterlab:bool=True) -> None: + +def main( + io_manager: Union[str, IOManager], saveall: bool = False, jupyterlab: bool = True +) -> None: """ Run main HSP2 program. Parameters @@ -30,36 +67,42 @@ def main(io_manager:Union[str, IOManager], saveall:bool=False, jupyterlab:bool=T Saves all calculated data ignoring SAVE tables. jupyterlab: bool, default=True Flag for specific output behavior for jupyter lab. - + Return ------------ None - + """ if isinstance(io_manager, str): hdf5_instance = HDF5(io_manager) io_manager = IOManager(hdf5_instance) - hdfname = io_manager._input.file_path - if not os.path.exists(hdfname): - raise FileNotFoundError(f'{hdfname} HDF5 File Not Found') + elif isinstance(io_manager, dict): + hdfname = '' + else: + hdfname = io_manager._input.file_path + if not os.path.exists(hdfname): + raise FileNotFoundError(f"{hdfname} HDF5 File Not Found") msg = messages() - msg(1, f'Processing started for file {hdfname}; saveall={saveall}') + msg(1, f"Processing started for file {hdfname}; saveall={saveall}") # read user control, parameters, states, and flags uci and map to local variables - uci_obj = io_manager.read_uci() + if isinstance(io_manager, dict): + uci_obj = convert_uci_InMem(io_manager) + else: + uci_obj = io_manager.read_uci() opseq = uci_obj.opseq ddlinks = uci_obj.ddlinks ddmasslinks = uci_obj.ddmasslinks ddext_sources = uci_obj.ddext_sources ddgener = uci_obj.ddgener uci = uci_obj.uci - siminfo = uci_obj.siminfo + siminfo = uci_obj.siminfo ftables = uci_obj.ftables specactions = uci_obj.specactions monthdata = uci_obj.monthdata - - start, stop = siminfo['start'], siminfo['stop'] + + start, stop = siminfo["start"], siminfo["stop"] copy_instances = {} gener_instances = {} @@ -71,190 +114,356 @@ def main(io_manager:Union[str, IOManager], saveall:bool=False, jupyterlab:bool=T state = init_state_dicts() state_siminfo_hsp2(uci_obj, siminfo) # Add support for dynamic functions to operate on STATE - # - Load any dynamic components if present, and store variables on objects - state_load_dynamics_hsp2(state, io_manager, siminfo) - # Iterate through all segments and add crucial paths to state + if not isinstance(io_manager, dict): + # - Load any dynamic components if present, and store variables on objects + state_load_dynamics_hsp2(state, io_manager, siminfo) + else: + state['state_step_hydr'] = 'disabled' + state['hsp2_local_py'] = False + state["state_step_om"] = "disabled" + # Iterate through all segments and add crucial paths to state # before loading dynamic components that may reference them state_init_hsp2(state, opseq, activities) # - finally stash specactions in state, not domain (segment) dependent so do it once - state['specactions'] = specactions # stash the specaction dict in state - om_init_state(state) # set up operational model specific state entries - specl_load_state(state, io_manager, siminfo) # traditional special actions - state_load_dynamics_om(state, io_manager, siminfo) # operational model for custom python + state["specactions"] = specactions # stash the specaction dict in state + om_init_state(state) # set up operational model specific state entries + specl_load_state(state, io_manager, siminfo) # traditional special actions + if not isinstance(io_manager, dict): + state_load_dynamics_om( + state, io_manager, siminfo + ) # operational model for custom python # finalize all dynamically loaded components and prepare to run the model state_om_model_run_prep(state, io_manager, siminfo) ####################################################################################### # main processing loop - msg(1, f'Simulation Start: {start}, Stop: {stop}') + msg(1, f"Simulation Start: {start}, Stop: {stop}") tscat = {} for _, operation, segment, delt in opseq.itertuples(): - msg(2, f'{operation} {segment} DELT(minutes): {delt}') - siminfo['delt'] = delt - siminfo['tindex'] = date_range(start, stop, freq=Minute(delt))[1:] - siminfo['steps'] = len(siminfo['tindex']) - - if operation == 'COPY': - copy_instances[segment] = activities[operation](io_manager, siminfo, ddext_sources[(operation,segment)]) - elif operation == 'GENER': + msg(2, f"{operation} {segment} DELT(minutes): {delt}") + siminfo["delt"] = delt + siminfo["tindex"] = date_range(start, stop, freq=Minute(delt))[1:] + siminfo["steps"] = len(siminfo["tindex"]) + + if operation == "COPY": + copy_instances[segment] = activities[operation]( + io_manager, siminfo, ddext_sources[(operation, segment)] + ) + elif operation == "GENER": try: - ts = get_timeseries(io_manager, ddext_sources[(operation, segment)], siminfo) - ts = get_gener_timeseries(ts, gener_instances, ddlinks[segment], ddmasslinks) - get_flows(io_manager, ts, {}, uci, segment, ddlinks, ddmasslinks, siminfo['steps'], msg) - gener_instances[segment] = activities[operation](segment, siminfo, copy_instances, gener_instances, ddlinks, ddmasslinks, ts, ddgener) + if not isinstance(io_manager, dict): + ts = get_timeseries( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + else: + ts = get_timeseries_InMem( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + ts = get_gener_timeseries( + ts, gener_instances, ddlinks[segment], ddmasslinks + ) + get_flows( + io_manager, + ts, + {}, + uci, + segment, + ddlinks, + ddmasslinks, + siminfo["steps"], + msg, + ) + gener_instances[segment] = activities[operation]( + segment, + siminfo, + copy_instances, + gener_instances, + ddlinks, + ddmasslinks, + ts, + ddgener, + ) except NotImplementedError as e: print(f"GENER '{segment}' may not function correctly. '{e}'") else: - # now conditionally execute all activity modules for the op, segment - ts = get_timeseries(io_manager,ddext_sources[(operation,segment)],siminfo) - ts = get_gener_timeseries(ts, gener_instances, ddlinks[segment],ddmasslinks) - flags = uci[(operation, 'GENERAL', segment)]['ACTIVITY'] - if operation == 'RCHRES': + if not isinstance(io_manager, dict): + ts = get_timeseries( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + else: + ts = get_timeseries_InMem( + io_manager, ddext_sources[(operation, segment)], siminfo + ) + ts = get_gener_timeseries( + ts, gener_instances, ddlinks[segment], ddmasslinks + ) + flags = uci[(operation, "GENERAL", segment)]["ACTIVITY"] + if operation == "RCHRES": # Add nutrient adsorption flags: - if flags['NUTRX'] == 1: - flags['TAMFG'] = uci[(operation, 'NUTRX', segment)]['FLAGS']['NH3FG'] - flags['ADNHFG'] = uci[(operation, 'NUTRX', segment)]['FLAGS']['ADNHFG'] - flags['PO4FG'] = uci[(operation, 'NUTRX', segment)]['FLAGS']['PO4FG'] - flags['ADPOFG'] = uci[(operation, 'NUTRX', segment)]['FLAGS']['ADPOFG'] - - get_flows(io_manager, ts, flags, uci, segment, ddlinks, ddmasslinks, siminfo['steps'], msg) + if flags["NUTRX"] == 1: + flags["TAMFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "NH3FG" + ] + flags["ADNHFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "ADNHFG" + ] + flags["PO4FG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "PO4FG" + ] + flags["ADPOFG"] = uci[(operation, "NUTRX", segment)]["FLAGS"][ + "ADPOFG" + ] + + get_flows( + io_manager, + ts, + flags, + uci, + segment, + ddlinks, + ddmasslinks, + siminfo["steps"], + msg, + ) for activity, function in activities[operation].items(): - if function == noop: #or not flags[activity]: + if function == noop: # or not flags[activity]: continue if (activity in flags) and (not flags[activity]): continue - if (activity == 'RQUAL') and (not flags['OXRX']) and (not flags['NUTRX']) and (not flags['PLANK']) and (not flags['PHCARB']): + if ( + (activity == "RQUAL") + and (not flags["OXRX"]) + and (not flags["NUTRX"]) + and (not flags["PLANK"]) + and (not flags["PHCARB"]) + ): continue - msg(3, f'{activity}') + msg(3, f"{activity}") # Set context for dynamic executables and special actions state_context_hsp2(state, operation, segment, activity) - - ui = uci[(operation, activity, segment)] # ui is a dictionary - if operation == 'PERLND' and activity == 'SEDMNT': + + ui = uci[(operation, activity, segment)] # ui is a dictionary + if operation == "PERLND" and activity == "SEDMNT": # special exception here to make CSNOFG available - ui['PARAMETERS']['CSNOFG'] = uci[(operation, 'PWATER', segment)]['PARAMETERS']['CSNOFG'] - if operation == 'PERLND' and activity == 'PSTEMP': + ui["PARAMETERS"]["CSNOFG"] = uci[(operation, "PWATER", segment)][ + "PARAMETERS" + ]["CSNOFG"] + if operation == "PERLND" and activity == "PSTEMP": # special exception here to make AIRTFG available - ui['PARAMETERS']['AIRTFG'] = flags['ATEMP'] - if operation == 'PERLND' and activity == 'PWTGAS': + ui["PARAMETERS"]["AIRTFG"] = flags["ATEMP"] + if operation == "PERLND" and activity == "PWTGAS": # special exception here to make CSNOFG available - ui['PARAMETERS']['CSNOFG'] = uci[(operation, 'PWATER', segment)]['PARAMETERS']['CSNOFG'] - if operation == 'RCHRES': - if not 'PARAMETERS' in ui: - ui['PARAMETERS'] = {} - ui['PARAMETERS']['NEXITS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['NEXITS'] - if activity == 'ADCALC': - ui['PARAMETERS']['ADFG'] = flags['ADCALC'] - ui['PARAMETERS']['KS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['KS'] - ui['PARAMETERS']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - ui['PARAMETERS']['ROS'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['ROS'] - nexits = uci[(operation, 'HYDR', segment)]['PARAMETERS']['NEXITS'] + ui["PARAMETERS"]["CSNOFG"] = uci[(operation, "PWATER", segment)][ + "PARAMETERS" + ]["CSNOFG"] + if operation == "RCHRES": + if not "PARAMETERS" in ui: + ui["PARAMETERS"] = {} + ui["PARAMETERS"]["NEXITS"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["NEXITS"] + if activity == "ADCALC": + ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] + ui["PARAMETERS"]["KS"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["KS"] + ui["PARAMETERS"]["VOL"] = uci[(operation, "HYDR", segment)][ + "STATES" + ]["VOL"] + ui["PARAMETERS"]["ROS"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["ROS"] + nexits = uci[(operation, "HYDR", segment)]["PARAMETERS"][ + "NEXITS" + ] for index in range(nexits): - ui['PARAMETERS']['OS' + str(index + 1)] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['OS'+ str(index + 1)] - if activity == 'HTRCH': - ui['PARAMETERS']['ADFG'] = flags['ADCALC'] - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] + ui["PARAMETERS"]["OS" + str(index + 1)] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["OS" + str(index + 1)] + if activity == "HTRCH": + ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - if activity == 'CONS': - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] - if activity == 'SEDTRN': - ui['PARAMETERS']['ADFG'] = flags['ADCALC'] - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] + if activity == "CONS": + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + if activity == "SEDTRN": + ui["PARAMETERS"]["ADFG"] = flags["ADCALC"] + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] # ui['STATES']['VOL'] = uci[(operation, 'HYDR', segment)]['STATES']['VOL'] - ui['PARAMETERS']['HTFG'] = flags['HTRCH'] - ui['PARAMETERS']['AUX3FG'] = 0 - if flags['HYDR']: - ui['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN'] - ui['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH'] - ui['PARAMETERS']['DB50'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DB50'] - ui['PARAMETERS']['AUX3FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX3FG'] - if activity == 'GQUAL': - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] - ui['PARAMETERS']['HTFG'] = flags['HTRCH'] - ui['PARAMETERS']['SEDFG'] = flags['SEDTRN'] + ui["PARAMETERS"]["HTFG"] = flags["HTRCH"] + ui["PARAMETERS"]["AUX3FG"] = 0 + if flags["HYDR"]: + ui["PARAMETERS"]["LEN"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["LEN"] + ui["PARAMETERS"]["DELTH"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DELTH"] + ui["PARAMETERS"]["DB50"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DB50"] + ui["PARAMETERS"]["AUX3FG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["AUX3FG"] + if activity == "GQUAL": + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + ui["PARAMETERS"]["HTFG"] = flags["HTRCH"] + ui["PARAMETERS"]["SEDFG"] = flags["SEDTRN"] # ui['PARAMETERS']['REAMFG'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['REAMFG'] - ui['PARAMETERS']['HYDRFG'] = flags['HYDR'] - if flags['HYDR']: - ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] - ui['PARAMETERS']['AUX1FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX1FG'] - ui['PARAMETERS']['AUX2FG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['AUX2FG'] - ui['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN'] - ui['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH'] - if flags['OXRX']: - ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] - ui['PARAMETERS']['CFOREA'] = uci[(operation, 'OXRX', segment)]['PARAMETERS']['CFOREA'] - if flags['SEDTRN']: - ui['PARAMETERS']['SSED1'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED1'] - ui['PARAMETERS']['SSED2'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED2'] - ui['PARAMETERS']['SSED3'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED3'] - if flags['HTRCH']: - ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'HTRCH', segment)]['PARAMETERS']['CFSAEX'] - elif flags['PLANK']: - if 'CFSAEX' in uci[(operation, 'PLANK', segment)]['PARAMETERS']: - ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'PLANK', segment)]['PARAMETERS']['CFSAEX'] - - if activity == 'RQUAL': + ui["PARAMETERS"]["HYDRFG"] = flags["HYDR"] + if flags["HYDR"]: + ui["PARAMETERS"]["LKFG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LKFG"] + ui["PARAMETERS"]["AUX1FG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["AUX1FG"] + ui["PARAMETERS"]["AUX2FG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["AUX2FG"] + ui["PARAMETERS"]["LEN"] = uci[(operation, "HYDR", segment)][ + "PARAMETERS" + ]["LEN"] + ui["PARAMETERS"]["DELTH"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DELTH"] + if flags["OXRX"]: + ui["PARAMETERS"]["LKFG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LKFG"] + ui["PARAMETERS"]["CFOREA"] = uci[ + (operation, "OXRX", segment) + ]["PARAMETERS"]["CFOREA"] + if flags["SEDTRN"]: + ui["PARAMETERS"]["SSED1"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED1"] + ui["PARAMETERS"]["SSED2"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED2"] + ui["PARAMETERS"]["SSED3"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED3"] + if flags["HTRCH"]: + ui["PARAMETERS"]["CFSAEX"] = uci[ + (operation, "HTRCH", segment) + ]["PARAMETERS"]["CFSAEX"] + elif flags["PLANK"]: + if ( + "CFSAEX" + in uci[(operation, "PLANK", segment)]["PARAMETERS"] + ): + ui["PARAMETERS"]["CFSAEX"] = uci[ + (operation, "PLANK", segment) + ]["PARAMETERS"]["CFSAEX"] + + if activity == "RQUAL": # RQUAL inputs: - ui['advectData'] = uci[(operation, 'ADCALC', segment)]['adcalcData'] - if flags['HYDR']: - ui['PARAMETERS']['LKFG'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LKFG'] - - ui['FLAGS']['HTFG'] = flags['HTRCH'] - ui['FLAGS']['SEDFG'] = flags['SEDTRN'] - ui['FLAGS']['GQFG'] = flags['GQUAL'] - ui['FLAGS']['OXFG'] = flags['OXFG'] - ui['FLAGS']['NUTFG'] = flags['NUTRX'] - ui['FLAGS']['PLKFG'] = flags['PLANK'] - ui['FLAGS']['PHFG'] = flags['PHCARB'] - if flags['CONS']: - if 'PARAMETERS' in uci[(operation, 'CONS', segment)]: - if 'NCONS' in uci[(operation, 'CONS', segment)]['PARAMETERS']: - ui['PARAMETERS']['NCONS'] = uci[(operation, 'CONS', segment)]['PARAMETERS']['NCONS'] + ui["advectData"] = uci[(operation, "ADCALC", segment)][ + "adcalcData" + ] + if flags["HYDR"]: + ui["PARAMETERS"]["LKFG"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LKFG"] + + ui["FLAGS"]["HTFG"] = flags["HTRCH"] + ui["FLAGS"]["SEDFG"] = flags["SEDTRN"] + ui["FLAGS"]["GQFG"] = flags["GQUAL"] + ui["FLAGS"]["OXFG"] = flags["OXFG"] + ui["FLAGS"]["NUTFG"] = flags["NUTRX"] + ui["FLAGS"]["PLKFG"] = flags["PLANK"] + ui["FLAGS"]["PHFG"] = flags["PHCARB"] + if flags["CONS"]: + if "PARAMETERS" in uci[(operation, "CONS", segment)]: + if ( + "NCONS" + in uci[(operation, "CONS", segment)]["PARAMETERS"] + ): + ui["PARAMETERS"]["NCONS"] = uci[ + (operation, "CONS", segment) + ]["PARAMETERS"]["NCONS"] # OXRX module inputs: - ui_oxrx = uci[(operation, 'OXRX', segment)] - - if flags['HYDR']: - ui_oxrx['PARAMETERS']['LEN'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['LEN'] - ui_oxrx['PARAMETERS']['DELTH'] = uci[(operation, 'HYDR', segment)]['PARAMETERS']['DELTH'] - - if flags['HTRCH']: - ui_oxrx['PARAMETERS']['ELEV'] = uci[(operation, 'HTRCH', segment)]['PARAMETERS']['ELEV'] - - if flags['SEDTRN']: - ui['PARAMETERS']['SSED1'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED1'] - ui['PARAMETERS']['SSED2'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED2'] - ui['PARAMETERS']['SSED3'] = uci[(operation, 'SEDTRN', segment)]['STATES']['SSED3'] + ui_oxrx = uci[(operation, "OXRX", segment)] + + if flags["HYDR"]: + ui_oxrx["PARAMETERS"]["LEN"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["LEN"] + ui_oxrx["PARAMETERS"]["DELTH"] = uci[ + (operation, "HYDR", segment) + ]["PARAMETERS"]["DELTH"] + + if flags["HTRCH"]: + ui_oxrx["PARAMETERS"]["ELEV"] = uci[ + (operation, "HTRCH", segment) + ]["PARAMETERS"]["ELEV"] + + if flags["SEDTRN"]: + ui["PARAMETERS"]["SSED1"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED1"] + ui["PARAMETERS"]["SSED2"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED2"] + ui["PARAMETERS"]["SSED3"] = uci[ + (operation, "SEDTRN", segment) + ]["STATES"]["SSED3"] # PLANK module inputs: - if flags['HTRCH']: - ui['PARAMETERS']['CFSAEX'] = uci[(operation, 'HTRCH', segment)]['PARAMETERS']['CFSAEX'] + if flags["HTRCH"]: + ui["PARAMETERS"]["CFSAEX"] = uci[ + (operation, "HTRCH", segment) + ]["PARAMETERS"]["CFSAEX"] # NUTRX, PLANK, PHCARB module inputs: - ui_nutrx = uci[(operation, 'NUTRX', segment)] - ui_plank = uci[(operation, 'PLANK', segment)] - ui_phcarb = uci[(operation, 'PHCARB', segment)] + ui_nutrx = uci[(operation, "NUTRX", segment)] + ui_plank = uci[(operation, "PLANK", segment)] + ui_phcarb = uci[(operation, "PHCARB", segment)] ############ calls activity function like snow() ############## - if operation not in ['COPY','GENER']: - if (activity == 'HYDR'): - errors, errmessages = function(io_manager, siminfo, ui, ts, ftables, state) - elif (activity == 'SEDTRN'): - errors, errmessages = function(io_manager, siminfo, ui, ts, state) - elif (activity != 'RQUAL'): + if operation not in ["COPY", "GENER"]: + if activity == "HYDR": + errors, errmessages = function( + io_manager, siminfo, ui, ts, ftables, state + ) + elif activity == "SEDTRN": + errors, errmessages = function( + io_manager, siminfo, ui, ts, state + ) + elif activity != "RQUAL": errors, errmessages = function(io_manager, siminfo, ui, ts) - else: - errors, errmessages = function(io_manager, siminfo, ui, ui_oxrx, ui_nutrx, ui_plank, ui_phcarb, ts, monthdata) + else: + errors, errmessages = function( + io_manager, + siminfo, + ui, + ui_oxrx, + ui_nutrx, + ui_plank, + ui_phcarb, + ts, + monthdata, + ) ############################################################### for errorcnt, errormsg in zip(errors, errmessages): if errorcnt > 0: - msg(4, f'Error count {errorcnt}: {errormsg}') + msg(4, f"Error count {errorcnt}: {errormsg}") # default to hourly output outstep = 2 @@ -262,165 +471,371 @@ def main(io_manager:Union[str, IOManager], saveall:bool=False, jupyterlab:bool=T outstep_nutrx = 2 outstep_plank = 2 outstep_phcarb = 2 - if 'BINOUT' in uci[(operation, 'GENERAL', segment)]: - if activity in uci[(operation, 'GENERAL', segment)]['BINOUT']: - outstep = uci[(operation, 'GENERAL', segment)]['BINOUT'][activity] - elif activity == 'RQUAL': - outstep_oxrx = uci[(operation, 'GENERAL', segment)]['BINOUT']['OXRX'] - outstep_nutrx = uci[(operation, 'GENERAL', segment)]['BINOUT']['NUTRX'] - outstep_plank = uci[(operation, 'GENERAL', segment)]['BINOUT']['PLANK'] - outstep_phcarb = uci[(operation, 'GENERAL', segment)]['BINOUT']['PHCARB'] - - if 'SAVE' in ui: - save_timeseries(io_manager,ts,ui['SAVE'],siminfo,saveall,operation,segment,activity,jupyterlab,outstep) - - if (activity == 'RQUAL'): - if 'SAVE' in ui_oxrx: save_timeseries(io_manager,ts,ui_oxrx['SAVE'],siminfo,saveall,operation,segment,'OXRX',jupyterlab,outstep_oxrx) - if 'SAVE' in ui_nutrx and flags['NUTRX'] == 1: save_timeseries(io_manager,ts,ui_nutrx['SAVE'],siminfo,saveall,operation,segment,'NUTRX',jupyterlab,outstep_nutrx) - if 'SAVE' in ui_plank and flags['PLANK'] == 1: save_timeseries(io_manager,ts,ui_plank['SAVE'],siminfo,saveall,operation,segment,'PLANK',jupyterlab,outstep_plank) - if 'SAVE' in ui_phcarb and flags['PHCARB'] == 1: save_timeseries(io_manager,ts,ui_phcarb['SAVE'],siminfo,saveall,operation,segment,'PHCARB',jupyterlab,outstep_phcarb) - - msglist = msg(1, 'Done', final=True) - - df = DataFrame(msglist, columns=['logfile']) - io_manager.write_log(df) - - if jupyterlab: - df = versions(['jupyterlab', 'notebook']) - io_manager.write_versioning(df) - print('\n\n', df) + if "BINOUT" in uci[(operation, "GENERAL", segment)]: + if activity in uci[(operation, "GENERAL", segment)]["BINOUT"]: + outstep = uci[(operation, "GENERAL", segment)]["BINOUT"][ + activity + ] + elif activity == "RQUAL": + outstep_oxrx = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "OXRX" + ] + outstep_nutrx = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "NUTRX" + ] + outstep_plank = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "PLANK" + ] + outstep_phcarb = uci[(operation, "GENERAL", segment)]["BINOUT"][ + "PHCARB" + ] + + if "SAVE" in ui: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui["SAVE"], + siminfo, + saveall, + operation, + segment, + activity, + jupyterlab, + outstep, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui["SAVE"], + siminfo, + saveall, + operation, + segment, + activity, + jupyterlab, + outstep, + ) + + if activity == "RQUAL": + if "SAVE" in ui_oxrx: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_oxrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "OXRX", + jupyterlab, + outstep_oxrx, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_oxrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "OXRX", + jupyterlab, + outstep_oxrx, + ) + if "SAVE" in ui_nutrx and flags["NUTRX"] == 1: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_nutrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "NUTRX", + jupyterlab, + outstep_nutrx, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_nutrx["SAVE"], + siminfo, + saveall, + operation, + segment, + "NUTRX", + jupyterlab, + outstep_nutrx, + ) + if "SAVE" in ui_plank and flags["PLANK"] == 1: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_plank["SAVE"], + siminfo, + saveall, + operation, + segment, + "PLANK", + jupyterlab, + outstep_plank, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_plank["SAVE"], + siminfo, + saveall, + operation, + segment, + "PLANK", + jupyterlab, + outstep_plank, + ) + if "SAVE" in ui_phcarb and flags["PHCARB"] == 1: + if not isinstance(io_manager, dict): + save_timeseries( + io_manager, + ts, + ui_phcarb["SAVE"], + siminfo, + saveall, + operation, + segment, + "PHCARB", + jupyterlab, + outstep_phcarb, + ) + else: + save_timeseries_InMem( + io_manager, + ts, + ui_phcarb["SAVE"], + siminfo, + saveall, + operation, + segment, + "PHCARB", + jupyterlab, + outstep_phcarb, + ) + + msglist = msg(1, "Done", final=True) + + df = DataFrame(msglist, columns=["logfile"]) + if isinstance(io_manager, dict): + io_manager["logfile"] = df + else: + io_manager.write_log(df) + if jupyterlab: + df = versions(["jupyterlab", "notebook"]) + io_manager.write_versioning(df) + print("\n\n", df) return + def messages(): - '''Closure routine; msg() prints messages to screen and run log''' + """Closure routine; msg() prints messages to screen and run log""" start = dt.now() mlist = [] + def msg(indent, message, final=False): now = dt.now() - m = str(now)[:22] + ' ' * indent + message + m = str(now)[:22] + " " * indent + message if final: - mn,sc = divmod((now-start).seconds, 60) - ms = (now-start).microseconds // 100_000 - m = '; '.join((m, f'Run time is about {mn:02}:{sc:02}.{ms} (mm:ss)')) + mn, sc = divmod((now - start).seconds, 60) + ms = (now - start).microseconds // 100_000 + m = "; ".join((m, f"Run time is about {mn:02}:{sc:02}.{ms} (mm:ss)")) print(m) mlist.append(m) return mlist + return msg -def get_flows(io_manager:SupportsReadTS, ts, flags, uci, segment, ddlinks, ddmasslinks, steps, msg): + +def get_flows( + io_manager: SupportsReadTS, + ts, + flags, + uci, + segment, + ddlinks, + ddmasslinks, + steps, + msg, +): # get inflows to this operation for x in ddlinks[segment]: - if x.SVOL != 'GENER': # gener already handled in get_gener_timeseries + if x.SVOL != "GENER": # gener already handled in get_gener_timeseries recs = [] - if x.MLNO == '': # Data from NETWORK part of Links table + if x.MLNO == "": # Data from NETWORK part of Links table rec = {} - rec['MFACTOR'] = x.MFACTOR - rec['SGRPN'] = x.SGRPN - rec['SMEMN'] = x.SMEMN - rec['SMEMSB1'] = x.SMEMSB1 - rec['SMEMSB2'] = x.SMEMSB2 - rec['TMEMN'] = x.TMEMN - rec['TMEMSB1'] = x.TMEMSB1 - rec['TMEMSB2'] = x.TMEMSB2 - rec['SVOL'] = x.SVOL + rec["MFACTOR"] = x.MFACTOR + rec["SGRPN"] = x.SGRPN + rec["SMEMN"] = x.SMEMN + rec["SMEMSB1"] = x.SMEMSB1 + rec["SMEMSB2"] = x.SMEMSB2 + rec["TMEMN"] = x.TMEMN + rec["TMEMSB1"] = x.TMEMSB1 + rec["TMEMSB2"] = x.TMEMSB2 + rec["SVOL"] = x.SVOL recs.append(rec) - else: # Data from SCHEMATIC part of Links table + else: # Data from SCHEMATIC part of Links table mldata = ddmasslinks[x.MLNO] for dat in mldata: - if dat.SMEMN != '': + if dat.SMEMN != "": rec = {} - rec['MFACTOR'] = dat.MFACTOR - rec['SGRPN'] = dat.SGRPN - rec['SMEMN'] = dat.SMEMN - rec['SMEMSB1'] = dat.SMEMSB1 - rec['SMEMSB2'] = dat.SMEMSB2 - rec['TMEMN'] = dat.TMEMN - rec['TMEMSB1'] = dat.TMEMSB1 - rec['TMEMSB2'] = dat.TMEMSB2 - rec['SVOL'] = dat.SVOL + rec["MFACTOR"] = dat.MFACTOR + rec["SGRPN"] = dat.SGRPN + rec["SMEMN"] = dat.SMEMN + rec["SMEMSB1"] = dat.SMEMSB1 + rec["SMEMSB2"] = dat.SMEMSB2 + rec["TMEMN"] = dat.TMEMN + rec["TMEMSB1"] = dat.TMEMSB1 + rec["TMEMSB2"] = dat.TMEMSB2 + rec["SVOL"] = dat.SVOL recs.append(rec) else: # this is the kind that needs to be expanded if dat.SGRPN == "ROFLOW" or dat.SGRPN == "OFLOW": - recs = expand_masslinks(flags,uci,dat,recs) + recs = expand_masslinks(flags, uci, dat, recs) for rec in recs: - mfactor = rec['MFACTOR'] - sgrpn = rec['SGRPN'] - smemn = rec['SMEMN'] - smemsb1 = rec['SMEMSB1'] - smemsb2 = rec['SMEMSB2'] - tmemn = rec['TMEMN'] - tmemsb1 = rec['TMEMSB1'] - tmemsb2 = rec['TMEMSB2'] - - if x.AFACTR != '': + mfactor = rec["MFACTOR"] + sgrpn = rec["SGRPN"] + smemn = rec["SMEMN"] + smemsb1 = rec["SMEMSB1"] + smemsb2 = rec["SMEMSB2"] + tmemn = rec["TMEMN"] + tmemsb1 = rec["TMEMSB1"] + tmemsb2 = rec["TMEMSB2"] + + if x.AFACTR != "": afactr = x.AFACTR factor = afactr * mfactor else: factor = mfactor # KLUDGE until remaining HSP2 modules are available. - if tmemn not in {'IVOL', 'ICON', 'IHEAT', 'ISED', 'ISED1', 'ISED2', 'ISED3', - 'IDQAL', 'ISQAL1', 'ISQAL2', 'ISQAL3', - 'OXIF', 'NUIF1', 'NUIF2', 'PKIF', 'PHIF', - 'ONE', 'TWO'}: + if tmemn not in { + "IVOL", + "ICON", + "IHEAT", + "ISED", + "ISED1", + "ISED2", + "ISED3", + "IDQAL", + "ISQAL1", + "ISQAL2", + "ISQAL3", + "OXIF", + "NUIF1", + "NUIF2", + "PKIF", + "PHIF", + "ONE", + "TWO", + }: continue - if (sgrpn == 'OFLOW' and smemn == 'OVOL') or (sgrpn == 'ROFLOW' and smemn == 'ROVOL'): - sgrpn = 'HYDR' - if (sgrpn == 'OFLOW' and smemn == 'OHEAT') or (sgrpn == 'ROFLOW' and smemn == 'ROHEAT'): - sgrpn = 'HTRCH' - if (sgrpn == 'OFLOW' and smemn == 'OSED') or (sgrpn == 'ROFLOW' and smemn == 'ROSED'): - sgrpn = 'SEDTRN' - if (sgrpn == 'OFLOW' and smemn == 'ODQAL') or (sgrpn == 'ROFLOW' and smemn == 'RODQAL'): - sgrpn = 'GQUAL' - if (sgrpn == 'OFLOW' and smemn == 'OSQAL') or (sgrpn == 'ROFLOW' and smemn == 'ROSQAL'): - sgrpn = 'GQUAL' - if (sgrpn == 'OFLOW' and smemn == 'OXCF2') or (sgrpn == 'ROFLOW' and smemn == 'OXCF1'): - sgrpn = 'OXRX' - if (sgrpn == 'OFLOW' and (smemn == 'NUCF9' or smemn == 'OSNH4' or smemn == 'OSPO4')) or (sgrpn == 'ROFLOW' and (smemn == 'NUCF1' or smemn == 'NUFCF2')): - sgrpn = 'NUTRX' - if (sgrpn == 'OFLOW' and smemn == 'PKCF2') or (sgrpn == 'ROFLOW' and smemn == 'PKCF1'): - sgrpn = 'PLANK' - if (sgrpn == 'OFLOW' and smemn == 'PHCF2') or (sgrpn == 'ROFLOW' and smemn == 'PHCF1'): - sgrpn = 'PHCARB' - - if tmemn == 'ISED' or tmemn == 'ISQAL': - tmemn = tmemn + tmemsb1 # need to add sand, silt, clay subscript - if (sgrpn == 'HYDR' and smemn == 'OVOL') or (sgrpn == 'HTRCH' and smemn == 'OHEAT'): - smemsb2 = '' - if sgrpn == 'GQUAL' and smemsb2 == '': - smemsb2 = '1' - - smemn, tmemn = expand_timeseries_names(sgrpn, smemn, smemsb1, smemsb2, tmemn, tmemsb1, tmemsb2) - - path = f'RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' - MFname = f'{x.SVOL}{x.SVOLNO}_MFACTOR' - AFname = f'{x.SVOL}{x.SVOLNO}_AFACTR' - data = f'{smemn}{smemsb1}{smemsb2}' - - data_frame = io_manager.read_ts(Category.RESULTS,x.SVOL,x.SVOLNO, sgrpn) + if (sgrpn == "OFLOW" and smemn == "OVOL") or ( + sgrpn == "ROFLOW" and smemn == "ROVOL" + ): + sgrpn = "HYDR" + if (sgrpn == "OFLOW" and smemn == "OHEAT") or ( + sgrpn == "ROFLOW" and smemn == "ROHEAT" + ): + sgrpn = "HTRCH" + if (sgrpn == "OFLOW" and smemn == "OSED") or ( + sgrpn == "ROFLOW" and smemn == "ROSED" + ): + sgrpn = "SEDTRN" + if (sgrpn == "OFLOW" and smemn == "ODQAL") or ( + sgrpn == "ROFLOW" and smemn == "RODQAL" + ): + sgrpn = "GQUAL" + if (sgrpn == "OFLOW" and smemn == "OSQAL") or ( + sgrpn == "ROFLOW" and smemn == "ROSQAL" + ): + sgrpn = "GQUAL" + if (sgrpn == "OFLOW" and smemn == "OXCF2") or ( + sgrpn == "ROFLOW" and smemn == "OXCF1" + ): + sgrpn = "OXRX" + if ( + sgrpn == "OFLOW" + and (smemn == "NUCF9" or smemn == "OSNH4" or smemn == "OSPO4") + ) or (sgrpn == "ROFLOW" and (smemn == "NUCF1" or smemn == "NUFCF2")): + sgrpn = "NUTRX" + if (sgrpn == "OFLOW" and smemn == "PKCF2") or ( + sgrpn == "ROFLOW" and smemn == "PKCF1" + ): + sgrpn = "PLANK" + if (sgrpn == "OFLOW" and smemn == "PHCF2") or ( + sgrpn == "ROFLOW" and smemn == "PHCF1" + ): + sgrpn = "PHCARB" + + if tmemn == "ISED" or tmemn == "ISQAL": + tmemn = tmemn + tmemsb1 # need to add sand, silt, clay subscript + if (sgrpn == "HYDR" and smemn == "OVOL") or ( + sgrpn == "HTRCH" and smemn == "OHEAT" + ): + smemsb2 = "" + if sgrpn == "GQUAL" and smemsb2 == "": + smemsb2 = "1" + + smemn, tmemn = expand_timeseries_names( + sgrpn, smemn, smemsb1, smemsb2, tmemn, tmemsb1, tmemsb2 + ) + + path = f"RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}" + MFname = f"{x.SVOL}{x.SVOLNO}_MFACTOR" + AFname = f"{x.SVOL}{x.SVOLNO}_AFACTR" + data = f"{smemn}{smemsb1}{smemsb2}" + + if not isinstance(io_manager, dict): + data_frame = io_manager.read_ts( + Category.RESULTS, x.SVOL, x.SVOLNO, sgrpn + ) + else: + path = f'/RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' + try: + data_frame = io_manager[path] + except KeyError: + data_frame = pd.DataFrame() + try: - if data in data_frame.columns: t = data_frame[data].astype(float64).to_numpy()[0:steps] - else: t = data_frame[smemn].astype(float64).to_numpy()[0:steps] + if data in data_frame.columns: + t = data_frame[data].astype(float64).to_numpy()[0:steps] + else: + t = data_frame[smemn].astype(float64).to_numpy()[0:steps] if MFname in ts and AFname in ts: t *= ts[MFname][:steps] * ts[AFname][0:steps] - msg(4, f'MFACTOR modified by timeseries {MFname}') - msg(4, f'AFACTR modified by timeseries {AFname}') + msg(4, f"MFACTOR modified by timeseries {MFname}") + msg(4, f"AFACTR modified by timeseries {AFname}") elif MFname in ts: t *= afactr * ts[MFname][0:steps] - msg(4, f'MFACTOR modified by timeseries {MFname}') + msg(4, f"MFACTOR modified by timeseries {MFname}") elif AFname in ts: t *= mfactor * ts[AFname][0:steps] - msg(4, f'AFACTR modified by timeseries {AFname}') + msg(4, f"AFACTR modified by timeseries {AFname}") else: t *= factor # if poht to iheat, imprecision in hspf conversion factor requires a slight adjustment - if (smemn == 'POHT' or smemn == 'SOHT') and tmemn == 'IHEAT': + if (smemn == "POHT" or smemn == "SOHT") and tmemn == "IHEAT": t *= 0.998553 - if (smemn == 'PODOXM' or smemn == 'SODOXM') and tmemn == 'OXIF1': + if (smemn == "PODOXM" or smemn == "SODOXM") and tmemn == "OXIF1": t *= 1.000565 # ??? ISSUE: can fetched data be at different frequency - don't know how to transform. @@ -430,27 +845,6 @@ def get_flows(io_manager:SupportsReadTS, ts, flags, uci, segment, ddlinks, ddmas ts[tmemn] = t except KeyError: - print('ERROR in FLOWS, cant resolve ', path + ' ' + smemn) + print("ERROR in FLOWS, cant resolve ", path + " " + smemn) return - -''' - - # This table defines the expansion to INFLOW, ROFLOW, OFLOW for RCHRES networks - d = [ - ['IVOL', 'ROVOL', 'OVOL', 'HYDRFG', 'HYDR'], - ['ICON', 'ROCON', 'OCON', 'CONSFG', 'CONS'], - ['IHEAT', 'ROHEAT', 'OHEAT', 'HTFG', 'HTRCH'], - ['ISED', 'ROSED', 'OSED', 'SEDFG', 'SEDTRN'], - ['IDQAL', 'RODQAL', 'ODQAL', 'GQALFG', 'GQUAL'], - ['ISQAL', 'ROSQAL', 'OSQAL', 'GQALFG', 'GQUAL'], - ['OXIF', 'OXCF1', 'OXCF2', 'OXFG', 'OXRX'], - ['NUIF1', 'NUCF1', 'NUCF1', 'NUTFG', 'NUTRX'], - ['NUIF2', 'NUCF2', 'NUCF9', 'NUTFG', 'NUTRX'], - ['PKIF', 'PKCF1', 'PKCH2', 'PLKFG', 'PLANK'], - ['PHIF', 'PHCF1', 'PHCF2', 'PHFG', 'PHCARB']] - df = pd.DataFrame(d, columns=['INFLOW', 'ROFLOW', 'OFLOW', 'Flag', 'Name']) - df.to_hdf(h2name, '/FLOWEXPANSION', format='t', data_columns=True) - - -''' From 6f96ff3313bb8ba728e1ab62e7d78d8b87c65307 Mon Sep 17 00:00:00 2001 From: Your Name Date: Mon, 8 Sep 2025 14:30:51 -0400 Subject: [PATCH 4/5] added prototype IOManagerDF to unify nonHDF and HDF code behind one interface --- src/hsp2/hsp2io/io.py | 55 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/src/hsp2/hsp2io/io.py b/src/hsp2/hsp2io/io.py index e41e0ecc..01ece829 100644 --- a/src/hsp2/hsp2io/io.py +++ b/src/hsp2/hsp2io/io.py @@ -2,6 +2,7 @@ from pandas.core.frame import DataFrame from hsp2.hsp2io.protocols import Category, SupportsReadUCI, SupportsReadTS, SupportsWriteTS, SupportsWriteLogging from typing import Union, List +from HSP2utilitiesInMem import convert_uci_InMem from hsp2.hsp2.uci import UCI @@ -103,6 +104,16 @@ def read_ts(self, data_frame = self._get_in_memory(category, operation, segment, activity) if data_frame is not None: return data_frame + data_frame = self.read_input(category, operation, segment, activity) + return data_frame + + # do the actual reading from the data store. This can be subclassed to support any manner of storage csv, wdm, etc + def read_input(self, + category:Category, + operation:Union[str,None]=None, + segment:Union[str,None]=None, + activity:Union[str,None]=None, + *args, **kwargs) -> pd.DataFrame: if category == Category.INPUTS: data_frame = self._input.read_ts(category, operation, segment, activity) key = (category, operation, segment, activity) @@ -111,7 +122,7 @@ def read_ts(self, if category == Category.RESULTS: return self._output.read_ts(category, operation, segment, activity) return pd.DataFrame - + def write_log(self, data_frame)-> None: if self._log: self._log.write_log(data_frame) @@ -128,3 +139,45 @@ def _get_in_memory(self, return self._in_memory[key].copy(deep=True) except KeyError: return None + +class IOManagerPandas(IOManager): + def __init__(self, + io_combined: Union[SupportsReadUCI, SupportsReadTS, SupportsWriteTS, None] = None, + uci: Union[SupportsReadUCI,None]=None, + input: Union[SupportsReadTS,None]=None, + output: Union[SupportsReadTS,SupportsWriteTS,None]=None, + log: Union[SupportsWriteLogging,None]=None,) -> None: + self._input = io_combined if input is None else input + self._output = io_combined if output is None else output + self._uci = io_combined if uci is None else uci + self._log = io_combined if log is None else log + self.uci_df = io_combined + # Essentially, we should be able to set the pandas df here and then we can fully support reads + self._in_memory = io_combined + + def read_uci(self, *args, **kwargs): + uci = convert_uci_InMem(args) + return uci + + def read_ts(self, + category:Category, + operation:Union[str,None]=None, + segment:Union[str,None]=None, + activity:Union[str,None]=None, + *args, **kwargs) -> pd.DataFrame: + data_frame = self._get_in_memory(category, operation, segment, activity) + return data_frame + + def _get_in_memory(self, + category:Category, + operation:Union[str,None]=None, + segment:Union[str,None]=None, + activity:Union[str,None]=None) -> Union[pd.DataFrame, None]: + key = (category, operation, segment, activity) + path = f'/RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' + # TODO: this ehavior is different than the base class - get with Paul Duda to see if it needs to be + try: + data_frame = io_manager[path] + except KeyError: + data_frame = pd.DataFrame() + return data_frame \ No newline at end of file From 3ade7e0537ddf48ccb90650cddff8a6d411e1c6d Mon Sep 17 00:00:00 2001 From: Burgholzer Date: Mon, 8 Sep 2025 14:48:56 -0400 Subject: [PATCH 5/5] added child class IOManagerDF for uniting under single interface --- src/hsp2/hsp2/main.py | 8 +++----- src/hsp2/hsp2io/io.py | 4 ++-- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/hsp2/hsp2/main.py b/src/hsp2/hsp2/main.py index 5564d2b4..04d78c56 100644 --- a/src/hsp2/hsp2/main.py +++ b/src/hsp2/hsp2/main.py @@ -47,7 +47,7 @@ versions, ) from hsp2.hsp2io.hdf import HDF5 -from hsp2.hsp2io.io import Category, IOManager, SupportsReadTS +from hsp2.hsp2io.io import Category, IOManager, IOManagerPandas, SupportsReadTS # special in-memory version of these from .HSP2utilitiesInMem import get_timeseries_InMem, save_timeseries_InMem @@ -78,6 +78,7 @@ def main( io_manager = IOManager(hdf5_instance) elif isinstance(io_manager, dict): hdfname = '' + io_manager = IOManagerPandas(io_manager) else: hdfname = io_manager._input.file_path if not os.path.exists(hdfname): @@ -87,10 +88,7 @@ def main( msg(1, f"Processing started for file {hdfname}; saveall={saveall}") # read user control, parameters, states, and flags uci and map to local variables - if isinstance(io_manager, dict): - uci_obj = convert_uci_InMem(io_manager) - else: - uci_obj = io_manager.read_uci() + uci_obj = io_manager.read_uci() opseq = uci_obj.opseq ddlinks = uci_obj.ddlinks ddmasslinks = uci_obj.ddmasslinks diff --git a/src/hsp2/hsp2io/io.py b/src/hsp2/hsp2io/io.py index 01ece829..f0466e8d 100644 --- a/src/hsp2/hsp2io/io.py +++ b/src/hsp2/hsp2io/io.py @@ -140,7 +140,7 @@ def _get_in_memory(self, except KeyError: return None -class IOManagerPandas(IOManager): +class IOManagerDF(IOManager): def __init__(self, io_combined: Union[SupportsReadUCI, SupportsReadTS, SupportsWriteTS, None] = None, uci: Union[SupportsReadUCI,None]=None, @@ -174,7 +174,7 @@ def _get_in_memory(self, segment:Union[str,None]=None, activity:Union[str,None]=None) -> Union[pd.DataFrame, None]: key = (category, operation, segment, activity) - path = f'/RESULTS/{x.SVOL}_{x.SVOLNO}/{sgrpn}' + path = f'/RESULTS/{operation}_{segment}/{activity}' # TODO: this ehavior is different than the base class - get with Paul Duda to see if it needs to be try: data_frame = io_manager[path]