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.