The vector fields $X$ that satisfy this condition are highly constrained, as there exist only a finite-dimensional family of $C^5$ solutions in a neighborhood of any given point in the plane. Here is how one can see this:
Using the standard coordinates, we have $g = dx^2 + dy^2$. Suppose that $X$ is defined on a simply-connected domain $D$ in the plane. Then, since $X$ is divergence-free, it can be written in the form
$$
X = u_y\,\frac{\partial}{\partial x} - u_x\,\frac{\partial}{\partial y}
$$
for some function $u$ on $D$. Define $\mathcal{L}_X^2g = \mathcal{L}_X(\mathcal{L}_X g)$ and, inductively $\mathcal{L}_X^{k+1}g = \mathcal{L}_X(\mathcal{L}^k_X g)$. As has been observed, we have $\mathrm{tr}_g(\mathcal{L}_X g) = 0$, since $X$ is divergence free. A direct computation gives
$$
\mathrm{tr}_g(\mathcal{L}^2_X g) = 2\left((u_{xx}-u_{yy})^2 + (2u_{xy})^2\right) = 8c^2
$$
for some constant $c\ge0$. If $c=0$, then either $X = a_0\,(y-y_0)\partial_x - a_0\,(x-x_0)\partial_y$ for some constants $a_0,a_1,a_2$, and the flow of $X$ is rotation around $(x_0,y_0)$ or $X = a_1\partial_x + a_2\partial_y$ and hence the flow of $X$ is simply translation.
Thus, from now on assume that $c>0$. Then the Hessian matrix of $u$ is of the form
$$
\begin{pmatrix} u_{xx} & u_{xy}\\u_{xy} & u_{yy}\end{pmatrix}
= \begin{pmatrix} h + c\,\cos v & c\,\sin v\\c\,\sin v & h-c\,\cos v\end{pmatrix}
$$
for some functions $h$ (unique) and $v$ (unique up to an additive constant that is an integral multiple of $2\pi$). The identities $(u_{xx})_y = (u_{xy})_x$ and $(u_{xy})_y = (u_{yy})_x$ imply that
$$
\mathrm{d}h = c\bigl(v_y\,\cos v-v_x\,\sin v\bigr)\,dx
+ c\bigl(v_x\,\cos v+v_y\,\sin v\bigr)\,dy,
$$
and, since $c\not=0$, the identity $\mathrm{d}(\mathrm{d}h) = 0$ gives
$$
(v_{xx}-v_{yy}+2v_xv_y)\,\cos v + (2v_{xy}-{v_x}^2+{v_y}^2)\,\sin v = 0.
$$
Thus, the Hessian of $v$ can be parametrized as
$$
\begin{pmatrix} v_{xx} & v_{xy}\\v_{xy} & v_{yy}\end{pmatrix}
= \begin{pmatrix}
q+p\,\sin v - v_xv_y & \tfrac12({v_x}^2-{v_y}^2)-p\,\cos v\\
\tfrac12({v_x}^2-{v_y}^2)-p\,\cos v & q-p\,\sin v + v_xv_y\end{pmatrix}
$$
for some functions $p$ and $q$. Now, a straightforward calculation yields that
$$
\mathrm{tr}_g(\mathcal{L}^3_X g) = 0
\quad\text{while}\quad
\mathrm{tr}_g(\mathcal{L}^4_X g) = 32\,c^4 - 8c^2\,(2h-u_xv_y+u_yv_x)^2.
$$
It follows that we must have $h = b + \tfrac12(u_xv_y-u_yv_x)$ for some constant $b$. Differentiating this equation and using the various formulae above allows one to solve for $p$ and $q$ in the form.
$$
p = \frac{P(\cos v, \sin v, v_x,v_y,u_x,u_y)}{\bigl(({u_x}^2-{u_y}^2)\cos v +2u_xu_y\sin v\bigr)}
$$
and
$$
q = \frac{Q(\cos v, \sin v, v_x,v_y,u_x,u_y)}{\bigl(({u_x}^2-{u_y}^2)\cos v +2u_xu_y\sin v\bigr)},
$$
where $P$ and $Q$ are certain polynomials in their arguments, which I won't write out here. Moreover, using these formulae, we find that
$$
\begin{align}
\mathcal{L}_X g
&= 2c\bigl(\sin v\,dx^2-2\,\cos v\,dx\,dy-\sin v\,dy^2\bigr)\\
\mathcal{L}^2_X g
&= 4c\bigl((c{+}b\,\cos v)\,dx^2+2b\sin v\,dx\,dy + (c{-}b\,\cos v)\,dx^2\bigr)
\end{align}
$$
and
$$
\mathcal{L}^3_X g = 4(c^2{-}b^2)\,\mathcal{L}_X g,
$$
from which it follows by induction that $\mathcal{L}^{k+2}_X g = 4(c^2{-}b^2)\,\mathcal{L}^k_X g$ for all $k\ge 1$, so that
$$
\mathrm{tr}_g(\mathcal{L}^{2k+1}_X g) = 0
\quad\text{and}\quad
\mathrm{tr}_g(\mathcal{L}^{2k+2}_X g) = 8c^2 \bigl(4(c^2{-}b^2)\bigr)^k
$$
for all $k\ge0$.
Moreover, because the equations express all of the second derivatives of $u$ and $v$ algebraically in terms of $u$, $v$ and their first derivatives, it follows that any solution $(u,v)$ that is $C^2$ is real-analytic. Thus, any $C^2$ solution $(u,v)$ of the equations above defines a vector field $X$ whose flow has singular values that depend only on $t$.
Now, it is easy to show that the equation $({u_x}^2-{u_y}^2)\cos v +2u_xu_y\sin v = 0$ cannot hold on any open domain. Moreover, given constants $(r_1,r_2,s,s_1,s_2)$ that satisfy $({r_1}^2-{r_2}^2)\cos s +2r_1r_2\sin s \not= 0$, and fixing constants $b$ and $c\not=0$, there is at most one solution on a neighborhood of $(x,y) = (0,0)$ of the above equations such that
$$
\bigl(u_x(0,0),u_y(0,0),v(0,0),v_x(0,0),v_y(0,0)\bigr) = (r_1,r_2,s,s_1,s_2).
$$
In fact, one finds that, in the case $b=0$, there exists such a local solution, as the PDE system is involutive (i.e., Frobenius). Unfortunately, there is good evidence that the 'generic' solution with $b=0$ cannot be written in terms of elementary functions. Sometimes, though, one can express a solution implicitly: As an example of a particular solution with $b=0$, one can show that when $x+i\,y$ is parametrized by a holomorphic parameter $w$ such that
$$
dx + i\,dy = \mathrm{e}^{w^2}\,dw,
$$
then, setting
$$
u_x + i\,u_y = \overline{w}\,\mathrm{e}^{w^2},
$$
one finds that, for all $t$, the time $t$ flow of the unimodular vector field $X = u_y\,\partial_y - u_x\,\partial_x$ has constant singular values. (This is most easily checked in the $w$ coordinate, where $X = 2i\left(w\,\partial_{\bar w} - \bar w\,\partial_w\right)$ and $g = \mathrm{e}^{w^2+\bar w^2}\,dw\circ d\bar w$.)
When $b\not=0$, it turns out that there are more restrictions on the initial conditions. In fact, an analysis shows that, when $b\not=0$, either $v$ is constant, or else the equation $(u_xv_x+u_yv_y)=4(c^2{-}b^2)$ must hold.
When $v$ is constant, the vector field $X$ is an affine vector field
whose flow is a $1$-parameter family of (unimodular) affine transformations of the plane, which is clearly a solution.
When $v$ is not constant and $b\not=0$, it turns out that one can integrate the equations explicitly. One finds that, up to a rigid motion in the plane, these solutions can be written in polar coordinates $(r,\theta)$ (where $g = dr^2 + r^2\,d\theta^2$) in the form
$$
X = a_1\left(r\,\frac{\partial}{\partial r}-2\theta\,\frac{\partial}{\partial\theta}\right) + (a_2 - 2b\log r)\,\frac{\partial}{\partial\theta}
$$
where $a_1, a_2, b$ are constants and $c^2=a_1^2+b^2$. Note that the case $a_1=a_2=0$ is, up to a constant multiple, the example provide by the OP.
