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$.