enhanced grid support for netcdf

Author: mjreno

flopy/discretization/grid.py

@@ -1185,12 +1185,14 @@ def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
         yul = None
         if os.path.exists(reffile):
             with open(reffile) as input:
+                print(f"Updating grid based on reference: {reffile}")
                 for line in input:
                     if len(line) > 1:
                         if line.strip()[0] != "#":
                             info = line.strip().split("#")[0].split()
                             if len(info) > 1:
                                 data = " ".join(info[1:]).strip("'").strip('"')
+                                print(f"Grid update on reference: {info}")
                                 if info[0] == "xll":
                                     self._xoff = float(data)
                                 elif info[0] == "yll":
@@ -1201,12 +1203,14 @@ def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
                                     yul = float(data)
                                 elif info[0] == "rotation":
                                     self._angrot = float(data)
-                                elif info[0] == "epsg":
+                                elif info[0] == "length_units":
+                                    self._units = data.lower()
+                                elif info[0].lower() == "epsg":
                                     self.epsg = int(data)
                                 elif info[0] == "proj4":
                                     self.crs = data
-                                elif info[0] == "start_date":
-                                    start_datetime = data
+                                else:
+                                    print("   ->warn: update not applied.")
 
             # model must be rotated first, before setting xoff and yoff
             # when xul and yul are provided.

flopy/discretization/modeltime.py

@@ -1,4 +1,5 @@
 import calendar
+import os
 from dataclasses import dataclass, field
 from datetime import datetime, timedelta
 from difflib import SequenceMatcher
@@ -802,3 +803,39 @@ def reverse(self) -> "ModelTime":
             self.start_datetime,
             self.steady_state[::-1] if self.steady_state is not None else None,
         )
+
+    def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
+        """read spatial reference info from the usgs.model.reference file
+        https://water.usgs.gov/ogw/policy/gw-model/modelers-setup.html"""
+        if os.path.exists(reffile):
+            start_date_time = ""
+            with open(reffile) as input:
+                print(f"Updating modeltime based on reference: {reffile}")
+                for line in input:
+                    if len(line) > 1:
+                        if line.strip()[0] != "#":
+                            info = line.strip().split("#")[0].split()
+                            if len(info) > 1:
+                                data = " ".join(info[1:]).strip("'").strip('"')
+                                print(f"ModelTime update on reference: {info}")
+                                if info[0].lower() == "time_units":
+                                    self.time_units = data.lower()
+                                elif info[0] == "start_date":
+                                    if len(start_date_time) > 0:
+                                        start_date_time = f"{data} {start_date_time}"
+                                    else:
+                                        start_date_time = data
+                                elif info[0] == "start_time":
+                                    if len(start_date_time) > 0:
+                                        start_date_time = f"{start_date_time} {data}"
+                                    else:
+                                        start_date_time = data
+                                else:
+                                    print("   ->warn: update not applied.")
+
+            if len(start_date_time) > 0:
+                self.start_datetime = start_date_time
+
+            return True
+        else:
+            return False

flopy/discretization/structuredgrid.py

@@ -1770,6 +1770,46 @@ def get_plottable_layer_array(self, a, layer):
         assert plotarray.shape == required_shape, msg
         return plotarray
 
