barycentric interpolation started

pyicon/pyicon_calc_xr.py

@@ -416,6 +416,77 @@ def xr_edges2edges_via_cell(ds_IcD, vn_e, scalar, dze='const'):
     out_vn_e = ()
     return out_vn_e
 
+## Mapping between cells and vertices
+
+def xr_calc_cell2vertex_coeff(ds_IcD):
+    raise(NotImplementedError)
+    dist_vector = ds_IcD.vert_cart_vec - ds_IcD.cell_cart_vec.isel(cell=ds_IcD.cells_of_vertex)
+    orientation = (dist_vector*ds_IcD.edge_prim_norm).sum(dim='cart')
+    dist_vector *= np.sign(orientation)
+    edge2cell_coeff_cc = (  ds_IcD.edge_prim_norm*ds_IcD.grid_sphere_radius
+                            * np.sqrt((dist_vector**2).sum(dim='cart'))
+                            / ds_IcD.dual_edge_length )
+    edge2cell_coeff_cc = edge2cell_coeff_cc.transpose('edge', 'nc', 'cart')
+    return edge2cell_coeff_cc
+
+def xr_cell2vertex(ds_IcD, da_b, cell2vertex_coeff_cc=None, fixed_vol_norm=None):
+    """ Interpolate cell variables on to vertices
+    Parameters
+    ----------
+    ds_IcD : xr.Dataset
+        A converted tgrid dataset (see pyic.convert_tgrid)
+
+    da_b : xr.DataArray
+        The data to be mapped onto cell vertices. This should be a scalar.
+
+    Returns
+    -------
+    ds_icd : xr.Dataset
+        A tgrid dataset compatible with pyicon functions
+
+    Notes
+    -----
+    This function uses hexagonal barycentric interpolation. Each vertex is
+    (normally) surrounded by 6 cell centres. The value at the vertex should
+    depend only upon these vertices. We express the position of the cell centre
+    as a sum of contributions from each point of the hexagon. We then use these
+    weights to interpolate. Hexagonal barycentric interpolant coefficients
+    aren't uniquely defined. Here, for the hexagon ABCDEF and point P, we
+    choose coefficients
+        x_A = [ABC]/[PBC]
+        y_B = [BCD]/[PCD]
+        z_C = [CDE]/[PDE]
+        u_D = [DEF]/[PEF]
+        v_E = [EFA]/[PFA]
+        w_F = [FAB]/[PAB]
+    where [xyz] represents the area of triangle xyz. We can then interpolate
+    a function f onto the point P using
+    f(P) = x_A f(A) + y_B f(B) + z_C f(C) + u_D f(D) + v_E f(E) + w_F f(F)
+
+    What isn't immediately clear to me is what happens at edge points where
+    there may be less valid centres surrounding the vertex? Perhaps we need to
+    do some masking following the calculation?
+
+    Furthermore, do we need to apply any volume normalisation? 
+    """
+    raise(NotImplementedError)
+    if cell2vertex_coeff_cc is None:
+        cell2vertex_coeff_cc = xr_calc_cell2vertex_coeff(ds_IcD)
+    ic0 = ds_IcD.adjacent_cell_of_edge.isel(nc=0).data
+    ic1 = ds_IcD.adjacent_cell_of_edge.isel(nc=1).data
+    ptp_vn = (
+        (
+            p_vn_c.isel(cell=ic0).rename({'cell': 'edge'})#.chunk(dict(edge=ic0.size))
+            * edge2cell_coeff_cc_t.isel(nc=0)
+        ).sum(dim='cart')                                            
+        +
+        (
+            p_vn_c.isel(cell=ic1).rename({'cell': 'edge'})#.chunk(dict(edge=ic0.size))
+            * edge2cell_coeff_cc_t.isel(nc=1)
+        ).sum(dim='cart') 
+    )
+    return ptp_vn
+
 ## Divergence
 
 def xr_calc_div_coeff(ds_IcD):