218 getting an extra layer when reading geotop from cache

nlmod/cache.py

@@ -140,7 +140,7 @@ def decorator(*args, cachedir=None, cachename=None, **kwargs):
                     f"module of function {func.__name__} recently modified, not using cache"
                 )
 
-            cached_ds = xr.open_dataset(fname_cache, mask_and_scale=False)
+            cached_ds = xr.open_dataset(fname_cache)
 
             if pickle_check:
                 # add netcdf hash to function arguments dic, see #66
@@ -188,7 +188,7 @@ def decorator(*args, cachedir=None, cachename=None, **kwargs):
                 result.to_netcdf(fname_cache)
 
             # add netcdf hash to function arguments dic, see #66
-            temp = xr.open_dataset(fname_cache, mask_and_scale=False)
+            temp = xr.open_dataset(fname_cache)
             func_args_dic["_nc_hash"] = dask.base.tokenize(temp)
             temp.close()
 
@@ -220,12 +220,6 @@ def _check_ds(ds, ds2):
         True if the two datasets have the same grid and time discretization.
     """
 
-    # first remove _FillValue from all coordinates
-    for coord in ds2.coords:
-        if coord in ds.coords:
-            if "_FillValue" in ds2[coord].attrs and "_FillValue" not in ds[coord].attrs:
-                del ds2[coord].attrs["_FillValue"]
-
     for coord in ds2.coords:
         if coord in ds.coords:
             try:
@@ -352,8 +346,8 @@ def _get_modification_time(func):
 
 def _update_docstring_and_signature(func):
     """Add function arguments 'cachedir' and 'cachename' to the docstring and signature
-    of a function. 
-    
+    of a function.
+
     The function arguments are added before the "Returns" header in the
     docstring. If the function has no Returns header in the docstring, the function
     arguments are not added to the docstring.

nlmod/dims/base.py

