---

This paper presents a nonparametric Bayesian method for estimating nonstationary covariance structures in big spatial datasets. The approach extends the Vecchia approximation and assumes conditional independence among ordered data points, leading to a sparse precision matrix and sparse Cholesky decomposition. This enables modeling an $n$-dimensional Gaussian process as a sequence of Bayesian linear regressions. Data ordering via maximum minimum distance improves model performance. Applying the grouping algorithm to ordered data removes weak dependencies and defines a block-sparse covariance structure, significantly reducing computational burden and enhancing accuracy. Simulations and real data analysis show that posterior samples from the proposed method yield narrower uncertainty intervals than those from ungrouped approaches.