We study the non-parametric covariance estimation of a stationary Gaussian field X observed on a lattice. To tackle this issue, we have introduced a model selection procedure in a previous paper. 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. Our contribution is threefold. We propose a data-driven algorithm for tuning the penalty function. Moreover, we extend the procedure to non-toroidal lattices. Finally, we study the numerical performances of this new method on simulated examples. These simulations suggest that Gaussian Markov random field selection is often a good alternative to variogram estimation.