uds.ugrid.to_netcdf() writes dataset that cannot be depth-sliced afterwards

[veenstrajelmer]

When writing+reading a subset of a dataset, the get_Dataset_atdepths plots an unexpected empty (all-nan) field. This happens because variables that are missing a _FillValue attribute in the encoding are set to _FillValue=nan upon xarray's .to_dataset(). The default fillvalue parsing that normally happens in get_Dataset_atdepths then does not work anymore.

MWE:

import dfm_tools as dfmt
import xugrid as xu
import matplotlib.pyplot as plt
plt.close("all")

file_nc = r'p:\11209731-002-nutrient-reduction-ta\runs_OSPAR\B05_waq_withDOM_2015\DFM_OUTPUT_DCSM-FM_0_5nm_waq\DCSM-FM_0_5nm_waq_0000_map.nc'
uds = xu.open_dataset(file_nc)

uds_sel = uds[['mesh2d_ucx', 'mesh2d_ucy',
               "mesh2d_s1", "mesh2d_bldepth", "mesh2d_sigmazdepth", "mesh2d_interface_z", "mesh2d_interface_sigma",
               'mesh2d_layer_sigma', 'mesh2d_layer_z', "mesh2d_layer_sigma_z", "mesh2d_interface_sigma_z", 
               'mesh2d_flowelem_bl']]

write_read = True
if write_read:
    file_temp = "test_sonia.nc"
    uds_sel.ugrid.to_netcdf(file_temp)
    # Re-reading subsetted file
    uds_sel_fromfile = dfmt.open_partitioned_dataset(file_temp)
    print("read_write finished")
else:
    uds_sel_fromfile = uds_sel

print("print fillvalues in original dataset")
ds = uds_sel
for varn in ds.obj.variables:
    if '_FillValue' in ds[varn].encoding:
        print(varn, ds[varn].encoding['_FillValue'])
print()
print("print fillvalues in write+read dataset")
ds = uds_sel_fromfile
for varn in ds.obj.variables:
    if '_FillValue' in ds[varn].encoding:
        print(varn, ds[varn].encoding['_FillValue'])

# depth extraction
uds_crop_sel_atdepth = dfmt.get_Dataset_atdepths(uds_sel, depths=0, reference='bedlevel')
uds_crop_sel_atdepth_fromfile = dfmt.get_Dataset_atdepths(uds_sel_fromfile, depths=0, reference='bedlevel')

# plotting
fig, ax = plt.subplots()
uds_crop_sel_atdepth.mesh2d_ucx.isel(time=0).ugrid.plot(cmap='jet')
fig, ax = plt.subplots()
uds_crop_sel_atdepth_fromfile.mesh2d_ucx.isel(time=0).ugrid.plot(cmap='jet') # blank plot in this case (only nans in this partition)

Expected (result without write+read):

Actual (with write+read):

[veenstrajelmer]

In the linked PR there is a workaround implemented in dfmt.open_partitioned_dataset that pops the _FillValue attrs that are nan. However, it is still inconvenient behaviour by xarray and would be good if solved in the future. Therefore we added test_remove_nan_fillvalue_attrs, to keep track if xarray fixes this themselves. Follow-up issue in Deltares/xugrid#176