Computational efficiency is at the forefront of many cutting edge spatial modeling techniques. Non-stationarity, on the other hand, has often been an unintentional feature of the approximations used in these spatial models. Deriving from the well known multivariate linear regression model, we propose a non-stationary and non-isotropic spatial model. In order to remain relevant with today’s massive datasets challenges, we apply the concept of nearest-neighbors to our normal-inverse-Wishart framework. The model, called Nearest-Neighbor Gaussian Process with Random Covariance matrices (NN-RCM) is developed for both univariate and multivariate spatial settings and allow for specific characteristics such as duplicate observations and missing data. In a first dive into nearest-neighbor models such as Nearest-neighbor Gaussian Processes (NNGP), we initially look at divide-and-conquer implementations. Ultimately, we devise a technique to use NNGP models on blocks of data to then aggregate the posterior inference without loss of information. The models are illustrated in a case study of albedo assessments over CONUS from the Geostationary Operational Environmental Satellites (GOES) East and West. First, the divide-and-conquer algorithm is used to output multi-day successive predictions of albedo assessments in a sequential updating scheme. By applying the distributed model on each day of data, one can concatenate the posterior parameters of the blocks of interest to output evolving predictions. Finally, we apply the bivariate NN-RCM model using each satellite as a source of information. The objective is to merge the albedo assessments while also quantifying the discrepancy between the sources.