fix(modelgrid): retain crs data from classic nam files

autotest/test_mfreadnam.py

@@ -1,7 +1,10 @@
 import pytest
 from autotest.conftest import get_example_data_path
 
-from flopy.utils.mfreadnam import get_entries_from_namefile
+from flopy.utils.mfreadnam import (
+    attribs_from_namfile_header,
+    get_entries_from_namefile,
+)
 
 _example_data_path = get_example_data_path()
 
@@ -39,3 +42,68 @@ def test_get_entries_from_namefile_mf2005(path):
     entry = entries[0]
     assert path.parent.name in entry[0]
     assert entry[1] == package
+
+
+@pytest.mark.parametrize(
+    "path,expected",
+    [
+        pytest.param(
+            None,
+            {
+                "crs": None,
+                "rotation": 0.0,
+                "xll": None,
+                "xul": None,
+                "yll": None,
+                "yul": None,
+            },
+            id="None",
+        ),
+        pytest.param(
+            _example_data_path / "freyberg" / "freyberg.nam",
+            {
+                "crs": None,
+                "rotation": 0.0,
+                "xll": None,
+                "xul": None,
+                "yll": None,
+                "yul": None,
+            },
+            id="freyberg",
+        ),
+        pytest.param(
+            _example_data_path
+            / "freyberg_multilayer_transient"
+            / "freyberg.nam",
+            {
+                "crs": "+proj=utm +zone=14 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
+                "rotation": 15.0,
+                "start_datetime": "1/1/2015",
+                "xll": None,
+                "xul": 619653.0,
+                "yll": None,
+                "yul": 3353277.0,
+            },
+            id="freyberg_multilayer_transient",
+        ),
+        pytest.param(
+            _example_data_path
+            / "mt3d_test"
+            / "mfnwt_mt3dusgs"
+            / "sft_crnkNic"
+            / "CrnkNic.nam",
+            {
+                "crs": "EPSG:26916",
+                "rotation": 0.0,
+                "start_datetime": "1-1-1970",
+                "xll": None,
+                "xul": 0.0,
+                "yll": None,
+                "yul": 15.0,
+            },
+            id="CrnkNic",
+        ),
+    ],
+)
+def test_attribs_from_namfile_header(path, expected):
+    assert attribs_from_namfile_header(path) == expected

autotest/test_modflow.py

@@ -7,6 +7,7 @@
 import numpy as np
 import pytest
 from autotest.conftest import get_example_data_path
+from modflow_devtools.misc import has_pkg
 from modflow_devtools.markers import excludes_platform, requires_exe
 
 from flopy.discretization import StructuredGrid
@@ -29,6 +30,8 @@
 from flopy.seawat import Seawat
 from flopy.utils import Util2d
 
+_example_data_path = get_example_data_path()
+
 
 @pytest.fixture
 def model_reference_path(example_data_path) -> Path:
@@ -76,6 +79,64 @@ def test_modflow_load(namfile, example_data_path):
     assert model.model_ws == str(mpath)
 
 
+@pytest.mark.parametrize(
+    "path,expected",
+    [
+        pytest.param(
+            _example_data_path / "freyberg" / "freyberg.nam",
+            {
+                "crs": None,
+                "epsg": None,
+                "angrot": 0.0,
+                "xoffset": 0.0,
+                "yoffset": 0.0,
+            },
+            id="freyberg",
+        ),
+        pytest.param(
+            _example_data_path
+            / "freyberg_multilayer_transient"
+            / "freyberg.nam",
+            {
+                "proj4": "+proj=utm +zone=14 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
+                "angrot": 15.0,
+                "xoffset": 622241.1904510253,
+                "yoffset": 3343617.741737109,
+            },
+            id="freyberg_multilayer_transient",
+        ),
+        pytest.param(
+            _example_data_path
+            / "mt3d_test"
+            / "mfnwt_mt3dusgs"
+            / "sft_crnkNic"
+            / "CrnkNic.nam",
+            {
+                "epsg": 26916,
+                "angrot": 0.0,
+                "xoffset": 0.0,
+                "yoffset": 0.0,
+            },
+            id="CrnkNic",
+        ),
+    ],
+)
+def test_modflow_load_modelgrid(path, expected):
+    """Check modelgrid metadata from NAM file."""
+    model = Modflow.load(path.name, model_ws=path.parent, load_only=[])
+    modelgrid = model.modelgrid
+    for key, expected_value in expected.items():
+        if key == "proj4" and has_pkg("pyproj"):
+            # skip since pyproj will usually restructure proj4 attribute
+            # otherwise continue test without pyproj, as it should be preserved
+            continue
+        modelgrid_value = getattr(modelgrid, key)
+        if isinstance(modelgrid_value, float):
+            assert modelgrid_value == pytest.approx(expected_value), key
+        else:
+            assert modelgrid_value == expected_value, key
+
+
 def test_modflow_load_when_nam_dne():
     with pytest.raises(OSError):
         Modflow.load("nonexistent.nam", check=False)
