improve getting layer/interface dimensions?

[veenstrajelmer]

The attributes "layer_dimension" and "interface_dimension" in the topology variable of Delft3DFM output cannot be found in the ugrid conventions. For Delft3D4 we manually add these here:

dfm_tools/dfm_tools/xugrid_helpers.py

Lines 477 to 488 in 00c3802

	# set vertical dimensions attr
	# TODO: would be more convenient to do within xu.Ugrid2d(): https://github.com/Deltares/xugrid/issues/195#issuecomment-2111841390
	grid_attrs = {"vertical_dimensions": ds.grid.attrs["vertical_dimensions"]}
	ds_temp["mesh2d"] = ds_temp["mesh2d"].assign_attrs(grid_attrs)
	# drop attrs pointing to the removed grid variable (topology is now in mesh2d)
	# TODO: this is not possible on the xu.UgridDataset directly
	for varn in ds_temp.data_vars:
	if "grid" in ds_temp[varn].attrs.keys():
	del ds_temp[varn].attrs["grid"]
	uds = xu.UgridDataset(ds_temp)

Note: a "vertical_dimensions" attribute is added, while dfm_tools needs "layer_dimension" and "interface_dimension".

Either way, would it be possible to extract the vertical dimensions for sigma/z/z-sigma in a more generic way? At the moment dfm_tools depends on them in dfmt.get_vertical_dimensions():

dfm_tools/dfm_tools/xugrid_helpers.py

Lines 40 to 43 in 00c3802

	gridname = uds.grid.name
	grid_info = uds.grid.to_dataset()[gridname]
	if hasattr(grid_info,'layer_dimension'):
	return grid_info.layer_dimension, grid_info.interface_dimension

sigma-layers
reconstruct_zw_zcc_fromsigma does not use dfmt.get_vertical_dimensions(), but depend on formula terms only. So no issue here. Testcase is the curvedbend demo dataset.

z-layers
reconstruct_zw_zcc_fromz uses dfmt.get_vertical_dimensions() and datasets do not have formula terms variables, but in theory the layering variables should have standard_name="altitude". However, for Grevelingen, this results in more variable names than only the layers:

import dfm_tools as dfmt
uds = dfmt.data.fm_grevelingen_map()
uds_alt = uds.filter_by_attrs(standard_name="altitude")
print(uds_alt.data_vars)

Gives:

Data variables:
    mesh2d_layer_z      (nmesh2d_layer) float64 288B dask.array<chunksize=(36,), meta=np.ndarray>
    mesh2d_interface_z  (nmesh2d_interface) float64 296B dask.array<chunksize=(37,), meta=np.ndarray>
    mesh2d_node_z       (nmesh2d_node) float64 185kB dask.array<chunksize=(23108,), meta=np.ndarray>
    mesh2d_flowelem_bl  (nmesh2d_face) float64 358kB dask.array<chunksize=(44796,), meta=np.ndarray>

z-sigma layers
reconstruct_zw_zcc_fromzsigma uses dfmt.get_vertical_dimensions(), but also uses formula terms (for interfaces). From these we can get the layer variables and from these we can get the layer dimensions. In principle.

filtering ideas
We can filter via formula terms (sigma and z-sigma) and standard_name="altitude" (fullgrid output and z-models). The sigma and z-sigma models can them be distinguished via the formula terms names (sigma vs sigma_z). The fullgrid output and the z-models can only be distinguished because the fullgrid output is set as coords and the z-coordinate variables are not. This seems not super robust, since this inconsistency can easily change in an update of Delft3DFM.

This codeblock shows in comments all the variables with a formula_terms attribute or a standard_name="altitude" attribute:

import dfm_tools as dfmt

