An alternative approach to using a matrix representation for the rotation operator is to use clifford algebra instead. You might know clifford algebra in the guise of Pauli and Dirac matrices, from quantum mechanics.
A boost in the $i$th direction is described fully by a "rotor", which takes the form
$$q = \exp(-\gamma_0 \gamma_i \phi/2) = \cosh \frac{\phi}{2} - \gamma_0 \gamma_i \sinh \frac{\phi}{2}$$
If you're used to thinking of the Dirac matrices as, well, matrices, then feel free to write $I \cosh \frac{\phi}{2}$ instead. Here, however, I eschew treating these objects as matrices--I only need their multiplication law to use their clifford algebra properties--and as such, I consider such a term to be "scalar" compared to the vector space formed by the $\gamma_\mu$.
Now, if you're familiar with quaternions, you might recognize what we're doing: indeed, we can proceed exactly as in quaternions (and a good exercise is to derive quaternions using the Pauli algebra). A vector $v$ that is a linear combination $v = v^\mu \gamma_\mu$ can be boosted by
$$R(v) = q v q^{-1} = \exp(-\gamma_0 \gamma_i \phi/2) v \exp(\gamma_0 \gamma_i \phi/2)$$
Multiply all this out, and you get the usual formula for a boost.
Now let's take your example, and let's multiply two rotors: one boosting along $e_x = \gamma_1$ and another boosting along $e_y = \gamma_2$:
$$\left(\cosh \frac{\phi}{2} - \gamma_0 \gamma_1 \sinh \frac{\phi}{2} \right) \left(\cosh \frac{\phi}{2} - \gamma_0 \gamma_2 \sinh \frac{\phi}{2} \right)$$
Remember the clifford multiplication rules: $\gamma_\mu \gamma_\nu = -\gamma_\nu \gamma_\mu$ if $\mu \neq \nu$. and $\gamma_\mu \gamma_\mu = \eta_{\mu \mu}$ (no sum). The gammas are associative, so we can manipulate the expression to get
$$\cosh^2 \frac{\phi}{2} - \left(\sinh\frac{\phi}{2} \cosh \frac{\phi}{2}\right) [\gamma_0 \gamma_1 + \gamma_0 \gamma_2] - \gamma_1 \gamma_2 \sinh^2 \frac{\phi}{2}$$
using $\eta = (+,-,-,-)$ convention.
The presence of the $\gamma_1 \gamma_2$ term tells us that there is a rotation introduced when we compose boosts in this way (at least, when those boosts don't use the same plane).
Does it make sense to talk about the boost velocity for an operation that has both boost terms and pure rotation terms? Remember, this is not merely a boost and then a rotation, nor a rotation and then a boost--they're both happening together as you apply this operator.
If you have a definition of boost velocity that applies in this case, though, I'd be happy to turn the crank and compute it.
I'm presenting my overall thought process as I went through the derivation. There might well be a simpler way to go about it. We start with
$$
ds^2 = dv^2 - v^2 du^2 = v^2 [ \frac{dv^2}{v^2} - du^2 ] = v^2 d ( \ln v - u ) d ( \ln v + u )
$$
Define
$$
u+\ln v = U , \qquad u - \ln v = V ~~\implies~~ u = \frac{1}{2}(U+V) , \qquad \ln v = \frac{1}{2} ( U - V )
$$
The metric is
$$
ds^2 = - e^U e^{-V} dU dV = d(e^U) d(e^{-V})
$$
Next, we define
$$
e^U = x + t , \qquad e^{-V} = x - t .
$$
The metric is now
$$
ds^2 = (dx+dt)(dx-dt) = dx^2 - dt^2.
$$
This is the final form of the metric you want. We can now eliminate $U$ and $V$ and directly determine the relationship between $(t,x)$ and $(u,v)$
$$
v^2 = e^{U} e^{-V} = (x+t)(x-t) = x^2-t^2 \implies v = \sqrt{x^2-t^2}
$$
and
$$
u = \frac{1}{2} ( U + V ) = \frac{1}{2} \ln(x+t) - \frac{1}{2} \ln (x-t) = \frac{1}{2} \ln \frac{x+t}{x-t} .
$$
Inversely,
$$
x = v \cosh u , \qquad t = v \sinh u .
$$
Best Answer
In the general case you want the Cartan-Karlhede algorithm. It is an algorithm for producing a complete set of classifying invariants for a metric, expressed as functions of the coordinates. Given the components of the metric $g$ in the coordinates $x_1, x_2, \ldots$, the algorithm produces a list \begin{align} \Lambda & = \Lambda(x_i) \\ \Psi_k & = \Psi_k(x_i) \quad k = 0,\ldots, 4 \\ R_{kj} & = R_{kj}(x_i) \quad k,j = 0,\ldots, 2 \\ \Lambda_{00} & = \Lambda_{00}(x_i)\\ & \;\; \vdots \end{align} where each quantity is defined in a way that is coordinate independent. (The names here are standard notation, but what each of them is, is a little beyond the scope of this answer.) This is in contrast to a quantity like $g_{00}$ whose value at a point depends on your coordinates. Of course, expressed as a function of coordinates, $\Lambda$ and the others may look very different in various coordinate systems, but at corresponding points, the value is the same.
Then if we have two metrics given, we can run the algorithm on both. If the metrics are really the same, but different coordinates, the invariants must agree. This gives a system of equations, \begin{align} \Lambda(x_i) & = \Lambda'(y_j) \\ \Psi_k(x_i) & = \Psi_k'(y_j) \\ & \;\; \vdots \end{align} which may or may not have a solution, $x_i = x_i(y_j)$. (For example if you do this with two Schwarzschild metrics in the standard coordinates, you find that it is necessary that $m = m'$.) If there is a solution, this is your change of coordinates.
There is a caveat to the preceding. It may be that not all the equations are independent. In $n$ dimensions we need $n$ equations but the algorithm may produce fewer independent equations. This happens precisely when there is a symmetry in the spacetime. Then there cannot a unique change of coordinates, because at least one coordinate is superfluous. Indeed for the case of flat metrics the entire system is just $0 = 0$.
In this case the algorithm only establishes that there exists a change of coordinates, but you have to look at some other invariant information to find a coordinate change. (You will not be able to find a unique change of coordinates because there are many.) One piece of such information is the Killing vectors.
(This particular case is amenable to the brute force method demonstrated in the other answers, but a more complicated metric in more than two dimensions is not.)