We propose a computationally efficient approach to construct nonstationary covariance structures for Gaussian processes (GPs) in point-referenced geostatistical models. Current methods that impose nonstationarity on the covariance functions often suffer from computational bottlenecks, causing researchers to choose less appropriate alternatives in many applications. A main contribution of this paper is the development of a well-defined spatial process with the nonstationary covariance using multiple yet simple directed acyclic graphs (DAGs), leading to computational efficiency, flexibility, and interpretability. Inspired by graph-based approaches, we use a “bag of DAGs,” each of which is chosen to represent a different possible dependence structure. Combined with a domain partitioning idea from the scalable GP literature, our methods enjoy efficient computations. We are motivated by spatiotemporal modelling of air pollutants in which a DAG represents a prevailing wind direction causing some associated covariance in the pollutants. We establish Bayesian hierarchical models embedding the nonstationary GP and analyze spatiotemporal variability of fine particulate matter (PM2.5).