To understand the local conditions, it's convenient to establish canonically associated local coordinate expressions for the quantities involved. Thus, let $(M^2,g,X)$ and $(N^2,h,Y)$ be as described and suppose that we want to test whether, for a given $p\in M$ and $q\in N$, there exists an open $p$-neighborhood $U\subset M$ and a local diffeomorphism $f:U\to N$ satisfying $f(p) = q$ with the desired properties. To simplify notation a little bit, let us fix orientations on $M$ and $N$ and require $f$ to be orientation-preserving. (We'll see what comes of this choice later.)
Then it is easy to show that there exist oriented, $p$-centered coordinates $(r,\theta):V\to (-\epsilon, \epsilon)\times (-\epsilon,\epsilon)$ on an open $p$-neighborhood $V\subset M$ and a smooth function $A:(-\epsilon,\epsilon)\to\mathbb{R}$ with $A(0)=0$ and $A'>0$ such that, on $V$, we have $g = \mathrm{d}r^2 + A'(r)^2\,\mathrm{d}\theta^2$ and $X = \partial/\partial\theta$. Similarly, there exist oriented, $q$-centered coordinates $(s,\phi):W\to (-\delta, \delta)\times (-\delta,\delta)$ on an open $q$-neighborhood $W\subset M$ and a smooth function $B:(-\delta,\delta)\to\mathbb{R}$ with $B(0)=0$ and $B'>0$ such that, on $W$, we have $h = \mathrm{d}s^2 + B'(s)^2\,\mathrm{d}\phi^2$ and $Y = \partial/\partial\phi$. Such adapted coordinates are locally unique.
[Note that, since $M$ is simply-connected (and, let's assume, connected, though the OP didn't include that condition), the functions $r$ and $\theta$ extend globally to $M$ uniquely so that $g = \mathrm{d}r^2 + a^2\,\mathrm{d}\theta^2$ and $X=\partial/\partial\theta$ where $a$ is a positive smooth function on $M$ satisfying $\mathrm{d}r\wedge\mathrm{d}a = 0$. The smooth mapping $(r,\theta):M\to\mathbb{R}^2$ is an immersion, but there is no reason to believe, under the given hypotheses, that it is an embedding, nor is it necessarily true that $a = A'(r)$ for some function $A:r(M)\to\mathbb{R}$.
The situation with $(N,h,Y)$ is similar.]
Supposing that an $f$ exists with all the specified properties, we can, by shrinking $\epsilon$, assume that $f(V)\subset W$ and hence, using the fact that $f_\ast(X) = Y$, conclude that
$$
f^*(s) = s\circ f = R(r)\quad\text{and}\quad f^*(\mathrm{d}\phi) = \mathrm{d}\theta + M(r)\,\mathrm{d}r.\tag1
$$
for some functions $R$ and $M$ on $(-\epsilon,\epsilon)$ with $R(0)=0$ and $R'>0$. This implies that, relative to the orthonormal coframings, we must have
$$
f^*\begin{pmatrix}\mathrm{d}s\\ B'(s)\,\mathrm{d}\phi\end{pmatrix}
=
\begin{pmatrix}R'(r) & 0\\ B'(R(r))M(r) & B'(R(r))/A'(r)\end{pmatrix}
\begin{pmatrix}\mathrm{d}r\\ A'(r)\,\mathrm{d}\theta\end{pmatrix}\tag2
$$
Now, the constancy of the singular values implies that, in particular, the determinant of the above coefficient matrix must be constant, i.e., that there must be a constant $c_2>0$ such that
$$
R'(r)B'(R(r))/A'(r) = c_2\,.\tag3
$$
Since $B(0) = R(0) = A(0) = 0$, we then integrate to get $B(R(r)) = c_2\,A(r)$. In particular, since $B$ is invertible, $R(r) = B^{-1}\bigl(c_2\,A(r)\bigr)$ for some positive constant $c_2$.
Now, the sum of the squares of the singular values of the coefficient matrix must be another constant $c_1 > 2c_2$ (so that the two constant singular values will be distinct) such that
$$
R'(r)^2 + B'(R(r))^2\,M(r)^2 + B'(R(r))^2/A'(r)^2 = c_1\,.\tag4
$$
Using the above formula for $R'(0) = c_2 A'(0)/B'(0)$, we see that, by taking $c_1$ sufficiently large, we can guarantee that the above equation for $M(r)$ has (two) real solutions on a neighborhood of $r=0$.
Conversely, for $c_2>0$ and $c_1 > 2c_2$ sufficiently large and $\epsilon>0$ sufficiently small, there will be functions $R(r)$ and $M(r)$ that satisfy the above equations (3) and (4) and the initial condition $R(0)=0$, and hence, via (1) and the initial condition $f(p)=q$, they will determine a unique mapping $f$ with the desired properties.
Thus, local mappings $f$ with constant, distinct singular values always exist carrying any desired point to any other. Moreover, it is clear that there is a 2-parameter family of such local mappings carrying any given point in the domain to any given point in the range.
The existence of a global such mapping $f:M\to N$ depends on the growth properties of the functions $A$ and $B$ and the validity of their domains. Little more can be said about this without more information or hypotheses about the functions $A$ and $B$.
