<jupyter_start><jupyter_text>Local grid refinement   *<NAME>, Artesia, 2021*This notebook shows an example of a vertex model create with the `nlmod` package.  Contents1. [Create model](create)2. [Local grid refinement](lgr)3. [Run Model](run)3. [Visualise](Visualise)<jupyter_code>import matplotlib.pyplot as plt
import flopy
import os
import geopandas as gpd
import logging

import nlmod
print(f'nlmod version: {nlmod.__version__}')

# toon informatie bij het aanroepen van functies
logging.basicConfig(level=logging.INFO)<jupyter_output>nlmod version: 0.2.0<jupyter_text>[1. Create model](TOC)Modflow 6 makes it possible to use locally refined grids. In NLmod you can use a shapefile and a number of levels to specify where and how much you want to use local grid refinement. Below we use a shapefile of the Planetenweg in IJmuiden and set the refinement levels at 2. This well create a grid with cells of 100x100m except at the Planetenweg where the cells will be refined to 25x25m. See also figures below.<jupyter_code># model settings vertex
model_ws = 'model3'
model_name = 'IJm_planeten'
refine_shp_fname = os.path.join(
    nlmod.NLMOD_DATADIR, 'shapes', 'planetenweg_ijmuiden')
levels = 2
extent = [95000., 105000., 494000., 500000.]
delr = 100.
delc = 100.
steady_state = False
steady_start = True
transient_timesteps = 5
perlen = 1.
start_time = '2015-1-1'
gridtype = 'vertex'
use_regis = True
regis_botm_layer = b'MSz1'
use_geotop = True
add_northsea = True
starting_head = 1.0
# create model time dataset
model_ds = nlmod.mdims.get_empty_model_ds(model_name, model_ws)

# add time discretisation
model_ds = nlmod.mdims.set_model_ds_time(model_ds,
                                         start_time,
                                         steady_state,
                                         steady_start=steady_start,
                                         transient_timesteps=transient_timesteps,
                                         perlen=perlen)

# create modflow packages
sim, gwf = nlmod.mfpackages.sim_tdis_gwf_ims_from_model_ds(model_ds)

extent, nrow, ncol = nlmod.read.regis.fit_extent_to_regis(extent,
                                                          delr,
                                                          delc)

# layer model
layer_model = nlmod.read.regis.get_combined_layer_models(extent, delr, delc,
                                                         use_regis=use_regis,
                                                         regis_botm_layer=regis_botm_layer,
                                                         use_geotop=use_geotop,
                                                         cachedir=model_ds.cachedir,
                                                         cachename='combined_layer_ds.nc')<jupyter_output>INFO:nlmod.mfpackages.mfpackages:creating modflow SIM, TDIS, GWF and IMS
INFO:nlmod.read.regis:redefining current extent: [95000.0, 105000.0, 494000.0, 500000.0], fit to regis raster
INFO:nlmod.read.regis:new extent is [95000.0, 105000.0, 494000.0, 500000.0] model has 60 rows and 100 columns
INFO:nlmod.cache:using cached data -> combined_layer_ds.nc<jupyter_text>[2. Local grid refinement](TOC)the code below applies a local grid refinement to the layer model. The local grid refinement is based on the shapefile 'planetenweg_ijmuiden.shp', which contains a line shape of the Planetenweg, and the levels, which is 2. This means that the model cells at the Planetenweg will get a size of 25 x 25m.<jupyter_code># use gridgen to create vertex grid
model_ds.attrs['gridgen_ws'] = os.path.join(model_ws, 'gridgen')
if not os.path.exists(model_ds.gridgen_ws):
    os.mkdir(model_ds.gridgen_ws)
    
gridprops = nlmod.mdims.create_vertex_grid(model_name, model_ds.gridgen_ws, gwf,
                                           [(refine_shp_fname, 'line', levels)],
                                           extent,
                                           layer_model.dims['layer'],
                                           nrow, ncol,
                                           delr, delc,
                                           cachedir=model_ds.cachedir, cachename='gridprops.pklz')

# add layer model to vertex grid
layer_model_unstr = nlmod.mdims.get_resampled_ml_layer_ds_vertex(raw_ds=layer_model,
                                                                 extent=extent,
                                                                 gridprops=gridprops)

# combine model time dataset with layer model dataset
model_ds = nlmod.mdims.update_model_ds_from_ml_layer_ds(model_ds,
                                                        layer_model_unstr,
                                                        gridtype,
                                                        keep_vars=['x', 'y'],
                                                        gridprops=gridprops,
                                                        add_northsea=add_northsea,
                                                        cachedir=model_ds.cachedir)

# Create discretization
disv = nlmod.mfpackages.disv_from_model_ds(model_ds, gwf, gridprops)
# create node property flow
nlmod.mfpackages.npf_from_model_ds(model_ds, gwf)

# Create the initial conditions package
nlmod.mfpackages.ic_from_model_ds(model_ds, gwf, starting_head=starting_head)

# Create the output control package
oc = nlmod.mfpackages.oc_from_model_ds(model_ds, gwf)
# voeg grote oppervlaktewaterlichamen toe
da_name = 'rws_oppwater'
rws_ds = nlmod.read.rws.get_surface_water(model_ds,
                                          da_name,
                                          gridprops=gridprops,
                                          cachedir=model_ds.cachedir,
                                          cachename=da_name+'.nc')
