From 1967697313960bb94027effb19e0d397714a2e36 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:02:10 +0200 Subject: [PATCH 1/9] add method to get time coordinate from gwf/gwt --- nlmod/dims/time.py | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py index 0d693844..0bc9a1af 100644 --- a/nlmod/dims/time.py +++ b/nlmod/dims/time.py @@ -9,6 +9,7 @@ import numpy as np import pandas as pd +from xarray import IndexVariable logger = logging.getLogger(__name__) @@ -145,7 +146,7 @@ def estimate_nstp( ---------- forcing : array-like Array with a forcing value for each stress period. Forcing can be - for example a pumping rate of a rainfall intensity. + for example a pumping rate or a rainfall intensity. perlen : float or array of floats (nper) An array of the stress period lengths. tsmult : float or array of floats (nper) @@ -216,3 +217,36 @@ def estimate_nstp( else: return nstp_ceiled + + +def ds_time_from_model(gwf): + """Get time index variable from model (gwf or gwt). + + Parameters + ---------- + gwf : flopy MFModel object + groundwater flow or groundwater transport model + + Returns + ------- + IndexVariable + time coordinate for xarray data-array or dataset + """ + + start_datetime = gwf.simulation.get_package("TDIS").start_date_time.data + + if start_datetime is not None: + dt = pd.to_timedelta( + np.cumsum(gwf.dimensions.simulation_time.get_perioddata()["perlen"]), "D" + ) + times = pd.Timestamp(start_datetime) + dt + + else: + times = np.cumsum(gwf.dimensions.simulation_time.get_perioddata()["perlen"]) + + time = IndexVariable(["time"], times) + time.attrs["time_units"] = gwf.dimensions.simulation_time.get_time_units() + if start_datetime is not None: + time.attrs["start"] = str(start_datetime) + + return time From 1c7be698f2eb9c97b1841d60b438f232517c1033 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:05:18 +0200 Subject: [PATCH 2/9] add budget output reader function --- nlmod/gwf/output.py | 89 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py index 689f9237..6ee0e74e 100644 --- a/nlmod/gwf/output.py +++ b/nlmod/gwf/output.py @@ -113,6 +113,95 @@ def get_heads_da(ds=None, gwf=None, fname_hds=None): return head_ar +def get_budget_da(text, ds=None, gwf=None, fname_cbc=None, kstpkper=None): + """Reads budget file given either a dataset or a groundwater flow object. + + Parameters + ---------- + text : str + record to get from budget file + ds : xarray.Dataset, optional + xarray dataset with model data. One of ds or gwf must be provided. + gwf : flopy ModflowGwf, optional + Flopy groundwaterflow object. One of ds or gwf must be provided. + fname_cbc : path, optional + specify the budget file to load, if not provided budget file will + be obtained from ds or gwf. + + Returns + ------- + q_ar : xarray.DataArray + budget data array. + """ + cbfobj = _get_cbf(ds=ds, gwf=gwf, fname_cbc=fname_cbc) + + q = cbfobj.get_data(text=text, kstpkper=kstpkper, full3D=True) + q = np.stack(q) + + if gwf is not None: + gridtype = gwf.modelgrid.grid_type + else: + gridtype = ds.gridtype + + if gridtype == "vertex": + q_ar = xr.DataArray( + data=q, + dims=("time", "layer", "icell2d"), + coords={}, + attrs={"units": "m3/d"}, + ) + + elif gridtype == "structured": + if gwf is not None: + delr = np.unique(gwf.modelgrid.delr).item() + delc = np.unique(gwf.modelgrid.delc).item() + extent = gwf.modelgrid.extent + x, y = get_xy_mid_structured(extent, delr, delc) + + else: + x = ds.x + y = ds.y + + q_ar = xr.DataArray( + data=q, + dims=("time", "layer", "y", "x"), + coords={ + "x": x, + "y": y, + }, + attrs={"units": "m3/d"}, + ) + else: + assert 0, "Gridtype not supported" + + # set layer and time coordinates + if gwf is not None: + q_ar.coords["layer"] = np.arange(gwf.modelgrid.nlay) + q_ar.coords["time"] = ds_time_from_model(gwf) + else: + q_ar.coords["layer"] = ds.layer + q_ar.coords["time"] = ds.time + + if ds is not None and "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0: + affine = get_affine(ds) + q_ar.rio.write_transform(affine, inplace=True) + + elif gwf is not None and gwf.modelgrid.angrot != 0.0: + attrs = dict( + delr=np.unique(gwf.modelgrid.delr).item(), + delc=np.unique(gwf.modelgrid.delc).item(), + xorigin=gwf.modelgrid.xoffset, + yorigin=gwf.modelgrid.yoffset, + angrot=gwf.modelgrid.angrot, + extent=gwf.modelgrid.extent, + ) + affine = get_affine(attrs) + q_ar.rio.write_transform(affine, inplace=True) + + q_ar.rio.write_crs("EPSG:28992", inplace=True) + + return q_ar + def _get_hds(ds=None, gwf=None, fname_hds=None): msg = "Load the heads using either the ds, gwf or fname_hds" From eb6bfaa7681b9c9b8ecd2c6f286e0addf97ff42a Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:05:35 +0200 Subject: [PATCH 3/9] rename get cellbudget file method --- nlmod/gwf/output.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py index 6ee0e74e..a29d5d32 100644 --- a/nlmod/gwf/output.py +++ b/nlmod/gwf/output.py @@ -218,7 +218,7 @@ def _get_hds(ds=None, gwf=None, fname_hds=None): return headobj -def _get_cbc(ds=None, gwf=None, fname_cbc=None): +def _get_cbf(ds=None, gwf=None, fname_cbc=None): msg = "Load the budgets using either the ds or the gwf" assert ((ds is not None) + (gwf is not None)) == 1, msg From 21a5604f546d8917753e0835ee84d21e0e6af38c Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:06:06 +0200 Subject: [PATCH 4/9] use new ds_time_from_model func to get time coordinates --- nlmod/gwf/output.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py index a29d5d32..b46a6b1f 100644 --- a/nlmod/gwf/output.py +++ b/nlmod/gwf/output.py @@ -10,6 +10,7 @@ from ..dims.grid import modelgrid_from_ds from ..dims.resample import get_affine, get_xy_mid_structured +from ..dims.time import ds_time_from_model logger = logging.getLogger(__name__) @@ -17,7 +18,7 @@ def get_heads_da(ds=None, gwf=None, fname_hds=None): """Reads heads file given either a dataset or a groundwater flow object. - Note: Calling this function with ds is currently prevered over calling it + Note: Calling this function with ds is currently preferred over calling it with gwf, because the layer and time coordinates can not be fully reconstructed from gwf. @@ -86,11 +87,12 @@ def get_heads_da(ds=None, gwf=None, fname_hds=None): else: assert 0, "Gridtype not supported" - if ds is not None: + # set layer and time coordinates + if gwf is not None: + head_ar.coords["layer"] = np.arange(gwf.modelgrid.nlay) + head_ar.coords["time"] = ds_time_from_model(gwf) + else: head_ar.coords["layer"] = ds.layer - - # TODO: temporarily only add time for when ds is passed because unable to - # exactly recreate ds.time from gwf. head_ar.coords["time"] = ds.time if ds is not None and "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0: From 648a9d24a2bc8dc434c6d1c8c3611fd4ad6daa1d Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:06:14 +0200 Subject: [PATCH 5/9] use new ds_time_from_model func to get time coordinates --- nlmod/gwt/output.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/nlmod/gwt/output.py b/nlmod/gwt/output.py index 2a6f08d3..dc3322d6 100644 --- a/nlmod/gwt/output.py +++ b/nlmod/gwt/output.py @@ -3,11 +3,11 @@ import flopy import numpy as np -import pandas as pd import xarray as xr from ..dims.layers import calculate_thickness from ..dims.resample import get_affine, get_xy_mid_structured +from ..dims.time import ds_time_from_model logger = logging.getLogger(__name__) @@ -97,20 +97,13 @@ def get_concentration_da(ds=None, gwt=None, fname_conc=None): else: assert 0, "Gridtype not supported" - if ds is not None: + # set layer and time coordinates + if gwt is not None: + conc_ar.coords["layer"] = np.arange(gwt.modelgrid.nlay) + conc_ar.coords["time"] = ds_time_from_model(gwt) + else: conc_ar.coords["layer"] = ds.layer - - # TODO: temporarily only add time for when ds is passed because unable to - # exactly recreate ds.time from gwt. - times = np.array( - [ - pd.Timestamp(ds.time.start) - + pd.Timedelta(t, unit=ds.time.time_units[0]) - for t in concobj.get_times() - ], - dtype=np.datetime64, - ) - conc_ar.coords["time"] = times + conc_ar.coords["time"] = ds.time if ds is not None and "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0: affine = get_affine(ds) From 3f1d03799f0f965deaa9215b3fa3433140a03c68 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:10:40 +0200 Subject: [PATCH 6/9] isort+black --- nlmod/dims/time.py | 2 +- nlmod/gwf/output.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py index 0bc9a1af..fb6e0c16 100644 --- a/nlmod/dims/time.py +++ b/nlmod/dims/time.py @@ -232,7 +232,7 @@ def ds_time_from_model(gwf): IndexVariable time coordinate for xarray data-array or dataset """ - + start_datetime = gwf.simulation.get_package("TDIS").start_date_time.data if start_datetime is not None: diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py index b46a6b1f..11ed5cd9 100644 --- a/nlmod/gwf/output.py +++ b/nlmod/gwf/output.py @@ -115,6 +115,7 @@ def get_heads_da(ds=None, gwf=None, fname_hds=None): return head_ar + def get_budget_da(text, ds=None, gwf=None, fname_cbc=None, kstpkper=None): """Reads budget file given either a dataset or a groundwater flow object. From bc843bc2382279fa67ebfcd6cf87bbb1e7c657e1 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:11:16 +0200 Subject: [PATCH 7/9] modpath bug is fixed, so set start_date_time in simulation --- nlmod/sim/sim.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/sim/sim.py b/nlmod/sim/sim.py index 75332164..e31a2ec9 100644 --- a/nlmod/sim/sim.py +++ b/nlmod/sim/sim.py @@ -182,7 +182,7 @@ def tdis(ds, sim, pname="tdis"): pname=pname, time_units=ds.time.time_units, nper=len(ds.time), - # start_date_time=ds.time.start, # disable until fix in modpath + start_date_time=pd.Timestamp(ds.time.start).isoformat(), perioddata=tdis_perioddata, ) From 60ab0f6757600020b87b787385cbe484b933e252 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:44:52 +0200 Subject: [PATCH 8/9] use time units from groundwater model object --- nlmod/dims/time.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py index fb6e0c16..6f435631 100644 --- a/nlmod/dims/time.py +++ b/nlmod/dims/time.py @@ -236,8 +236,10 @@ def ds_time_from_model(gwf): start_datetime = gwf.simulation.get_package("TDIS").start_date_time.data if start_datetime is not None: + time_units = gwf.dimensions.simulation_time.get_time_units() dt = pd.to_timedelta( - np.cumsum(gwf.dimensions.simulation_time.get_perioddata()["perlen"]), "D" + np.cumsum(gwf.dimensions.simulation_time.get_perioddata()["perlen"]), + time_units, ) times = pd.Timestamp(start_datetime) + dt From c143a598034398ab7f348b8cd2459aaca4db23c2 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Mon, 15 May 2023 15:48:11 +0200 Subject: [PATCH 9/9] rename cbf to cbc --- nlmod/gwf/output.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/nlmod/gwf/output.py b/nlmod/gwf/output.py index 11ed5cd9..3551a19b 100644 --- a/nlmod/gwf/output.py +++ b/nlmod/gwf/output.py @@ -136,9 +136,9 @@ def get_budget_da(text, ds=None, gwf=None, fname_cbc=None, kstpkper=None): q_ar : xarray.DataArray budget data array. """ - cbfobj = _get_cbf(ds=ds, gwf=gwf, fname_cbc=fname_cbc) + cbcobj = _get_cbc(ds=ds, gwf=gwf, fname_cbc=fname_cbc) - q = cbfobj.get_data(text=text, kstpkper=kstpkper, full3D=True) + q = cbcobj.get_data(text=text, kstpkper=kstpkper, full3D=True) q = np.stack(q) if gwf is not None: @@ -221,18 +221,18 @@ def _get_hds(ds=None, gwf=None, fname_hds=None): return headobj -def _get_cbf(ds=None, gwf=None, fname_cbc=None): +def _get_cbc(ds=None, gwf=None, fname_cbc=None): msg = "Load the budgets using either the ds or the gwf" assert ((ds is not None) + (gwf is not None)) == 1, msg if fname_cbc is None: if ds is None: - cbf = gwf.output.budget() + cbc = gwf.output.budget() else: fname_cbc = os.path.join(ds.model_ws, ds.model_name + ".cbc") if fname_cbc is not None: - cbf = flopy.utils.CellBudgetFile(fname_cbc) - return cbf + cbc = flopy.utils.CellBudgetFile(fname_cbc) + return cbc def get_gwl_from_wet_cells(head, layer="layer", botm=None):