Staked intelligent metasurface (SIM) based techniques are developed to perform two-dimensional (2D) direction-of-arrival (DOA) estimation. In contrast to the conventional designs, an advanced SIM in front of the receiving array automatically performs the 2D discrete Fourier transform (DFT) as the incident waves propagate through it. To arrange for the SIM to carry out this task, we design a gradient descent algorithm for iteratively updating the phase shift of each meta-atom in the SIM to minimize the fitting error between the SIM's response and the 2D DFT matrix. To further improve the DOA estimation accuracy, we configure the phase shifts in the input layer of SIM to generate a set of 2D DFT matrices having orthogonal spatial frequency bins. Extensive numerical simulations verify the capability of a well-trained SIM to perform 2D DFT. Specifically, it is demonstrated that the SIM having an optical computational speed achieves an MSE of $10^{-4}$ in 2D DOA estimation.