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 <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mi>β</mml:mi> <mml:mo stretchy="false">(</mml:mo> <mml:mi mathvariant="bold-italic">x</mml:mi> <mml:mo stretchy="false">)</mml:mo> <mml:mo>=</mml:mo> <mml:munderover> <mml:mo>∑</mml:mo> <mml:mrow> <mml:mi>k</mml:mi> <mml:mo>=</mml:mo> <mml:mn>1</mml:mn> </mml:mrow> <mml:mi>K</mml:mi> </mml:munderover> <mml:msub> <mml:mi>β</mml:mi> <mml:mi>k</mml:mi> </mml:msub> <mml:mo> </mml:mo> <mml:msup> <mml:mtext>e</mml:mtext> <mml:mrow> <mml:mtext>i</mml:mtext> <mml:msub> <mml:mi>ω</mml:mi> <mml:mi>k</mml:mi> </mml:msub> <mml:mi mathvariant="bold-italic">x</mml:mi> </mml:mrow> </mml:msup> </mml:math> approximating the velocity field, with randomly sampled frequencies <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:msub> <mml:mi>ω</mml:mi> <mml:mi>k</mml:mi> </mml:msub> </mml:math> and amplitudes <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:msub> <mml:mi>β</mml:mi> <mml:mi>k</mml:mi> </mml:msub> </mml:math> trained to minimize a loss function. We include a physically motivated divergence penalty <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mrow> <mml:mo stretchy="false">|</mml:mo> </mml:mrow> <mml:mi mathvariant="normal">∇</mml:mi> <mml:mo>⋅</mml:mo> <mml:mi>β</mml:mi> <mml:mo stretchy="false">(</mml:mo> <mml:mi mathvariant="bold-italic">x</mml:mi> <mml:mo stretchy="false">)</mml:mo> <mml:msup> <mml:mrow> <mml:mo stretchy="false">|</mml:mo> </mml:mrow> <mml:mn>2</mml:mn> </mml:msup> </mml:math> , as well as a penalty on the Sobolev norm of <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mi>β</mml:mi> </mml:math> . 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.