Answer modified on 1 February 2020:
It's not true 'locally' in the sense that non-affine $f$'s satisfying this system of PDE can be constructed on some open sets in $\mathbb{R}^2$. This first order, determined PDE system is hyperbolic, so there are many local solutions. However, it turns out (see below) that all $C^3$ solutions with domain equal to $\mathbb{R}^2$ are affine. (The proof I give below does not work for solutions of lower regularity.)
Let $D\subset\mathbb{R}^2$ be a $1$-connected open domain on which there exists a $C^3$ mapping $f:D\to\mathbb{R}^2$ whose differential $\mathrm{d}f$ has constant, distinct singular values $0<\sigma_1<\sigma_2$. Because $D$ is simply connected, one can choose an orthonormal frame field $E_1,E_2$ on $D$ such that, at each point $p\in D$, the image vectors $F_i(p) = \mathrm{d}f\bigl(E_i(p)\bigr)$ are orthogonal and satisfy $|F_i(p)|=\sigma_i$.
Let $\omega = (\omega_1,\omega_2)$ be the dual coframing on $D$, which is of regularity type $C^2$. The $1$-forms $\eta_i = \sigma_i\,\omega_i$ for $i=1,2$ have the property that $(\eta_1)^2+(\eta_2)^2$, being the $f$-pullback of the flat metric on $\mathbb{R}^2$, must also be a flat metric.
Let $\omega_{12}$ be the connection $1$-form associated to the coframing $\omega$, i.e., it satisfies the structure equations
$$
\mathrm{d}\omega_1 = -\omega_{12}\wedge\omega_2
\qquad\text{and}\qquad
\mathrm{d}\omega_2 = \omega_{12}\wedge\omega_1\,.\tag1
$$
Write $\omega_{12} = -\kappa_1\,\omega_1 + \kappa_2\,\omega_2$. The function $\kappa_i$ is the curvature of the $E_i$-integral curve. Since $\omega_{12}$ is $C^1$, so are the functions $\kappa_i$. A straightforward computation shows that the $1$-form $\eta_{12}$ that satisfies the corresponding structure equations
$$
\mathrm{d}\eta_1 = -\eta_{12}\wedge\eta_2
\qquad\text{and}\qquad
\mathrm{d}\eta_2 = \eta_{12}\wedge\eta_1\,.\tag2
$$
is given by
$$
\eta_{12} = -(\sigma_1/\sigma_2)\,\kappa_1\omega_1
+ (\sigma_2/\sigma_1)\,\kappa_2\omega_2\,.
$$
Since $\sigma_1\not=\sigma_2$, the conditions $\mathrm{d}\omega_{12} = \mathrm{d}\eta_{12}=0$ (which hold because the domain metric and the $f$-pullback of the range metric are both flat) are equivalent to
$$
0 = \mathrm{d}(\kappa_i\,\omega_i) = \bigl(\mathrm{d}\kappa_i - {\kappa_i}^2\,\omega_{3-i}\bigr)\wedge\omega_i\,\qquad i = 1,2.\tag3
$$
Proposition: If $D = \mathbb{R}^2$, then $\kappa_1 \equiv \kappa_2 \equiv 0$, and $f$ is an affine map.
Proof: Suppose that, say, $\kappa_1$ be nonzero at some point $p\in\mathbb{R}^2$ and consider the value of $\kappa_1$ along the $E_2$ integral curve through $p$, which, since $E_2$ has unit length, is necessarily complete. Let $p(s)$ be the flow of $E_2$ by time $s$ starting at $p = p(0)$. Then (3) implies that the function $\lambda(s) = \kappa_1\bigl(p(s)\bigr)$ satisfies $\lambda'(s) = \lambda(s)^2$. Consequently,
$$
\kappa_1\bigl(p(s)\bigr) = \frac{\kappa_1\bigl(p(0)\bigr)}{1-\kappa_1\bigl(p(0)\bigr)s}.
$$
Hence $\kappa_1$ cannot be continuous along this integral curve, which is a contradiction. Thus, $\kappa_1$ and, similarly, $\kappa_2$ must vanish identically when $D = \mathbb{R}^2$. In particular, $\mathrm{d}\omega_i = 0$, from which one easily concludes that $f$ is affine. QED
More interesting, locally, is what happens near a point where $\kappa_1\kappa_2\not=0$. (There is a similar analysis when one of $\kappa_i$ vanishes identically that can safely be left to the reader, but see the note at the end.) One might as well assume that $\kappa_1\kappa_2$ is nowhere vanishing on $D$. Then one can write
$$
\kappa_1\,\omega_1 = \mathrm{d}u
\qquad\text{and}\qquad
\kappa_2\,\omega_2 = \mathrm{d}v
$$
for some $C^2$ functions $u$ and $v$ on $D$, uniquely defined up to additive constants.
Writing $\omega_1 = p\,\mathrm{d}u$ and $\omega_2 = q\,\mathrm{d}v$ for some non-vanishing functions $p$ and $q$, one finds that the structure equations (1), with $\omega_{12} = -\mathrm{d}u + \mathrm{d}v$, yield the equations
$$
p_v = - q \qquad\text{and}\qquad
q_u = -p.
$$
In particular, note that $p_v$ is $C^1$ and $p_{uv}-p = 0$.
Conversely, if $p$ be any nonvanishing $C^2$ function on a domain $D'$ in the $uv$-plane that satisfies the hyperbolic equation $p_{uv}-p=0$ and is such that $p_v$ is also nonvanishing on $D'$, then the $1$-forms
$$
\omega_1 = p\,\mathrm{d}u,\quad \omega_2 = -p_v\,\mathrm{d}v,\qquad
\omega_{12} = -\mathrm{d}u+\mathrm{d}v\tag4
$$
satisfy the structure equations of a flat metric, and so do
$$
\eta_1 = \sigma_1\,p\,\mathrm{d}u,\quad \eta_2 = -\sigma_2\,p_v\,\mathrm{d}v,\qquad
\eta_{12} = -(\sigma_1/\sigma_2)\,\mathrm{d}u+(\sigma_2/\sigma_1)\,\mathrm{d}v.\tag5
$$
Indeed, one now sees that the $1$-forms
$$
\begin{aligned}
\alpha_1 &= \cos(u{-}v)\,p\,\mathrm{d}u +\sin(u{-}v)\,p_v\,\mathrm{d}v\\
\alpha_2 &= \sin(u{-}v)\,p\,\mathrm{d}u -\cos(u{-}v)\,p_v\,\mathrm{d}v
\end{aligned}
$$
are closed, and therefore can be written in the form $\alpha_i = \mathrm{d}x_i$ for some $C^3$ functions $x_i$ on $D'$.
$$
(\mathrm{d}x_1)^2 + (\mathrm{d}x_2)^2 = (\alpha_1)^2 + (\alpha_2)^2 = (\omega_1)^2 + (\omega_2)^2
$$
and, hence, they define a $C^3$ submersion $x = (x_1,x_2):D'\to\mathbb{R}^2$ that pulls back the standard flat metric on $\mathbb{R}^2$ to the metric $(\omega_1)^2 + (\omega_2)^2$ on $D'$.
Similarly, setting $\rho = \sigma_2/\sigma_1$ and
$$
\begin{aligned}
\beta_1 &= \cos(u/\rho{-}\rho v)\,\sigma_1\,p\,\mathrm{d}u
+\sin(u/\rho{-}\rho v)\,\sigma_2\,p_v\,\mathrm{d}v\\
\beta_2 &= \sin(u/\rho{-}\rho v)\,\sigma_1\,p\,\mathrm{d}u
-\cos(u/\rho{-}\rho v)\,\sigma_2\,p_v\,\mathrm{d}v,
\end{aligned}
$$
one finds that $\mathrm{d}\beta_i = 0$ and hence there exist $C^3$ functions $y_i$ on $D'$ such that $\beta_i = \mathrm{d}y_i$. Set $y = (y_1,y_2)$.
Restricting to a subdomain $D''\subset D'$ on which $x$ is $1$-to-$1$ onto its image $D = x(D'')$ yields a domain on which $x^{-1}:D\to D''$ is a $C^3$ diffeomorphism. Now set $f = y\circ x^{-1}:D\to\mathbb{R}^2$, and one has a $C^3$ solution of the original PDE system.
This completely determines the structure of the 'generic' local $C^3$ solutions.
The case when one of the $\kappa_i$, say, $\kappa_1$, vanishes identically (so that the corresponding integral curves are straight lines) and the other is nonvanishing can easily be reduced to the normal form
$$
\omega_1 = \mathrm{d}u,\qquad \omega_2 = \bigl(p(v)-u\bigr)\,\mathrm{d}v,\qquad
\omega_{12} = \mathrm{d}v\tag6
$$
where, now, $p$ is a $C^2$ function of $v$, and the rest of the analysis goes through essentially unchanged.
