Parsing of 'c' value as float instead of timedelta[ns]

nlmod/read/regis.py

@@ -112,10 +112,10 @@ def get_regis(
         included in the model. the Default is "AKc" which is the bottom
         layer of regis. call nlmod.read.regis.get_layer_names() to get a list
         of regis names.
-    variables : tuple, optional
-        a tuple of the variables to keep from the regis Dataset. Possible
-        entries in the list are 'top', 'botm', 'kD', 'c', 'kh', 'kv', 'sdh' and
-        'sdv'. The default is ("top", "botm", "kh", "kv").
+    variables : tuple or list, optional
+        The variables to keep from the regis Dataset. Possible entries in the list are
+        'top', 'botm', 'kD', 'c', 'kh', 'kv', 'sdh' and 'sdv'. The default is
+        ("top", "botm", "kh", "kv").
     remove_nan_layers : bool, optional
         When True, layers that do not occur in the requested extent (layers that contain
         only NaN values for the botm array) are removed. The default is True.
@@ -135,7 +135,7 @@ def get_regis(
     regis_ds : xarray dataset
         dataset with regis data projected on the modelgrid.
     """
-    ds = xr.open_dataset(REGIS_URL, decode_coords="all")
+    ds = xr.open_dataset(REGIS_URL, decode_times=False, decode_coords="all")
     if "crs" in ds.coords:
         # remove the crs coordinate, as rioxarray does not recognise the crs
         # and we set the crs at the end of this method by hand
@@ -167,7 +167,11 @@ def get_regis(
     ds = ds.rename_vars({"bottom": "botm"})
 
     # slice data vars
-    if variables is not None:
+    if variables is None:
+        variables = list(ds.data_vars)
+    else:
+        if isinstance(variables, str):
+            variables = [variables]
         if probabilities:
             variables = variables + ("sdh", "sdv")
         ds = ds[list(variables)]
@@ -180,12 +184,20 @@ def get_regis(
 
     if remove_nan_layers:
         # only keep layers with at least one active cell
-        ds = ds.sel(layer=~(np.isnan(ds["botm"])).all(ds["botm"].dims[1:]))
+        if "botm" in ds:
+            mask_layer = ~(np.isnan(ds["botm"])).all(ds["botm"].dims[1:])
+        else:
+            var = variables[0]
+            mask_layer = ~(np.isnan(ds[var])).all(ds[var].dims[1:])
+            for var in variables[1:]:
+                mask_layer = mask_layer | ~(np.isnan(ds[var])).all(ds[var].dims[1:])
+        ds = ds.sel(layer=mask_layer)
+
         if len(ds.layer) == 0:
             msg = "No data found. Please supply valid extent in the Netherlands in RD-coordinates"
             raise (Exception(msg))
 
-    if drop_layer_dim_from_top:
+    if drop_layer_dim_from_top and "botm" in ds and "top" in ds:
         ds = remove_layer_dim_from_top(ds)
 
     ds.attrs["gridtype"] = "structured"
@@ -198,6 +210,10 @@ def get_regis(
             ds[datavar].attrs["units"] = "mNAP"
         elif datavar in ["kh", "kv"]:
             ds[datavar].attrs["units"] = "m/day"
+        elif datavar in ["c"]:
+            ds[datavar].attrs["units"] = "day"
+        elif datavar in ["kD"]:
+            ds[datavar].attrs["units"] = "m2/day"
 
     # set the crs to dutch rd-coordinates
     ds.rio.write_crs(28992, inplace=True)

tests/test_002_regis_geotop.py

@@ -1,5 +1,5 @@
 import matplotlib.pyplot as plt
-
+import numpy as np
 import nlmod
 
 
@@ -20,6 +20,18 @@ def test_get_regis_botm_layer_BEk1(
     assert regis_ds.layer.values[-1] == botm_layer
 
 
+def test_get_regis_only_c(extent=[98700.0, 99000.0, 489500.0, 489700.0]):
+    regis_ds = nlmod.read.regis.get_regis(extent, variables="c")
+    assert np.all([x == "c" for x in regis_ds.data_vars])
+    assert regis_ds.sizes["layer"] == 8
+
+
+def test_get_regis_only_c_and_kd(extent=[98700.0, 99000.0, 489500.0, 489700.0]):
+    regis_ds = nlmod.read.regis.get_regis(extent, variables=["c", "kD"])
+    assert np.all([x in ["c", "kD"] for x in regis_ds.data_vars])
+    assert regis_ds.sizes["layer"] == 18
+
+
 def test_get_geotop(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
     geotop_ds = nlmod.read.geotop.get_geotop(extent)
     line = [(extent[0], extent[2]), (extent[1], extent[3])]