The method of characteristics is a bit strange here because the equation is underdetermined, so one can't expect to be able to specify a solution by fixing initial data for $u$ and $v$ along a surface in $xyz$ space. It would be more efficient to put the equation in normal form and use the integration method appropriate to that normal form.
It actually helps to consider a slightly more general form of a linear, first-order PDE for $u$ and $v$ as functions of $x$, $y$, and $z$ as follows:
$$
A \frac{\partial u}{\partial x} + B \frac{\partial u}{\partial y} + C \frac{\partial u}{\partial z} + D \frac{\partial v}{\partial x} + E \frac{\partial v}{\partial y} + F \frac{\partial v}{\partial z}
+ G + Hu+Iv = 0,\tag0
$$
where $A,B,C,D,E,F,G,H,I$ are given functions of $(x,y,z)$. The key to understanding these equations is to understand the properties of the two vector fields
$$
U = A\,\frac{\partial\phantom{u}}{\partial x} +B\,\frac{\partial\phantom{u}}{\partial y} + C\, \frac{\partial\phantom{u}}{\partial z}
\quad\text{and}\quad
V = D\,\frac{\partial\phantom{u}}{\partial x} +E\,\frac{\partial\phantom{u}}{\partial y} + F\, \frac{\partial\phantom{u}}{\partial z}.
$$
In the generic situation, $U$ and $V$ will be linearly independent in the domain $\mathcal{D}$ in $xyz$-space that is being considered, and, moreover, their Lie bracket $W = [U,V]$ will be linearly independent from $U$ and $V$ in $\mathcal{D}$ as well. In this case, it is not hard to show that, by making a change of variables in $xyz$-space together with a change of variables of the form $(u,v)\mapsto (Pu+Qv+M,Ru+Sv+N)$, where $P,Q,R,S,M,N$ are functions of $xyz$ with $PS-RQ\not=0$, the equation can, locally, be put in the normal form
$$
u_x - v_y - x \,v_z = 0.\tag1
$$
In this case, a general solution is constructed by choosing $v(x,y,z)$ arbitrarily and then integrating $v_y+x\,v_z$ with respect to $x$ to find $u$.
In the special case that $U$ and $V$ are linearly independent in the domain $\mathcal{D}$ in $xyz$-space that is being considered, and their Lie bracket $W = [U,V]$ is a linear combination of $U$ and $V$ in $\mathcal{D}$, the local normal form of the equation under coordinate changes as above simpifies to
$$
u_x - v_y = 0,\tag2
$$
which has, as, general solution, $(u,v) = (P_y,P_x)$, where $P$ is any function of $xyz$.
Finally, if $U$ and $V$ do not simultaneously vanish but they are linearly dependent everywhere in $\mathcal{D}$, then, using the above change of coordinates, one can put the equation in the normal form
$$
u_x - P\,v = 0\tag3
$$
where $P$ is a function of $xyz$. On the open set where $P$ is nonvanishing, the normal form can be simplified to $u_x-v=0$, while, if $P$ vanishes identically, the normal form simplifies to $u_x=0$.
In each case, finding the normal form is a standard ODE problem in 3 variables in the domain $\mathcal{D}$. Any invocation of the method of characteristics would have to work in the case of these three or four normal forms, and its implementation would probably involve solving the ODE that put the equations in one of the above normal forms anyway.
Of course, there are intermediate cases, such as when $U$ and $V$ are linearly dependent along a subset of $\mathcal{D}$ or $U$, $V$, and $W = [U,V]$ are linearly dependent along a subset of $\mathcal{D}$, etc.
