added node-to-faces conversion (#644)

* added node-to-faces

* added testcase

* updated whatsnew

dfm_tools/xugrid_helpers.py

@@ -20,6 +20,7 @@
  "open_dataset_curvilinear",
  "open_dataset_delft3d4",
  "uda_edges_to_faces",
+ "uda_nodes_to_faces",
  "uda_interfaces_to_centers",
  "add_network_cellinfo",
  "enrich_rst_with_map",
@@ -462,6 +463,57 @@ def uda_edges_to_faces(uda_edge : xu.UgridDataArray) -> xu.UgridDataArray:
  return uda_face
 
 
+def uda_nodes_to_faces(uda_node : xu.UgridDataArray) -> xu.UgridDataArray:
+ """
+ Interpolates a ugrid variable (xu.DataArray) with an node dimension to the faces by averaging the 3/4 nodes around each face.
+ Since node variables are mostly defined on interfaces, it also interpolates from interfaces to layers
+
+ Parameters
+ ----------
+ uda_node : xu.UgridDataArray
+ DESCRIPTION.
+
+ Raises
+ ------
+ KeyError
+ DESCRIPTION.
+
+ Returns
+ -------
+ uda_face : xu.UgridDataArray
+ DESCRIPTION.
+
+ """
+
+ dimn_faces = uda_node.grid.face_dimension
+ dimn_maxfn = 'nMax_face_nodes' #arbitrary dimname that is reduced anyway
+ dimn_nodes = uda_node.grid.node_dimension
+ fill_value = uda_node.grid.fill_value
+
+ if dimn_nodes not in uda_node.sizes:
+ raise KeyError(f'varname "{uda_node.name}" does not have an node dimension ({dimn_nodes})')
+
+ # construct indexing array
+ data_fnc = xr.DataArray(uda_node.grid.face_node_connectivity,dims=(dimn_faces,dimn_maxfn))
+ data_fnc_validbool = data_fnc!=fill_value
+ data_fnc = data_fnc.where(data_fnc_validbool,-1)
+
+ print('node-to-face interpolation: ',end='')
+ dtstart = dt.datetime.now()
+ # for each face, select all corresponding node values (this takes some time)
+ uda_face_allnodes = uda_node.isel({dimn_nodes:data_fnc})
+ # replace nonexistent nodes with nan
+ uda_face_allnodes = uda_face_allnodes.where(data_fnc_validbool) #replace all values for fillvalue nodes (-1) with nan
+ # average node values per face
+ uda_face = uda_face_allnodes.mean(dim=dimn_maxfn,keep_attrs=True)
+ #update attrs from node to face
+ face_attrs = {'location': 'face', 'cell_methods': f'{dimn_faces}: mean'}
+ uda_face = uda_face.assign_attrs(face_attrs)
+ print(f'{(dt.datetime.now()-dtstart).total_seconds():.2f} sec')
+
+ return uda_face
+
+
 def uda_interfaces_to_centers(uda_int : xu.UgridDataArray) -> xu.UgridDataArray:
  dimn_layer, dimn_interface = get_vertical_dimensions(uda_int)
 

docs/whats-new.md

@@ -1,3 +1,9 @@
+## UNRELEASED
+
+### Feat
+- interpolation of node variables to faces with `dfmt.uda_nodes_to_faces()` by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#644](https://github.com/Deltares/dfm_tools/pull/644)
+
+
 ## 0.16.0 (2023-11-03)
 
 ### Feat

tests/test_xugrid_helpers.py

@@ -56,6 +56,17 @@ def test_open_2Dnetwork_with_1Dtopology():
  assert node_x.attrs['standard_name'] == "longitude"
 
 
+@pytest.mark.unittest
+def test_uda_edges_to_faces():
+ file_nc = dfmt.data.fm_grevelingen_map(return_filepath=True) #zlayer
+
+ uds = xu.open_dataset(file_nc.replace('0*','0002')) #partition 0002 of grevelingen contains both triangles as squares
+
+ uda_edge = uds['mesh2d_vicwwu'].isel(time=0, nmesh2d_interface=0)
+ uda_face = dfmt.uda_edges_to_faces(uda_edge)
+ assert uda_face.dims == (uds.grid.face_dimension,)
+
+
 @pytest.mark.systemtest
 def test_uda_edges_to_faces_interfaces_to_centers():
  file_nc = dfmt.data.fm_grevelingen_map(return_filepath=True) #zlayer
@@ -92,6 +103,17 @@ def test_uda_edges_to_faces_interfaces_to_centers():
  assert (np.abs(uda_face_sel.to_numpy()-uda_face_sel_expected)<1e-6).all()
 
 
+@pytest.mark.unittest
+def test_uda_nodes_to_faces():
+ file_nc = dfmt.data.fm_grevelingen_net(return_filepath=True)
+ uds = dfmt.open_partitioned_dataset(file_nc)
+
+ uda_node = uds.mesh2d_node_z
+ uda_face = dfmt.uda_nodes_to_faces(uda_node)
+
+ assert uda_face.dims == (uds.grid.face_dimension,)
+
+
 @pytest.mark.requireslocaldata
 @pytest.mark.unittest
 def test_enrich_rst_with_map():