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
$$
\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](https://en.wikipedia.org/wiki/Laguerre_polynomials#Recursive_definition,_closed_form,_and_generating_function), 
$$
\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](http://www.hep.anl.gov/czachos/a.pdf),  Thomas L. Curtright, David B. Fairlie, & Cosmas K. Zachos, *A Concise Treatise on Quantum Mechanics in Phase Space*, World Scientific, 2014.