Spatial generalized linear mixed models (SGLMMs) are popular and flexible models for spatial non-Gaussian data. They are useful for spatial interpolations as well as for fitting regression models that account for spatial dependence, and are commonly used in many disciplines such as epidemiology, atmospheric science, and sociology. Inference for SGLMMs is typically carried out under the Bayesian framework. Maximum likelihood inference is also available but computational issues often make it problematic, especially when high-dimensional spatial data are involved. Here we provide a computationally efficient projection-based maximum likelihood approach for routinely fitting SGLMMs. Our methodology is very general and applies to both discrete-domain (Gaussian Markov random field) as well as continuous-domain (Gaussian process) spatial models.