We study the nonparametric covariance estimation of a stationary Gaussian field X observed on a lattice. To tackle this issue, a neighborhood selection procedure has been recently introduced. This procedure amounts to selecting a neighborhood m by a penalization method and estimating the covariance of X in the space of Gaussian Markov random fields (GMRFs) with neighborhood m. Such a strategy is shown to satisfy oracle inequalities as well as minimax adaptive properties. However, it suffers several drawbacks which make the method difficult to apply in practice. On the one hand, the penalty depends on some unknown quantities. On the other hand, the procedure is only defined for toroidal lattices. The present contribution is threefold. A data-driven algorithm is proposed for tuning the penalty function. Moreover, the procedure is extended to non-toroidal lattices. Finally, numerical study illustrate the performances of the method on simulated examples. These simulations suggest that Gaussian Markov random field selection is often a good alternative to variogram estimation.