Abstract: A highly scalable method for (non-)Gaussian random fields estimation is proposed. In particular, a novel (a) symmetric weight function based on nearest neighbors for the method of maximum weighted composite likelihood based on pairs (WCLP) is studied. The new weight function allows estimating massive (up to millions) spatial datasets and improves the statistical efficiency of the WCLP method using symmetric weights based on distances, as shown in the numerical examples. As an application of the proposed method, the estimation of a novel non-Gaussian random field named Tukey-hh random field that has flexible marginal distributions (possibly skewed and/or heavy-tailed) is considered. In an extensive simulation study the statistical efficiency of the proposed nearest neighbors WCLP method with respect to the WCLP method using weights based on distances is explored when estimating the parameters of the Tukey-hh random field. In the Gaussian case the proposed method is compared with the Vecchia approximation from computational and statistical viewpoints. Finally, the effectiveness of the proposed methodology is illustrated by estimating a large dataset of mean temperatures in South -America. The proposed methodology has been implemented in an open-source package for the R statistical environment.