It has taken me a while to find time to write a more comprehensive answer to the above question. It turns out that for general dimension $N$, the overdetermined PDE system involved is not involutive, but, depending on the assumptions one makes, it can be prolonged to an involutive system, and the conclusion is that, for suitably generic $D$ (where `suitable' depends on the eigenvalue matrix $\Lambda(x)$), the space of $f$ for which a local $g$ exists as above depends on $N^2$ functions of $1$ variable.
First, though, I should make explicit my assumptions: The OP didn't explicitly say this, but I will assume that $D$ is a constant diagonal matrix. (Since the OP never wrote $D(x)$ anywhere, I think that this is a reasonable assumption.). I will also assume that $\Lambda(x)$ is invertible for all $x$ in the domain of $f$, in particular that $\Lambda(x) = \mathrm{diag}\bigl(\lambda_1(x),\ldots,\lambda_N(x)\bigr)$, where the $\lambda_i$ are nonvanishing and $\lambda_1 < \cdots < \lambda_N$. I will also assume that $D = \mathrm{diag}\bigl(D_1,\ldots,D_N\bigr)$ has distinct eigenvalues and that certain inequalities (to be described below) also hold on the eigenvalues of $D\Lambda(x)$. (Some assumptions about $D$ need to be made in order to make the problem interesting. After all, if $D$ were the identity matrix, then one could always solve the problem by taking $g=f$.)
I'm now going to describe a method of constructing pairs $f$ and $g$ of mappings $f,g:X\to\mathbb{R}^N$ (where $X\subset\mathbb{R}^N$ is an open domain) whose Jacobians satisfy the above relations for a given constant matrix $D$. Then, assuming that $D$ is sufficiently 'generic' with respect to $\Lambda(x)$, I will prove that this method constructs all such pairs $(f,g)$, locally. (N.B.: Just to avoid any confusion, let me remark that I regard $\mathbb{R}^N$ as the set of column vectors of height $N$ with entries in $\mathbb{R}$.)
Now consider an `auxilliary' coordinate space $\mathbb{R}^N$, with coordinates $u = \textstyle{\left(\matrix{u^1\\\vdots\\u^N}\right)}$, and the equations
$$
\mathrm{d}X = T(u)\,\mathrm{d}u,\qquad \mathrm{d}F = T(u)\mu(u)\,\mathrm{d}u,\qquad \mathrm{d}G = T(u)D\mu(u)\,\mathrm{d}u
\tag1
$$
for smooth mappings $X,F,G:U\to\mathbb{R}^N$ and $T:U\to\mathrm{GL}(N,\mathbb{R})$,
with $\mu(u) = \mathrm{diag}\bigl(\mu_1(u^1),\ldots,\mu_N(u^N)\bigr)$ and where each $\mu_i$ is a nonvanishing smooth function on an open interval $U_i=(a_i,b_i)\subset\mathbb{R}$, with $\mu_i(U_i)$ to the left of $\mu_j(U_j)$ in $\mathbb{R}$ when $1\le i < j\le N$. Set $U= U_1\times U_2\times\cdots\times U_N\subset\mathbb{R}^N$. The idea is to choose a matrix-valued function $T:U\to \mathrm{GL}(N,\mathbb{R})$ so that it satisfies the first-order, linear system of PDE
$$
0 = \mathrm{d}\bigl(T(u)\,\mathrm{d}u\bigr),\qquad
0 = \mathrm{d}\bigl(T(u)\mu(u)\,\mathrm{d}u\bigr),\qquad
0 = \mathrm{d}\bigl(T(u)D\mu(u)\,\mathrm{d}u\bigr),
\tag2
$$
for then the desired $X,F,G:U\to\mathbb{R}^N$ satisfying (1) will exist (by Poincaré's Lemma).
Moreover, because $T$ is invertible, the map $X:U\to\mathbb{R}^N$ will be (locally) invertible (by the Inverse Function Theorem). Given any point of $U$, we can, by shrinking $U$ to a neighborhood of that point, ensure that $X:U\to X(U)\subset\mathbb{R}^N$ is one-to-one, so that $X^{-1}:X(U)\to U$ exists and is smooth. Then, by the Chain Rule, the functions $f = F\circ X^{-1}$ and $g = G\circ X^{-1}$, defined on the open set $X(U)\subset\mathbb{R}^N$, will satisfy the desired relations, with $\Lambda = \mu\circ X^{-1}$ and $S = T\circ X^{-1}$.
Of course, it may not be clear how to choose $T$ so that (2) will be satisfied. Because of the way $\mu$ was defined, (2) is equivalent to
$$
0 = \mathrm{d}\bigl(T(u)\bigr)\wedge\mathrm{d}u = \mathrm{d}\bigl(T(u)\bigr)\wedge \mu(u)\,\mathrm{d}u
= \mathrm{d}\bigl(T(u)\bigr)\wedge D\mu(u)\,\mathrm{d}u,
\tag{2'}
$$
which, depending on the values of $\mu(u)$ and $D$, can be as many as $N^3{-}N^2$ independent linear equations on the $N^2\cdot N = N^3$ first derivatives of the $n^2$ components of $T$. Fortuitously, it turns out that there is an $N^2$-parameter family of solutions of these equations that does not depend on how $\mu(u)$ and $D$ are specified: Set
$$
\mathrm{d}\bigl(T(u)\bigr) = P(u)\,\mathrm{diag}(\mathrm{d}u^1,\ldots,\mathrm{d}u^N),
\tag3
$$
where $P:U\to M_n(\mathbb{R})$ is any $N$-by-$N$ matrix-valued function on $U$.
Since $\mathrm{diag}(\mathrm{d}u^1,\ldots,\mathrm{d}u^N)\wedge\mathrm{d}u$ vanishes identically and since all diagonal matrices commute, it follows that the formula (3) does indeed satisfy (2'). In other words, the $i$th column of $T$ should depend on $u^i$ only. Thus, if we take $T$ of the form
$$
T(u) = \bigl(\,T_1(u^1)\ \ T_2(u^2)\ \ \cdots\ \ T_N(u^N)\,\bigr)
\tag 4
$$
where $T_i:U_i\to\mathbb{R}^N$ is smooth and $T_1(u^1)\wedge T_2(u^2)\wedge\cdots\wedge T_N(u^N)\not=0$ for all $u^i\in U_i$, then we have a solution of (2), which, in turn, yields a solution $(X,F,G)$ of (1). Note that choosing $T$ amounts to choosing $N$ functions of each variable $u^I$, i.e., $N^2$ functions of one variable all together.
Remark: The alert reader might note that we also get to choose the functions $\mu_i$ for $1\le i\le N$, which might lead one to think that the solutions $(f,g)$ obtained in this way actually depend on $N^2+N$ functions of one variable, but reparametrization in each $u^i$ variable separately will `wash out' of the formulae $(f,g) = (F\circ X^{-1}, G\circ X^{-1})$, so one should subtract "$N$ functions of one variable" worth of variability in the final formula, restoring the true "$N^2$ functions of 1 variable" generality.
Finally, let me address the question of whether this ansatz captures all of the local solutions $(f,g)$: The result is that, when $\Lambda$ and $D$ are sufficiently generic (i.e., satisfy some simple algebraic inequalities at every point), then every such solution $(f,g)$ is locally of the above form for some choice of $T$ satisfying (4). The proof of this result is not hard, but it takes a bit of fiddling and some care with notation. If there is interest, I'll put in the proof when I have some more time. (The key point is that, when $\mu$ and $D$ are sufficiently generic, every solution of (2') is of the form (3).)
