diff --git a/docs/source/conf.py b/docs/source/conf.py index 2b609afc..6ff2bb15 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -10,11 +10,12 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # +from nlmod import __version__ import os import sys sys.path.insert(0, os.path.abspath('.')) -from nlmod import __version__ + # -- Project information ----------------------------------------------------- diff --git a/nlmod/mdims/mbase.py b/nlmod/mdims/mbase.py index 1cd0592d..058e4056 100644 --- a/nlmod/mdims/mbase.py +++ b/nlmod/mdims/mbase.py @@ -6,7 +6,8 @@ from .. import util -def get_empty_model_ds(model_name, model_ws, mfversion="mf6", + +def get_empty_model_ds(model_name, model_ws, mfversion="mf6", exe_name=None): """ get an empty model dataset @@ -34,7 +35,7 @@ def get_empty_model_ds(model_name, model_ws, mfversion="mf6", model_ds.attrs['model_name'] = model_name model_ds.attrs['mfversion'] = mfversion model_ds.attrs['model_dataset_created_on'] = dt.datetime.now().strftime("%Y%m%d_%H:%M:%S") - + if exe_name is None: exe_name = os.path.join(os.path.dirname(__file__), '..', '..', 'bin', model_ds.mfversion) @@ -42,16 +43,14 @@ def get_empty_model_ds(model_name, model_ws, mfversion="mf6", # if working on Windows add .exe extension if sys.platform.startswith('win'): exe_name += ".exe" - + model_ds.attrs["exe_name"] = exe_name - - # add some directories + + # add some directories if model_ws is not None: figdir, cachedir = util.get_model_dirs(model_ws) model_ds.attrs['model_ws'] = model_ws model_ds.attrs['figdir'] = figdir model_ds.attrs['cachedir'] = cachedir - - return model_ds diff --git a/nlmod/mdims/mtime.py b/nlmod/mdims/mtime.py index ade1ad68..356b083e 100644 --- a/nlmod/mdims/mtime.py +++ b/nlmod/mdims/mtime.py @@ -14,7 +14,7 @@ def set_model_ds_time(model_ds, start_time, steady_state, steady_start=False, time_units='DAYS', transient_timesteps=0, - perlen=1.0, steady_perlen=3650, + perlen=1.0, steady_perlen=3650, nstp=1, tsmult=1.0): """Set timing for a model dataset. diff --git a/nlmod/mdims/resample.py b/nlmod/mdims/resample.py index 01bef40f..680f16c7 100644 --- a/nlmod/mdims/resample.py +++ b/nlmod/mdims/resample.py @@ -179,7 +179,7 @@ def resample_dataarray_to_structured_grid(da_in, extent=None, delr=None, delc=No Also flips the y-coordinates to make them descending instead of ascending. This makes it easier to export array to flopy. In other words, make sure that both lines of code create the same plot:: - + da_in['top'].sel(layer=b'Hlc').plot() plt.imshow(da_in['top'].sel(layer=b'Hlc').data) @@ -635,7 +635,7 @@ def fillnan_dataarray_unstructured_grid(xar_in, gridprops=None, def resample_unstr_2d_da_to_struc_2d_da(da_in, model_ds, cellsize=25, - method ='nearest'): + method='nearest'): """resample a 2d dataarray (xarray) from an unstructured grid to a new dataaraay from a structured grid. @@ -657,22 +657,22 @@ def resample_unstr_2d_da_to_struc_2d_da(da_in, model_ds, cellsize=25, """ points_unstr = np.array([model_ds.x.values, model_ds.y.values]).T - modelgrid_x = np.arange(model_ds.x.values.min(), - model_ds.x.values.max(), + modelgrid_x = np.arange(model_ds.x.values.min(), + model_ds.x.values.max(), cellsize) - modelgrid_y = np.arange(model_ds.y.values.max(), - model_ds.y.values.min(), + modelgrid_y = np.arange(model_ds.y.values.max(), + model_ds.y.values.min(), -cellsize) mg = np.meshgrid(modelgrid_x, modelgrid_y) - points = np.vstack((mg[0].ravel(), mg[1].ravel())).T - + points = np.vstack((mg[0].ravel(), mg[1].ravel())).T + arr_out_1d = griddata(points_unstr, da_in.values, points, method=method) - arr_out2d = arr_out_1d.reshape(len(modelgrid_y), + arr_out2d = arr_out_1d.reshape(len(modelgrid_y), len(modelgrid_x)) - + da_out = xr.DataArray(arr_out2d, - dims=('y','x'), + dims=('y', 'x'), coords={'y': modelgrid_y, 'x': modelgrid_x}) - + return da_out diff --git a/nlmod/mfpackages/surface_water.py b/nlmod/mfpackages/surface_water.py index 4eb2601f..28b0ab35 100644 --- a/nlmod/mfpackages/surface_water.py +++ b/nlmod/mfpackages/surface_water.py @@ -287,7 +287,7 @@ def estimate_polygon_length(gdf): def distribute_cond_over_lays(cond, cellid, rivbot, laytop, laybot, idomain=None, kh=None, stage=None): - if isinstance(rivbot, np.ndarray) or isinstance(rivbot, xr.DataArray): + if isinstance(rivbot, (np.ndarray, xr.DataArray)): rivbot = float(rivbot[cellid]) if len(laybot.shape) == 3: # the grid is structured grid @@ -435,5 +435,5 @@ def build_spd(celldata, pkg, model_ds): spd.append([cid, stage, cond, rbot] + auxlist) elif pkg in ["DRN", "GHB"]: spd.append([cid, stage, cond] + auxlist) - + return spd diff --git a/nlmod/read/geotop.py b/nlmod/read/geotop.py index 726e40d2..885b17a7 100644 --- a/nlmod/read/geotop.py +++ b/nlmod/read/geotop.py @@ -71,7 +71,7 @@ def get_geotop_dataset(extent, delr, delc, geotop_ds = mdims.get_resampled_ml_layer_ds_struc(raw_ds=geotop_ds_raw, extent=extent, delr=delr, delc=delc) - + for datavar in geotop_ds: geotop_ds[datavar].attrs['source'] = 'Geotop' geotop_ds[datavar].attrs['url'] = geotop_url @@ -98,8 +98,6 @@ def get_geotop_raw_within_extent(extent, url): slices geotop netcdf. """ - - geotop_ds_raw = xr.open_dataset(url) # set x and y dimensions to cell center diff --git a/nlmod/read/jarkus.py b/nlmod/read/jarkus.py index b477f4de..ada55536 100644 --- a/nlmod/read/jarkus.py +++ b/nlmod/read/jarkus.py @@ -184,12 +184,11 @@ def bathymetry_to_model_dataset(model_ds, model_ds_out['bathymetry'] = xr.where( model_ds['northsea'], da_bathymetry, np.nan) - + for datavar in model_ds_out: model_ds_out[datavar].attrs['source'] = 'Jarkus' model_ds_out[datavar].attrs['url'] = url model_ds_out[datavar].attrs['source'] = dt.datetime.now().strftime('%Y%m%d') - return model_ds_out diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py index 5ac42bc1..08a333a8 100644 --- a/nlmod/read/knmi.py +++ b/nlmod/read/knmi.py @@ -190,7 +190,7 @@ def add_knmi_to_model_dataset(model_ds, elif model_ds.gridtype == 'unstructured': model_ds_out['recharge'].loc[loc_sel.index, :] = model_recharge.values - + for datavar in model_ds_out: model_ds_out[datavar].attrs['source'] = 'KNMI' model_ds_out[datavar].attrs['date'] = dt.datetime.now().strftime('%Y%m%d') diff --git a/nlmod/read/regis.py b/nlmod/read/regis.py index f24b9682..76cb818b 100644 --- a/nlmod/read/regis.py +++ b/nlmod/read/regis.py @@ -80,8 +80,7 @@ def get_layer_models(extent, delr, delc, def get_combined_layer_models(extent, delr, delc, regis_botm_layer=b'AKc', use_regis=True, use_geotop=True, - remove_nan_layers=True, - cachedir=None, use_cache=False): + remove_nan_layers=True): """combine layer models into a single layer model. Possibilities so far include: @@ -110,12 +109,6 @@ def get_combined_layer_models(extent, delr, delc, if True regis and geotop layers with only nans are removed from the model. if False nan layers are kept which might be usefull if you want to keep some layers that exist in other models. The default is True. - cachedir : str - directory to store cached values, if None a temporary directory is - used. default is None - use_cache : bool, optional - if True the cached resampled regis dataset is used. - The default is False. Raises ------ @@ -216,7 +209,7 @@ def get_regis_dataset(extent, delr, delc, botm_layer=b'AKc'): regis_ds = mdims.get_resampled_ml_layer_ds_struc(raw_ds=regis_ds_raw2, extent=extent, delr=delr, delc=delc) - + for datavar in regis_ds: regis_ds[datavar].attrs['source'] = 'REGIS' regis_ds[datavar].attrs['url'] = regis_url @@ -332,8 +325,6 @@ def add_geotop_to_regis_hlc(regis_ds, geotop_ds, regis_geotop_ds[key].attrs['regis_url'] = regis_ds[key].url regis_geotop_ds[key].attrs['geotop_url'] = geotop_ds[key].url regis_geotop_ds[key].attrs['date'] = dt.datetime.now().strftime('%Y%m%d') - - return regis_geotop_ds diff --git a/nlmod/read/rws.py b/nlmod/read/rws.py index 45560240..68c7b66b 100644 --- a/nlmod/read/rws.py +++ b/nlmod/read/rws.py @@ -94,7 +94,7 @@ def surface_water_to_model_dataset(model_ds, modelgrid, da_name): area = xr.zeros_like(model_ds['top']) cond = xr.zeros_like(model_ds['top']) peil = xr.zeros_like(model_ds['top']) - for i, row in gdf.iterrows(): + for _, row in gdf.iterrows(): area_pol = mdims.polygon_to_area(modelgrid, row['geometry'], xr.ones_like(model_ds['top']), model_ds.gridtype) @@ -106,7 +106,7 @@ def surface_water_to_model_dataset(model_ds, modelgrid, da_name): model_ds_out[f'{da_name}_area'] = area model_ds_out[f'{da_name}_cond'] = cond model_ds_out[f'{da_name}_peil'] = peil - + for datavar in model_ds_out: model_ds_out[datavar].attrs['source'] = 'RWS' model_ds_out[datavar].attrs['date'] = dt.datetime.now().strftime('%Y%m%d') diff --git a/nlmod/util.py b/nlmod/util.py index 4a9b0290..6e986928 100644 --- a/nlmod/util.py +++ b/nlmod/util.py @@ -19,20 +19,16 @@ from shapely.geometry import box from shutil import copyfile -import __main__ - - logger = logging.getLogger(__name__) - def write_and_run_model(gwf, model_ds, write_model_ds=True, nb_path=None): """ write modflow files and run the model. 2 extra options: 1. write the model dataset to cache 2. copy the modelscript (typically a Jupyter Notebook) to the model workspace with a timestamp. - + Parameters ---------- @@ -49,29 +45,26 @@ def write_and_run_model(gwf, model_ds, write_model_ds=True, of a Jupyter Notebook from within the notebook itself. """ - - if not nb_path is None: + + if nb_path is not None: new_nb_fname = f'{dt.datetime.now().strftime("%Y%m%d")}' + os.path.split(nb_path)[-1] - dst = os.path.join(model_ds.model_ws, new_nb_fname) + dst = os.path.join(model_ds.model_ws, new_nb_fname) logger.info(f'write script {new_nb_fname} to model workspace') copyfile(nb_path, dst) - - + if write_model_ds: logger.info('write model dataset to cache') - model_ds.to_netcdf(os.path.join(model_ds.attrs['cachedir'], + model_ds.to_netcdf(os.path.join(model_ds.attrs['cachedir'], 'full_model_ds.nc')) model_ds.attrs['model_dataset_written_to_disk_on'] = dt.datetime.now().strftime("%Y%m%d_%H:%M:%S") - + logger.info('write modflow files to model workspace') gwf.simulation.write_simulation() model_ds.attrs['model_data_written_to_disk_on'] = dt.datetime.now().strftime("%Y%m%d_%H:%M:%S") - + logger.info('run model') assert gwf.simulation.run_simulation()[0] model_ds.attrs['model_ran_on'] = dt.datetime.now().strftime("%Y%m%d_%H:%M:%S") - - def get_model_dirs(model_ws): @@ -105,7 +98,6 @@ def get_model_dirs(model_ws): os.mkdir(cachedir) return figdir, cachedir - def get_model_ds_empty(model_ds): @@ -160,7 +152,7 @@ def check_delr_delc_extent(dic, model_ds): if 'extent' in dic.keys(): key_check = (np.array(dic['extent']) == np.array( model_ds.attrs['extent'])).all() - + if key_check: logger.info('extent of current grid is the same as cached grid') else: @@ -698,6 +690,7 @@ def getmfexes(pth='.', version='', pltfrm=None): return + def add_heads_to_model_ds(model_ds, fname_hds=None): """reads the heads from a modflow .hds file and returns an xarray DataArray. @@ -717,10 +710,8 @@ def add_heads_to_model_ds(model_ds, fname_hds=None): if fname_hds is None: fname_hds = os.path.join(model_ds.model_ws, model_ds.model_name + '.hds') - + head_filled = get_heads_array(fname_hds, gridtype=model_ds.gridtype) - - if model_ds.gridtype == 'unstructured': head_ar = xr.DataArray(data=head_filled[:, :, :], @@ -728,7 +719,7 @@ def add_heads_to_model_ds(model_ds, fname_hds=None): coords={'cid': model_ds.cid, 'layer': model_ds.layer, 'time': model_ds.time}) - elif model_ds.gridtype =='structured': + elif model_ds.gridtype == 'structured': head_ar = xr.DataArray(data=head_filled, dims=('time', 'layer', 'y', 'x'), coords={'x': model_ds.x, @@ -738,6 +729,7 @@ def add_heads_to_model_ds(model_ds, fname_hds=None): return head_ar + def get_heads_array(fname_hds, gridtype='structured', fill_nans=True): """reads the heads from a modflow .hds file and returns a numpy array. @@ -767,7 +759,7 @@ def get_heads_array(fname_hds, gridtype='structured', if gridtype == 'unstructured': head_filled = np.ones((head.shape[0], head.shape[1], head.shape[3])) * np.nan - + for t in range(head.shape[0]): for lay in range(head.shape[1] - 1, -1, -1): head_filled[t][lay] = head[t][lay][0] @@ -777,7 +769,7 @@ def get_heads_array(fname_hds, gridtype='structured', head_filled[t][lay + 1], head_filled[t][lay]) - elif gridtype =='structured': + elif gridtype == 'structured': head_filled = np.zeros_like(head) for t in range(head.shape[0]): for lay in range(head.shape[1] - 1, -1, -1): @@ -789,7 +781,7 @@ def get_heads_array(fname_hds, gridtype='structured', head_filled[t][lay]) else: raise ValueError('wrong gridtype') - + return head_filled @@ -811,5 +803,3 @@ def download_mfbinaries(binpath=None, version='6.0'): pltfrm = get_platform(None) # Download and unpack mf6 exes getmfexes(pth=binpath, version=version, pltfrm=pltfrm) - - diff --git a/tests/test_001_model.py b/tests/test_001_model.py index 72b43bc0..65b22c77 100644 --- a/tests/test_001_model.py +++ b/tests/test_001_model.py @@ -97,21 +97,14 @@ """ -import datetime as dt import os import tempfile - -import flopy -import geopandas as gpd import nlmod import pytest import xarray as xr -import test_002_regis_geotop - tmpdir = tempfile.gettempdir() - tst_model_dir = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'data') @@ -149,7 +142,7 @@ def test_model_ds_time_transient(tmpdir, modelname='test'): def test_create_seamodel_grid_only_without_northsea(tmpdir): model_ds = test_model_ds_time_transient(tmpdir) extent = [95000., 105000., 494000., 500000.] - extent, nrow, ncol = nlmod.read.regis.fit_extent_to_regis(extent, 100, 100) + extent, _, _ = nlmod.read.regis.fit_extent_to_regis(extent, 100, 100) regis_geotop_ds = nlmod.read.regis.get_layer_models(extent, 100., 100., use_regis=True, use_geotop=True) @@ -184,7 +177,7 @@ def test_create_small_model_grid_only(tmpdir): 'x', 'y'], gridtype='structured') - sim, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds) + _, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds) # Create discretization nlmod.mfpackages.dis_from_model_ds(model_ds, gwf) @@ -241,7 +234,7 @@ def test_create_sea_model(tmpdir): model_ds = xr.open_dataset(os.path.join(tst_model_dir, 'basic_sea_model.nc')) # create modflow packages - sim, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds) + _, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds) # Create discretization nlmod.mfpackages.dis_from_model_ds(model_ds, gwf) @@ -260,14 +253,14 @@ def test_create_sea_model(tmpdir): model_ds = nlmod.read.rws.surface_water_to_model_dataset(model_ds, gwf.modelgrid, da_name) - ghb = nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name) + nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name) # surface level drain model_ds = nlmod.read.ahn.get_ahn_dataset(model_ds) - drn = nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) + nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) # add constant head cells at model boundaries - chd = nlmod.mfpackages.chd_at_model_edge_from_model_ds( + nlmod.mfpackages.chd_at_model_edge_from_model_ds( model_ds, gwf, head='starting_head') # add knmi recharge to the model datasets @@ -332,14 +325,14 @@ def test_create_sea_model_perlen_list(tmpdir): model_ds = nlmod.read.rws.surface_water_to_model_dataset(model_ds, gwf.modelgrid, da_name) - ghb = nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name) + nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name) # surface level drain model_ds = nlmod.read.ahn.get_ahn_dataset(model_ds) - drn = nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) + nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) # add constant head cells at model boundaries - chd = nlmod.mfpackages.chd_at_model_edge_from_model_ds( + nlmod.mfpackages.chd_at_model_edge_from_model_ds( model_ds, gwf, head='starting_head') # add knmi recharge to the model datasets @@ -396,14 +389,14 @@ def test_create_sea_model_perlen_14(tmpdir): model_ds = nlmod.read.rws.surface_water_to_model_dataset(model_ds, gwf.modelgrid, da_name) - ghb = nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name) + nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name) # surface level drain model_ds = nlmod.read.ahn.get_ahn_dataset(model_ds) - drn = nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) + nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf) # add constant head cells at model boundaries - chd = nlmod.mfpackages.chd_at_model_edge_from_model_ds( + nlmod.mfpackages.chd_at_model_edge_from_model_ds( model_ds, gwf, head='starting_head') # add knmi recharge to the model datasets diff --git a/tests/test_007_run_notebooks.py b/tests/test_007_run_notebooks.py index a23dff88..a9a7eda0 100644 --- a/tests/test_007_run_notebooks.py +++ b/tests/test_007_run_notebooks.py @@ -43,10 +43,12 @@ def test_run_notebook_03_local_grid_refinement(): def test_run_notebook_04_modifying_layermodels(): _run_notebook(nbdir, '04_modifying_layermodels.ipynb') + @pytest.mark.notebooks def test_run_notebook_05_caching(): _run_notebook(nbdir, '05_caching.ipynb') - + + @pytest.mark.notebooks def test_run_notebook_06_compare_layermodels(): - _run_notebook(nbdir, '06_compare_layermodels.ipynb') \ No newline at end of file + _run_notebook(nbdir, '06_compare_layermodels.ipynb')