Modeling the spatial correlation of ground motion residuals, caused by coherent contributions from source, path, and site, can provide valuable loss and hazard information, as well as a more realistic depiction of ground motion intensities. The U.S. Geological Survey (USGS) software package, ShakeMap, utilizes a deterministic empirical approach to estimate median ground shaking in conjunction with observed seismic data. ShakeMap-based shaking estimates are used in concert with loss estimation algorithms to estimate fatalities and economic losses after significant seismic events around the globe. Incorporating the spatial correlation of ground motion residuals has been shown to improve seismic loss estimates. In particular, Park, Bazzuro, and Baker (Applications of Statistics and Probability in Civil Engineering, 2007) investigated computing spatially correlated random fields of residuals. However, for large scale ShakeMap grids, computational requirements of the method are prohibitive. In this work, a memory efficient algorithm is developed to compute the random fields and implemented using the ShakeMap framework. This new, iterative parallel algorithm is based on decay properties of an associated ground motion correlation function and is shown to significantly reduce computational requirements associated with adding spatial variability to the ShakeMap g round motion estimates. Further, we demonstrate and quantify the impact of adding peak ground motion spatial variability on resulting earthquake loss estimates.