initial work on face to node

test/test_topological_agg.py

@@ -7,7 +7,6 @@
 
 current_path = Path(os.path.dirname(os.path.realpath(__file__)))
 
-
 ds_path = current_path / 'meshfiles' / "mpas" / "QU" / 'oQU480.231010.nc'
 
 AGGS = ["topological_mean",
@@ -42,9 +41,12 @@ def test_node_to_edge_aggs():
 
 
 def test_edge_to_face_aggs():
-    grid_path = "/Users/aaronzedwick/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc"
-    data_path = "/Users/aaronzedwick/uxarray/test/meshfiles/ugrid/quad-hexagon/random-edge-data.nc"
+    grid_path = '/Users/aaronzedwick/uxarray/test/meshfiles/mpas/QU/mesh.QU.1920km.151026.nc'
+
+    uxds = ux.open_dataset(grid_path, grid_path)
+
+    uxds = uxds['latCell'].subset.nearest_neighbor(k=3, center_coord=[0, 0])
 
-    uxds = ux.open_dataset(grid_path, data_path)
+    uxda_edge_face_agg = uxds.topological_mean(destination="node")
 
-    test = uxds["random_data_edge"].topological_mean("face")
+    print(uxda_edge_face_agg)

uxarray/conventions/ugrid.py

@@ -224,6 +224,12 @@
 
 N_NODES_PER_FACE_DIMS = ["n_face"]
 
+N_FACES_PER_NODE_ATTRS = {
+    "cf_role": "n_faces_per_node",
+    "long name": "Number of faces per node",
+}
+
+N_FACES_PER_NODE_DIMS = ["n_node"]
 
 CONNECTIVITY_NAMES = [
     "face_node_connectivity",

uxarray/core/aggregation.py

@@ -56,15 +56,15 @@ def _uxda_grid_aggregate(uxda, destination, aggregation, **kwargs):
 
     elif uxda._face_centered():
         # aggregation of a face-centered data variable
-        raise NotImplementedError(
-            "Aggregation of face-centered data variables is not yet supported."
-        )
-        # if destination == "node":
-        #     pass
-        # elif destination == "edge":
-        #     pass
-        # else:
-        #     raise ValueError("TODO: ")
+        # raise NotImplementedError(
+        #     "Aggregation of face-centered data variables is not yet supported."
+        # )
+        if destination == "node":
+            return _face_to_node_aggregation(uxda, aggregation, kwargs)
+        elif destination == "edge":
+            pass
+        else:
+            raise ValueError("TODO: ")
 
     else:
         raise ValueError(
@@ -221,7 +221,7 @@ def _edge_to_face_aggregation(uxda, aggregation, aggregation_func_kwargs):
     """Applies an Edge to Face Topological aggregation."""
     if not uxda._edge_centered():
         raise ValueError(
-            f"Data Variable must be mapped to the corner nodes of each face, with dimension "
+            f"Data Variable must be mapped to the edge centers of each face, with dimension "
             f"{uxda.uxgrid.n_edge}."
         )
 
@@ -244,3 +244,64 @@ def _edge_to_face_aggregation(uxda, aggregation, aggregation_func_kwargs):
         dims=uxda.dims,
         name=uxda.name,
     ).rename({"n_edge": "n_face"})
+
+
+def _apply_face_to_node_aggregation_numpy(
+    uxda, aggregation_func, aggregation_func_kwargs
+):
+    """Applies an Edge to Face Topological aggregation on a Numpy array."""
+    data = uxda.values
+    node_face_conn = uxda.uxgrid.node_face_connectivity.values
+    n_nodes_per_face = uxda.uxgrid.n_faces_per_node.values
+
+    (
+        change_ind,
+        n_nodes_per_face_sorted_ind,
+        element_sizes,
+        size_counts,
+    ) = get_face_node_partitions(n_nodes_per_face)
+
+    result = np.empty(shape=(data.shape[:-1]) + (uxda.uxgrid.n_face,))
+
+    for e, start, end in zip(element_sizes, change_ind[:-1], change_ind[1:]):
+        face_inds = n_nodes_per_face_sorted_ind[start:end]
+        face_nodes_par = node_face_conn[face_inds, 0:e]
+
+        # apply aggregation function to current node face partition
+        aggregation_par = aggregation_func(
+            data[..., face_nodes_par], axis=-1, **aggregation_func_kwargs
+        )
+
+        # store current aggregation
+        result[..., face_inds] = aggregation_par
+
+    return result
+
+
+def _face_to_node_aggregation(uxda, aggregation, aggregation_func_kwargs):
+    """Applies an Edge to Face Topological aggregation."""
+    if not uxda._face_centered():
+        raise ValueError(
+            f"Data Variable must be mapped to the face centers of each face, with dimension "
+            f"{uxda.uxgrid.n_face}."
+        )
+
+    if isinstance(uxda.data, np.ndarray):
+        # apply aggregation using numpy
+        aggregated_var = _apply_face_to_node_aggregation_numpy(
+            uxda, NUMPY_AGGREGATIONS[aggregation], aggregation_func_kwargs
+        )
+    elif isinstance(uxda.data, da.Array):
+        # apply aggregation on dask array, TODO:
+        aggregated_var = _apply_face_to_node_aggregation_numpy(
+            uxda, NUMPY_AGGREGATIONS[aggregation], aggregation_func_kwargs
+        )
+    else:
+        raise ValueError
+
+    return uxarray.core.dataarray.UxDataArray(
+        uxgrid=uxda.uxgrid,
+        data=aggregated_var,
+        dims=uxda.dims,
+        name=uxda.name,
+    ).rename({"n_face": "n_node"})

uxarray/grid/connectivity.py

@@ -158,6 +158,42 @@ def _build_n_nodes_per_face(face_nodes, n_face, n_max_face_nodes):
     return n_nodes_per_face
 
 
+def _populate_n_faces_per_node(grid):
+    """Constructs the connectivity variable (``n_nodes_per_face``) and stores
+    it within the internal (``Grid._ds``) and through the attribute
+    (``Grid.n_nodes_per_face``)."""
+
+    n_faces_per_node = _build_n_nodes_per_face(
+        grid.node_face_connectivity.values, grid.n_node, grid.n_max_node_faces
+    )
+
+    if n_faces_per_node.ndim == 0:
+        # convert scalar value into a [1, 1] array
+        n_faces_per_node = np.expand_dims(n_faces_per_node, 0)
+
+    # add to internal dataset
+    grid._ds["n_faces_per_node"] = xr.DataArray(
+        data=n_faces_per_node,
+        dims=ugrid.N_FACES_PER_NODE_DIMS,
+        attrs=ugrid.N_NODES_PER_FACE_ATTRS,
+    )
+
+
+@njit(cache=True)
+def _build_n_faces_per_node(node_faces, n_nodes, n_max_node_faces):
+    """Constructs ``n_nodes_per_face``, which contains the number of non-fill-
+    value nodes for each face in ``face_node_connectivity``"""
+
+    # padding to shape [n_face, n_max_face_nodes + 1]
+    closed = np.ones((n_nodes, n_max_node_faces + 1), dtype=INT_DTYPE) * INT_FILL_VALUE
+
+    closed[:, :-1] = node_faces.copy()
+
+    n_faces_per_node = np.argmax(closed == INT_FILL_VALUE, axis=1)
+
+    return n_faces_per_node
+
+
 def _populate_edge_node_connectivity(grid):
     """Constructs the UGRID connectivity variable (``edge_node_connectivity``)
     and stores it within the internal (``Grid._ds``) and through the attribute

uxarray/grid/grid.py

@@ -54,7 +54,7 @@
     _populate_n_nodes_per_face,
     _populate_node_face_connectivity,
     _populate_edge_face_connectivity,
-    _populate_face_face_connectivity,
+    _populate_face_face_connectivity, _populate_n_faces_per_node,
 )
 
 from uxarray.grid.geometry import (
@@ -848,6 +848,23 @@ def n_nodes_per_face(self, value):
         assert isinstance(value, xr.DataArray)
         self._ds["n_nodes_per_face"] = value
 
+    @property
+    def n_faces_per_node(self) -> xr.DataArray:
+        """The number of faces that surround each node.
+
+        Dimensions: ``(n_face, )``
+        """
+        if "n_faces_per_node" not in self._ds:
+            _populate_n_faces_per_node(self)
+
+        return self._ds["n_faces_per_node"]
+
+    @n_faces_per_node.setter
+    def n_faces_per_node(self, value):
+        """Setter for ``n_faces_per_node``"""
+        assert isinstance(value, xr.DataArray)
+        self._ds["n_faces_per_node"] = value
+
     @property
     def node_lon(self) -> xr.DataArray:
         """Longitude of each node in degrees.