netcdf attr, coordinate updates

flopy4/mf6/netcdf.py

@@ -24,22 +24,39 @@
 from flopy4.version import __version__
 
 
-def _cf_var_attrs(dims: list[str], mesh: str | None, grid) -> dict:
+def _cf_var_attrs(dims: list[str], mesh: str | None, grid, layer: int | None = None) -> dict:
     """Return {"attrs": {...}, "encoding": {...}} for a NetCDF data variable.
 
     In mesh context x/y are abstract row/col indices, not geographic — only nmesh_face
     vars get coordinates linking. In structured context x/y ARE geographic.
+
+    layer is the 1-based layer number for a per-layer variable, or None
+    otherwise. z_l{layer} is appended to coordinates only when set -- the
+    only case it shares nmesh_face with the variable (CF-1.13 5.2).
     """
     attrs: dict[str, str] = {}
     encoding: dict[str, object] = {}
     has_crs = grid is not None and getattr(grid, "crs", None) is not None
+    # grid_mapping is CF-required only on variables spanning the full 2D
+    # horizontal spatial extent (CF-1.13 5.6); 1D dimension-definition arrays
+    # such as dis_delr/dis_delc are not georeferenced fields and must not
+    # carry it, matching MF6's own ncvar_gridmap (DisNCStructured.f90).
+    has_full_extent = ("x" in dims and "y" in dims) or "nmesh_face" in dims
 
-    if has_crs:
+    if has_crs and has_full_extent:
         attrs["grid_mapping"] = "projection"
     if "nmesh_face" in dims:
-        attrs["coordinates"] = "mesh_face_x mesh_face_y"
+        if layer:
+            attrs["coordinates"] = f"mesh_face_x mesh_face_y z_l{layer}"
+        else:
+            attrs["coordinates"] = "mesh_face_x mesh_face_y"
         attrs["mesh"] = "mesh"
         attrs["location"] = "face"
+    elif "layer" in dims:
+        # Structured per-layer fields link to the z (cell center elevation)
+        # coordinate, matching MF6's own DisNCStructured.f90 -- non-layered
+        # fields (e.g. dis_top, a y/x-only 2D array) do not.
+        attrs["coordinates"] = "z"
 
     return {"attrs": attrs, "encoding": encoding}
 
