Spatial heteroskedasticity has been observed in many spatial data applications such as air pollution and vegetation. We propose a model, the volatility modulated moving average, to account for changing variances across space. This stochastic process is driven by Gaussian noise and involves a stochastic volatility field. It is conditionally non-stationary but unconditionally stationary: a useful property for theory and practice. We develop a discrete convolution algorithm as well as a two-step moments-matching estimation method for simulation and inference respectively. These are tested via simulation experiments and the consistency of the estimators is proved under suitable double asymptotics. To illustrate the advantages that this model has over the usual Gaussian moving average or process convolution, sea surface temperature anomaly data from the International Research Institute for Climate and Society are analysed.