In this paper, we investigate the massive multiple-input multiple-output orthogonal frequency division multiplexing channel estimation for low-earth-orbit satellite communication systems. First, we use the angle-delay domain channel to characterize the space-frequency domain channel. Then, we show that the asymptotic minimum mean square error (MMSE) of the channel estimation can be minimized if the array response vectors of the user terminals (UTs) that use the same pilot are orthogonal. Inspired by this, we design an efficient graph-based pilot allocation strategy to enhance the channel estimation performance. To reduce the computational complexity, we devise a novel two-stage channel estimation (TSCE) approach, in which the received signals at the satellite are manipulated with per-subcarrier space domain processing followed by per-user frequency domain processing. Moreover, the space domain processing of each UT is shown to be identical for all the subcarriers, and an asymptotically optimal vector for the per-subcarrier space domain linear processing is derived. The frequency domain processing can be efficiently implemented by means of the fast Toeplitz system solver. Simulation results show that the proposed TSCE approach can achieve a near performance to the MMSE estimation with much lower complexity.