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):V\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.
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)$, 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$.
Thus, for $c_2>0$ and $c_1 > 2c_2$ sufficiently large, there will be functions $R(r)$ and $M(r)$ that satisfy the above equations, and hence they will determine the desired mapping $f$.
Thus, local solutions always exist carrying any desired point to any other. Moreover, it is clear that there is a 2-parameter family of local solutions carrying any given point in the domain to any given point in the range.
The existence of a global solution $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.