I believe the correct Wigner function for the eigenstates is half yours, so take it to be \begin{equation} W_n(x,p) = \frac{(-1)^n}{\hbar \pi} e^{- z/2} L_n(z ), \end{equation} where $z=4 H/\hbar\omega $.
You know that, since the resolution of the identity must be the resolution of the identity, $$ \sum_n W_n= \frac{1}{2\pi \hbar}=1/h, $$ which, indeed, holds (trivially checkable) by dint of the standard generating function of the Laguerre polynomials, $$ \sum_n t^n L_n(z)= \frac{e^{-tz/(1-t)}} {1-t} ~~. $$
Your sum then readily collapses to $$ \sum_n \frac{\bar{n}_\mathrm{th}^n}{(1+\bar{n}_\mathrm{th})^n} W_n(x,p)= \frac{e^{-z/2}}{\pi \hbar} \sum_n \left (\frac{- \bar{n}_\mathrm{th}}{1+ \bar{n}_\mathrm{th}} \right )^n L_n (z) = \frac{ (1+\bar{n}_\mathrm{th})}{\pi \hbar (1+2\bar{n}_\mathrm{th})} ~ e^{-z / 2 (1+2 \bar{n}_\mathrm{th}) } , $$ a gaussian in x and p with the requisite widths.
These are the basic maneuvers in phase-space quantization, Thomas L. Curtright, David B. Fairlie, & Cosmas K. Zachos, A Concise Treatise on Quantum Mechanics in Phase Space, World Scientific, 2014.