The most conceptually simple approach is to select an interpolation scheme, interpolate through the data points that you have, then compute the derivative of the interpolant. This is a natural generalization of the typical finite-difference differentiation formulas.
While you don't mention your exact requirements (in particular: do you need the derivatives at the mesh points, or within the mesh triangles?), the simplest approach is to interpolate the function values linearly in each triangle, and estimate the derivative as the derivative of the linear interpolant; the interpolant will have a piecewise constant derivative, and be discontinuous at the mesh points and edges, where you can take linear combinations of adjacent interpolants.
You can also look at other inteprolation methods, like radial basis functions.
Also note that many software packages support 2d interpolation themselves, in which case it might not be necessary to do the interpolation step yourself. If a library can interpolate a function, but doesn't calculate the interpolant's derivatives, you can recover the derivatives by using a finite-difference approximation on the interpolant.