discretization creates dataset

Author: mjreno

.docs/Notebooks/netcdf01_tutorial.py

@@ -62,7 +62,8 @@
 # For the purposes of this tutorial, the specifics of this simulation
 # other than it is a candidate for NetCDF input are not a focus. It
 # is a NetCDF input candidate because it defines a candidate model type
-# (`GWF6`) with packages that support NetCDF input parameters.
+# (`GWF6`) with a structured discretization and packages that support
+# NetCDF input parameters.
 #
 # A NetCDF dataset will be created from array data in the `DIS`, `NPF` and
 # `GHBG` packages. Data will be copied from the package objects into dataset
@@ -274,10 +275,14 @@ def _data_shape(shape):
     return dataset
 
 
+# ## Create simulation workspace
+
 # create temporary directories
 temp_dir = TemporaryDirectory()
 workspace = Path(temp_dir.name)
 
+# ## Write and run baseline simulation
+
 # run the non-netcdf simulation
 sim = create_sim(ws=workspace)
 sim.write_simulation()
@@ -300,9 +305,14 @@ def _data_shape(shape):
     print(fh.read())
 
 # ## Create dataset
+#
+# Create the base xarray dataset from the modelgrid object. This
+# will add required dimensions and coordinate variables to the
+# dataset according to the grid specification. Modeltime is needed
+# to for timeseries support.
 
 # create the dataset
-ds = xr.Dataset()
+ds = gwf.modelgrid.dataset(modeltime=gwf.modeltime)
 
 # ## Access model NetCDF attributes
 #
@@ -322,27 +332,15 @@ def _data_shape(shape):
 for a in nc_info["attrs"]:
     ds.attrs[a] = nc_info["attrs"][a]
 
-# ## Define dimensions relevant to NetCDF input file
+# ## Map dataset dimension names to values
 
 # define dimensional info
-dis = gwf.modelgrid
-xoff = dis.xoffset
-yoff = dis.yoffset
-x = xoff + dis.xycenters[0]
-y = yoff + dis.xycenters[1]
-z = [float(x) for x in range(1, dis.nlay + 1)]
-nstp = sum(gwf.modeltime.nstp)
-time = gwf.modeltime.tslen
-nlay = dis.nlay
-nrow = dis.nrow
-ncol = dis.ncol
-dimmap = {"time": nstp, "z": nlay, "y": nrow, "x": ncol}
-
-# ## Create dataset dimensions
-
-# create dataset coordinate vars
-var_d = {"time": (["time"], time), "z": (["z"], z), "y": (["y"], y), "x": (["x"], x)}
-ds = ds.assign(var_d)
+dimmap = {
+    "time": sum(gwf.modeltime.nstp),
+    "z": gwf.modelgrid.nlay,
+    "y": gwf.modelgrid.nrow,
+    "x": gwf.modelgrid.ncol,
+}
 
 # ## Access package NetCDF attributes
 #
@@ -406,7 +404,6 @@ def _data_shape(shape):
 with open(workspace / "netcdf" / "uzf01.dis", "r") as fh:
     print(fh.read())
 
-
 # ## Update MODFLOW 6 package input file
 #
 # Follow the same process as above for the `NPF` package.

flopy/discretization/grid.py

@@ -1260,6 +1260,20 @@ def write_shapefile(self, filename="grid.shp", crs=None, prjfile=None, **kwargs)
         )
         return
 
+    def dataset(self, mesh=None):
+        """
+        Method to generate baseline Xarray dataset
+
+        Parameters
+        ----------
+        mesh
+
+        Returns
+        -------
+            Xarray dataset
+        """
+        raise NotImplementedError("dataset must be defined in the child class")
+
     # initialize grid from a grb file
     @classmethod
     def from_binary_grid_file(cls, file_path, verbose=False):

flopy/discretization/structuredgrid.py

@@ -1675,6 +1675,69 @@ def get_plottable_layer_array(self, a, layer):
         assert plotarray.shape == required_shape, msg
         return plotarray
 
+    def dataset(self, modeltime=None, mesh=None):
+        import xarray as xr
+
+        lenunits = {0: "m", 1: "ft", 2: "m", 3: "m"}
+
+        ds = xr.Dataset()
+
+        x = self.xoffset + self.xycenters[0]
+        y = self.yoffset + self.xycenters[1]
+        z = [float(x) for x in range(1, self.nlay + 1)]
+
+        # set coordinate var bounds
+        x_bnds = []
+        xv = self.xoffset + self.xyedges[0]
+        for idx, val in enumerate(xv):
+            if idx + 1 < len(xv):
+                bnd = []
+                bnd.append(xv[idx])
+                bnd.append(xv[idx + 1])
+                x_bnds.append(bnd)
+
+        y_bnds = []
+        yv = self.yoffset + self.xyedges[1]
+        for idx, val in enumerate(yv):
+            if idx + 1 < len(yv):
+                bnd = []
+                bnd.append(yv[idx + 1])
+                bnd.append(yv[idx])
+                y_bnds.append(bnd)
+
+        # create dataset coordinate vars
+        var_d = {
+            "time": (["time"], modeltime.totim),
+            "z": (["z"], z),
+            "y": (["y"], y),
+            "x": (["x"], x),
+        }
+        ds = ds.assign(var_d)
+
+        # create bound vars
+        var_d = {"x_bnds": (["x", "bnd"], x_bnds), "y_bnds": (["y", "bnd"], y_bnds)}
+        ds = ds.assign(var_d)
+
+        ds["time"].attrs["calendar"] = "standard"
+        ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
+        ds["time"].attrs["axis"] = "T"
+        ds["time"].attrs["standard_name"] = "time"
+        ds["time"].attrs["long_name"] = "time"
+        ds["z"].attrs["units"] = "layer"
+        ds["z"].attrs["long_name"] = "layer number"
+        ds["y"].attrs["units"] = lenunits[self.lenuni]
+        ds["y"].attrs["axis"] = "Y"
+        ds["y"].attrs["standard_name"] = "projection_y_coordinate"
+        ds["y"].attrs["long_name"] = "Northing"
+        ds["y"].attrs["bounds"] = "y_bnds"
+        ds["x"].attrs["units"] = lenunits[self.lenuni]
+        ds["x"].attrs["axis"] = "X"
+        ds["x"].attrs["standard_name"] = "projection_x_coordinate"
+        ds["x"].attrs["long_name"] = "Easting"
+        ds["x"].attrs["bounds"] = "x_bnds"
+
+        return ds
+
     def _set_structured_iverts(self):
         """
         Build a list of the vertices that define each model cell and the x, y