I did something similar for reporting stress in a linear elasticity finite element program, which approximates stress as a constant quantity across a linear tetrahedral element. It is often desirable to plot stress as a continuous quantity, because a smoothly varying quantity is more believable. I figured this was a good compromise. This approach should be valid for whatever quantity you wish to interpolate.
In short, I interpolated the stress at each of the vertices, while still honoring the elemental average values.
For each element
e, approximate the value at each node
a corresponding to that element using your basis functions
\sigma_e = \sum_a N_a\sigma_a d\Gamma
Integrate over the element.
\int_e\sigma_e d\Gamma = \int_e\ \sum_a N_a\sigma_a d\Gamma
\sigma_eV_e = \sum_a (\int_e N_a d\Gamma)\sigma_a
Therefore, you end up with an underdetermined system of
e equations and
v variables, where
e is the number of elements and
v is the number of vertices. An underdetermined system results in a least squares approximation that honors each of your equations. Technically, the nodes at the boundary are a linear extrapolation of that data, but the average quantities are still honored.
Here is a plot of the original "real" data (stress is a symmetric 3x3 matrix -- this is pressure which is a scalar reduction of that quantity) that is constant at each element.
Here is a plot of the interpolated "nodal" data that varies linearly between vertices.
EDIT: I only just saw that you wanted to take directional derivatives. However, I think you can do this as you now have a linear approximation of the quantity inside the element, which is certainly differentiable.