how to interpret the 2D FFT
Yes, due to the separability of the kernel :

$$e^{-j \left(\frac{2\pi}{N_1} n_1 k_1 + \frac{2\pi}{N_2} n_2 k_2 \right) } = e^{-j \frac{2\pi}{N_1} n_1 k_1} \cdot e^{-j \frac{2\pi}{N_2} n_2 k_2}$$

the 2D-DFT sum 
$$X[k_1,k_2] = \sum_{n_1} \sum_{n_2} f[n_1,n_2] e^{-j \frac{2\pi}{N_1} n_1 k_1} e^{-j \frac{2\pi}{N_2} n_2 k_2} $$

can be implemented using row-column (or column-row) decompositions. 

$$X[k_1,k_2] = \sum_{n_1} \left( \sum_{n_2} f[n_1,n_2] e^{-j \frac{2\pi}{N_2} n_2 k_2} \right) e^{-j \frac{2\pi}{N_1} n_1 k_1} $$

The following MATLAB/OCTAVE code shows how to apply :

N = 8;            % length of columns-rows 
x = randn(N,N);   % row data

X = fft(x,N) ;    % 1D-fft along columns of x[n1,n2]
X = fft(X.',N).' ;% 1D-fft along rows of intermediate X[k1,k2]