Our main contribution is to introduce a novel method for Gaussian process (GP) modeling of massive datasets. The key idea is to build an ensemble of independent GPs that use the same hyperparameters but distribute the entire training dataset among themselves. This is motivated by our observation that estimates of the GP hyperparameters change negligibly as the size of the training data exceeds a certain level, which can be found in a systematic way. For inference, the predictions from all GPs in the ensemble are pooled to efficiently exploit the entire training dataset for prediction. We name our modeling approach globally approximate Gaussian process (GAGP), which, unlike most largescale supervised learners such as neural networks and trees, is easy to fit and can interpret the model behavior. These features make it particularly useful in engineering design with big data. We use analytical examples to demonstrate that GAGP achieves very high predictive power that matches or exceeds that of state-of-the-art machine learning methods. We illustrate the application of GAGP in engineering design with a problem on data-driven metamaterials design where it is used to link reduced-dimension geometrical descriptors of unit cells and their properties. Searching for new unit cell designs with desired properties is then accomplished by employing GAGP in inverse optimization.