The short answer is "The operators that are diagonal in the basis of coherent states are those that have a smooth Glauber–Sudarshan P transform".
The Wigner-Weyl transform beautifully maps all kinds of operators on Hilbert space (density matrices or otherwise, coherent-state-diagonal or not) to real-valued functions on phase space. However, this doesn't directly tell us which ones are diagonal since all "nice" operators can be represented in this way. (See here for some discussion of regularity conditions.)
But we might take a naive stab at what an operator would look like if it were diagonal in the basis of coherent states $\vert \alpha \rangle$, where $\alpha = (x,p)$ is a variable ranging over phase space. First, note that we can write the position operator in the position basis as
$$
\hat{X}=\int \mathrm{d}x \, x \vert x \rangle \langle x \vert.
$$
So, given an arbitrary real-valued (classical) function $f(\alpha)$ on phase-space (e.g. position, momentum, or energy) consider the following method of defining a corresponding operator:
$$
\hat{\Omega}[f] := \int \mathrm{d}\alpha\, f(\alpha) \vert \alpha \rangle \langle \alpha \vert.
$$
(The integral on the right should be normalized by a factor of $1/2 \pi$ but I'm going to ignore normalization in everything I have to say.) Intuitively, this operator acts on a state vector by "measuring it" in the wavepacket basis and weighting the corresponding components by the value of the function at that point in phase space. In fact, one can check that $\hat{\Omega}[X] = \hat{X}$, where $X(\alpha) = X(x,p) = x$ is the position function on phase space.
How does this relate to the Wigner-Weyl transform? Well, the W transform of $\hat{\Omega}[f]$ is
$$
\begin{align}
W\left[\hat{\Omega}[f]\right](x,p) &= \int \mathrm{d}\Delta x \, e^{i p \Delta x} \langle x + \frac{\Delta x}{2} \vert \hat{\Omega}[f] \vert x - \frac{\Delta x}{2} \rangle \\
&= \int \mathrm{d}\Delta x \int \mathrm{d}\alpha\, e^{i p \Delta x} f(\alpha) \langle x + \frac{\Delta x}{2} \vert \alpha \rangle \langle \alpha \vert x - \frac{\Delta x}{2} \rangle.\\
&= \int \mathrm{d}\alpha\, e^{-\left[\alpha -(x,p)\right]^2} f(\alpha).\\
&= (f \circ g) (x,p)
\end{align}
$$
where $\circ$ denotes the convolution and $g(\alpha) = e^{-\alpha^2}$ is a Gaussian kernel. Another way to write this is to use the convolution theorem:
$$
\tilde{W}\left[\hat{\Omega}[f]\right](\xi) = \tilde{f}(\xi) \tilde{g}(\xi)
$$
where $\tilde{\phantom{g}}$ denotes the Fourier transform, which makes $\tilde{g}(\xi) = e^{-\xi^2/4}$, where and $\tilde{W}$ is the characteristic function associated with $W$. [Usually the "symplectic Fourier transform" is used for the characteristic function, but that's differs from the above only by reparameterizing as $\xi=(\xi_1,\xi_2)\to(-\xi_2,\xi_1)$. Note also that the factors like the 4 in $\tilde{g}(\xi) = e^{-\xi^2/4}$ is dependent on the choice of Fourier convention, and whether $\hbar = 1/2,1,$ or $2$.]
So we see that the W transform of a "diagonal" operator like $\hat{\Omega}[f]$ is just $f$ smoothed by a convolution with the Gaussian kernel $g$. But that's exactly what makes $f$ the Glauber-Sudarshan P transform of an operator $\hat{\Omega}[f]$. Indeed, the definition of $\hat{\Omega}[f]$ given above is the same as the defining equation for a P transform for a density matrix, only generalized to operators that aren't necessarily density matrices.
Now, all smooth operators on the Hilbert space will have W transforms associated with them, but not all such W transforms can be de-convolved to a smooth P transform. In particular, if we move to Fourier space, de-convolution is just multiplication by $e^{+\xi^2/4}$ which, in general, will produce a function on which the inverse Fourier transform does not converge. There are ways to try construct a P transform in this case, but they are highly singular, involving infinite derivatives of Dirac delta functions. (See Wikipedia, references therein, and especially Bonifacio et al. and Sudarshan.)
