Small improvements

nlmod/dims/grid.py

@@ -309,7 +309,7 @@ def modelgrid_to_ds(mg=None, grbfile=None):
         Dataset containing grid information
     """
     if mg is None and grbfile is not None:
-        mg = flopy.utils.MfGrdFile(grbfile).modelgrid
+        mg = flopy.mf6.utils.MfGrdFile(grbfile).modelgrid
     elif mg is None and grbfile is None:
         raise ValueError("Either 'mg' or 'grbfile' should be specified!")
     if mg.grid_type == "structured":
@@ -2124,3 +2124,59 @@ def get_affine(ds, sx=None, sy=None):
         xoff = attrs["extent"][0]
         yoff = attrs["extent"][3]
         return Affine.translation(xoff, yoff) * Affine.scale(sx, sy)
+
+
+def _shoelace_formula(x, y):
+    """Calculate the area of a polygon using the shoelace formula.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        x-coordinates of the polygon.
+    y : np.ndarray
+        y-coordinates of the polygon.
+
+    Returns
+    -------
+    area : float
+        area of the polygon.
+    """
+    x = x - np.min(x)
+    y = y - np.min(y)
+    return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
+
+
+def get_area(ds):
+    """Calculate the area of each cell in the model grid.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        model dataset.
+
+    Returns
+    -------
+    ds : xr.Dataset
+        model dataset with an area variable
+    """
+    if ds.gridtype == "structured":
+        area = xr.DataArray(
+            np.outer(get_delc(ds), get_delr(ds)),
+            dims=("y", "x"),
+            coords={"y": ds["y"], "x": ds["x"]},
+        )
+    elif ds.gridtype == "vertex":
+        area = np.zeros(ds["icell2d"].size)
+        for icell2d in ds["icell2d"]:
+            area[icell2d] = _shoelace_formula(
+                ds["xv"][ds["icvert"].isel(icell2d=icell2d)],
+                ds["yv"][ds["icvert"].isel(icell2d=icell2d)],
+            )
+        area = xr.DataArray(
+            area,
+            dims=("icell2d"),
+            coords={"icell2d": ds["icell2d"]},
+        )
+    else:
+        raise ValueError("function only support structured or vertex gridtypes")
+    return area

nlmod/dims/time.py

@@ -247,7 +247,9 @@ def set_ds_time(
         # parse start
         if isinstance(start, (int, np.integer, float)):
             if isinstance(time[0], (int, np.integer, float, str)):
-                raise TypeError("Make sure 'start' or 'time' argument is a valid TimeStamp")
+                raise TypeError(
+                    "Make sure 'start' or 'time' argument is a valid TimeStamp"
+                )
             start = time[0] - pd.to_timedelta(start, "D")
         elif isinstance(start, str):
             start = pd.Timestamp(start)
@@ -278,13 +280,14 @@ def set_ds_time(
         elif isinstance(time[0], cftime.datetime):
             start = _pd_timestamp_to_cftime(start)
         else:
-            msg = (
-                f"Cannot process 'time' argument. Datatype -> {type(time)} not understood."
-            )
+            msg = f"Cannot process 'time' argument. Datatype -> {type(time)} not understood."
             raise TypeError(msg)
     except (OutOfBoundsDatetime, OutOfBoundsTimedelta) as e:
-        msg = "cannot convert 'start' and 'time' to pandas datetime, use cftime types for 'start' and 'time'"
-        raise type(e)(msg)
+        msg = (
+            "cannot convert 'start' and 'time' to pandas datetime, use cftime "
+            "types for 'start' and 'time'"
+        )
+        raise type(e)(msg) from e
 
     if time[0] <= start:
         msg = (
@@ -735,13 +738,13 @@ def dataframe_to_flopy_timeseries(
     )
 
 
-def ds_time_to_pandas_index(ds, include_start=True):
+def ds_time_to_pandas_index(da_or_ds, include_start=True):
     """Convert xarray time index to pandas datetime index.
 
     Parameters
     ----------
-    ds : xarray.Dataset
-        dataset with time index
+    da_or_ds : xarray.Dataset or xarray.DataArray
+        dataset with time index or time data array
     include_start : bool, optional
         include the start time in the index, by default True
 
@@ -750,13 +753,22 @@ def ds_time_to_pandas_index(ds, include_start=True):
     pd.DatetimeIndex
         pandas datetime index
     """
+    if isinstance(da_or_ds, xr.Dataset):
+        da = da_or_ds["time"]
+    elif isinstance(da_or_ds, xr.DataArray):
+        da = da_or_ds
+    else:
+        raise ValueError(
+            "Either pass model dataset or time dataarray (e.g. `ds.time`)."
+        )
+
     if include_start:
-        if ds.time.dtype.kind == "M": # "M" is a numpy datetime-type
-            return ds.time.to_index().insert(0, pd.Timestamp(ds.time.start))
-        elif ds.time.dtype.kind == "O":
-            start = _pd_timestamp_to_cftime(pd.Timestamp(ds.time.start))
-            return ds.time.to_index().insert(0, start)
-        elif ds.time.dtype.kind in ["i", "f"]:
-            return ds.time.to_index().insert(0, 0)
+        if da.dtype.kind == "M":  # "M" is a numpy datetime-type
+            return da.to_index().insert(0, pd.Timestamp(da.start))
+        elif da.dtype.kind == "O":
+            start = _pd_timestamp_to_cftime(pd.Timestamp(da.start))
+            return da.to_index().insert(0, start)
+        elif da.dtype.kind in ["i", "f"]:
+            return da.to_index().insert(0, 0)
     else:
-        return ds.time.to_index()
+        return da.to_index()

nlmod/plot/plot.py

@@ -46,7 +46,7 @@ def modelgrid(ds, ax=None, **kwargs):
     if ax is None:
         _, ax = plt.subplots(figsize=(10, 10))
         ax.set_aspect("auto")
-    modelgrid = modelgrid_from_ds(ds)
+    modelgrid = modelgrid_from_ds(ds, rotated=kwargs.get("rotated", False))
     extent = None if ax.get_autoscale_on() else ax.axis()
     modelgrid.plot(ax=ax, **kwargs)
     if extent is not None:

nlmod/plot/plotutil.py

@@ -70,10 +70,14 @@ def add_background_map(ax, crs=EPSG_28992, map_provider="nlmaps.standaard", **kw
     """
     import contextily as ctx
 
-    if isinstance(crs, (str, int)):
+    if isinstance(crs, str):
         import pyproj
 
-        proj = pyproj.Proj(crs)
+        proj = pyproj.CRS.from_string(crs)
+    elif isinstance(crs, int):
+        import pyproj
+
+        proj = pyproj.CRS.from_epsg(crs)
 
     providers = _list_contextily_providers()
     ctx.add_basemap(ax, source=providers[map_provider], crs=proj.srs, **kwargs)