the value of the function f is known on each grid point, to approximate the gradient of the function f, i.e. (df/dx, df/dy) which is the good choice, central difference? or one-sided difference? or some more complicate form?

(1). It depends on the nature of this function. (2). In principle, the central difference is more accurate than the one-sided difference, based on the Taylor series expansion. (2). Most numerical analysis books have a chapter on the numerical interpolation. And if you are dealing with test data, with some random distributions, then you will have to create a new approxmate smooth function first, in order to derive the general trend or gradient.

For finite differnece method,the central difference is usually used to calculate the gradient on account of the accuracy.Meanwhile,for finite volume method,it is better to use Gauss integral theorem over a control unit to calculate the gradient than central differnce.For example,when you calculate the derivative of velocity on the cell interface for viscous term in the N-S equations,the central differnce only takes account of one direction,i.e.you deal with the thin layer N-S eq or PNS,but the Gauss integral considers all direction,i.e.you deal with the full N-S eq.