There is an excellent explanation of this in Chapter 4 of L. N. Trefethen's Approximation Theory and Approximation Practice (henceforth ATAP; the first 6 chapters are available for free online). I will try to summarize it in your notation, but I recommend reading his explanation.
Chebyshev interpolation
Interpolation at the Chebyshev points gives your $\Pi^{CG}_n f$, which can also be written in terms of the first $n$ Chebyshev polynomials:
$$\Pi^{CG}_n f = \sum_{k=0}^n c_k T_k(x).$$
The question then becomes:
How are the coefficients $c_k$ of the degree-$n$ Chebyshev interpolant related to the coefficients $\hat{f}_k$ of the infinite Chebyshev series?
Aliasing
It turns out that on the Chebyshev grid with $n+1$ points, the following Chebyshev polynomials all take exactly the same values:
$$T_{2jn\pm m}(x), \ \ \ \ j = 0, 1, 2 \dots$$
The above statement is valid as long as $0\le m \le n$ and where the subscript is positive (of course). See Thm. 4.1 of ATAP. Thus each $T_k$ for $k > n$ is "aliased" to appear like some corresponding $T_k$ with $k\le n$.
Relation between the Chebyshev interpolant and the Chebyshev series
Once you know about aliasing, the connection between the coefficients $c_k$ and $\hat{f}_k$ is immediate. Since the Chebyshev interpolant and the infinite Chebyshev series must be equal at the Chebyshev points, we have
$$
\sum_{k=0}^n c_k T_k(x_j) = \sum_{k=0}^\infty \hat{f}_k T_k(x_j),
$$
where $x_j$ is any of the Chebyshev points.
But from the aliasing result, we know that each term with $k>n$ in the right-hand-side sum can be rewritten in terms of $T_k(x_j)$ for some $k\le n$. The result is Thm. 4.2 of ATAP, which gives the formulae (for the degree-$n$ interpolant coefficients)
\begin{align}
c_0 - f_0 & = \sum_{j=1}^\infty f_{2jn} \\
c_n - f_n & = \sum_{j=1}^\infty f_{(2j+1)n} \\
c_k - f_k & = \sum_{j=1}^\infty (f_{2jn+k} + f_{2jn-k}) & 1\le k \le n-1.
\end{align}
So the difference between the coefficients of the Chebyshev interpolant and the Chebyshev projection is given by certain sums of coefficients of the remainder of the Chebyshev series.
One of the major themes of ATAP is that the Chebyshev interpolant is almost as good as the Chebyshev projection for approximation purposes; see Theorems 7.2, 8.2, and 16.1 therein. The great virtue of the interpolant is that it is much easier to compute.
Best Answer
Evaluating polynomials of arbitrarily large degree in a Chebyshev basis is practical, and provably numerically stable, using a barycentric interpolation formula. In this case, extended precision isn't needed, even for order 1,000,000 polynomials. See the first section of this paper and the references, or here (Myth #2) for more details. I'll summarize briefly. Let's say you have a polynomial $f$ in a Chebyshev basis, and you know its values at the Chebyshev nodes $$ f_j = f(x_j) $$ $$ x_j = \cos\left(\frac{\pi j}{N}\right), \;\; 0\leq j\leq N $$
Then for any $x$ which isn't one of the Chebyshev nodes $x_j$, we have $$ f(x) = \frac{\displaystyle \sum_{j=0}^N \frac{w_j}{x-x_j}f_j}{\displaystyle \sum_{j=0}^N \frac{w_j}{x-x_j}}, $$ where $$ w_j = \left\{ \begin{array}{cc} (-1)^j/2, & j=0\text{ or }j=N,\\ (-1)^j, & \text{otherwise} \end{array} \right. $$
I believe a naive implementation of the above for various values of $x$ provides a stable algorithm that does not suffer the numerical difficulties encountered when trying to sum up the Chebyshev polynomials directly. The key is to work with the representation of the function by its values, not by its coefficients in a Chebyshev basis.