support structured lat/lon

Author: mjreno

.docs/Notebooks/mf6_netcdf01_tutorial.py

@@ -164,10 +164,10 @@
 
 # ## Create NetCDF based simulation method 2
 #
-# In this method we will set the FloPy `netcdf` argument to `nofile`
-# when `write_simulation()` is called. As such, FloPy will not generate
-# the NetCDF file automatically. We will manage the NetCDF file
-# generation ourselves in method 2.
+# In this method we will set the FloPy `netcdf` argument to `` (empty
+# string) when `write_simulation()` is called. As such, FloPy will
+# not generate the NetCDF file automatically. We will manage the NetCDF
+# file generation ourselves in method 2.
 #
 # Reset the simulation path and set the `GWF` name file `nc_filerecord`
 # attribute to the name of the intended input NetCDF file. Display
@@ -188,7 +188,7 @@
 gwf.name_file.nc_filerecord = "uzf01.structured.nc"
 # write simulation with ASCII inputs tagged for NetCDF
 # but do not create NetCDF file
-sim.write_simulation(netcdf="nofile")
+sim.write_simulation(netcdf="")
 
 # ## Show name file with NetCDF input configured
 
@@ -223,7 +223,7 @@
 # details related to creating the variable, including dimensions,
 # shape and data type.
 #
-# The dictionaries are available via the `netcdf_info()` call. Their
+# The dictionaries are available via the `netcdf_meta()` call. Their
 # content also varies depending on the desired type of dataset (i.e.
 # `structured` or `layered`). In this step we will access the model
 # NetCDF metadata and use it to update dataset scoped attributes.
@@ -238,12 +238,12 @@
 # Examples of this approach (with package objects) are shown below.
 
 # get model netcdf info
-nc_info = gwf.netcdf_info()
-pprint(nc_info)
+nc_meta = gwf.netcdf_meta()
+pprint(nc_meta)
 
 # update dataset directly with required attributes
-for a in nc_info["attrs"]:
-    ds.attrs[a] = nc_info["attrs"][a]
+for a in nc_meta["attrs"]:
+    ds.attrs[a] = nc_meta["attrs"][a]
 
 # ## Update the dataset with supported `DIS` arrays
 #
@@ -260,16 +260,16 @@
 # ## Access `NPF` package NetCDF attributes
 #
 # Access package scoped NetCDF details by storing the dictionary returned
-# from `netcdf_info()`.
+# from `netcdf_meta()`.
 #
 # The contents of the info dictionary are shown and then, in the following
 # step, the dictionary and the dataset are passed to a helper routine that
 # create the intended array variables.
 
 # get npf package netcdf info
 npf = gwf.get_package("npf")
-nc_info = npf.netcdf_info()
-pprint(nc_info)
+nc_meta = npf.netcdf_meta()
+pprint(nc_meta)
 
 # ## Update package NetCDF metadata dictionary and dataset
 #
@@ -284,9 +284,9 @@
 # values.
 
 # update dataset with `NPF` arrays
-nc_info["k"]["varname"] = "npf_k_updated"
-nc_info["k"]["attrs"]["standard_name"] = "soil_hydraulic_conductivity_at_saturation"
-ds = npf.update_dataset(ds, netcdf_info=nc_info)
+nc_meta["k"]["varname"] = "npf_k_updated"
+nc_meta["k"]["attrs"]["standard_name"] = "soil_hydraulic_conductivity_at_saturation"
+ds = npf.update_dataset(ds, netcdf_meta=nc_meta)
 
 # ## Show dataset `NPF K` parameter with updates
 
@@ -335,7 +335,7 @@
 # ## Method 3: DIY with xarray
 #
 # The above method still uses FloPy objects to update the dataset arrays
-# to values consistent with the state of the objects. The `netcdf_info`
+# to values consistent with the state of the objects. The `netcdf_meta`
 # dictionary is intended to supported creation of the dataset without
 # an existing simulation defined. The base dataset can be defined with
 # `modelgrid` and `modeltime` objects, while the full package netcdf
@@ -345,16 +345,16 @@
 
 # ## Demonstrate static call on MFPackage (structured dataset):
 
