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 by Mehta et al. (2014b) 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 GP (laGP) method allow accurate full-scale emulation. The new methodological contributions, which involve a multilevel global/local modeling approach, and a setwise 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-world satellites: the Hubble space telescope (HST) and the gravity recovery and climate experiment (GRACE).