We introduce a novel type of approximation spaces for functions with values in a nonlinear manifold. The discrete functions are constructed by piecewise polynomial interpolation in a Euclidean embedding space, and then projected pointwise onto the manifold. We show optimal interpolation error bounds with respect to Lebesgue and Sobolev norms. Additionally, we show similar bounds for the test functions, i.e., variations of discrete functions. Combining these results with a nonlinear Céa lemma, we prove optimal L2 and H1 discretization error bounds for harmonic maps from a planar domain into a smooth manifold. All these error bounds are also verified numerically.