Geophysical potential-field data generally are observed at scattered sampling points and are contaminated by measurement and preprocessing errors. Analysis and further treatment require a representation of the underlying potential fields on a regular grid. In the presence of noise, an approximation of the data is more appropriate than an exact interpolation. Analysis of the smoothness of potential fields indicates an approximation by band-limited functions. The resulting least-squares problem can be formulated as a linear system of equations in which the system matrix is of block-Toeplitz type. This special structure allows a computationally attractive and robust solution via the fast Fourier transform and the conjugate gradient algorithm. On this basis, we present a gridding method that incorporates additional knowledge about the physical properties of potential fields and the statistics of the data. Experimental results for synthetic and real-world gravity data demonstrate good performance for highly scattered noisy data.