  We provide a computationally and statistically efficient method for
estimating the parameters of a stochastic covariance model observed on a
regular spatial grid in any number of dimensions. Our proposed method, which we
call the Debiased Spatial Whittle likelihood, makes important corrections to
the well-known Whittle likelihood to account for large sources of bias caused
by boundary effects and aliasing. We generalise the approach to flexibly allow
for significant volumes of missing data including those with lower-dimensional
substructure, and for irregular sampling boundaries. We build a theoretical
framework under relatively weak assumptions which ensures consistency and
asymptotic normality in numerous practical settings including missing data and
non-Gaussian processes. We also extend our consistency results to multivariate
processes. We provide detailed implementation guidelines which ensure the
estimation procedure can be conducted in O(n log n) operations, where n is the
number of points of the encapsulating rectangular grid, thus keeping the
computational scalability of Fourier and Whittle-based methods for large data
sets. We validate our procedure over a range of simulated and real-world
settings, and compare with state-of-the-art alternatives, demonstrating the
enduring practical appeal of Fourier-based methods, provided they are corrected
by the procedures developed in this paper.
