The answer is already 'no' when $n=2$, which implies that the answer is 'no' for all $n\ge 2$.
To see why, note that the local orthogonal diffeomorphisms that you describe depend essentially only on arbitrary functions of one variable: The reparametrizations $\nu_i$ and conformal diffeomorphisms of domains in $\mathbb{R}^n$, which either depend on constants (if $n>2$) or two functions of one variable (if $n=2$).
Thus, if we construct were to construct orthogonal diffeomorphisms $u:\mathbb{R}^2\to\mathbb{R}^2$ that essentially depend locally on a function of $2$ variables, then the generic such orthogonal diffeomorphism cannot be written in the above form.
To do this, let $\theta:\mathbb{R}^2\to\mathbb{R^2}$ be any smooth function that has compact support contained in the first quadrant (i.e., $\theta(x^1,x^2)=0$ unless $0<x^i<L$ for some $L>0$). Consider two functions $f_1>0$, and $f_2>0$ satisfying the differential equations
$$
\mathrm{d}\bigl(f_1( \cos\theta\,\mathrm{d}x^1+\sin\theta\,\mathrm{d}x^2) \bigr)
= \mathrm{d}\bigl(f_2( -\sin\theta\,\mathrm{d}x^1+\cos\theta\,\mathrm{d}x^2) \bigr) = 0
$$
and the initial conditions $f_1(x^1,0) = 1$ and $f_2(0,x^2) = 1$. These are linear first order PDE with noncharacteristic initial conditions for $f_1$ and $f_2$ of the form
$$
\sin\theta\,\frac{\partial f^1}{\partial x^1}
- \cos\theta\,\frac{\partial f^1}{\partial x^2} - A(x^1,x^2)\,\,f^1 = 0,\qquad f_1(x^1,0) = 1.
$$
and
$$
\cos\theta\,\frac{\partial f^2}{\partial x^1}
+ \sin\theta\,\frac{\partial f^2}{\partial x^2} - B(x^1,x^2)\,\,f^2 = 0,\qquad f_2(0,x^2) = 1.
$$
Looking at the characteristic curves (which are the coordinate curves outside of the square $S$ defined by $0<x^i<L$), one sees that there exist unique solutions $f^1$ and $f^2$ to these equations with initial conditions, and the homogeneous linearity implies that these two solutions are nowhere vanishing (and hence everywhere positive). In fact, outside of the half-slab $H_1$ defined by $0<x^1<L$ and $0<x^2$, one has that $f^1$ is identically $1$, while $f^2$ is identically $1$ outside of the half-slab $H_2$ defined by $0<x^2<L$ and $0<x^1$.
Now, since $\mathbb{R}^2$ is simply-connected, it follows that there exist functions $u^i$ on $\mathbb{R^2}$ such that
$$
\mathrm{d}u^1 = f_1( \cos\theta\,\mathrm{d}x^1+\sin\theta\,\mathrm{d}x^2)\quad\text{and}\quad
\mathrm{d}u^2 = f_2( -\sin\theta\,\mathrm{d}x^1+\cos\theta\,\mathrm{d}x^2).
$$
Thus, $u^1$ and $u^2$ are orthogonal coordinates on $\mathbb{R}^2$, and one finds that $(u^1,u^2):\mathbb{R}^2\to\mathbb{R}^2$ is a global diffeomorphism that is the identity outside of the union of the two half-slabs $H_1\cup H_2$.
Since $\theta$ was essentially arbitrary, subject to the condition of being compactly supported, it follows that these (global) orthogonal diffeomorphisms depend essentially on one arbitrary function of two variables and hence cannot be written, even locally, in the suggested form.
Once one has such a thing in two variables, by taking products, one can arrange to construct examples in arbitrarily high dimension that cannot be written in the proposed form.
Addendum: Concerning the higher dimensional case: In general, if $u:\mathbb{R}^n\to\mathbb{R}^n$ is an orthogonal coordinate system on $\mathbb{R}^n$, the $1$-forms
$$
\omega^i = \frac{\mathrm{d}u^i}{|\mathrm{d}u^i|}
$$
are orthonormal and integrable, i.e.,
$$
\mathrm{d}x^T{\circ}\mathrm{d}x = {\omega_1}^2 + \cdots + {\omega_n}^2
$$
and $\mathrm{d}\omega_i\wedge\omega_i = 0$. In particular, the orthogonality implies that, when we write the matrix equation
$$
\omega = A\,\mathrm{d}x
$$
(i.e., $\omega_i = A_{ij}\,\mathrm{d}x^i$), then $A$ takes values in the orthgonal matrices.
Conversely, given a mapping $A:\mathbb{R}^n\to\mathrm{O}(n)$ that satisfies the differential equations $\mathrm{d}\omega_i\wedge\omega_i = 0$, where $\omega_i = A_{ij}\,\mathrm{d}x^j$, one can then look for functions $f_i>0$ on $\mathbb{R}^n$ that satisfy $\mathrm{d}(f_i\omega_i) = 0$ (no sum on $i$). These will always exist locally, and, if the $\omega_i$ are not too 'wild', globally as well. Then there will exist functions $u_i:\mathbb{R}^n\to\mathbb{R}$ such that $\mathrm{d}u^i = f_i\,\omega_i$, and this will furnish an orthogonal coordinate system.
Thus, the essential problem is what one can say about the question of how many mappings $A:\mathbb{R}^n\to\mathrm{O}(n)$ satisfy the differential equations $\mathrm{d}\omega_i\wedge\omega_i = 0$. Here is where Cartan's results are relevant: What he proves is this: If $\Sigma\subset \mathbb{R}^n$ is a surface (i.e., $2$-dimensional) that is analytic and $A:\Sigma\to\mathrm{O}(n)$ is an analytic mapping such that the corresponding $1$-forms $\omega_i = A_{ij}\,\mathrm{d}x^j$ satisfy $\omega_i\wedge\omega_j$ is nonvanishing on $\Sigma$ for $i\not=j$, then there is an open $\Sigma$-neighborhood $U\subset\mathbb{R}^n$ to which $A$ extends as $A:U\to\mathrm{O}(n)$ so that the extended $1$-forms $\omega_i = A_{ij}\,\mathrm{d}x^j$ on $U$ satisfy the conditions $\mathrm{d}\omega_i\wedge\omega_i = 0$. Moreover, such an extension is germ-unique in the sense that any other extension $A':U'\to\mathrm{O}(n)$ of $A:\Sigma\to\mathrm{O}(n)$ agrees with $A:U\to\mathrm{O}(n)$ on some open $\Sigma$-neighborhood $V\subset U\cap U'$.
Thus, the 'space' of local orthogonal coordinates in $n$-dimensions depends on $\tfrac12n(n{-}1)$ functions of two variables in Cartan's sense, and the 'generic' such system is 'irreducible' in the sense that it is not a product of lower dimensional ones.
I think that, by choosing $\Sigma$ close to a generic $2$-plane in $\mathbb{R}^n$ and choosing $A:\Sigma\to\mathrm{O}(n)$ so that it stays very close to the identity (or even approaches the identity asymptotically as you leave every compact set in $\Sigma$) one might even be able to construct (or at least show the existence of) global orthogonal coordinates on $\mathbb{R}^n$ that depend on $\tfrac12n(n{-}1)$ functions of two variables and are 'irreducible' in any reasonable sense.
