How about doing Monte Carlo integration to compute the pushforward of the Weyl integration formula? Write a maximal torus for $G_2$ as $(\mathbb{R}/2 \pi \mathbb{Z})^2$ with coordinates $\alpha$ and $\beta$ corresponding to the long and short simple roots respectively. (Conveniently, in $G_2$ the root lattice and the weight lattice are the same.) Generate $10^4$ (say) points on this torus uniformly at random. For each of them, compute the trace of your $7$ dim rep: $$t:=2 \cos (\alpha+\beta) + 2 \cos(\alpha+2 \beta) + 2 \cos \beta +1.$$ (Check before using!) Divide $[-5,7]$ into $10^2$ buckets (say) according to the value of $t$ and sort your pairs $(\alpha, \beta)$ into these buckets. According to the Weyl integration formula, the volume of conjugacy class $(\alpha, \beta)$ is proportional to $$\mu := \sin^2 \left( \frac{\alpha}{2} \right) \sin^2 \left( \frac{\alpha+\beta}{2} \right) \sin^2 \left( \frac{\alpha+2\beta}{2} \right) \sin^2 \left( \frac{\alpha+3\beta}{2} \right) \sin^2 \left( \frac{2\alpha+3\beta}{2} \right) \sin^2 \left( \frac{\beta}{2} \right) d \alpha d \beta.$$ (Definitely check before using!) So sum up the values of $\mu$ in each bucket and plot the results.