Over the past five years, Nearest Neighbor Gaussian Processes (NNGP) arose as a computationally scalable method forspatial statistical models, but remain hampered by problems caused by the behavior of Markov Chain Monte-Carlo (MCMC)algorithms. Several approaches allow to alleviate those issues but they restrict the flexibility of the original model.This work keeps the ``jack of all trades" basic model and tackles its MCMC weak points with several strategies. Therobustness and efficiency of high-level parameters estimation is boosted using interweaving strategies. Lower-leveloperations are parallelized using Chromatic Sampling. Efficient Hamiltonian methods are developed for NNGP models.In a second time, the versatility of the NNGP model is used in order to tackle nonstationary modeling. An originalparametrization and model architecture are proposed in order to ease model interpretation and selection while capturingcomplex nonstationarity patterns. An innovative MCMC strategy based on Hamiltonian methods and Nested Interweaving isproposed.