The following grow longer than intended, but you do not need to read everything.

For simplicity we will, in the following, ignore the constant part and, as the signal is the sum of a x and y part (the data values are the outer sum), treat x and y data separately.

Fourier transforms describes periodic functions. To sample such a signal, you must choose the sample rate so that you sample at least one whole period. This implies that you sample in each dimension so that n+1 sample would replicate the first sample.

Now consider your signal. In x direction you have 2 frequencies, namely: 1/4 and 7/16. The GCD of these is: 1/16. Therefore in the x direction you should sample 16 times (note we start at 0 and ignore the constant part):

Table[2*Cos[x Pi/2] + 6*Sin[x Pi 7/8], {x, 0, 16 - 1}]

In the y direction we have a frequency of 5/16. We must sample so that the n+1 sample replicates the first one. Or in other words, the angle must be a multiple of 2Pi. This happens the first time for n+1== 16. Therefore, we sample at:

Table[4*Sin[y Pi 5/8], {y, 0, 16 - 1}]

The Fourier transforms of our data are:

fdatx = Fourier[datx];
fdaty = Fourier[daty];

Now, as Fourier transform can be look at as a base transformation, we need the new base functions. Here comes a small trick. The base functions should be Exp[+/- 2Pi I f x], with pos. and neg. frequencies. However for the points where we need these functions: x==0,1,2.. these have the same values as (note the minus sign and Sqrt[n] come due to the definition in MMA):

base = Table[-Exp[2 Pi I f i/16 ], {i, 0, 15}] / Sqrt[16]

Additionally they are in the same order as Fourier returns them. This order is: DC-value, lowest neg. freq., second lowest neg. freq.,.. either single highest freq. or 2 highest freq., dependingon even/odd number of data points,... second highest pos. freq.,... lowest pos. freq (note other definitions exchange neg/pos). That is really the case can be shown by explicitly constructing the base and compare them:

base0 = Join[{1}, Table[-Exp[2 Pi I  i x/16 ], {i, 7}], {-Exp[Pi I x]}, 
  Table[-Exp[-2 Pi I  i x/16 ], {i, 7, 1, -1}]]/Sqrt[16];
(base == base0) /. x -> Range[0, 15]
(* True *)

We can now assemble the functions for the x and y data:

funx[x_] = fdatx . base // Chop;
funy[x_] = fdaty . base // Chop;

and compare them to the original data:

(funx[x] /. x -> Range[0, 15] // Chop) == datx
(* True *)

(funy[x] /. x -> Range[0, 15] // Chop) == daty
(* True *)