# file_nc = dfmt.data.fm_curvedbend_map(return_filepath=True)
# ['mesh2d_layer_sigma', 'mesh2d_interface_sigma']
# ['mesh2d_node_z', 'mesh2d_flowelem_bl']
# sigma-layer model, computing zw/zcc (fullgrid) values and treat as fullgrid model from here
# file_nc = dfmt.data.fm_grevelingen_map(return_filepath=True)
# []
# ['mesh2d_layer_z', 'mesh2d_interface_z', 'mesh2d_node_z', 'mesh2d_flowelem_bl']
# z-layer model, computing zw/zcc (fullgrid) values and treat as fullgrid model from here
# file_nc = dfmt.data.fm_westernscheldt_map(return_filepath=True)
# ['mesh2d_layer_sigma_z', 'mesh2d_interface_sigma_z']
# ['mesh2d_node_z', 'mesh2d_layer_z', 'mesh2d_interface_z', 'mesh2d_flowelem_bl']
# zsigma-layer model, computing zw/zcc (fullgrid) values and treat as fullgrid model from here
# file_nc = r'p:\1204257-dcsmzuno\2006-2012\3D-DCSM-FM\A18b_ntsu1\DFM_OUTPUT_DCSM-FM_0_5nm\DCSM-FM_0_5nm_0000_map.nc' # szigma fullgrid
# []
# ['mesh2d_node_z', 'mesh2d_flowelem_bl', 'mesh2d_flowelem_zcc', 'mesh2d_flowelem_zw']
# zw/zcc (fullgrid) values already present in Dataset in variable mesh2d_flowelem_zw
# file_nc = r'p:\archivedprojects\11206813-006-kpp2021_rmm-2d\C_Work\31_RMM_FMmodel\computations\model_setup\run_207\results\RMM_dflowfm_0000_map.nc' # 2D model
# []
# ['mesh2d_node_z', 'mesh2d_flowelem_bl']
# KeyError: 'layers present, but unknown layertype, expected one of variables: mesh2d_flowelem_zw, mesh2d_layer_sigma, mesh2d_layer_z' >> makes sense since 2D
file_nc = r'p:\archivedprojects\11203379-005-mwra-updated-bem\03_model\02_final\A72_ntsu0_kzlb2\DFM_OUTPUT_MB_02\MB_02_0000_map.nc' # zlayer
# []
# ['mesh2d_node_z', 'mesh2d_layer_z', 'mesh2d_interface_z', 'mesh2d_flowelem_bl']
# z-layer model, computing zw/zcc (fullgrid) values and treat as fullgrid model from here
uds = dfmt.open_partitioned_dataset(file_nc)

# uds = dfmt.data.d3d_curvedbend_trim()
# ['SIG_LYR', 'SIG_INTF']
# []
# KeyError: 'layers present, but unknown layertype, expected one of variables: mesh2d_flowelem_zw, mesh2d_layer_sigma, mesh2d_layer_z'
# uds = dfmt.data.d3d_westernscheldt_trim()
# ['SIG_LYR', 'SIG_INTF']
# []
# KeyError: 'layers present, but unknown layertype, expected one of variables: mesh2d_flowelem_zw, mesh2d_layer_sigma, mesh2d_layer_z'

# names of variables containing attribute "formula_terms", only present in z-sigma and sigma models
print(list(uds.reset_coords().filter_by_attrs(formula_terms=lambda v: v is not None).variables))

# z-layer models have int/cen depths in variable with standard_name=altitude
# however, also mesh2d_node_z and mesh2d_flowelem_bl have standard_name=altitude
# we need reset_coords() since fullgrid output (mesh2d_flowelem_zcc and mesh2d_flowelem_zw) are set as coords
# we use data_vars to avoid getting node_x/node_y/etc
print(list(uds.reset_coords().filter_by_attrs(standard_name="altitude").data_vars))
# inconveniently enough, both fullgrid output and z-layer models are both to be identified with standard_name="altitude"
# there seems no way to distinguish between them

dfmt.reconstruct_zw_zcc(uds)

Furthermore, in dfmt.reconstruct_zw_zcc() we already deduce the layertype (is not called in case of 2D models). This is not a super generic function, but the above analysis also shows why. However, since SIG_LYR and SIG_INTF are index variables, the Delft3D4 models are not seen as sigma. This can be fixed with the updated filter above. However, next it fails on the fact that "interface" or "layer" is not part of the varname, this can be resolved by checking these in the long_name instead since that seems to be shared. This is probably more generic, albeit probably not according to cf-conventions (there is no way to distinguish interfaces from sigma in the conventions).

additional checks

• does all the above also works for three delft3d4 mapfiles from notebook? >> they both have sigma layers defined by formula_terms attrs. So with this more generic approach, these files also work
• fullgrid output: reconstruct function is not called, but do we need dfmt.get_vertical_dimensions() for something else?
• more generic: are there any other locations in the dfm_tools code where we access these attrs?
• run updated methods on all ugrid testdatasets in the postprocessing example script
• check for hisfiles or are these always fullgrid? (also check delft3d4 mapfiles)
• If all else fails: at least add support for layers for delft3d4, this can be relatively easily done since they have formula_terms