@@ -209,12 +226,22 @@ def from_model(
             raise ValueError("model must have a 'name' attribute")
 
         modeltype = model.__class__.__name__.lower()
-        attrs = {"title": f"{model.name.upper()} model input"}
+        # Matches MF6's own NCModel.f90 title convention (model-type-specific
+        # fragment + " array input" for the pre-solve/validate-mode INPUT
+        # file this class always produces); models MF6 itself doesn't
+        # support for NetCDF export (anything but GWF/GWT/GWE) fall back to
+        # a generic title.
+        _title_fragment = {"gwf": "hydraulic head", "gwt": "concentration", "gwe": "temperature"}
+        if modeltype in _title_fragment:
+            title = f"{model.name.upper()} {_title_fragment[modeltype]} array input"
+        else:
+            title = f"{model.name.upper()} model input"
+        attrs = {"title": title}
         packages = []
         distype = None
 
         if netcdf_format == NetCDFFormat.LAYERED_MESH:
-            attrs["mesh"] = NetCDFFormat.LAYERED_MESH.value
+            attrs["modflow_mesh"] = NetCDFFormat.LAYERED_MESH.value
 
         # Resolve nper from time arg, simulation tdis, or model's data dims.
         if time is not None:
@@ -260,11 +287,12 @@ def from_model(
                     val = getattr(package, f.name)
                     if val is None:
                         continue
-                    arr = np.asarray(val, dtype=np.float64)
+                    _dtype = _PKG_DTYPE_MAP.get(to_field_type(f.type), np.float64)
+                    arr = np.asarray(val, dtype=_dtype)
                     # Only broadcast scalars to full grid for nodes-shaped fields
                     shape_meta = f.metadata.get("shape", ())
                     if "nodes" in shape_meta and arr.size < _nodes:
-                        arr = np.full(_nodes, float(arr.ravel()[0]))
+                        arr = np.full(_nodes, arr.ravel()[0], dtype=_dtype)
                     p["params"].append({"name": f.name, "data": arr})
                 elif f.metadata.get("fk"):
                     # dynamically named arrays (RCHA's aux): one param per auxiliary name
@@ -285,6 +313,14 @@ def from_model(
                                 "data": dense(named_periods, _nper, carry_forward=False),
                             }
                         )
+                elif (
+                    f.metadata.get("block") == "period" and f.metadata.get("reader") == "readarray"
+                ):
+                    val = getattr(package, f.name)
+                    if val is None:
+                        continue
+                    _dtype = _PKG_DTYPE_MAP.get(to_field_type(f.type), np.float64)
+                    p["params"].append({"name": f.name, "data": np.asarray(val, dtype=_dtype)})
                 else:
                     # period arrays: {kper: array}, periods not given filled.
                     # Time-array series references stay in the package file.
@@ -345,12 +381,12 @@ def to_xarray(self) -> xr.Dataset:
         meta = self.model_dump(by_alias=True)
 
         if self._grid is not None and self._time is not None:  # type: ignore
-            conventions = "CF-1.11"  # type: ignore
-            if meta["attrs"]["mesh"] is not None:
+            conventions = "CF-1.13"  # type: ignore
+            if meta["attrs"]["modflow_mesh"] is not None:
                 conventions = f"{conventions} UGRID-1.0"
             _fmt = (
                 NetCDFFormat.LAYERED_MESH
-                if meta["attrs"]["mesh"] is not None
+                if meta["attrs"]["modflow_mesh"] is not None
                 else NetCDFFormat.STRUCTURED
             )
             dss.append(self._grid.to_xarray(netcdf_format=_fmt, modeltime=self._time))
@@ -413,7 +449,13 @@ def validate_attrs(cls, v: dict[str, str]) -> dict[str, str]:
         """
         validate model (dataset) scoped attributes dictionary
         """
+        # title is free-text (e.g. "GWFMODEL hydraulic head array input") and must
+        # keep its original casing, matching MF6 -- unlike the other keys/values
+        # here, which are lowercase identifiers by convention.
+        title = v.get("title")
         v = lower(v)
+        if title is not None:
+            v["title"] = title
         return v
 
     @staticmethod
@@ -431,7 +473,9 @@ def _backfill_meta(meta: dict, context: dict, verbose: bool = True) -> dict:
 
         _packages = []
         for pkg in _meta["packages"]:
-            pkgctx = {"mesh": _meta["attrs"]["mesh"]} if "mesh" in _meta["attrs"] else {}
+            pkgctx = (
+                {"mesh": _meta["attrs"]["modflow_mesh"]} if "modflow_mesh" in _meta["attrs"] else {}
+            )
             pkgctx["modelname"] = _meta["modelname"]
             pkgctx["gridtype"] = _meta["gridtype"]
             pkgctx |= context
@@ -443,17 +487,17 @@ def _backfill_meta(meta: dict, context: dict, verbose: bool = True) -> dict:
 
 class NetCDFModelAttrs(BaseModel):
     # order of params dictates when data added to info dict
-    mesh: str | None = Field(default=None)
+    modflow_mesh: str | None = Field(default=None)
     modflow_grid: str = Field()
     modflow_model: str = Field()
 
     model_config = ConfigDict(extra="allow")
 
-    @field_validator("mesh", mode="before")
+    @field_validator("modflow_mesh", mode="before")
     @classmethod
-    def validate_mesh(cls, v: str | None, info: ValidationInfo) -> str | None:
+    def validate_modflow_mesh(cls, v: str | None, info: ValidationInfo) -> str | None:
         """
-        validate model mesh attribute
+        validate model modflow_mesh attribute
         """
         if v is not None:
             if v.lower() != "layered":
@@ -728,6 +772,7 @@ def to_xarray(self) -> xr.Dataset:
             [str(d) for d in ds[varname].dims],
             mesh,
             self._context.get("grid"),
+            layer=meta["attrs"].get("layer"),
         )
         ds[varname].attrs.update(cf["attrs"])
         ds[varname].encoding.update(cf["encoding"])

flopy4/mf6/utils/crs.py

