Addressing modern problems in environmental science often requires data collected on large, heterogeneous domains. Gaussian process models are useful to account for spatial correlation, which assist in prediction accuracy as well as correctly define errors for regression coefficients. Second order non-stationary models are often needed to capture the complexity of spatial relationships between environmental variables. However, the size of these datasets (n>10,000) often prohibits formation of the covariance matrix due to memory or computing time limits. We propose a new computational method that allows evolutionary spectrum models to be applied data on a partial grid with potentially missing values and to irregularly spaced data where the locations are approximated with a fine grid. We generalize the use of stochastic score approximations, circulant embedding, and pre-conditioned conjugate gradient to this non-stationary case and demonstrate how matrix-free kriging can be performed in this model. We apply our method to predict a fine grid of average daily temperatures across the US for several days in 2018.