Often in spatial regression problems, the covariates could be high-dimensional and have a non-linear relationship with the response. Furthermore, the functional relationship between the response and the covariates are often smoother than the spatial correlation. We propose a Gaussian spatial additive model on regular lattice where the large scale effects of the spatial covariates are modeled by smooth functions and the small scale spatial variability is modeled using a random field on the lattice. In order to facilitate linearity selection, we use penalized thin plate spline basis functions and derive a penalized h-likelihood method for simultaneous non-linearity selection and spatial adjustments. We derive novel estimating equations for estimating the precision parameters based on the profile h-likelihood. We demonstrate our method using Arsenic contamination data from Bangladesh.