Generating model datasets

This notebook contains two workflows to generate a model dataset and add layer information to it:

• Get data (download or read from file)

• Regrid (convert data to match your desired grid)

or

• Create grid (define your desired grid)

• Add data to dataset (add data based on this grid)

from shapely.geometry import LineString

import nlmod

nlmod.util.get_color_logger("INFO")
nlmod.show_versions()

Python version     : 3.11.14
NumPy version      : 2.4.4
Xarray version     : 2026.4.0
Matplotlib version : 3.10.9
Flopy version      : 3.10.0

nlmod version      : 0.11.3dev

Get REGIS data

extent = [100000, 101000, 400000, 401000]

# downlaod regis
regis = nlmod.read.regis.download_regis(extent)

f, ax = nlmod.plot.get_map(extent, figsize=5)
nlmod.plot.data_array(regis["top"], edgecolor="k");

Define some properties of the grid

We choose a resolution of the calculation grid (delr and delc) larger than the resolution of regis (100 meter).

delr = 200.0
delc = 200.0

polygon = LineString([(extent[0], extent[2]), (extent[1], extent[3])]).buffer(10)
refinement_features = [([polygon], "polygon", 3)]

Regrid data

create a model dataset from regis and refine.

ds = nlmod.to_model_ds(regis, delr=delr, delc=delc)
ds = nlmod.grid.refine(ds, model_ws="01_layer_generation", refinement_features=refinement_features)

INFO:nlmod.dims.base.to_model_ds:resample layer model data to structured modelgrid
INFO:nlmod.dims.layers.get_kh_kv:kv and kh both undefined in layer HLc
INFO:nlmod.dims.layers._fill_var:Filling 153 values in active cells of kh by multipying kv with an anisotropy of 10
INFO:nlmod.dims.layers._fill_var:Filling 304 values in active cells of kv by dividing kh by an anisotropy of 10
INFO:nlmod.dims.layers._fill_var:Filling 22 values in active cells of kh with a value of 1.0 m/day
INFO:nlmod.dims.layers._fill_var:Filling 22 values in active cells of kv with a value of 0.1 m/day
INFO:nlmod.dims.grid.refine:create vertex grid using gridgen
INFO:nlmod.dims.grid.ds_to_gridprops:resample model Dataset to vertex modelgrid

When we plot the model top, we see that all information has a resolution of 200 meter, also in the smaller cells. The original resolution of the regis data is 100 meter however. So information is lost in the process.

f, ax = nlmod.plot.get_map(extent, figsize=5)
nlmod.plot.data_array(ds["top"], ds=ds, edgecolor="k");

Add data to existing grid

We can also first create a grid first, and then warp the information from regis to this grid.

ds = nlmod.get_ds(extent, delr=delr, delc=delc)
ds = nlmod.grid.refine(ds, model_ws="01_layer_generation", refinement_features=refinement_features)
ds = nlmod.grid.update_ds_from_layer_ds(ds, regis, method="average")

INFO:nlmod.dims.base.to_model_ds:resample layer model data to structured modelgrid
INFO:nlmod.dims.grid.refine:create vertex grid using gridgen
INFO:nlmod.dims.grid.ds_to_gridprops:resample model Dataset to vertex modelgrid
INFO:nlmod.dims.layers.get_kh_kv:kv and kh both undefined in layer HLc
INFO:nlmod.dims.layers._fill_var:Filling 1934 values in active cells of kh by multipying kv with an anisotropy of 10.0
INFO:nlmod.dims.layers._fill_var:Filling 4034 values in active cells of kv by dividing kh by an anisotropy of 10.0
INFO:nlmod.dims.layers._fill_var:Filling 308 values in active cells of kh with a value of 1.0
INFO:nlmod.dims.layers._fill_var:Filling 308 values in active cells of kv with a value of 0.1

When we plot the model top, we now see that in the cells with an equal resolution to that of Regis (or smaller) no information is lost.

f, ax = nlmod.plot.get_map(extent, figsize=5)
nlmod.plot.data_array(ds["top"], ds=ds, edgecolor="k");