@@ -472,7 +472,7 @@ def _get_vertex_grid_ds(
     for i in range(ncpl):
         icvert[i, : cell2d[i][3]] = cell2d[i][4:]
     ds["icvert"] = ("icell2d", "icv"), icvert
-    ds["icvert"].attrs["_FillValue"] = nodata
+    ds["icvert"].attrs["nodata"] = nodata
 
     if crs is not None:
         ds.rio.set_crs(crs)

nlmod/dims/grid.py

@@ -169,7 +169,7 @@ def modelgrid_to_vertex_ds(mg, ds, nodata=-1):
     for i in range(mg.ncpl):
         icvert[i, : cell2d[i][3]] = cell2d[i][4:]
     ds["icvert"] = ("icell2d", "icv"), icvert
-    ds["icvert"].attrs["_FillValue"] = nodata
+    ds["icvert"].attrs["nodata"] = nodata
     return ds
 
 
@@ -278,7 +278,7 @@ def gridprops_to_vertex_ds(gridprops, ds, nodata=-1):
     for i in range(gridprops["ncpl"]):
         icvert[i, : cell2d[i][3]] = cell2d[i][4:]
     ds["icvert"] = ("icell2d", "icv"), icvert
-    ds["icvert"].attrs["_FillValue"] = nodata
+    ds["icvert"].attrs["nodata"] = nodata
     return ds
 
 
@@ -300,8 +300,8 @@ def get_cell2d_from_ds(ds):
     x = ds["x"].data
     y = ds["y"].data
     icvert = ds["icvert"].data
-    if "_FillValue" in ds["icvert"].attrs:
-        nodata = ds["icvert"].attrs["_FillValue"]
+    if "nodata" in ds["icvert"].attrs:
+        nodata = ds["icvert"].attrs["nodata"]
     else:
         nodata = -1
         icvert = icvert.copy()
@@ -924,7 +924,7 @@ def da_to_reclist(
     else:
         if first_active_layer:
             fal = get_first_active_layer(ds)
-            cellids = np.where((mask.squeeze()) & (fal != fal.attrs["_FillValue"]))
+            cellids = np.where((mask.squeeze()) & (fal != fal.attrs["nodata"]))
             layers = col_to_list(fal, ds, cellids)
         elif only_active_cells:
             cellids = np.where((mask) & (ds["idomain"][layer] == 1))
@@ -1128,7 +1128,6 @@ def gdf_to_da(
     da = util.get_da_from_da_ds(ds, dims=ds.top.dims, data=fill_value)
     for ind, row in gdf_agg.iterrows():
         da.values[ind] = row[column]
-    da.attrs["_FillValue"] = fill_value
     return da
 
 

nlmod/dims/layers.py

@@ -1102,7 +1102,7 @@ def get_first_active_layer_from_idomain(idomain, nodata=-999):
             i,
             first_active_layer,
         )
-    first_active_layer.attrs["_FillValue"] = nodata
+    first_active_layer.attrs["nodata"] = nodata
     return first_active_layer
 
 
@@ -1132,7 +1132,7 @@ def get_last_active_layer_from_idomain(idomain, nodata=-999):
             i,
             last_active_layer,
         )
-    last_active_layer.attrs["_FillValue"] = nodata
+    last_active_layer.attrs["nodata"] = nodata
     return last_active_layer
 
 

nlmod/gis.py

@@ -366,7 +366,7 @@ def ds_to_ugrid_nc_file(
     # direction. Flopy specifies them in clockwise direction, so we need to
     # reverse the direction.
     data = np.flip(ds[face_node_connectivity].data, 1)
-    nodata = ds[face_node_connectivity].attrs.get("_FillValue")
+    nodata = ds[face_node_connectivity].attrs.get("nodata")
     if nodata is not None:
         # move the nodata values from the first columns to the last
         data_new = np.full(data.shape, nodata)

nlmod/gwf/horizontal_flow_barrier.py

@@ -240,7 +240,7 @@ def polygon_to_hfb(
     else:
         # find connections
         icvert = ds["icvert"].data
-        nodata = ds["icvert"].attrs["_FillValue"]
+        nodata = ds["icvert"].attrs["nodata"]
 
         edges = []
         for icell2d in range(icvert.shape[0]):

nlmod/plot/plotutil.py

@@ -14,8 +14,8 @@ def get_patches(ds, rotated=False):
         affine = get_affine_mod_to_world(ds)
         xy[:, 0], xy[:, 1] = affine * (xy[:, 0], xy[:, 1])
     icvert = ds["icvert"].data
-    if "_FillValue" in ds["icvert"].attrs:
-        nodata = ds["icvert"].attrs["_FillValue"]
+    if "nodata" in ds["icvert"].attrs:
+        nodata = ds["icvert"].attrs["nodata"]
     else:
         nodata = -1
         icvert = icvert.copy()

nlmod/read/boundaries.py

@@ -24,7 +24,7 @@ def get_municipalities(source="cbs", drop_water=True, **kwargs):
             gdf = gdf[gdf["water"] == "NEE"]
         gdf = gdf.set_index("gemeentenaam")
     else:
-        raise (Exception(f"Unknown source: {source}"))
+        raise ValueError(f"Unknown source: {source}")
     return gdf
 
 

nlmod/read/knmi.py

@@ -171,7 +171,7 @@ def get_locations_vertex(ds):
     """
     # get active locations
     fal = get_first_active_layer(ds)
-    icell2d_active = np.where(fal != fal.attrs["_FillValue"])[0]
+    icell2d_active = np.where(fal != fal.attrs["nodata"])[0]
 
     # create dataframe from active locations
     x = ds["x"].sel(icell2d=icell2d_active)
@@ -206,7 +206,7 @@ def get_locations_structured(ds):
 
     # store x and y mids in locations of active cells
     fal = get_first_active_layer(ds)
-    rows, columns = np.where(fal != fal.attrs["_FillValue"])
+    rows, columns = np.where(fal != fal.attrs["nodata"])
     x = np.array([ds["x"].data[col] for col in columns])
     y = np.array([ds["y"].data[row] for row in rows])
     if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:

nlmod/read/regis.py

@@ -175,8 +175,6 @@ def get_regis(
             ds[datavar].attrs["units"] = "mNAP"
         elif datavar in ["kh", "kv"]:
             ds[datavar].attrs["units"] = "m/day"
-        # set _FillValue to NaN, otherise problems with caching will arise
-        ds[datavar].encoding["_FillValue"] = np.NaN
 
     # set the crs to dutch rd-coordinates
     ds.rio.set_crs(28992)

nlmod/read/rws.py

@@ -154,7 +154,7 @@ def add_northsea(ds, cachedir=None):
 
     # fill top, bot, kh, kv at sea cells
     fal = dims.get_first_active_layer(ds)
-    fill_mask = (fal == fal.attrs["_FillValue"]) * ds["northsea"]
+    fill_mask = (fal == fal.attrs["nodata"]) * ds["northsea"]
     ds = dims.fill_top_bot_kh_kv_at_mask(ds, fill_mask)
 
     # add bathymetry noordzee

tests/test_004_northsea.py

@@ -39,7 +39,7 @@ def test_fill_top_bot_kh_kv_seamodel():
     ds.update(nlmod.read.rws.get_northsea(ds))
 
     fal = nlmod.layers.get_first_active_layer(ds)
-    fill_mask = (fal == fal._FillValue) * ds["northsea"]
+    fill_mask = (fal == fal.nodata) * ds["northsea"]
     nlmod.layers.fill_top_bot_kh_kv_at_mask(ds, fill_mask)
 
 
@@ -49,7 +49,7 @@ def test_fill_top_bot_kh_kv_nosea():
     ds.update(nlmod.read.rws.get_northsea(ds))
 
     fal = nlmod.layers.get_first_active_layer(ds)
-    fill_mask = (fal == fal._FillValue) * ds["northsea"]
+    fill_mask = (fal == fal.nodata) * ds["northsea"]
     nlmod.layers.fill_top_bot_kh_kv_at_mask(ds, fill_mask)
 
 
@@ -78,6 +78,6 @@ def test_add_bathymetrie_to_top_bot_kh_kv_seamodel():
     ds.update(nlmod.read.jarkus.get_bathymetry(ds, ds["northsea"]))
 
     fal = nlmod.layers.get_first_active_layer(ds)
-    fill_mask = (fal == fal._FillValue) * ds["northsea"]
+    fill_mask = (fal == fal.nodata) * ds["northsea"]
 
     nlmod.read.jarkus.add_bathymetry_to_top_bot_kh_kv(ds, ds["bathymetry"], fill_mask)