We investigate the use of spatial interpolation methods for reconstructing the horizontal near-surface wind field given a sparse set of measurements. In particular, random Fourier features is compared with a set of benchmark methods including kriging and inverse distance weighting. Random Fourier features is a linear model beta(x)= Sigma(K)(k-1) beta(k) e(i omega kx) approximating the velocity field, with randomly sampled frequencies omega(k) and amplitudes beta(k) trained to minimize a loss function. We include a physically motivated divergence penalty vertical bar del. beta(x)vertical bar(2), as well as a penalty on the Sobolev norm of beta. We derive a bound on the generalization error and a sampling density that minimizes the bound. We then devise an adaptive Metropolis-Hastings algorithm for sampling the frequencies of the optimal distribution. In our experiments, our random Fourier features model outperforms the benchmark models.