model_ds.update(rws_ds)
ghb = nlmod.mfpackages.ghb_from_model_ds(model_ds, gwf, da_name)

# surface level drain
ahn_ds = nlmod.read.ahn.get_ahn(model_ds, gridprops=gridprops, 
                                cachedir=model_ds.cachedir, cachename='ahn.nc')
model_ds.update(ahn_ds)

drn = nlmod.mfpackages.surface_drain_from_model_ds(model_ds, gwf)


# add constant head cells at model boundaries
model_ds.update(nlmod.mfpackages.constant_head.get_chd_at_model_edge(model_ds, model_ds['idomain']))    
chd = nlmod.mfpackages.chd_from_model_ds(model_ds, gwf, head='starting_head')
# add knmi recharge to the model datasets
knmi_ds = nlmod.read.knmi.get_recharge(model_ds, cachedir=model_ds.cachedir, cachename='recharge')
model_ds.update(knmi_ds)

# create recharge package
rch = nlmod.mfpackages.rch_from_model_ds(model_ds, gwf)
model_ds<jupyter_output><empty_output><jupyter_text>[3. Write and Run](TOC)<jupyter_code>nlmod.util.write_and_run_model(gwf, model_ds, write_model_ds=True, nb_path='03_local_grid_refinement.ipynb')<jupyter_output>INFO:nlmod.util:write script 2021121303_local_grid_refinement.ipynb to model workspace
INFO:nlmod.util:write model dataset to cache
INFO:nlmod.util:write modflow files to model workspace<jupyter_text>[4. Visualise](TOC)Using the `model_ds` and `gwf` variables it is quite easy to visualise model data. Below the modelgrid together with the surface water is shown.<jupyter_code>plan_weg_gdf = gpd.read_file(refine_shp_fname+'.shp')

# plot modelgrid
ax = nlmod.visualise.plots.plot_modelgrid(model_ds, gwf)
plan_weg_gdf.plot(ax=ax, color='r', label='Planetenweg')
ax.legend()
ax.get_figure().savefig(os.path.join(model_ds.figdir, 'mgrid_swater.png'), bbox_inches='tight')

# plot zoomed modelgrid
ax = nlmod.visualise.plots.plot_modelgrid(model_ds, gwf)
ax.set_title('Planetenweg')
plan_weg_gdf.plot(ax=ax, color='r', label='Planetenweg')
ax.set_xlim(100000, 103000)
ax.set_ylim(495000, 497500)
ax.legend()
ax.get_figure().savefig(os.path.join(model_ds.figdir, 'Planetenweg.png'), bbox_inches='tight')<jupyter_output><empty_output><jupyter_text>The model dataset of a vertex model differs from a structured model dataset. The data is stored relative to the cell-id instead of the row and column number. Therefore the model dataset has the dimension cid instead of the dimensions x and y.<jupyter_code>model_ds<jupyter_output><empty_output><jupyter_text>To plot the same rasters as for the previous model we can use the `nlmod.visualise.plots.plot_vertex_array()` function.<jupyter_code>model_ds['vertices'] = nlmod.mdims.get_vertices(model_ds, gridprops=gridprops)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 11))
nlmod.visualise.plots.plot_vertex_array(model_ds['ahn'], model_ds['vertices'], ax=axes[0][0])
nlmod.visualise.plots.plot_vertex_array(model_ds['bot'][0], model_ds['vertices'], ax=axes[0][1])
nlmod.visualise.plots.plot_vertex_array(model_ds['idomain'][0], model_ds['vertices'], ax=axes[1][0])
nlmod.visualise.plots.plot_vertex_array(model_ds['chd'][0], model_ds['vertices'], ax=axes[1][1])
for axes1 in axes:
    for ax in axes1:
        ax.axis('scaled')

fig.savefig(os.path.join(model_ds.figdir, 'ahn_bot_idom_chd.png'), bbox_inches='tight')

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 11))
nlmod.visualise.plots.plot_vertex_array(model_ds['bathymetry'], model_ds['vertices'], ax=axes[0][0])
nlmod.visualise.plots.plot_vertex_array(model_ds['northsea'], model_ds['vertices'], ax=axes[0][1])
nlmod.visualise.plots.plot_vertex_array(model_ds['kh'][1], model_ds['vertices'], ax=axes[1][0])
nlmod.visualise.plots.plot_vertex_array(model_ds['recharge'][:, 0], model_ds['vertices'], ax=axes[1][1])

for axes1 in axes:
    for ax in axes1:
        ax.axis('scaled')
fig.savefig(os.path.join(model_ds.figdir, 'bath_nsea_kh_rch.png'), bbox_inches='tight')
model_ds<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0950067197353458, 'fraction_numerical': 0.02129639201902202, 'mean_word_length': 2.841938046068308, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 60, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '649590', 'license': 'MIT', 'n_tokens_mistral': 3245, 'n_tokens_neox': 3053, 'n_words': 651, 'path': 'examples/03_local_grid_refinement.ipynb', 'repo_name': 'ArtesiaWater/nlmod', 'size': '9678'}