Large spatial data sets are now ubiquitous in environmental science. Fine spatial sampling or many observations across large domains provides a wealth of information and can often address new scientific questions.  The richness and scale of large datasets, however,  often reveal heterogeneity in spatial processes that add more complexity to a statistical analysis.  Our new approach is to estimate spatially varying covariance parameters in a local manner but then encode these into a sparse Markov random field model for a global representation. This strategy makes it possible to estimate and then simulate (unconditional) non-stationary Gaussian processes.  This approach is illustrated for the emulation of surface temperature fields from an ensemble of climate model experiments (Community Earth System Model Large Ensemble) and showcases efficient computation using parallel methods and sparse matrices.  Current methods in spatial statistics inherit the foundational work in nonparametric regression and splines that was pioneered by Grace Wahba and others. This talk will also