@@ -551,12 +612,12 @@ def test_read_usgs_model_reference(function_tmpdir, model_reference_path):
 
 
 def mf2005_model_namfiles():
-    path = get_example_data_path() / "mf2005_test"
+    path = _example_data_path / "mf2005_test"
     return [str(p) for p in path.glob("*.nam")]
 
 
 def parameters_model_namfiles():
-    path = get_example_data_path() / "parameters"
+    path = _example_data_path / "parameters"
     skip = ["twrip.nam", "twrip_upw.nam"]  # TODO: why do these fail?
     return [str(p) for p in path.glob("*.nam") if p.name not in skip]
 
@@ -744,15 +805,15 @@ def test_mflist_add_record():
     np.testing.assert_array_equal(wel.stress_period_data[1], check1)
 
 
-__mf2005_test_path = get_example_data_path() / "mf2005_test"
-__mf2005_namfiles = [
-    Path(__mf2005_test_path) / f
-    for f in __mf2005_test_path.rglob("*")
+_mf2005_test_path = _example_data_path / "mf2005_test"
+_mf2005_namfiles = [
+    Path(_mf2005_test_path) / f
+    for f in _mf2005_test_path.rglob("*")
     if f.suffix == ".nam"
 ]
 
 
-@pytest.mark.parametrize("namfile", __mf2005_namfiles)
+@pytest.mark.parametrize("namfile", _mf2005_namfiles)
 def test_checker_on_load(namfile):
     # load all of the models in the mf2005_test folder
     # model level checks are performed by default on load()
@@ -770,7 +831,7 @@ def test_checker_on_load(namfile):
 
 @pytest.mark.parametrize("str_path", [True, False])
 def test_manual_check(function_tmpdir, str_path):
-    namfile_path = __mf2005_namfiles[0]
+    namfile_path = _mf2005_namfiles[0]
     summary_path = function_tmpdir / "summary"
     model = Modflow.load(namfile_path, model_ws=namfile_path.parent)
     model.change_model_ws(function_tmpdir)

flopy/discretization/grid.py

@@ -6,6 +6,13 @@
 
 import numpy as np
 
+try:
+    import pyproj
+
+    HAS_PYPROJ = True
+except ImportError:
+    HAS_PYPROJ = False
+
 from ..utils import geometry
 from ..utils.crs import get_crs
 from ..utils.gridutil import get_lni
@@ -204,15 +211,15 @@ def __init__(
             raise TypeError(f"unhandled keywords: {kwargs}")
         if prjfile is not None:
             self.prjfile = prjfile
-        try:
+        if HAS_PYPROJ:
             self._crs = get_crs(**get_crs_args)
-        except ImportError:
-            # provide some support without pyproj by retaining 'epsg' integer
-            if getattr(self, "_epsg", None) is None:
-                epsg = _get_epsg_from_crs_or_proj4(crs, self.proj4)
-                if epsg is not None:
-                    self.epsg = epsg
-
+        elif crs is not None:
+            # provide some support without pyproj
+            if isinstance(crs, str) and self.proj4 is None:
+                self._proj4 = crs
+            if self.epsg is None:
+                if epsg := _get_epsg_from_crs_or_proj4(crs, self.proj4):
+                    self._epsg = epsg
         self._prjfile = prjfile
         self._xoff = xoff
         self._yoff = yoff
@@ -304,9 +311,9 @@ def crs(self, crs):
         if crs is None:
             self._crs = None
             return
-        try:
+        if HAS_PYPROJ:
             self._crs = get_crs(crs=crs)
-        except ImportError:
+        else:
             warnings.warn(
                 "cannot set 'crs' property without pyproj; "
                 "try setting 'epsg' or 'proj4' instead",
@@ -336,11 +343,8 @@ def epsg(self, epsg):
             raise ValueError("epsg property must be an int or None")
         self._epsg = epsg
         # If crs was previously unset, use EPSG code
-        if self._crs is None and epsg is not None:
-            try:
-                self._crs = get_crs(crs=epsg)
-            except ImportError:
-                pass
+        if HAS_PYPROJ and self._crs is None and epsg is not None:
+            self._crs = get_crs(crs=epsg)
 
     @property
     def proj4(self):
@@ -364,11 +368,8 @@ def proj4(self, proj4):
             raise ValueError("proj4 property must be a str or None")
         self._proj4 = proj4
         # If crs was previously unset, use lossy PROJ string
-        if self._crs is None and proj4 is not None:
-            try:
-                self._crs = get_crs(crs=proj4)
-            except ImportError:
-                pass
+        if HAS_PYPROJ and self._crs is None and proj4 is not None:
+            self._crs = get_crs(crs=proj4)
 
     @property
     def prj(self):
@@ -402,10 +403,10 @@ def prjfile(self, prjfile):
             raise ValueError("prjfile property must be str, PathLike or None")
         self._prjfile = prjfile
         # If crs was previously unset, use .prj file input
-        if self._crs is None:
+        if HAS_PYPROJ and self._crs is None:
             try:
                 self._crs = get_crs(prjfile=prjfile)
-            except (ImportError, FileNotFoundError):
+            except FileNotFoundError:
                 pass
 
     @property
@@ -1007,14 +1008,14 @@ def set_coord_info(
         # Handle deprecated projection kwargs; warnings are raised in crs.py
         get_crs_args = {"crs": crs, "prjfile": prjfile}
         if "epsg" in kwargs:
-            self.epsg = get_crs_args["epsg"] = kwargs.pop("epsg")
+            self._epsg = get_crs_args["epsg"] = kwargs.pop("epsg")
         if "proj4" in kwargs:
-            self.proj4 = get_crs_args["proj4"] = kwargs.pop("proj4")
+            self._proj4 = get_crs_args["proj4"] = kwargs.pop("proj4")
         if kwargs:
             raise TypeError(f"unhandled keywords: {kwargs}")
-        try:
+        if HAS_PYPROJ:
             new_crs = get_crs(**get_crs_args)
-        except ImportError:
+        else:
             new_crs = None
             # provide some support without pyproj by retaining 'epsg' integer
             if getattr(self, "_epsg", None) is None:

flopy/mbase.py

@@ -410,6 +410,12 @@ def __init__(
         self._crs = kwargs.pop("crs", None)
         self._start_datetime = kwargs.pop("start_datetime", "1-1-1970")
 
+        if kwargs:
+            warnings.warn(
+                f"unhandled keywords: {kwargs}",
+                category=UserWarning,
+            )
+
         # build model discretization objects
         self._modelgrid = Grid(
             crs=self._crs,

flopy/modflow/mf.py

@@ -278,9 +278,14 @@ def modelgrid(self):
             ibound = self.bas6.ibound.array
         else:
             ibound = None
-
+        # take the first non-None entry
+        crs = (
+            self._modelgrid.crs
+            or self._modelgrid.proj4
+            or self._modelgrid.epsg
+        )
         common_kwargs = {
-            "crs": self._modelgrid.crs or self._modelgrid.epsg,
+            "crs": crs,
             "xoff": self._modelgrid.xoffset,
             "yoff": self._modelgrid.yoffset,
             "angrot": self._modelgrid.angrot,

flopy/utils/mfreadnam.py

@@ -205,7 +205,17 @@ def parsenamefile(namfilename, packages, verbose=True):
 
 
 def attribs_from_namfile_header(namefile):
-    # check for reference info in the nam file header
+    """Return spatial and temporal reference info from the nam header.
+
+    Parameters
+    ----------
+    namefile : str, PathLike or None
+        Path to NAM file to read.
+
+    Returns
+    -------
+    dict
+    """
     defaults = {
         "xll": None,
         "yll": None,
@@ -255,7 +265,6 @@ def attribs_from_namfile_header(namefile):
             except:
                 print(f"   could not parse rotation in {namefile}")
         elif "proj4_str" in item.lower():
-            # deprecated, use "crs" instead
             try:
                 proj4 = ":".join(item.split(":")[1:]).strip()
                 if proj4.lower() == "none":
@@ -277,6 +286,9 @@ def attribs_from_namfile_header(namefile):
                 defaults["start_datetime"] = start_datetime
             except:
                 print(f"   could not parse start in {namefile}")
+    if "proj4_str" in defaults and defaults["crs"] is None:
+        # handle deprecated keyword, use "crs" instead
+        defaults["crs"] = defaults.pop("proj4_str")
     return defaults