In this paper, we present a novel approach to geostatistical filtering which tackles two challenges encountered when applying this method to complex spatial datasets: modeling the non-stationarity of the data while still being able to work with large datasets. The approach is based on a finite element approximation of Gaussian random fields expressed as an expansion of the eigenfunctions of a Laplace--Beltrami operator defined to account for local anisotropies. The numerical approximation of the resulting random fields using a finite element approach is then leveraged to solve the scalability issue through a matrix-free approach. Finally, two cases of application of this approach, on simulated and real seismic data are presented.