Ocean buoy data in the form of high frequency multivariate time series are routinely recorded at many locations in the world's oceans. Such data can be used to characterise the ocean wavefield, which is important for numerous socio-economic and scientific reasons. This characterisation is typically achieved by modelling the frequency-direction spectrum, which decomposes spatiotemporal variability by both frequency and direction. State-of-the-art methods for estimating the parameters of such models do not make use of the full spatiotemporal content of the buoy observations due to unnecessary assumptions and smoothing steps. We explain how the multivariate debiased Whittle likelihood can be used to jointly estimate all parameters of such frequency-direction spectra directly from the recorded time series. When applied to North Sea buoy data, debiased Whittle likelihood inference reveals smooth evolution of spectral parameters over time. We discuss challenging practical issues including model misspecification, and provide guidelines for future application of the method.