diff --git a/src/troute-bmi/src/troute_nwm_bmi/troute_bmi.py b/src/troute-bmi/src/troute_nwm_bmi/troute_bmi.py index 08b1b423d..d46d39d1f 100644 --- a/src/troute-bmi/src/troute_nwm_bmi/troute_bmi.py +++ b/src/troute-bmi/src/troute_nwm_bmi/troute_bmi.py @@ -9,13 +9,30 @@ _VAR_NAME_UNITS_MAP = { 'land_surface_water_source__volume_flow_rate': ['streamflow_cms', 'm3 s-1'], + 'channel_exit_water_x-section__volume_flow_rate': ['streamflow_cms', 'm3 s-1'], + 'channel_water_flow__speed': ['streamflow_ms', 'm s-1'], + 'channel_water__mean_depth': ['streamflow_m', 'm'], + 'lake_water~incoming__volume_flow_rate': ['waterbody_cms', 'm3 s-1'], + 'lake_water~outgoing__volume_flow_rate': ['waterbody_cms', 'm3 s-1'], + 'lake_surface__elevation': ['waterbody_m', 'm'], } -_OUTPUT_VAR_NAMES = [] +_OUTPUT_VAR_NAMES = [ + "channel_water__id", + "channel_exit_water_x-section__volume_flow_rate", + "channel_water_flow__speed", + "channel_water__mean_depth", + "lake_water__id", + "lake_water~incoming__volume_flow_rate", + "lake_water~outgoing__volume_flow_rate", + "lake_surface__elevation" +] _INPUT_VAR_NAMES = [ "land_surface_water_source__id", + "land_surface_water_source__id" + _COUNT_SUFFIX, "land_surface_water_source__volume_flow_rate", + "land_surface_water_source__volume_flow_rate" + _COUNT_SUFFIX, "upstream_id", ] @@ -33,10 +50,11 @@ def __init__(self): "land_surface_water_source__volume_flow_rate" + _COUNT_SUFFIX: np.zeros(1, dtype=np.int64), "upstream_id": np.zeros(0, dtype=int), } + for var in _OUTPUT_VAR_NAMES: + self._values[var] = np.zeros(0, np.int32) # dtype subject to change after running model self._var_loc = "node" self._var_grid_id = 0 self._time_units = "s" - self._start_time = 0.0 self._end_time = np.finfo("d").max def initialize(self, bmi_cfg_file): @@ -95,17 +113,17 @@ def get_value_ptr(self, var_name: str): def get_start_time(self): """Start time of model.""" - return self._start_time + return 0.0 def get_end_time(self): """End time of model.""" - return self._end_time + return self._model.ngen_dt * (self._model.nts - 1) def get_current_time(self): return self._model.time def get_time_step(self): - return self._model.dt + return self._model.ngen_dt def get_time_units(self): return self._time_units @@ -113,7 +131,7 @@ def get_time_units(self): def finalize(self): """Finalize model.""" if self._model is not None: - self._model.run(self._values) + self._model.finalize(self._values) self._model = None # BMI functions that are not being used yet... @@ -124,11 +142,7 @@ def update_frac(self, time_frac: float): time_frac : float Fraction fo a time step. """ - time_step = self.get_time_step() - self._model.dt = int(time_frac * time_step) - if self._model.dt > 0: - self.update() - self._model.dt = time_step + raise NotImplementedError("update_frac") def get_var_type(self, var_name): """Data type of variable. diff --git a/src/troute-bmi/src/troute_nwm_bmi/troute_model.py b/src/troute-bmi/src/troute_nwm_bmi/troute_model.py index f61c05711..ded72c1c4 100644 --- a/src/troute-bmi/src/troute_nwm_bmi/troute_model.py +++ b/src/troute-bmi/src/troute_nwm_bmi/troute_model.py @@ -1,4 +1,5 @@ """Basic Model Interface backing model for NGEN t-route.""" +from __future__ import annotations import logging import time import yaml @@ -21,8 +22,6 @@ class Model: - dt: int - def __init__(self, config_file: str): self._main_start_time = time.time() self._time = 0.0 @@ -33,7 +32,18 @@ def __init__(self, config_file: str): log_level_set(self.log_parameters) - self.dt = int(self.forcing_parameters["dt"]) + output_type = self.stream_output.get("stream_output_type", None) + if output_type and (output_type != ".nc"): + error = "The stream output type can only be NetCDF. Current type: " + str(output_type) + LOG.error(error) + raise RuntimeError(error) + if self.compute_parameters.get("parallel_compute_method", "bmi") != "bmi": + LOG.warning( + "The parallel_compute_method is set to " + + '"' + str(self.compute_parameters.get("parallel_compute_method")) + '"' + + " in the configuration. To improve processing speed during catchment calculations," + + ' this will be set to "' + "bmi" '".' + ) LOG.info("Creating network of type " + self.supernetwork_parameters.get("network_type")) network_start_time = time.time() @@ -66,20 +76,20 @@ def __init__(self, config_file: str): ) else: raise Exception("Supernetwork network type must be HYFeaturesNetwork or NHDNetwork") + self.start_time = self._network.t0 + self._network.assemble_coastal_coupling_data() network_creation_time = time.time() - network_start_time # Data data assimilation LOG.debug("Creating DataAssimilation object") forcing_start_time = time.time() - da_run = {} + run_sets = self._create_run_sets() if self.data_assimilation_parameters: - run_set = { - "nts": self.nts, - "final_timestamp": self.t0 + timedelta(seconds=self.nts * self.dt) - } - da_sets = hnu.build_da_sets(self.data_assimilation_parameters, [run_set], self._network.t0) - if da_sets: - da_run = da_sets[0] + self._da_sets = hnu.build_da_sets(self.data_assimilation_parameters, run_sets, self._network.t0) + else: + self._da_sets = [{} for _ in run_sets] + self._da_index = 1 + self._data_assimilation = DataAssimilation( network=self._network, data_assimilation_parameters=self.data_assimilation_parameters, @@ -87,7 +97,7 @@ def __init__(self, config_file: str): waterbody_parameters=self.waterbody_parameters, from_files=True, value_dict=None, - da_run=da_run, + da_run=self._da_sets[0], ) forcing_time = time.time() - forcing_start_time @@ -111,14 +121,20 @@ def update(self, bmi_values: dict): step_time = self._network.t0 + timedelta(seconds=self.time) timestamp = step_time.strftime("%Y%m%d%H%M") self._df_data[timestamp] = np.array(qlat_values) - self._time += self.dt + self._time += self.ngen_dt self._timings["forcing_time"] += time.time() - start + if len(self._df_data) >= self.max_timestep_buffer: + self.run(bmi_values) + self._df_data.clear() def run(self, bmi_values: dict): - nts = self.nts - qts_subdivisions = self.forcing_parameters.get('qts_subdivisions', 12) - + nts = len(self._df_data) * self.qts_subdivisions + run_params = { + "t0": self._network.t0, + "dt": self.dt, + "nts": nts, + } LOG.debug("Assembling forcing dataframe") forcing_start_time = time.time() ## setup the qlats dataframe from the update() data @@ -136,6 +152,10 @@ def run(self, bmi_values: dict): else: flowveldepth_interorder = {} + usgs_df = self._data_assimilation.usgs_df + if not usgs_df.empty: + usgs_df = usgs_df.loc[:,self._network.t0:] + LOG.debug("Starting routing function") route_start_time = time.time() run_results, self._subnetwork = nwm_routing.nwm_route( @@ -143,19 +163,19 @@ def run(self, bmi_values: dict): upstream_connections=self._network.reverse_network, waterbodies_in_connections=self._network.waterbody_connections, reaches_bytw=self._network._reaches_by_tw, - parallel_compute_method=self.compute_parameters.get("parallel_compute_method", "serial"), + parallel_compute_method="bmi", compute_kernel=self.compute_parameters.get("compute_kernel"), subnetwork_target_size=self.compute_parameters.get('subnetwork_target_size'), cpu_pool=self.cpu_pool, - t0=self.t0, + t0=self._network.t0, dt=self.dt, nts=nts, - qts_subdivisions=qts_subdivisions, + qts_subdivisions=self.qts_subdivisions, independent_networks=self._network.independent_networks, param_df=self._network.dataframe, q0=self._network.q0, qlats=qlats, - usgs_df=self._data_assimilation.usgs_df, + usgs_df=usgs_df, lastobs_df=self._data_assimilation.lastobs_df, reservoir_usgs_df=self._data_assimilation.reservoir_usgs_df, reservoir_usgs_param_df=self._data_assimilation.reservoir_usgs_param_df, @@ -185,20 +205,18 @@ def run(self, bmi_values: dict): ) self._timings["route_time"] = time.time() - route_start_time - # create initial conditions for next loop iteration + # # create initial conditions for next loop iteration self._network.new_q0(run_results) self._network.update_waterbody_water_elevation() - # update reservoir parameters and lastobs_df + # # update reservoir parameters and lastobs_df self._data_assimilation.update_after_compute(run_results, self.dt * nts) LOG.debug("Generating output") output_start_time = time.time() - run_params = { - "t0": self.t0, - "dt": self.dt, - "nts": nts, - } + + # Note: After creating the output file the first run, all subsequent writes will append the results. + # TODO: Allow the results to be written into the middle of an existing file, e.g., NGEN loads a previous step and needs to re-write nwm_output_generator( run=run_params, results=run_results, @@ -207,7 +225,7 @@ def run(self, bmi_values: dict): parity_parameters=self.parity_parameters, restart_parameters=self.restart_parameters, parity_set={}, - qts_subdivisions=qts_subdivisions, + qts_subdivisions=self.qts_subdivisions, return_courant=self.compute_parameters.get("return_courant", False), cpu_pool=self.cpu_pool, waterbodies_df=self._network.waterbody_dataframe, @@ -219,9 +237,49 @@ def run(self, bmi_values: dict): link_lake_crosswalk=self._network.link_lake_crosswalk, nexus_dict=self._network.nexus_dict, poi_crosswalk=self._network.poi_nex_dict or {}, + filename_t0=self.start_time, ) + + self._network.new_t0(self.dt, nts) + if self._da_index < len(self._da_sets): + self._data_assimilation.update_for_next_loop( + self._network, + self._da_sets[self._da_index] + ) + self._da_index += 1 + + # compute BMI outputs + qvd_columns = pd.MultiIndex.from_product( + [range(nts), ["q", "v", "d"]] + ).to_flat_index() + flowveldepth = pd.concat( + [pd.DataFrame(r[1], index=r[0], columns=qvd_columns) for r in run_results], + copy=False, + ) + bmi_values["channel_exit_water_x-section__volume_flow_rate"] = flowveldepth.iloc[:,-3].to_numpy() + bmi_values["channel_water_flow__speed"] = flowveldepth.iloc[:,-2].to_numpy() + bmi_values["channel_water__mean_depth"] = flowveldepth.iloc[:,-1].to_numpy() + bmi_values["channel_water__id"] = flowveldepth.index.to_numpy() + + i_columns = pd.MultiIndex.from_product( + [range(int(nts)), ["i"]] + ).to_flat_index() + wbdy = pd.concat( + [pd.DataFrame(r[6], index=r[0], columns=i_columns) for r in run_results], + copy=False, + ) + + wbdy_id = self._network.waterbody_dataframe.index.values + bmi_values["lake_water__id"] = wbdy_id + bmi_values["lake_water~incoming__volume_flow_rate"] = wbdy.loc[wbdy_id].iloc[:,-1] + bmi_values["lake_water~outgoing__volume_flow_rate"] = flowveldepth.loc[wbdy_id].iloc[:,-3] + bmi_values["lake_surface__elevation"] = flowveldepth.loc[wbdy_id].iloc[:,-1] + self._timings["output_time"] = time.time() - output_start_time + def finalize(self, bmi_values: dict): + if len(self._df_data) > 0: + self.run(bmi_values) if self.show_timing: self._log_times() @@ -243,7 +301,7 @@ def log_parameters(self) -> dict: @property def compute_parameters(self) -> dict: - return self._config.get("compute_parameters", {}) + return self._config["compute_parameters"] @property def network_topology_parameters(self) -> dict: @@ -253,6 +311,10 @@ def network_topology_parameters(self) -> dict: def output_parameters(self) -> dict: return self._config.get("output_parameters", {}) + @property + def stream_output(self) -> dict: + return self.output_parameters.get("stream_output", {}) + @property def preprocessing_parameters(self) -> dict: return self.network_topology_parameters.get("preprocessing_parameters", {}) @@ -267,7 +329,7 @@ def supernetwork_parameters(self) -> dict: @property def forcing_parameters(self) -> dict: - return self.compute_parameters.get("forcing_parameters", {}) + return self.compute_parameters["forcing_parameters"] @property def restart_parameters(self) -> dict: @@ -285,6 +347,10 @@ def data_assimilation_parameters(self) -> dict: def parity_parameters(self) -> dict: return self.output_parameters.get("wrf_hydro_parity_check", {}) + @property + def max_timestep_buffer(self) -> int: + return self.bmi_parameters.get("max_timestep_buffer", 1000) + @property def show_timing(self): return bool(self.log_parameters.get("showtiming")) @@ -303,15 +369,42 @@ def time(self) -> float: return self._time @property - def t0(self) -> datetime: - return self._network.t0 + def qts_subdivisions(self) -> int: + return self.forcing_parameters["qts_subdivisions"] + + @property + def dt(self) -> int: + return self.forcing_parameters["dt"] + + @property + def ngen_dt(self) -> int: + return int(self.dt * self.qts_subdivisions) + + def _create_run_sets(self): + nts = self.nts + ngen_dt = self.ngen_dt + max_buffer = self.max_timestep_buffer + run_sets = [] + run_start = 0 + while run_start < nts: + run_nts = max_buffer + if run_start + run_nts > nts: + run_nts = nts - run_start + run_t0 = self.start_time + timedelta(seconds=run_start * ngen_dt) + run_end = run_t0 + timedelta(seconds=(run_nts - 1) * ngen_dt) + run_sets.append({ + "nts": run_nts * self.qts_subdivisions, + "final_timestamp": run_end + }) + run_start += max_buffer + return run_sets def _log_times(self): def sec_and_per(title, key: str): seconds = round(self._timings[key], 2) percent = round(self._timings[key] / process_time * 100, 2) LOG.info(f"{title}: {seconds} secs, {percent} %") - process_time = time.time() - self._main_start_time + process_time = sum(self._timings.values()) LOG.debug(f"Processes complete in {process_time} seconds.") LOG.info('************ TIMING SUMMARY ************') LOG.info('----------------------------------------') diff --git a/src/troute-config/troute/config/bmi_parameters.py b/src/troute-config/troute/config/bmi_parameters.py index b8a15a32e..e59255ae2 100644 --- a/src/troute-config/troute/config/bmi_parameters.py +++ b/src/troute-config/troute/config/bmi_parameters.py @@ -46,4 +46,8 @@ class BMIParameters(BaseModel): 'hydroseq', 'hl_uri', ] - ) \ No newline at end of file + ) + max_timestep_buffer: Optional[int] = 1000 + """ + The maximum number of timestamps the BMI can hold in memory before it forces a run, writes the outputs, and clears its memory. + """ diff --git a/src/troute-network/troute/AbstractNetwork.py b/src/troute-network/troute/AbstractNetwork.py index 830be578e..37785c7ac 100644 --- a/src/troute-network/troute/AbstractNetwork.py +++ b/src/troute-network/troute/AbstractNetwork.py @@ -135,7 +135,9 @@ def assemble_forcings(self, run,): "lateral inflow DataFrame creation complete in %s seconds." \ % (time.time() - start_time) ) + self.assemble_coastal_coupling_data() + def assemble_coastal_coupling_data(self): #--------------------------------------------------------------------- # Assemble coastal coupling data [WIP] #--------------------------------------------------------------------- diff --git a/src/troute-network/troute/nhd_io.py b/src/troute-network/troute/nhd_io.py index 57665dee1..808b220fe 100644 --- a/src/troute-network/troute/nhd_io.py +++ b/src/troute-network/troute/nhd_io.py @@ -2062,9 +2062,12 @@ def write_single_waterbody_netcdf( t0, dt, nts, - time_index + time_index, + filename_t0: datetime = None ): - netcdfname = 'troute_lakeout_' + t0.strftime('%Y%m%d%H%M') + '.nc' + if filename_t0 is None: + filename_t0 = t0 + netcdfname = 'troute_lakeout_' + filename_t0.strftime('%Y%m%d%H%M') + '.nc' # array of simulation time wbdy_time = np.array([t0 + timedelta(seconds = int(time_index[i] + 1) * dt) for i in range(len(time_index))]) #time_index @@ -2095,252 +2098,275 @@ def write_single_waterbody_netcdf( # open netCDF4 Dataset in write mode with netCDF4.Dataset( filename = str(wbdy_filepath) + '/' + netcdfname, - mode = 'w', + mode = 'a', format = "NETCDF4" ) as f: - # =========== DIMENSIONS =============== - _ = f.createDimension("time", num_timesteps) - _ = f.createDimension("feature_id", len(wbdy_feature_id)) - _ = f.createDimension("reference_time", 1) - _ = f.createDimension("latitude", len(wbdy_feature_id)) - _ = f.createDimension("longitude", len(wbdy_feature_id)) + # Get the current length of the time dimension. + # If this is 0, assume that the file was just created and variable need to be created within. + # If this is greater than zero, assume we're appending and use the time index for where to start new writes + time_len = len(f.dimensions["time"]) if "time" in f.dimensions else 0 # =========== time VARIABLE =============== - TIME = f.createVariable( - varname = "time", - datatype = 'int32', - dimensions = ("time",), - ) - TIME[:] = date2num( + if "time" in f.variables: + TIME = f.variables["time"] + else: + _ = f.createDimension("time", None) + TIME = f.createVariable( + varname = "time", + datatype = 'int32', + dimensions = ("time",), + ) + f['time'].setncatts( + { + 'long_name': 'valid output time', + 'standard_name': 'time', + 'valid_min': date2num( + t0, + units = "minutes since 1970-01-01 00:00:00 UTC", + calendar = "gregorian" + ), + 'valid_max': date2num( + wbdy_time[num_timesteps-1], + units = "minutes since 1970-01-01 00:00:00 UTC", + calendar = "gregorian" + ), + 'units': 'minutes since 1970-01-01T00:00:00+00:00', + 'calendar': 'proleptic_gregorian' + } + ) + TIME[time_len:] = date2num( wbdy_time, units = "minutes since 1970-01-01 00:00:00 UTC", calendar = "gregorian" ) - f['time'].setncatts( - { - 'long_name': 'valid output time', - 'standard_name': 'time', - 'valid_min': date2num( - t0, - units = "minutes since 1970-01-01 00:00:00 UTC", - calendar = "gregorian" - ), - 'valid_max': date2num( - wbdy_time[num_timesteps-1], - units = "minutes since 1970-01-01 00:00:00 UTC", - calendar = "gregorian" - ), - 'units': 'minutes since 1970-01-01T00:00:00+00:00', - 'calendar': 'proleptic_gregorian' - } - ) # =========== reference_time VARIABLE =============== - REF_TIME = f.createVariable( - varname = "reference_time", - datatype = 'int32', - dimensions = ("reference_time",), - ) - REF_TIME[:] = date2num( - t0, - units = "minutes since 1970-01-01 00:00:00 UTC", - calendar = "gregorian" - ) - f['reference_time'].setncatts( - { - 'long_name': 'model initialization time', - 'standard_name': 'forecast_reference_time', - 'units': 'minutes since 1970-01-01T00:00:00+00:00', - 'calendar': 'proleptic_gregorian' - } - ) + if "reference_time" not in f.variables: + _ = f.createDimension("reference_time", 1) + REF_TIME = f.createVariable( + varname = "reference_time", + datatype = 'int32', + dimensions = ("reference_time",), + ) + REF_TIME[:] = date2num( + t0, + units = "minutes since 1970-01-01 00:00:00 UTC", + calendar = "gregorian" + ) + f['reference_time'].setncatts( + { + 'long_name': 'model initialization time', + 'standard_name': 'forecast_reference_time', + 'units': 'minutes since 1970-01-01T00:00:00+00:00', + 'calendar': 'proleptic_gregorian' + } + ) # =========== feature_id VARIABLE =============== - FEATURE_ID = f.createVariable( - varname = "feature_id", - datatype = 'int64', - dimensions = ("feature_id",), - ) - FEATURE_ID[:] = wbdy_feature_id - f['feature_id'].setncatts( - { - 'long_name': 'Lake ComID', - 'comment': '', - 'cf_role:': 'timeseries_id' - } - ) + if "feature_id" not in f.variables: + _ = f.createDimension("feature_id", len(wbdy_feature_id)) + FEATURE_ID = f.createVariable( + varname = "feature_id", + datatype = 'int64', + dimensions = ("feature_id",), + ) + FEATURE_ID[:] = wbdy_feature_id + f['feature_id'].setncatts( + { + 'long_name': 'Lake ComID', + 'comment': '', + 'cf_role:': 'timeseries_id' + } + ) # =========== latitude VARIABLE =============== - LATITUDE = f.createVariable( - varname = "latitude", - datatype = 'f4', - dimensions = ("latitude",), - fill_value = np.nan - ) - if not waterbodies_df.empty: - LATITUDE[:] = waterbodies_df.lat.tolist() - else: - LATITUDE[:] = np.nan - f['latitude'].setncatts( - { - 'long_name': 'Lake latitude', - 'standard_name': 'latitude', - 'units': 'degrees_north' - } - ) - + if "latitude" not in f.variables: + _ = f.createDimension("latitude", len(wbdy_feature_id)) + LATITUDE = f.createVariable( + varname = "latitude", + datatype = 'f4', + dimensions = ("latitude",), + fill_value = np.nan + ) + if not waterbodies_df.empty: + LATITUDE[:] = waterbodies_df.lat.tolist() + else: + LATITUDE[:] = np.nan + f['latitude'].setncatts( + { + 'long_name': 'Lake latitude', + 'standard_name': 'latitude', + 'units': 'degrees_north' + } + ) + + # =========== longitude VARIABLE =============== - LONGITUDE = f.createVariable( - varname = "longitude", - datatype = 'f4', - dimensions = ("longitude",), - fill_value = np.nan - ) - if not waterbodies_df.empty: - LONGITUDE[:] = waterbodies_df.lon.tolist() - else: - LONGITUDE[:] = np.nan - f['longitude'].setncatts( - { - 'long_name': 'Lake longitude', - 'standard_name': 'longitude', - 'units': 'degrees_east' - } - ) - - # =========== reservoir type VARIABLE =============== - res_type = f.createVariable( - varname = "reservoir_type", - datatype = "i4", - dimensions = ("feature_id") - ) - res_type[:] = wbdy_type - f['reservoir_type'].setncatts( - { - 'coordinates': 'latitude longitude', - 'long_name': 'reservoir_type', - 'flag_values': [1, 2, 3, 4], - 'flag_meanings': 'Level_pool USGS-persistence USACE-persistence RFC-forecasts' - } - ) - - # =========== crs VARIABLE =============== - crs = f.createVariable( - varname = "crs", - datatype = "S1", - dimensions = () + if "longitude" not in f.variables: + _ = f.createDimension("longitude", len(wbdy_feature_id)) + LONGITUDE = f.createVariable( + varname = "longitude", + datatype = 'f4', + dimensions = ("longitude",), + fill_value = np.nan + ) + if not waterbodies_df.empty: + LONGITUDE[:] = waterbodies_df.lon.tolist() + else: + LONGITUDE[:] = np.nan + f['longitude'].setncatts( + { + 'long_name': 'Lake longitude', + 'standard_name': 'longitude', + 'units': 'degrees_east' + } + ) + + # =========== reservoir type VARIABLE =============== + if "reservoir_type" not in f.variables: + res_type = f.createVariable( + varname = "reservoir_type", + datatype = "i4", + dimensions = ("feature_id") + ) + res_type[:] = wbdy_type + f['reservoir_type'].setncatts( + { + 'coordinates': 'latitude longitude', + 'long_name': 'reservoir_type', + 'flag_values': [1, 2, 3, 4], + 'flag_meanings': 'Level_pool USGS-persistence USACE-persistence RFC-forecasts' + } + ) + + # =========== crs VARIABLE =============== + if crs not in f.variables: + crs = f.createVariable( + varname = "crs", + datatype = "S1", + dimensions = () + ) + if not waterbodies_df.empty: + crs[:] = np.array(waterbodies_df['crs'].iat[0],dtype = '|S1') + else: + crs[:] = np.array(['n'],dtype = '|S1') + f['crs'].setncatts( + { + 'transform_name': 'latitude longitude', + 'grid_mapping_name': 'latitude longitude', + 'esri_pe_string': 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision', + 'spatial_ref': 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision', + 'long_name': 'CRS definition', + 'longitude_of_prime_meridian': 0.0, + '_CoordinateAxes': 'latitude longitude', + 'semi_major_axis': 6378137.0, + 'semi_minor_axis': 6356752.5, + 'inverse_flattening': 298.25723 + } ) - if not waterbodies_df.empty: - crs[:] = np.array(waterbodies_df['crs'].iat[0],dtype = '|S1') - else: - crs[:] = np.array(['n'],dtype = '|S1') - f['crs'].setncatts( - { - 'transform_name': 'latitude longitude', - 'grid_mapping_name': 'latitude longitude', - 'esri_pe_string': 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision', - 'spatial_ref': 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision', - 'long_name': 'CRS definition', - 'longitude_of_prime_meridian': 0.0, - '_CoordinateAxes': 'latitude longitude', - 'semi_major_axis': 6378137.0, - 'semi_minor_axis': 6356752.5, - 'inverse_flattening': 298.25723 - } - ) - # =========== inflow VARIABLE =============== - inflow = f.createVariable( - varname = "inflow", - datatype = "f4", - dimensions = ("time","latitude","longitude","feature_id"), - fill_value = -999900 - ) - inflow[:] = max_i_vals - f['inflow'].setncatts( - { - 'long_name': 'Lake Inflow', - 'units': 'm3 s-1', - 'grid_mapping': 'crs', - 'valid_range': [-1000000,1000000], - 'coordinates': 'latitude longitude', - 'add_offset': 0.0, - 'scale_factor': 0.01, - 'missing_value': -999900 - } - ) + # =========== inflow VARIABLE =============== + if "inflow" in f.variables: + inflow = f.variables["inflow"] + else: + inflow = f.createVariable( + varname = "inflow", + datatype = "f4", + dimensions = ("time","latitude","longitude","feature_id"), + fill_value = -999900 + ) + f['inflow'].setncatts( + { + 'long_name': 'Lake Inflow', + 'units': 'm3 s-1', + 'grid_mapping': 'crs', + 'valid_range': [-1000000,1000000], + 'coordinates': 'latitude longitude', + 'add_offset': 0.0, + 'scale_factor': 0.01, + 'missing_value': -999900 + } + ) + inflow[time_len:, :, :, :] = max_i_vals - # =========== outflow VARIABLE =============== - outflow = f.createVariable( + # =========== outflow VARIABLE =============== + if "outflow" in f.variables: + outflow = f.variables["outflow"] + else: + outflow = f.createVariable( varname = "outflow", datatype = "f4", dimensions = ("time","latitude","longitude","feature_id"), fill_value = -999900 #np.nan ) - outflow[:] = max_q_vals - f['outflow'].setncatts( - { - 'long_name': 'Lake Outflow', - 'units': 'm3 s-1', - 'grid_mapping': 'crs', - 'valid_range': [-1000000,1000000], - 'coordinates': 'latitude longitude', - 'add_offset': 0.0, - 'scale_factor': 0.01, - 'missing_value': -999900 - } - ) + f['outflow'].setncatts( + { + 'long_name': 'Lake Outflow', + 'units': 'm3 s-1', + 'grid_mapping': 'crs', + 'valid_range': [-1000000,1000000], + 'coordinates': 'latitude longitude', + 'add_offset': 0.0, + 'scale_factor': 0.01, + 'missing_value': -999900 + } + ) + outflow[time_len:, :, :, :] = max_q_vals - # =========== depth VARIABLE =============== - depth = f.createVariable( + # =========== depth VARIABLE =============== + if "water_sfc_elev" in f.variables: + depth = f.variables["water_sfc_elev"] + else: + depth = f.createVariable( varname = "water_sfc_elev", datatype = "f4", dimensions = ("time","latitude","longitude","feature_id"), fill_value = np.nan ) - depth[:] = max_d_vals - f['water_sfc_elev'].setncatts( - { - 'long_name': 'Water Surface Elevation', - 'units': 'm', - 'comment': 'If reservoir_type = 4, water_sfc_elev is invalid because this value corresponds only to level pool', - 'coordinates': 'latitude longitude' - } - ) - - # =========== GLOBAL ATTRIBUTES =============== - if not waterbodies_df.empty: - f.setncatts( - { - 'TITLE': 'OUTPUT FROM T-ROUTE', - 'featureType': 'timeSeries', - 'proj4': '+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@', - 'model_initialization_time': t0.strftime('%Y-%m-%d_%H:%M:%S'), - 'station_dimension': 'lake_id', - 'model_total_valid_times': nts, - 'Conventions': 'CF-1.6', - 'code_version': '', - 'model_output_type': 'reservoir', - 'model_configuration': '' - } - ) - else: - f.setncatts( + f['water_sfc_elev'].setncatts( { - 'TITLE': 'OUTPUT FROM T-ROUTE', - 'basin_characteristics': 'There are no waterbodies in this basin', - 'featureType': 'timeSeries', - 'proj4': '+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@', - 'model_initialization_time': t0.strftime('%Y-%m-%d_%H:%M:%S'), - 'station_dimension': 'lake_id', - 'model_total_valid_times': nts, - 'Conventions': 'CF-1.6', - 'code_version': '', - 'model_output_type': 'reservoir', - 'model_configuration': '' + 'long_name': 'Water Surface Elevation', + 'units': 'm', + 'comment': 'If reservoir_type = 4, water_sfc_elev is invalid because this value corresponds only to level pool', + 'coordinates': 'latitude longitude' } ) + depth[time_len:, :, :, :] = max_d_vals + + # =========== GLOBAL ATTRIBUTES =============== + if time_len == 0: + if not waterbodies_df.empty: + f.setncatts( + { + 'TITLE': 'OUTPUT FROM T-ROUTE', + 'featureType': 'timeSeries', + 'proj4': '+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@', + 'model_initialization_time': t0.strftime('%Y-%m-%d_%H:%M:%S'), + 'station_dimension': 'lake_id', + 'model_total_valid_times': nts, + 'Conventions': 'CF-1.6', + 'code_version': '', + 'model_output_type': 'reservoir', + 'model_configuration': '' + } + ) + else: + f.setncatts( + { + 'TITLE': 'OUTPUT FROM T-ROUTE', + 'basin_characteristics': 'There are no waterbodies in this basin', + 'featureType': 'timeSeries', + 'proj4': '+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@', + 'model_initialization_time': t0.strftime('%Y-%m-%d_%H:%M:%S'), + 'station_dimension': 'lake_id', + 'model_total_valid_times': nts, + 'Conventions': 'CF-1.6', + 'code_version': '', + 'model_output_type': 'reservoir', + 'model_configuration': '' + } + ) def write_flowveldepth_csv_pkl(stream_output_directory, file_name, flow, velocity, depth, nudge_df, timestamps, @@ -2381,148 +2407,157 @@ def write_flowveldepth_netcdf(stream_output_directory, file_name, # Open netCDF4 Dataset in write mode with netCDF4.Dataset( filename=f"{stream_output_directory}/{file_name}", - mode='w', + mode='a', # append mode (creates if it doesn't exist) format='NETCDF4' ) as ncfile: + # assumes the feature IDs will be the same between instances # ============ DIMENSIONS =================== - _ = ncfile.createDimension('feature_id', len(flow)) - _ = ncfile.createDimension('time', len(timestamps)) max_str_len = max(len(str(x)) for x in flow.index.get_level_values('Type')) - _ = ncfile.createDimension('type_strlen', max_str_len) - #_ = ncfile.createDimension('gage', gage) - #_ = ncfile.createDimension('nudge_timestep', nudge_timesteps) # Add dimension for nudge time steps - - # =========== time VARIABLE =============== - TIME = ncfile.createVariable( - varname = "time", - datatype = 'float64', - dimensions = ("time",), - fill_value = -9999.0 - ) - TIME[:] = timestamps - ncfile['time'].setncatts( - { - 'long_name': 'valid output time', - 'standard_name': 'time', - 'units': f'seconds since {t0.strftime("%Y-%m-%d %H:%M:%S")}', - 'missing_value': -9999.0 - #'calendar': 'proleptic_gregorian' - } - ) - - # # =========== reference_time VARIABLE =============== - # REF_TIME = ncfile.createVariable( - # varname = "reference_time", - # datatype = 'S1', - # dimensions = ("reference_time",), - # ) - # REF_TIME[:] = t0.strftime('%Y-%m-%d_%H:%M:%S') - # ncfile['reference_time'].setncatts( - # { - # 'long_name': 'reference time', - # 'standard_name': 'reference_time', - # } - # ) + # Get the current length of the time dimension. + # If this is 0, assume that the file was just created and variable need to be created within. + # If this is greater than zero, assume we're appending and use the time index for where to start new writes + time_len = len(ncfile.dimensions["time"]) if "time" in ncfile.dimensions else 0 # =========== feature_id VARIABLE =============== - FEATURE_ID = ncfile.createVariable( - varname = "feature_id", - datatype = 'int64', - dimensions = ("feature_id",), - ) - FEATURE_ID[:] = flow.index.get_level_values('featureID') - ncfile['feature_id'].setncatts( - { - 'long_name': 'Segment ID', - } - ) + if "feature_id" not in ncfile.variables: + _ = ncfile.createDimension('type_strlen', max_str_len) + _ = ncfile.createDimension('feature_id', len(flow)) + FEATURE_ID = ncfile.createVariable( + varname = "feature_id", + datatype = 'int64', + dimensions = ("feature_id",), + ) + ncfile['feature_id'].setncatts( + { + 'long_name': 'Segment ID', + } + ) + FEATURE_ID[:] = flow.index.get_level_values('featureID') + + # =========== time VARIABLE =============== + if "time" in ncfile.variables: + TIME = ncfile.variables["time"] + else: + _ = ncfile.createDimension('time', None) + TIME = ncfile.createVariable( + varname = "time", + datatype = 'float64', + dimensions = ("time",), + fill_value = -9999.0 + ) + ncfile['time'].setncatts( + { + 'long_name': 'valid output time', + 'standard_name': 'time', + 'units': f'seconds since {t0.strftime("%Y-%m-%d %H:%M:%S")}', + 'missing_value': -9999.0 + #'calendar': 'proleptic_gregorian' + } + ) + TIME[time_len:] = timestamps + # =========== type VARIABLE =============== - TYPE = ncfile.createVariable( - varname="type", - datatype=f'S{max_str_len}', - dimensions=("feature_id",), - ) - TYPE[:] = np.array(flow.index.get_level_values('Type').astype(str).tolist(), dtype=f'S{max_str_len}') - ncfile['type'].setncatts( - { - 'long_name': 'Type', - } - ) + if "type" not in ncfile.variables: + TYPE = ncfile.createVariable( + varname="type", + datatype=f'S{max_str_len}', + dimensions=("feature_id",), + ) + ncfile['type'].setncatts( + { + 'long_name': 'Type', + } + ) + TYPE[:] = np.array(flow.index.get_level_values('Type').astype(str).tolist(), dtype=f'S{max_str_len}') - # =========== flow VARIABLE =============== - flow_var = ncfile.createVariable( - varname = "flow", - datatype = "f4", - dimensions = ("feature_id", "time"), - fill_value = -9999.0 + # =========== flow VARIABLE =============== + if "flow" in ncfile.variables: + flow_var = ncfile.variables["flow"] + else: + flow_var = ncfile.createVariable( + varname = "flow", + datatype = "f4", + dimensions = ("feature_id", "time"), + fill_value = -9999.0 ) + ncfile['flow'].setncatts( + { + 'long_name': 'Flow', + 'units': 'm3 s-1', + 'missing_value': -9999.0 + } + ) + flow_var[:, time_len:] = flow.to_numpy(dtype=np.float32) - flow_var[:] = flow.to_numpy(dtype=np.float32) - ncfile['flow'].setncatts( - { - 'long_name': 'Flow', - 'units': 'm3 s-1', - 'missing_value': -9999.0 - } - ) + # =========== velocity VARIABLE =============== + if "velocity" in ncfile.variables: + velocity_var = ncfile.variables["velocity"] + else: + velocity_var = ncfile.createVariable( + varname = "velocity", + datatype = "f4", + dimensions = ("feature_id", "time"), + fill_value = -9999.0 + ) + ncfile['velocity'].setncatts( + { + 'long_name': 'Velocity', + 'units': 'm/s', + 'missing_value': -9999.0 + } + ) + velocity_var[:, time_len:] = velocity.to_numpy(dtype=np.float32) - # =========== velocity VARIABLE =============== - velocity_var = ncfile.createVariable( - varname = "velocity", - datatype = "f4", - dimensions = ("feature_id", "time"), - fill_value = -9999.0 - ) - velocity_var[:] = velocity.to_numpy(dtype=np.float32) - ncfile['velocity'].setncatts( - { - 'long_name': 'Velocity', - 'units': 'm/s', - 'missing_value': -9999.0 - } - ) + # =========== depth VARIABLE =============== + if "depth" in ncfile.variables: + depth_var = ncfile.variables["depth"] + else: + depth_var = ncfile.createVariable( + varname = "depth", + datatype = "f4", + dimensions = ("feature_id", "time"), + fill_value = -9999.0 + ) + ncfile['depth'].setncatts( + { + 'long_name': 'Depth', + 'units': 'm', + 'missing_value': -9999.0 + } + ) + depth_var[:, time_len:] = depth.to_numpy(dtype=np.float32) - # =========== depth VARIABLE =============== - depth_var = ncfile.createVariable( - varname = "depth", - datatype = "f4", - dimensions = ("feature_id", "time"), - fill_value = -9999.0 - ) - depth_var[:] = depth.to_numpy(dtype=np.float32) - ncfile['depth'].setncatts( - { - 'long_name': 'Depth', - 'units': 'm', - 'missing_value': -9999.0 - } - ) - - # =========== nudge VARIABLE =============== - nudge = ncfile.createVariable( - varname = "nudge", - datatype = "f4", - dimensions = ("feature_id", "time"), - fill_value = -9999.0 - ) - nudge[:] = nudge_df.to_numpy(dtype=np.float32) - ncfile['nudge'].setncatts( - { - 'long_name': 'Streamflow Nudge Value', - 'units': 'm3 s-1', - 'missing_value': -9999.0 - } - ) + # =========== nudge VARIABLE =============== + if "nudge" in ncfile.variables: + nudge = ncfile.variables["nudge"] + else: + nudge = ncfile.createVariable( + varname = "nudge", + datatype = "f4", + dimensions = ("feature_id", "time"), + fill_value = -9999.0 + ) + ncfile['nudge'].setncatts( + { + 'long_name': 'Streamflow Nudge Value', + 'units': 'm3 s-1', + 'missing_value': -9999.0 + } + ) + nudge[:, time_len:] = nudge_df.to_numpy(dtype=np.float32) # =========== GLOBAL ATTRIBUTES =============== - ncfile.setncatts( - { - 'TITLE': 'OUTPUT FROM T-ROUTE', - 'file_reference_time': t0.strftime('%Y-%m-%d_%H:%M:%S'), - 'code_version': '', - } - ) + if time_len == 0: + ncfile.setncatts( + { + 'TITLE': 'OUTPUT FROM T-ROUTE', + 'file_reference_time': t0.strftime('%Y-%m-%d_%H:%M:%S'), + 'code_version': '', + } + ) + + def stream_output_mask_reader(stream_output_mask): if not stream_output_mask: return {} @@ -2648,7 +2683,8 @@ def write_flowveldepth( cpu_pool = 1, poi_crosswalk = None, nexus_dict= None, - ): + file_name_time=None, +): ''' Write the results of flowveldepth and nudge to netcdf- break. Arguments @@ -2667,8 +2703,8 @@ def write_flowveldepth( # timesteps = list(timesteps) n_timesteps = flowveldepth.shape[1]//3 ts = stream_output_internal_frequency//(dt//60) - ind = [i for i in range(ts-1,n_timesteps,ts)] - timestamps_sec = [(i+1)*dt for i in ind] + ind = [i for i in range(ts-1, n_timesteps, ts)] + timestamps_sec = [(i+1) * dt for i in ind] flow = flowveldepth.iloc[:,0::3].iloc[:,ind] velocity = flowveldepth.iloc[:,1::3].iloc[:,ind] @@ -2683,7 +2719,8 @@ def write_flowveldepth( pd.set_option('future.no_silent_downcasting', True) empty_df = pd.DataFrame(index=empty_ids, columns=nudge_df.columns).fillna(-9999.0) nudge_df = pd.concat([nudge_df, empty_df]).loc[flowveldepth.index] - file_name_time = t0 + if file_name_time is None: + file_name_time = t0 jobs = [] if stream_output_timediff > 0: @@ -2695,12 +2732,15 @@ def write_flowveldepth( for _ in range(num_files): filename = 'troute_output_' + file_name_time.strftime('%Y%m%d%H%M') + stream_output_type - args = (stream_output_directory,filename, - flow.iloc[:,0:ts_per_file], - velocity.iloc[:,0:ts_per_file], - depth.iloc[:,0:ts_per_file], - nudge_df.iloc[:,0:ts_per_file], - timestamps_sec[0:ts_per_file],t0) + args = ( + stream_output_directory,filename, + flow.iloc[:,0:ts_per_file], + velocity.iloc[:,0:ts_per_file], + depth.iloc[:,0:ts_per_file], + nudge_df.iloc[:,0:ts_per_file], + timestamps_sec[0:ts_per_file], + t0 + ) if stream_output_type == '.nc': if cpu_pool > 1 & num_files > 1: jobs.append(delayed(write_flowveldepth_netcdf)(*args)) diff --git a/src/troute-nwm/src/nwm_routing/output.py b/src/troute-nwm/src/nwm_routing/output.py index c9e59b9c7..6ef7eae26 100644 --- a/src/troute-nwm/src/nwm_routing/output.py +++ b/src/troute-nwm/src/nwm_routing/output.py @@ -133,7 +133,8 @@ def nwm_output_generator( link_lake_crosswalk = None, nexus_dict = None, poi_crosswalk = None, - logFileName='NONE' + logFileName='NONE', + filename_t0=None, ): dt = run.get("dt") @@ -302,7 +303,8 @@ def nwm_output_generator( cpu_pool = cpu_pool, poi_crosswalk = poi_crosswalk, nexus_dict= nexus_dict, - ) + file_name_time=filename_t0, + ) nTimeBins = int( len(flowveldepth.columns)/3) fCalc = int(dt/60) @@ -353,7 +355,8 @@ def nwm_output_generator( t0, dt, nts, - time_index + time_index, + filename_t0=filename_t0, ) elif waterbodies_df.empty: #there are no waterbodies in the basin. #this portion of the code is used to create lakeout nc files with zero values. @@ -373,7 +376,8 @@ def nwm_output_generator( t0, dt, nts, - time_index + time_index, + filename_t0=filename_t0, ) LOG.warning("lakeout_output specified with no waterbodies in the basin.") end_time = time.time() diff --git a/src/troute-routing/troute/routing/compute.py b/src/troute-routing/troute/routing/compute.py index 55f529596..3a9caf86d 100644 --- a/src/troute-routing/troute/routing/compute.py +++ b/src/troute-routing/troute/routing/compute.py @@ -1673,12 +1673,31 @@ def compute_nhd_routing_v02( reservoir_usace_prev_persisted_flow, reservoir_usace_persistence_update_time, reservoir_usace_persistence_index, + reservoir_rfc_df_sub, + reservoir_rfc_totalCounts, + reservoir_rfc_file, + reservoir_rfc_use_forecast, + reservoir_rfc_timeseries_idx, + reservoir_rfc_update_time, + reservoir_rfc_da_timestep, + reservoir_rfc_persist_days, + gl_df_sub, + gl_parm_lake_id_sub, + gl_param_flows_sub, + gl_param_time_sub, + gl_param_update_time_sub, + gl_climatology_df_sub, waterbody_types_df_sub, - ) = _prep_reservoir_da_dataframes( + ) = _prep_reservoir_da_dataframes( reservoir_usgs_df, reservoir_usgs_param_df, reservoir_usace_df, reservoir_usace_param_df, + reservoir_rfc_df, + reservoir_rfc_param_df, + great_lakes_df, + great_lakes_param_df, + great_lakes_climatology_df, waterbody_types_df_sub, t0, from_files, @@ -1725,6 +1744,25 @@ def compute_nhd_routing_v02( reservoir_usace_prev_persisted_flow.astype("float32"), reservoir_usace_persistence_update_time.astype("float32"), reservoir_usace_persistence_index.astype("float32"), + # RFC Reservoir DA data + reservoir_rfc_df_sub.values.astype("float32"), + reservoir_rfc_df_sub.index.values.astype("int32"), + reservoir_rfc_totalCounts.astype("int32"), + reservoir_rfc_file, + reservoir_rfc_use_forecast.astype("int32"), + reservoir_rfc_timeseries_idx.astype("int32"), + reservoir_rfc_update_time.astype("float32"), + reservoir_rfc_da_timestep.astype("int32"), + reservoir_rfc_persist_days.astype("int32"), + # Great Lakes DA data + gl_df_sub.lake_id.values.astype("int32"), + gl_df_sub.time.values.astype("int32"), + gl_df_sub.Discharge.values.astype("float32"), + gl_parm_lake_id_sub.astype("int32"), + gl_param_flows_sub.astype("float32"), + gl_param_time_sub.astype("int32"), + gl_param_update_time_sub.astype("int32"), + gl_climatology_df_sub.values.astype("float32"), { us: fvd for us, fvd in flowveldepth_interorder.items() @@ -1732,6 +1770,7 @@ def compute_nhd_routing_v02( }, assume_short_ts, return_courant, + from_files=from_files, ) )