Locally, this is always possible. Constructing such a coordinate system is equivalent to solving a first-order hyperbolic PDE system for two unknowns of two variables, so it always has local smooth solutions.
Here is a sketch of how this can be proved: You want to find a $g$-orthonormal coframing $g = {\omega_1}^2+{\omega_2}^2$ and a function $\lambda$ so that $\mathrm{e}^\lambda\omega_1$ and $\mathrm{e}^{-\lambda}\omega_2$ are both closed. Then they can be written as $\omega_1 = e^{-\lambda}\,\mathrm{d}x$ and $\omega_2 = e^{\lambda}\,\mathrm{d}y$ for some functions $x$ and $y$ and these are the coordinates you want.
Now, if you let $g = {\eta_1}^2+{\eta_2}^2$ be any local orthonormal coframing, then, up to a possible change of sign of $\omega_2$ (which won't affect the argument), the general orthonormal coframing of $g$ can be written as
$$
\omega_1 = \cos\theta\,\eta_1 -\sin\theta\,\eta_2
\quad\text{and}\quad
\omega_2 = \sin\theta\,\eta_1 +\cos\theta\,\eta_2\tag1
$$
for some arbitrary function $\theta$. Now consider the equations
$$
\mathrm{d}(e^\lambda\,\omega_1) = \mathrm{d}(e^{-\lambda}\,\omega_2) = 0.\tag2
$$
These constitute two first-order, quasi-linear equations for the two unknown function $\theta$ and $\lambda$. It is easy to show that they constitute a hyperbolic system. (The characteristics of a given solution are given by $\omega_1\pm\omega_2=0$.)
In the smooth category, the local existence of smooth solutions of such systems is well-known. The non-characteristic initial value problem (always locally solvable) is specified as follows: First choose an embedded smooth curve $C$ in the domain of the coframing $\eta$ and choose a smooth function $\theta$ along $C$ subject to the open condition that $\omega_1+\omega_2$ and $\omega_1-\omega_2$ do not vanish when pulled back to $C$. Then choose an arbitrary function $\lambda$ along $C$. Then there will be an open neighborhood $U$ of $C$ to which $\theta$ and $\lambda$ extend uniquely to $U$ satisfying $(2)$.