+    def latlon(self):
+        try:
+            import warnings
+
+            from pyproj import Proj
+
+            epsg = None
+            if self.crs is not None:
+                epsg = self.crs.to_epsg()
+
+            proj = Proj(
+                f"EPSG:{epsg}",
+            )
+
+            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)
+
+            return np.array(lats), np.array(lons)
+
+        except Exception as e:
+            warnings.warn(
+                f"Cannot create coordinates from CRS: {e}",
+                UserWarning,
+            )
+
+            return None, None
+
     def dataset(self, modeltime=None, mesh=None, configuration=None):
         """
         modeltime : FloPy ModelTime object
@@ -1803,7 +1843,9 @@ def _layered_mesh_dataset(self, ds, modeltime=None, configuration=None):
         }
         ds = ds.assign(var_d)
         ds["time"].attrs["calendar"] = "standard"
-        ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
+        ds["time"].attrs["units"] = (
+            f"{modeltime.time_units} since {modeltime.start_datetime}"
+        )
         ds["time"].attrs["axis"] = "T"
         ds["time"].attrs["standard_name"] = "time"
         ds["time"].attrs["long_name"] = "time"
@@ -1955,7 +1997,9 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
         ds = ds.assign(var_d)
 
         ds["time"].attrs["calendar"] = "standard"
-        ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
+        ds["time"].attrs["units"] = (
+            f"{modeltime.time_units} since {modeltime.start_datetime}"
+        )
         ds["time"].attrs["axis"] = "T"
         ds["time"].attrs["standard_name"] = "time"
         ds["time"].attrs["long_name"] = "time"
@@ -1980,46 +2024,13 @@ def _structured_dataset(self, ds, modeltime=None, configuration=None):
             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,
-                    )
+        if latlon_cfg:
+            lats = configuration["latitude"]
+            lons = configuration["longitude"]
+        elif self.crs is not None:
+            lats, lons = self.latlon()
 
+        if lats is not None and lons is not None:
             # create coordinate vars
             var_d = {
                 "lat": (["y", "x"], lats.reshape(yc.size, xc.size)),

flopy/discretization/vertexgrid.py

@@ -631,7 +631,9 @@ def dataset(self, modeltime=None, mesh=None, configuration=None):
         }
         ds = ds.assign(var_d)
         ds["time"].attrs["calendar"] = "standard"
-        ds["time"].attrs["units"] = f"days since {modeltime.start_datetime}"
+        ds["time"].attrs["units"] = (
+            f"{modeltime.time_units} since {modeltime.start_datetime}"
+        )
         ds["time"].attrs["axis"] = "T"
         ds["time"].attrs["standard_name"] = "time"
         ds["time"].attrs["long_name"] = "time"

flopy/mf6/mfmodel.py

@@ -458,14 +458,19 @@ def modelgrid(self):
                 if idomain is None:
                     force_resync = True
                     idomain = self._resolve_idomain(idomain, botm)
+                crs = self._modelgrid.crs
+                if crs is None and hasattr(dis, "crs"):
+                    crs = dis.crs.get_data()
+                    if crs is not None:
+                        crs = crs[0][1]
                 self._modelgrid = StructuredGrid(
                     delc=dis.delc.array,
                     delr=dis.delr.array,
                     top=dis.top.array,
                     botm=botm,
                     idomain=idomain,
                     lenuni=dis.length_units.array,
-                    crs=self._modelgrid.crs,
+                    crs=crs,
                     xoff=self._modelgrid.xoffset,
                     yoff=self._modelgrid.yoffset,
                     angrot=self._modelgrid.angrot,
@@ -496,14 +501,19 @@ def modelgrid(self):
                 if idomain is None:
                     force_resync = True
                     idomain = self._resolve_idomain(idomain, botm)
+                crs = self._modelgrid.crs
+                if crs is None and hasattr(dis, "crs"):
+                    crs = dis.crs.get_data()
+                    if crs is not None:
+                        crs = crs[0][1]
                 self._modelgrid = VertexGrid(
                     vertices=dis.vertices.array,
                     cell2d=dis.cell2d.array,
                     top=dis.top.array,
                     botm=botm,
                     idomain=idomain,
                     lenuni=dis.length_units.array,
-                    crs=self._modelgrid.crs,
+                    crs=crs,
                     xoff=self._modelgrid.xoffset,
                     yoff=self._modelgrid.yoffset,
                     angrot=self._modelgrid.angrot,
@@ -525,6 +535,11 @@ def modelgrid(self):
             idomain = dis.idomain.array
             if idomain is None:
                 idomain = np.ones(dis.nodes.array, dtype=int)
+            crs = self._modelgrid.crs
+            if crs is None and hasattr(dis, "crs"):
+                crs = dis.crs.get_data()
+                if crs is not None:
+                    crs = crs[0][1]
             if cell2d is None:
                 if (
                     self.simulation.simulation_data.verbosity_level.value
@@ -557,7 +572,7 @@ def modelgrid(self):
                 idomain=idomain,
                 lenuni=dis.length_units.array,
                 ncpl=ncpl,
-                crs=self._modelgrid.crs,
+                crs=crs,
                 xoff=self._modelgrid.xoffset,
                 yoff=self._modelgrid.yoffset,
                 angrot=self._modelgrid.angrot,
@@ -590,14 +605,19 @@ def modelgrid(self):
                 if idomain is None:
                     force_resync = True
                     idomain = self._resolve_idomain(idomain, botm)
+                crs = self._modelgrid.crs
+                if crs is None and hasattr(dis, "crs"):
+                    crs = dis.crs.get_data()
+                    if crs is not None:
+                        crs = crs[0][1]
                 self._modelgrid = VertexGrid(
                     vertices=dis.vertices.array,
                     cell1d=dis.cell1d.array,
                     top=None,
                     botm=botm,
                     idomain=idomain,
                     lenuni=dis.length_units.array,
-                    crs=self._modelgrid.crs,
+                    crs=crs,
                     xoff=self._modelgrid.xoffset,
                     yoff=self._modelgrid.yoffset,
                     angrot=self._modelgrid.angrot,
@@ -628,14 +648,19 @@ def modelgrid(self):
                 if idomain is None:
                     force_resync = True
                     idomain = self._resolve_idomain(idomain, botm)
+                crs = self._modelgrid.crs
+                if crs is None and hasattr(dis, "crs"):
+                    crs = dis.crs.get_data()
+                    if crs is not None:
+                        crs = crs[0][1]
                 self._modelgrid = StructuredGrid(
                     delc=dis.delc.array,
                     delr=dis.delr.array,
                     top=None,
                     botm=botm,
                     idomain=idomain,
                     lenuni=dis.length_units.array,
-                    crs=self._modelgrid.crs,
+                    crs=crs,
                     xoff=self._modelgrid.xoffset,
                     yoff=self._modelgrid.yoffset,
                     angrot=self._modelgrid.angrot,
@@ -666,14 +691,19 @@ def modelgrid(self):
                 if idomain is None:
                     force_resync = True
                     idomain = self._resolve_idomain(idomain, botm)
+                crs = self._modelgrid.crs
+                if crs is None and hasattr(dis, "crs"):
+                    crs = dis.crs.get_data()
+                    if crs is not None:
+                        crs = crs[0][1]
                 self._modelgrid = VertexGrid(
                     vertices=dis.vertices.array,
                     cell2d=dis.cell2d.array,
                     top=None,
                     botm=botm,
                     idomain=idomain,
                     lenuni=dis.length_units.array,
-                    crs=self._modelgrid.crs,
+                    crs=crs,
                     xoff=self._modelgrid.xoffset,
                     yoff=self._modelgrid.yoffset,
                     angrot=self._modelgrid.angrot,
@@ -1361,13 +1391,6 @@ def write(
             if write_netcdf:
                 # set data storage to write ascii for netcdf
                 pp._set_netcdf_storage()
-            if (
-                pp.package_type.startswith("dis")
-                and hasattr(pp, "crs")
-            ):
-                crs = pp.crs.get_data()
-                if crs is not None and self.modelgrid.crs is None:
-                    self.modelgrid.crs = crs[0][1]
             if (
                 self.simulation_data.verbosity_level.value
                 >= VerbosityLevel.normal.value
@@ -2366,7 +2389,7 @@ def _write_netcdf(self, mesh=None):
         nc_meta["attrs"]["history"] = f"first created {timestamp}"
         if mesh is None:
             nc_meta["attrs"]["Conventions"] = "CF-1.11"
-        elif mesh.upper() is "LAYERED":
+        elif mesh.upper() == "LAYERED":
             nc_meta["attrs"]["Conventions"] = "CF-1.11 UGRID-1.0"
 
         ds = self.update_dataset(