feat(netcdf): support writing simulation

[mjreno]

NetCDF simulation write support:

• netcdf argument supported for write_simulation()
• write and load support for READARRAYGRID packages including layered aux variables
• simulation load not yet supported for NetCDF input files
• write argument can create structured, layered or no NetCDF file (ASCII inputs updated in all cases)
• model and package support for interactively creating dataset

Non-interactive example:

# load ascii simulation 
sim = flopy.mf6.MFSimulation.load(sim_ws=ascii_ws)
# set simulation path and write simulation
sim.set_sim_path(netcdf_ws)
sim.write_simulation(netcdf="structured")

Internally, the dataset is generated by the model modelgrid object. Package variables and data are added by mf6 package objects. User facing functions are intended, however, to support different approaches to building the dataset.

Create a base dataset:
ds = gwf.modelgrid.dataset(modeltime=gwf.modeltime)

Add package variables with data:

dis = gwf.get_package("dis")
ds = dis.update_dataset(ds)

To retrieve information about package netcdf variables without requesting the package to update the dataset:

nc_meta = dis.netcdf_meta()

{'botm': {'attrs': {'_FillValue': 9.96920996838687e+36,
                    'longname': 'cell bottom elevation',
                    'modflow_input': 'UZF01/DIS/BOTM'},
          'netcdf_shape': ['z', 'y', 'x'],
          'varname': 'dis_botm',
          'xarray_type': <class 'numpy.float64'>},
 'delc': {'attrs': {'_FillValue': 9.96920996838687e+36,
                    'longname': 'spacing along a column',
                    'modflow_input': 'UZF01/DIS/DELC'},
          'netcdf_shape': ['y'],
          'varname': 'dis_delc',
          'xarray_type': <class 'numpy.float64'>},
 'delr': {'attrs': {'_FillValue': 9.96920996838687e+36,
                    'longname': 'spacing along a row',
                    'modflow_input': 'UZF01/DIS/DELR'},
          'netcdf_shape': ['x'],
          'varname': 'dis_delr',
          'xarray_type': <class 'numpy.float64'>},
 'idomain': {'attrs': {'_FillValue': np.int32(-2147483647),
                       'longname': 'idomain existence array',
                       'modflow_input': 'UZF01/DIS/IDOMAIN'},
             'netcdf_shape': ['z', 'y', 'x'],
             'varname': 'dis_idomain',
             'xarray_type': <class 'numpy.int32'>},
 'top': {'attrs': {'_FillValue': 9.96920996838687e+36,
                   'longname': 'cell top elevation',
                   'modflow_input': 'UZF01/DIS/TOP'},
         'netcdf_shape': ['y', 'x'],
         'varname': 'dis_top',
         'xarray_type': <class 'numpy.float64'>}}

Or, if a simulation doesn't yet exist:

import flopy
netcdf_meta = flopy.mf6.mfpackage.MFPackage.netcdf_package("GWF", "DIS")

In this case the dictionary contains very similar information to the one above but that is not always true. This approach offers the ability to construct a MODFLOW 6 compliant NetCDF input file with help from modelgrid / modeltime and FloPy static model and package functions.

Layered Mesh (UGRID) versions of these examples also are supported.

[codecov]

Codecov Report

❌ Patch coverage is 13.64162% with 747 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.2%. Comparing base (556c088) to head (5ea9e0a).
⚠️ Report is 75 commits behind head on develop.

Files with missing lines	Patch %	Lines
flopy/mf6/mfpackage.py	6.9%	229 Missing ⚠️
flopy/discretization/structuredgrid.py	2.1%	183 Missing ⚠️
flopy/mf6/mfmodel.py	14.4%	136 Missing ⚠️
flopy/discretization/vertexgrid.py	1.1%	83 Missing ⚠️
flopy/mf6/data/mfdataarray.py	39.2%	51 Missing ⚠️
flopy/discretization/modeltime.py	7.1%	26 Missing ⚠️
flopy/mf6/data/mfdatastorage.py	39.3%	20 Missing ⚠️
flopy/mf6/data/mffileaccess.py	33.3%	8 Missing ⚠️
flopy/mf6/utils/codegen/filters.py	0.0%	6 Missing ⚠️
flopy/mf6/data/mfstructure.py	81.8%	2 Missing ⚠️
... and 3 more

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #2564      +/-   ##
===========================================
+ Coverage     55.5%    72.2%   +16.6%     
===========================================
  Files          644      667      +23     
  Lines       124135   130721    +6586     
===========================================
+ Hits         68947    94422   +25475     
+ Misses       55188    36299   -18889     

Files with missing lines	Coverage Δ
flopy/mf6/mfbase.py	74.7% <100.0%> (-12.9%)	⬇️
flopy/mf6/mfsimbase.py	62.1% <100.0%> (-13.3%)	⬇️
flopy/utils/datautil.py	65.0% <100.0%> (-4.8%)	⬇️
flopy/discretization/grid.py	75.7% <87.5%> (-0.2%)	⬇️
flopy/discretization/unstructuredgrid.py	75.0% <50.0%> (-6.5%)	⬇️
flopy/mf6/utils/codegen/__init__.py	0.0% <0.0%> (ø)
flopy/mf6/data/mfstructure.py	74.1% <81.8%> (+0.3%)	⬆️
flopy/mf6/utils/codegen/filters.py	0.0% <0.0%> (ø)
flopy/mf6/data/mffileaccess.py	70.9% <33.3%> (-5.1%)	⬇️
flopy/mf6/data/mfdatastorage.py	69.8% <39.3%> (-4.1%)	⬇️
... and 6 more

... and 542 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[wpbonelli]

I think the pattern here is great — iteratively building a dataset from select packages.

A similar way to achieve the same thing with xarray could be a to_dataset() method on the packages. Then you can do xr.merge on any number of them to create a combined dataset. And I think this gets rid of the need for the "template" datasets with the right size and dummy data currently being created for dimension-consuming packages — as long as the dimension-defining package is in the merge, the others can reference its dims.