Methods for fitting continuous spatial autoregressive moving average (ARMA) models to unequally spaced observations in two dimensions are reviewed and extended. These are models with rational two-dimensional spectra. Assuming Gaussian input noise and observational errors, maximum likelihood methods are used to estimate the ARMA parameters and the regression coefficients of the deterministic trend. When the number of observations is too large for exact maximum likelihood estimation, approximate maximum likelihood estimation is used based on nearest neighbors. Comparisons of nearest-neighbor methods with exact likelihood methods are presented. Predictions of the height of the field at unobserved points can be calculated with confidence intervals.