@@ -242,44 +242,12 @@ def create_sim(ws):
242242 return sim
243243
244244
245- # ## Create helper function to update dataset
246- #
247- # This function updates an xarray dataset to add variables described
248- # in a FloPy provided dictionary.
249- #
250- # The dimmap variable relates NetCDF dimension names to a value.
251-
252-
253- # A subroutine that can update an xarray dataset with package
254- # netcdf information stored in a dict
255- def add_netcdf_vars (dataset , nc_info , dimmap ):
256- def _data_shape (shape ):
257- dims_l = []
258- for d in shape :
259- dims_l .append (dimmap [d ])
260-
261- return dims_l
262-
263- for v in nc_info :
264- varname = nc_info [v ]["varname" ]
265- data = np .full (
266- _data_shape (nc_info [v ]["netcdf_shape" ]),
267- nc_info [v ]["attrs" ]["_FillValue" ],
268- dtype = nc_info [v ]["xarray_type" ],
269- )
270- var_d = {varname : (nc_info [v ]["netcdf_shape" ], data )}
271- dataset = dataset .assign (var_d )
272- for a in nc_info [v ]["attrs" ]:
273- dataset [varname ].attrs [a ] = nc_info [v ]["attrs" ][a ]
274-
275- return dataset
276-
277-
278245# ## Create simulation workspace
279246
280247# create temporary directories
281- temp_dir = TemporaryDirectory ()
282- workspace = Path (temp_dir .name )
248+ # temp_dir = TemporaryDirectory()
249+ # workspace = Path(temp_dir.name)
250+ workspace = Path ("./working" )
283251
284252# ## Write and run baseline simulation
285253
@@ -323,49 +291,39 @@ def _data_shape(shape):
323291# First, retrieve and store the netcdf info dictionary and display
324292# its contents. Then, in the following step, update the dataset with
325293# the model scoped attributes defined in the dictionary.
294+ #
295+ # These 2 operations can also be accomplised by calling `update_dataset()`
296+ # on the model object. Analogous functions for the package are shown
297+ # below.
326298
327299# get model netcdf info
328300nc_info = gwf .netcdf_info ()
329301pprint (nc_info )
330302
331- # update dataset with required attributes
303+ # update dataset directly with required attributes
332304for a in nc_info ["attrs" ]:
333305 ds .attrs [a ] = nc_info ["attrs" ][a ]
334306
335- # ## Map dataset dimension names to values
336-
337- # define dimensional info
338- dimmap = {
339- "time" : sum (gwf .modeltime .nstp ),
340- "z" : gwf .modelgrid .nlay ,
341- "y" : gwf .modelgrid .nrow ,
342- "x" : gwf .modelgrid .ncol ,
343- }
344-
345- # ## Access package NetCDF attributes
346- #
347- # Access package scoped NetCDF details by storing the dictionary returned
348- # from `netcdf_info()`. We need to set package variable attributes that are
349- # stored in the package netcdf info dict, but we also need other information
350- # that is relevant to creating the variables themselves.
307+ # ## Update the dataset with supported `DIS` arrays
351308#
352- # The contents of the info dictionary are shown and then, in the following
353- # step, the dictionary and the dataset are passed to a helper routine that
354- # create the intended array variables.
309+ # Add NetCDF supported data arrays in package to dataset. Internally, this call
310+ # uses a `netcdf_info()` package dictionary to determine candidate variables
311+ # and relevant information about them. Alternatively, this dictionary can
312+ # be directly accessed, updated, and passed to the `update_dataset()` function.
313+ # That workflow will be demonstrated in the `NPF` package update which follows.
355314
356- # get dis package netcdf info
315+ # update dataset with `DIS` arrays
357316dis = gwf .get_package ("dis" )
358- nc_info = dis .netcdf_info ()
359- pprint (nc_info )
360-
361- # create dis dataset variables
362- ds = add_netcdf_vars (ds , nc_info , dimmap )
317+ ds = dis .update_dataset (ds )
363318
364319# ## Update array data
365320#
366321# We have created dataset array variables for the package but they do not yet
367322# define the expected input data for MODFLOW 6. We will take advantage of the
368323# existing simulation objects and update the dataset.
324+ #
325+ # Default dataset variable names are defined in the package `netcdf_info()`
326+ # dictionary.
369327
370328# update dataset from dis arrays
371329ds ["dis_delr" ].values = dis .delr .get_data ()
@@ -378,7 +336,7 @@ def _data_shape(shape):
378336#
379337# MODFLOW 6 input data for the package is now in the dataset. Once the NetCDF
380338# file is generated, we need to configure MODFLOW 6 so that it looks to that
381- # file for the package array input. The ASCII file will no longer defined the
339+ # file for the package array input. The ASCII file will no longer define the
382340# arrays- instead the array names will be followed by the NETCDF keyword.
383341#
384342# We will simply overwrite the entire MODFLOW 6 `DIS` package input file with the
@@ -404,21 +362,46 @@ def _data_shape(shape):
404362with open (workspace / "netcdf" / "uzf01.dis" , "r" ) as fh :
405363 print (fh .read ())
406364
407- # ## Update MODFLOW 6 package input file
365+ # ## Access `NPF` package NetCDF attributes
408366#
409- # Follow the same process as above for the `NPF` package.
367+ # Access package scoped NetCDF details by storing the dictionary returned
368+ # from `netcdf_info()`. We need to set package variable attributes that are
369+ # stored in the package netcdf info dict, but we also need other information
370+ # that is relevant to creating the variables themselves.
371+ #
372+ # The contents of the info dictionary are shown and then, in the following
373+ # step, the dictionary and the dataset are passed to a helper routine that
374+ # create the intended array variables.
410375
411376# get npf package netcdf info
412377npf = gwf .get_package ("npf" )
413378nc_info = npf .netcdf_info ()
414379pprint (nc_info )
415380
416- # create npf dataset variables
417- ds = add_netcdf_vars (ds , nc_info , dimmap )
381+ # ## Update package `netcdf_info` dictionary and dataset
382+ #
383+ # Here we replace the default name for the `NPF K` input parameter and add
384+ # the `standard_name` attribute to it's attribute dictionary. The dictionary
385+ # is then passed to the `update_dataset()` function. Note the udpated name
386+ # is used in the subsequent block when updating the array values.
387+
388+ # update dataset with `NPF` arrays
389+ nc_info ["k" ]["varname" ] = "npf_k_updated"
390+ nc_info ["k" ]["attrs" ]["standard_name" ] = "soil_hydraulic_conductivity_at_saturation"
391+ ds = npf .update_dataset (ds , netcdf_info = nc_info )
392+
393+ # ## Update array data
418394
419395# update dataset from npf arrays
420396ds ["npf_icelltype" ].values = npf .icelltype .get_data ()
421- ds ["npf_k" ].values = npf .k .get_data ()
397+ ds ["npf_k_updated" ].values = npf .k .get_data ()
398+
399+ # ## Show dataset `NPF K` parameter with udpates
400+
401+ # print dataset npf k variable
402+ print (ds ["npf_k_updated" ])
403+
404+ # ## Update MODFLOW 6 package input file
422405
423406# rewrite mf6 npf input to read from netcdf
424407with open (workspace / "netcdf" / "uzf01.npf" , "w" ) as f :
@@ -431,21 +414,13 @@ def _data_shape(shape):
431414with open (workspace / "netcdf" / "uzf01.npf" , "r" ) as fh :
432415 print (fh .read ())
433416
434- # ## Update MODFLOW 6 package input file
435- #
436- # Follow the same process as above for the `GHBG` package. The difference is
437- # that this is PERIOD input and therefore stored as timeseries data in the
438- # NetCDF file. As NETCDF timeseries are defined in terms of total number of
439- # simulation steps, care must be taken in the translation of FloPy period
440- # data to the timeseries.
417+ # ## Update the dataset with supported `GHBG` arrays
441418
442- # get ghbg package netcdf info
419+ # update dataset with 'GHBG' arrays
443420ghbg = gwf .get_package ("ghbg_0" )
444- nc_info = ghbg .netcdf_info ()
445- pprint (nc_info )
421+ ds = ghbg .update_dataset (ds )
446422
447- # create ghbg dataset variables
448- ds = add_netcdf_vars (ds , nc_info , dimmap )
423+ # ## Update array data
449424
450425# update bhead netcdf array from flopy perioddata
451426# timeseries step index is first of stress period
@@ -459,17 +434,7 @@ def _data_shape(shape):
459434 istp = sum (gwf .modeltime .nstp [0 :p ])
460435 ds ["ghbg_0_cond" ].values [istp ] = ghbg .cond .get_data ()[p ]
461436
462- # ## Display generated dataset
463-
464- # show the dataset
465- print (ds )
466-
467- # ## Export generated dataset to NetCDF
468-
469- # write dataset to netcdf
470- ds .to_netcdf (
471- workspace / "netcdf/uzf01.structured.nc" , format = "NETCDF4" , engine = "netcdf4"
472- )
437+ # ## Update MODFLOW 6 package input file
473438
474439# rewrite mf6 ghbg input to read from netcdf
475440with open (workspace / "netcdf/uzf01.ghbg" , "w" ) as f :
@@ -486,6 +451,18 @@ def _data_shape(shape):
486451with open (workspace / "netcdf" / "uzf01.ghbg" , "r" ) as fh :
487452 print (fh .read ())
488453
454+ # ## Display generated dataset
455+
456+ # show the dataset
457+ print (ds )
458+
459+ # ## Export generated dataset to NetCDF
460+
461+ # write dataset to netcdf
462+ ds .to_netcdf (
463+ workspace / "netcdf/uzf01.structured.nc" , format = "NETCDF4" , engine = "netcdf4"
464+ )
465+
489466# ## Run MODFLOW 6 simulation with NetCDF input
490467#
491468# The simulation generated by this tutorial should be runnable by
0 commit comments