Based on a quantitative version of the inverse function theorem and an appropriate saddle-point formulation, we derive a quasi-optimal error estimate for the finite element approximation of harmonic maps into spheres with a nodal discretization of the unit-length constraint. The error analysis is based on an equivalent formulation of the constrained discrete problem as a mixed finite element scheme. The resulting estimate holds under natural regularity requirements and appropriate geometric stability conditions on solutions. Extensions to other target manifolds, including boundaries of ellipsoids, are discussed.