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.