Here is a sketch of a possible proof $L^{\infty}$ (uniform convergence) approximation of continuous maps to $S^3 = SU(2)$. (Actually, I need my maps to be a little better than continuous, because I need them not to be space filling curves.) As Andre observes, this also implies the result for $S^2$.
I suspect that someone
with a stronger background in topology could deal with the space
filling curves. I can't decide whether or not I think we should be able to push $L^{\infty}$ to $C^{\infty}$.
We will describe $S^3$ as $\{ (w,x,y,z) : w^2+x^2+y^2+z^2 = 1 \}$. We'll also think of this as the unit quaternions, writing $w+xi+yj+zk$. Stripped of the sophisticated language, the question is the following: We have a smooth map $\gamma: S^1 \to S^3$. Write $\gamma(\theta) = (w(\theta), x(\theta), y(\theta), z(\theta))$. Can we approximate $\gamma$ by maps $\gamma_1(\theta) = (w_1(\theta), x_1(\theta), y_1(\theta), z_1(\theta))$ where $w_1$, $x_1$, $y_1$ and $z_1$ are given by Fourier series of finite length, and the identity
$$w_1^2 + x_1^2 + y_1^2 + z_1^2 = 1$$
must hold identically.
Since $\gamma$ is not a space filling curve, it misses some point of $S^3$ and we will assume without loss of generality that it misses the quaternion $-1$. We can define the logarithm on $SU(2) \setminus \{ -1 \}$: The logarithm of $w+xi+yj+zk$ is the unique quaternion $pi+qj+rk$ with $\exp(pi+qj+rk) = w+xi+yj+zk$ and $p^2+q^2+r^2 < \pi^2$. For $a \in \mathbb{R}$, and $w+xi+yz+zk \neq -1$, we will define $(w+xi+yz+zk)^a$ to be $\exp(a \log(w+xi+yz+zk))$.
Choose a large integer $n$. For $0 \leq s \leq n-1$, define $a_s : S_1 \to \mathbb{R}$ to be the piecewise linear function which is $1/2$ on $[2 \pi (s-1/2)/n, 2 \pi (s+1/2)/n]$, interpolates linearly between $0$ and $1/2$ on $[2 \pi (s-3/2)/n, 2 \pi (s-1/2)/n]$ and $[2 \pi (s+1/2)/n, 2 \pi (s+3/2)/n]$m and is $0$ on the rest of $S^1$. Define $g_s(\theta) = \gamma(2 \pi s/n)^{a_s(\theta)}$ and define $g(\theta) = \prod_{s=1}^n g_s(\theta)$. As $n \to \infty$, the function $g(\theta)$ uniformly approaches $\gamma(\theta)$. Therefore, it is enough to choose $n$ so that $g$ is within $\epsilon/2$ of $\gamma$, and then approximate each $g_s$ within $\epsilon/(2n)$ by finite Fourier series.
If we wanted to do $C^{\infty}$ approximation rather than $L^{\infty}$, we would of course want the bump functions $a_s$ to be smooth rather than $C^{\infty}$. I expect this isn't hard, but it's not what I want to focus on.
In short, we have reduced to approximating maps of the form $\theta \mapsto q^{a(\theta)}$ where $q$ is a fixed quaternion not equal to $-1$. We may conjugate $q$ to be of the form $a+bi$, so we need to approximate $w(\theta) + x(\theta) i$ by $w_1(\theta) + x_1(\theta) i + y_1(\theta) j + z_1(\theta) k$.
Take $w_1$ and $x_1$ to be finite Fourier series approximations of $w$ and $x$, say within $\epsilon^2/100$. Multiplying by a scalar near $1$, we can further arrange that $w_1^2+x_1^2 \in (1-\epsilon^2/20, 1-\epsilon^2/10)$. Define $u(\theta) = 1-w_1(\theta)^2-x_1(\theta)^2$.
We are now reduced to finding finite Fourier series $y_1(\theta)$ and $z_1(\theta)$ with $y_1^2 + z_1^2 = u$. We no longer need to impose that $y_1$ and $z_1$ are close to $y=z=0$, since this is automatic from the identity $y_1^2 + z_1^2 = u < \epsilon^2/10$.
Lemma Let $u: S^1 \to \mathbb{R}_{>0}$ be a finite Fourier series. Then there are finite Fourier series $y_1$ and $z_1$ with $u=y_1^2+z_1^2$.
Proof: We will write $u$ as $\sum u_m q^m$, where $q=e^{i \theta}$. The condition that $u$ is real valued means that $u_{-m}=\overline{u_m}$.
Considering $u(q)$ as a function of a complex variable $q$, we have $u(\overline{q}^{-1}) = \overline{u(q)}$. So, if $\alpha$ is a zero of $u$, so is $\overline{\alpha}^{-1}$. Also, none of these zeroes can be on the unit circle, since $u$ takes positive values. Therefore we can write
$$u(q) = c \prod (1-q/\alpha_i )(1-q^{-1}/\overline{\alpha_i})$$
for some positive constant $c$ and some list of zeroes $\alpha_i$.
Define $v(q) = \sqrt{c} \prod (1-q/\alpha_i)$; then we have $v(e^{i
\theta}) \overline{v(e^{i \theta})} = u(\theta)$. Take $y_1$ and
$z_1$ to be the real and imaginary parts of $v$. $\square$
Note that this proof gives us no control over the sizes of $y'_1$ and $z'_1$, which is why I can't decide whether we should be able to get convergence in $C^k$ for $k \geq 1$.
We can also approximate maps $S^1 \to SO(3)$. Multiplying $\gamma$ by $\begin{pmatrix} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix}$ if necessary, we may assume that $\gamma$ has a square root $\sqrt{\gamma}: S^1 \to SU(2)$. Approximate $\sqrt{\gamma}$ by polynomials and then observe that the squaring map $SU(2) \to SO(3)$ is given by polynomials (I wrote them down here.)
I think this argument might actually be extendible to all semisimple
compact Lie groups. Here is a crude sketch.
The bump function trick let's us focus on approximating maps $g$ that land
in the exponential of a one dimensional subspace $\mathfrak{s}$ of the Lie algebra.
Let $T$ be the closure of $\exp(\mathfrak{s})$, so $T$ is a torus; let
its Lie algebra be $\mathfrak{t}$. Multiplying by a cocharacter $S^1 \to T$, we may assume that $g: S^1 \to T$ is null homotopic and thus has a logarithm $\log g: S^1 \to \mathfrak{t}$. Let $\mathfrak{r}_1$, $\mathfrak{r}_2$, ...., $\mathfrak{r}_m$ in $\mathfrak{t}$ be the one dimensional subspaces spanned by the simple roots: we can write $\log g$ as a linear combination $\log g = \sum h_i$ of maps $h_i: S^1 \to \mathfrak{r}_i$ and can exponetiate these to maps $\exp(h_i)$ from $S^1$ to one dimensional torii $R_i$. Each $R_i$ sits in an $SU(2)$ or $SO(3)$, so we can approximate $\exp(h_i)$ by the above argument.