-netcdf_info = flopy.mf6.mfpackage.MFPackage.netcdf_package(
+netcdf_meta = flopy.mf6.mfpackage.MFPackage.netcdf_package(
     mtype="GWF",
     ptype="NPF",
     auxiliary=["CONCENTRATION"],
 )
-pprint(netcdf_info)
+pprint(netcdf_meta)
 
 # ## Demonstrate static call on MFPackage (layered dataset):
 
-netcdf_info = flopy.mf6.mfpackage.MFPackage.netcdf_package(
+netcdf_meta = flopy.mf6.mfpackage.MFPackage.netcdf_package(
     mtype="GWF", ptype="NPF", mesh="LAYERED", auxiliary=["CONCENTRATION"], nlay=2
 )
-pprint(netcdf_info)
+pprint(netcdf_meta)

autotest/regression/test_model_netcdf.py

@@ -87,7 +87,7 @@ def check_netcdf(path, mobj, mesh=None):
 @pytest.mark.regression
 def test_uzf01_model_scope_nofile(function_tmpdir, example_data_path):
     sim_name = "uzf01"
-    netcdf = "nofile"
+    netcdf = ""
     fname = f"{sim_name}.structured.nc"
     data_path_base = example_data_path / "mf6" / "netcdf"
     ws = function_tmpdir / sim_name
@@ -108,7 +108,7 @@ def test_uzf01_model_scope_nofile(function_tmpdir, example_data_path):
 @pytest.mark.regression
 def test_uzf02_model_scope_nofile(function_tmpdir, example_data_path):
     sim_name = "uzf02"
-    netcdf = "nofile"
+    netcdf = ""
     fname = f"{sim_name}.input.nc"  # default
     data_path_base = example_data_path / "mf6" / "netcdf"
     ws = function_tmpdir / sim_name
@@ -230,7 +230,7 @@ def test_uzf02_sim_scope_fname(function_tmpdir, example_data_path):
 @pytest.mark.regression
 def test_uzf01_model_scope_nomesh(function_tmpdir, example_data_path):
     sim_name = "uzf01"
-    netcdf = "nofile"
+    netcdf = ""
     fname = f"{sim_name}.structured.nc"
     data_path_base = example_data_path / "mf6" / "netcdf"
     ws = function_tmpdir / sim_name
@@ -258,7 +258,7 @@ def test_uzf01_model_scope_nomesh(function_tmpdir, example_data_path):
 @pytest.mark.regression
 def test_uzf01_model_scope_mesh(function_tmpdir, example_data_path):
     sim_name = "uzf01"
-    netcdf = "nofile"
+    netcdf = ""
     mesh = "layered"
     fname = f"{sim_name}.layered.nc"
     data_path_base = example_data_path / "mf6" / "netcdf"
@@ -287,7 +287,7 @@ def test_uzf01_model_scope_mesh(function_tmpdir, example_data_path):
 @pytest.mark.regression
 def test_uzf02_model_scope(function_tmpdir, example_data_path):
     sim_name = "uzf02"
-    netcdf = "nofile"
+    netcdf = ""
     mesh = "layered"
     fname = f"{sim_name}.layered.nc"
     data_path_base = example_data_path / "mf6" / "netcdf"
@@ -335,11 +335,11 @@ def test_uzf01_pkg_scope(function_tmpdir, example_data_path):
     ds = gwf.modelgrid.dataset(modeltime=gwf.modeltime)
 
     # get model netcdf info
-    nc_info = gwf.netcdf_info()
+    nc_meta = gwf.netcdf_meta()
 
     # update dataset directly with required attributes
-    for a in nc_info["attrs"]:
-        ds.attrs[a] = nc_info["attrs"][a]
+    for a in nc_meta["attrs"]:
+        ds.attrs[a] = nc_meta["attrs"][a]
 
     # add all packages and update data
     for p in gwf.packagelist:
@@ -373,25 +373,25 @@ def test_uzf01_pkg_scope_modify(function_tmpdir, example_data_path):
     ds = gwf.modelgrid.dataset(modeltime=gwf.modeltime)
 
     # get model netcdf info
-    nc_info = gwf.netcdf_info()
+    nc_meta = gwf.netcdf_meta()
 
     # update dataset directly with required attributes
-    for a in nc_info["attrs"]:
-        ds.attrs[a] = nc_info["attrs"][a]
+    for a in nc_meta["attrs"]:
+        ds.attrs[a] = nc_meta["attrs"][a]
 
     # update dataset with `DIS` arrays
     dis = gwf.get_package("dis")
     ds = dis.update_dataset(ds)
 
     # get npf package netcdf info
     npf = gwf.get_package("npf")
-    nc_info = npf.netcdf_info()
+    nc_meta = npf.netcdf_meta()
 
     # update dataset with `NPF` arrays
     # change k varname and add attribute
-    nc_info["k"]["varname"] = "npf_k_updated"
-    nc_info["k"]["attrs"]["standard_name"] = "soil_hydraulic_conductivity_at_saturation"
-    ds = npf.update_dataset(ds, netcdf_info=nc_info)
+    nc_meta["k"]["varname"] = "npf_k_updated"
+    nc_meta["k"]["attrs"]["standard_name"] = "soil_hydraulic_conductivity_at_saturation"
+    ds = npf.update_dataset(ds, netcdf_meta=nc_meta)
 
     # ic
     ic = gwf.get_package("ic")

flopy/discretization/grid.py

@@ -1261,15 +1261,15 @@ def write_shapefile(self, filename="grid.shp", crs=None, prjfile=None, **kwargs)
         )
         return
 
-    def dataset(self, modeltime=None, mesh=None, encoding=None):
+    def dataset(self, modeltime=None, mesh=None, configuration=None):
         """
         Method to generate baseline xarray dataset
 
         Parameters
         ----------
         modeltime : FloPy ModelTime object
         mesh : mesh type
-        encoding : variable encoding dictionary
+        configuration : configuration dictionary
 
         Returns
         -------

flopy/discretization/structuredgrid.py

@@ -1676,12 +1676,12 @@ def get_plottable_layer_array(self, a, layer):
         assert plotarray.shape == required_shape, msg
         return plotarray
 
-    def dataset(self, modeltime=None, mesh=None, encoding=None):
+    def dataset(self, modeltime=None, mesh=None, configuration=None):
         """
         modeltime : FloPy ModelTime object
         mesh : mesh type
                valid mesh types are "layered" or None
-        encoding : variable encoding dictionary
+        configuration : configuration dictionary
         """
         from ..utils import import_optional_dependency
 
@@ -1694,11 +1694,11 @@ def dataset(self, modeltime=None, mesh=None, encoding=None):
         ds.attrs["modflow_grid"] = "STRUCTURED"
 
         if mesh and mesh.upper() == "LAYERED":
-            return self._layered_mesh_dataset(ds, modeltime, encoding)
+            return self._layered_mesh_dataset(ds, modeltime, configuration)
         elif mesh is None:
-            return self._structured_dataset(ds, modeltime, encoding)
+            return self._structured_dataset(ds, modeltime, configuration)
 
-    def _layered_mesh_dataset(self, ds, modeltime=None, encoding=None):
+    def _layered_mesh_dataset(self, ds, modeltime=None, configuration=None):
         FILLNA_INT32 = np.int32(-2147483647)
         FILLNA_DBL = 9.96920996838687e36
         lenunits = {0: "m", 1: "ft", 2: "m", 3: "m"}
@@ -1801,29 +1801,31 @@ def _layered_mesh_dataset(self, ds, modeltime=None, encoding=None):
         ds["mesh_face_nodes"].attrs["_FillValue"] = FILLNA_INT32
         ds["mesh_face_nodes"].attrs["start_index"] = np.int32(1)
 
-        if encoding is not None and "wkt" in encoding and encoding["wkt"] is not None:
-            ds = ds.assign({"projection": ([], np.int32(1))})
-            # wkt override to existing crs
-            ds["projection"].attrs["wkt"] = encoding["wkt"]
+        wkt_configured = (
+            configuration is not None
+            and "wkt" in configuration
+            and configuration["wkt"] is not None
+        )
+
+        if wkt_configured or self.crs is not None:
             ds["mesh_node_x"].attrs["grid_mapping"] = "projection"
             ds["mesh_node_y"].attrs["grid_mapping"] = "projection"
             ds["mesh_face_x"].attrs["grid_mapping"] = "projection"
             ds["mesh_face_y"].attrs["grid_mapping"] = "projection"
-        elif self.crs is not None:
             ds = ds.assign({"projection": ([], np.int32(1))})
-            ds["projection"].attrs["wkt"] = self.crs.to_wkt()
-            ds["mesh_node_x"].attrs["grid_mapping"] = "projection"
-            ds["mesh_node_y"].attrs["grid_mapping"] = "projection"
-            ds["mesh_face_x"].attrs["grid_mapping"] = "projection"
-            ds["mesh_face_y"].attrs["grid_mapping"] = "projection"
+            if wkt_configured:
+                # wkt override to existing crs
+                ds["projection"].attrs["wkt"] = configuration["wkt"]
+            else:
+                ds["projection"].attrs["wkt"] = self.crs.to_wkt()
 
         return ds
 
-    def _structured_dataset(self, ds, modeltime=None, encoding=None):
+    def _structured_dataset(self, ds, modeltime=None, configuration=None):
         lenunits = {0: "m", 1: "ft", 2: "m", 3: "m"}
 
-        x = self.xoffset + self.xycenters[0]
-        y = self.yoffset + self.xycenters[1]
+        xc = self.xoffset + self.xycenters[0]
+        yc = self.yoffset + self.xycenters[1]
         z = [float(x) for x in range(1, self.nlay + 1)]
 
         # set coordinate var bounds
@@ -1849,8 +1851,8 @@ def _structured_dataset(self, ds, modeltime=None, encoding=None):
         var_d = {
             "time": (["time"], modeltime.totim),
             "z": (["z"], z),
-            "y": (["y"], y),
-            "x": (["x"], x),
+            "y": (["y"], yc),
+            "x": (["x"], xc),
         }
         ds = ds.assign(var_d)
 
@@ -1876,17 +1878,78 @@ def _structured_dataset(self, ds, modeltime=None, encoding=None):
         ds["x"].attrs["long_name"] = "Easting"
         ds["x"].attrs["bounds"] = "x_bnds"
 
-        if encoding is not None and "wkt" in encoding and encoding["wkt"] is not None:
-            ds = ds.assign({"projection": ([], np.int32(1))})
-            # wkt override to existing crs
-            ds["projection"].attrs["crs_wkt"] = encoding["wkt"]
+        latlon_cfg = (
+            configuration is not None
+            and "latitude" in configuration
+            and configuration["latitude"] is not None
+            and "longitude" in configuration
+            and configuration["longitude"] is not None
+        )
+
+        if latlon_cfg or self.crs is not None:
+            if latlon_cfg:
+                lats = configuration["latitude"]
+                lons = configuration["longitude"]
+            else:
+                try:
+                    import warnings
+
+                    from pyproj import Proj
+
+                    epsg_code = self.crs.to_epsg(min_confidence=90)
+                    proj = Proj(
+                        f"EPSG:{epsg_code}",
+                    )
+
+                    lats = []
+                    lons = []
+                    x_local = []
+                    y_local = []
+                    for y in self.xycenters[1]:
+                        for x in self.xycenters[0]:
+                            x_local.append(x)
+                            y_local.append(y)
+
+                    x_global, y_global = self.get_coords(x_local, y_local)
+
+                    for i, x in enumerate(x_global):
+                        lon, lat = proj(x, y_global[i], inverse=True)
+                        lats.append(lat)
+                        lons.append(lon)
+
+                    lats = np.array(lats)
+                    lons = np.array(lons)
+
+                except Exception as e:
+                    warnings.warn(
+                        f"Cannot create coordinates from CRS: {e}",
+                        UserWarning,
+                    )
+
+            # create coordinate vars
+            var_d = {
+                "lat": (["y", "x"], lats.reshape(yc.size, xc.size)),
+                "lon": (["y", "x"], lons.reshape(yc.size, xc.size)),
+            }
+            ds = ds.assign(var_d)
+
+            # set coordinate attributes
+            ds["lat"].attrs["units"] = "degrees_north"
+            ds["lat"].attrs["standard_name"] = "latitude"
+            ds["lat"].attrs["long_name"] = "latitude"
+            ds["lon"].attrs["units"] = "degrees_east"
+            ds["lon"].attrs["standard_name"] = "longitude"
+            ds["lon"].attrs["long_name"] = "longitude"
+
+        elif (
+            configuration is not None
+            and "wkt" in configuration
+            and configuration["wkt"] is not None
+        ):
             ds["x"].attrs["grid_mapping"] = "projection"
             ds["y"].attrs["grid_mapping"] = "projection"
-        elif self.crs is not None:
             ds = ds.assign({"projection": ([], np.int32(1))})
-            ds["projection"].attrs["crs_wkt"] = self.crs.to_wkt()
-            ds["x"].attrs["grid_mapping"] = "projection"
-            ds["y"].attrs["grid_mapping"] = "projection"
+            ds["projection"].attrs["crs_wkt"] = configuration["wkt"]
 
         return ds