This is more of an explanation about why non-isotropic stationary processes aren't often (or ever) really dealt with explicitly than a positive answer.

For a finite Borel measure $F$ on $\mathbb{R}^m$, write (loosely)

$F(k,dk) := F((k_1,k_1+dk_1) \times \dots \times (k_m,k_m+dk_m))$.

Bochner’s theorem essentially states that $B$ is a wide-sense stationary (WSS) covariance function iff

$B(x) = \int_{\mathbb{R}^m} F_B(k,dk) \exp(2\pi i \langle k, x \rangle)$

for some finite Borel measure $F_B$ called the spectral measure.

In one dimension a mild (absolute continuity) restriction and the Radon-Nikodym theorem allow us to write

$f_B(k) := F_B((-\infty, k)) \Rightarrow B(x) =\int_\mathbb{R} F_B(k,dk) \exp(2 \pi i k x) = \int_\mathbb{R} dk f'_B(k) \exp(2 \pi i k x)$,

which implies that $\Rightarrow \hat B = f'_B$, where the Fourier transform is indicated.

In particular, given the (Fourier transform of the WSS) covariance, the spectral measure is also in hand. **In higher dimensions, however, the situation is more difficult**. Although we can still write $B(x)$ along the lines above, actual computations are hard. To drive this point home, notice that

$d\nu_g(k) = dx \lvert \det \nabla_k g \rvert \overset{g \in C^1(\mathbb{R}^m,\mathbb{R}^m)}{\Leftrightarrow} \nu_g = \lambda \circ g$,

(where $\lambda$ is Lebesgue measure) does not provide a solution for the inverse problem of obtaining the diffeomorphism $g$ from $\nu_g$.

Isotropy makes things much easier, as you've noticed.