Antwon(theswan),
As you say, let $f \in \mathbb{P_{2n+1}}$; then $f$ factors by long division into $f(x) = q(x) \phi_{n+1} + r(x)$. But then $f(x)$ is of degree $2n+1$, and $\phi_{n+1}$ is of degree $n+1$ and $q(x),r(x)$ are of degree at most $n$.
Now we take the $x_j$ to be the zeroes of $\phi_{n+1}$, so you have $f(x_j) = r(x_j)$ for $j = 0, \dots, n$.
I think you made it at least this far since you understand the factoring of the integral and understand the orthogonality principle behind the first term of the factorization being $0$.
Now the second integral is where the "clever" factoring comes into play.
We have:
$\int^b_a w(x)f(x)dx =^* \int^b_a w(x)r(x) =^{**} \color{red}{\sum^n_{j=0} P_j r(x_j)} = \sum^n_{j=0} \int^b_a w(x) L_{n,j}(x)r(x_j)$
where the first equivalence ($*$) comes from the orthogonality (which you recognize) $\langle\phi_{n+1},q_n\rangle = 0$ and the exactness of ($**$), which pertains to your question, comes from the Newton-Cotes formulas, which tell us that ($**$) is exact for any set of roots ($x_j$) for polynomials of degree at most $n$. And since we have chosen our roots wisely, i.e. the ($x_j$) of $\phi_{n+1}$, we can use what we deduced earlier, that $f(x_j) = r(x_j)$, which leads to the result.
I think your hang-up was due to the fact that the second integral in red is of degree at most $n$. For this reason, Newton-Cotes can be applied to draw exactness. This concept should jive with the idea of Gauss 'improving' the exactness of quadrature rules for polynomials up to $2n+1$ from the Newton-Cotes formulas, and this proof that you cite is using the fact that we already know exactness holds for one of these polynomials up to degree $n$.
NOTE: I highlighted the formula in red because I am having difficulty deciphering your notation. Explicitly, $\int^b_a w(x)f(x) dx \approx \sum^n_{j=0} P_j f(x_j)$ will be exact, if $P_j = \int^b_a w(x) p_j(x) dx$ for $j = 0, \dots, n$ and the $p_j$ are the Lagrange interpolating polynomials.
