This work is devoted to present the Legendre space-time spectral method for two-dimensional (2D) Sobolev equations. Considering the asymmetry of the first-order differential operator, the Legendre-tau-Galerkin method is employed in time discretization and its multi-interval form is also investigated. In the theoretical analysis, rigorous proof of the stability and $L^2(\mathrm{\Sigma })$-error estimates is given for the fully discrete schemes in both single-interval and multi-interval forms. Being different from the general Legendre-Galerkin method, we specifically take the Fourier-like basis functions in space to save the computing time and memory in the algorithm of the proposed method. Numerical experiments were included to confirm that our method attains exponential convergence in both time and space and that the multi-interval form can achieve improved numerical results compared with the single interval form.