Obtaining accurate estimates of satellite drag coefficients in low Earth orbit is a crucial component in positioning and collision avoidance. Simulators can produce accurate estimates, but their computational expense is much too large for real-time application. A pilot study showed that Gaussian process (GP) surrogate models could accurately emulate simulations. However, cubic runtime for training GPs means that they could only be applied to a narrow range of input configurations to achieve the desired level of accuracy. In this paper we show how extensions to the local approximate Gaussian Process (laGP) method allow accurate full-scale emulation. The new methodological contributions, which involve a multi-level global/local modeling approach, and a set-wise approach to local subset selection, are shown to perform well in benchmark and synthetic data settings. We conclude by demonstrating that our method achieves the desired level of accuracy, besting simpler viable (i.e., computationally tractable) global and local modeling approaches, when trained on seventy thousand core hours of drag simulations for two real satellites: the HST and the GRACE.