Abstract We measure the three-dimensional power spectrum (P3D) of the transmitted flux in the Lyman-α (Ly-α) forest using the complete extended Baryon Oscillation Spectroscopic Survey data release 16 (eBOSS DR16). This sample consists of ∼ 205,000 quasar spectra in the redshift range 2 ≤ z ≤ 4 at an effective redshift z = 2.334. We propose a pair-count spectral estimator in configuration space, weighting each pair by exp (ik · r), for wave vector k and pixel pair separation r, effectively measuring the anisotropic power spectrum without the need for fast Fourier transforms. This accounts for the window matrix in a tractable way, avoiding artifacts found in Fourier-transform based power spectrum estimators due to the sparse sampling transverse to the line-of-sight of Ly-α skewers. We extensively test our pipeline on two sets of mocks: (i) idealized Gaussian random fields with a sparse sampling of Ly-α skewers, and (ii) log-normal LyaCoLoRe mocks including realistic noise levels, the eBOSS survey geometry and contaminants. On eBOSS DR16 data, the Kaiser formula with a non-linear correction term obtained from hydrodynamic simulations yields a good fit to the power spectrum data in the range (0.02 ≤ k ≤ 0.35) h Mpc−1 at the 1 − 2σ level with a covariance matrix derived from LyaCoLoRe mocks. We demonstrate a promising new approach for full-shape cosmological analyses of Ly-α forest data from cosmological surveys such as eBOSS, the currently observing Dark Energy Spectroscopic Instrument and future surveys such as the Prime Focus Spectrograph, WEAVE-QSO and 4MOST.