Your solution for $2$ dimensions works for any dimension. You simply solve for co-ordinates $x^1, \dots, x^n$ such that the metric tensor $g_{ij}dx^idx^j$ satisfies
$$
g_{11} = g_{22},\ g_{12} = 0
$$
and set $f = x^1$, $g = x^2$. This is effectively parameterized isothermal co-ordinates, so your title is justified. You should be able to get solutions near a point by fixing a foliation of $2$-dimensional surfaces and solving for isothermal co-ordinates along each surface in such a way that the solution depends smoothly on the other $n-2$ co-ordinates.
ADDED: Here's a more explicit description in dimension 3: Start with co-ordinates $y^1$, $y^2$, and $y^3$ (where the level sets of $y^3$ give the foliation.) and a metric $g = g_{ij}dy^idy^j$. We want to solve for functions $\phi^1$ and $\phi^2$ such that if
$$
y^1 = \phi^1(x^1,x^2,x^3),\ y^2 = \phi^2(x^1,x^2,x^3),\ y^3 = x^3
$$
then the metric written with respect to $x^1, x^2, x^3$, $g = \hat{g}_{ij}dx^idx^j$ satisfies
$$
\hat{g}_{11} = \hat{g}_{22} \text{ and }\hat{g}_{12} = 0.
$$
A straightforward calculation shows that this is equivalent to
$$
g_{ab}\partial_1\phi^a\partial_1\phi^b = g_{ab}\partial_1\phi^a\partial_1\phi^b\text{ and }
g_{ab}\partial_1\phi^a\partial_2\phi^b = 0,
$$
where we sum over $1 \le a, b \le 2$. This is the same system of PDE's corresponding to finding co-ordinates $x^1, x^2, x^3$ such that $x^1, x^2$ are isothermal co-ordinates on any level set of $x^3$. In particular, since it involves only differentiation in the $x^1$ and $x^2$ directions and is elliptic, it can be solved for each value of the "parameter" $x^3$ using the standard approach for constructing isothermal co-ordinates. You can then verify that $\Phi = dx^1 + i dx^2$ satisfies $\langle d\Phi,d\Phi\rangle_g = 0$.
CORRECTION: Robert Bryant points out that the wrong problem is being solved above. The correct PDE is
$$
0 = \langle\partial\Phi,\partial\Phi\rangle_g = g^{ij}\partial_i\Phi\partial_j\Phi.
$$
Since $\Phi$ is complex-valued, this consists of 2 real equations for two real-valued functions. The linearization can be shown to be elliptic, if $\partial\Phi$ has maximal rank. Therefore, local solutions exist. However, the construction above of isothermal co-ordinates on each leaf of a foliation does not work.
CORRECTION to the CORRECTION: The linearized system is not elliptic. It is, however, generically hypoelliptic subelliptic. Given this, I believe that local solutions exist, but the only way I can see how to prove this is by proving smooth tame estimates for a solution to the hypoelliptic subelliptic PDE and using the Nash-Moser implicit function theorem.