@@ -0,0 +1,149 @@
+from typing import Any
+
+import numpy as np
+from pyproj import CRS
+
+# CF-standard numeric grid_mapping parameters shared across all three
+# supported projected grid_mapping_name values.
+_SHARED_CF_PARAMS = (
+    "longitude_of_central_meridian",
+    "latitude_of_projection_origin",
+    "false_easting",
+    "false_northing",
+    "semi_major_axis",
+    "inverse_flattening",
+)
+
+
+def cf_grid_mapping_params(crs: CRS) -> tuple[dict[str, Any], str | None]:
+    """
+    Return CF-standard numeric grid_mapping parameters for *crs*, in
+    addition to wkt/crs_wkt: some netCDF readers -- notably ArcGIS's
+    classic netCDF connector -- position a projected grid using these
+    individual attributes rather than parsing wkt/crs_wkt text, and
+    default every parameter to 0 if absent.
+
+    Covers transverse_mercator, lambert_conformal_conic (2SP only -- CF
+    has no scale_factor attribute for this method, so a 1SP CRS cannot
+    be exactly represented), and albers_conical_equal_area. Always
+    derive these from the original, unwrapped CRS: pyproj.CRS.to_cf()
+    does not "see through" a wrapped DerivedProjectedCRS.
+
+    Parameters
+    ----------
+    crs : pyproj.CRS
+        The (unwrapped) CRS to extract parameters from.
+
+    Returns
+    -------
+    tuple[dict, str | None]
+        The parameters (empty if the grid_mapping_name is unsupported or
+        unrecognized), and a warning message if the grid_mapping_name is
+        recognized but unsupported (currently: 1SP lambert_conformal_conic
+        variants), else None.
+    """
+    cf = crs.to_cf()
+    gmn = cf.get("grid_mapping_name")
+
+    if gmn == "transverse_mercator":
+        keys = ("scale_factor_at_central_meridian", *_SHARED_CF_PARAMS)
+    elif gmn in ("lambert_conformal_conic", "albers_conical_equal_area"):
+        method = crs.coordinate_operation.method_name if crs.coordinate_operation else ""
+        if "1SP" in method:
+            return {}, (
+                "1SP lambert_conformal_conic variants are not supported for "
+                "CF grid_mapping parameter extraction (CF has no "
+                "scale_factor attribute for this method); some netCDF "
+                "readers may not correctly position this grid."
+            )
+        keys = ("standard_parallel", *_SHARED_CF_PARAMS)
+    else:
+        return {}, None
+
+    params: dict[str, Any] = {}
+    for k in keys:
+        if k in cf:
+            params[k] = np.asarray(cf[k]) if k == "standard_parallel" else cf[k]
+    return params, None
+
+
+def wrap_rotated_crs(
+    crs: CRS,
+    xorigin: float,
+    yorigin: float,
+    angrot: float,
+) -> CRS | None:
+    """
+    Wrap a projected CRS in a derived CRS encoding MF6 grid rotation.
+
+    Builds a DerivedProjectedCRS over *crs* with an EPSG:9624 (affine
+    parametric transformation) deriving conversion parameterized from
+    xorigin/yorigin/angrot (degrees), matching MF6's DisNCStructured.f90
+    CRS_WKT rotation-wrapping exactly. Per ISO 19111, a deriving
+    conversion is directed base->derived, so the parameters encode the
+    world->local (inverse) rotation; a CRS-aware consumer applies the
+    inverse to resolve true position from local coordinates.
+
+    Parameters
+    ----------
+    crs : pyproj.CRS
+        The base CRS to wrap. Must be a projected CRS.
+    xorigin, yorigin : float
+        Grid origin in *crs*.
+    angrot : float
+        Grid rotation angle in degrees.
+
+    Returns
+    -------
+    pyproj.CRS | None
+        The wrapped CRS, or None if *crs* is not a projected CRS.
+    """
+    if not crs.is_projected:
+        return None
+
+    ang = np.radians(angrot)
+    a0 = -(xorigin * np.cos(ang) + yorigin * np.sin(ang))
+    a1 = np.cos(ang)
+    a2 = np.sin(ang)
+    b0 = xorigin * np.sin(ang) - yorigin * np.cos(ang)
+    b1 = -np.sin(ang)
+    b2 = np.cos(ang)
+
+    derived = {
+        "type": "DerivedProjectedCRS",
+        "name": "MODFLOW 6 rotated grid CRS",
+        "base_crs": crs.to_json_dict(),
+        "conversion": {
+            "name": "MODFLOW 6 grid rotation",
+            "method": {
+                "name": "Affine parametric transformation",
+                "id": {"authority": "EPSG", "code": 9624},
+            },
+            "parameters": [
+                {"name": "A0", "value": a0, "unit": "metre"},
+                {"name": "A1", "value": a1, "unit": "unity"},
+                {"name": "A2", "value": a2, "unit": "unity"},
+                {"name": "B0", "value": b0, "unit": "metre"},
+                {"name": "B1", "value": b1, "unit": "unity"},
+                {"name": "B2", "value": b2, "unit": "unity"},
+            ],
+        },
+        "coordinate_system": {
+            "subtype": "Cartesian",
+            "axis": [
+                {
+                    "name": "Easting",
+                    "abbreviation": "X",
+                    "direction": "east",
+                    "unit": "metre",
+                },
+                {
+                    "name": "Northing",
+                    "abbreviation": "Y",
+                    "direction": "north",
+                    "unit": "metre",
+                },
+            ],
+        },
+    }
+    return CRS.from_json_dict(derived)