GaussianProcessRegressor with geometric anisotropy

Is it possible to define a kernel for the GaussianProcessRegressor which uses "geometric anisotropy"?

I know it's possible to define anisotropy with some of the existing kernels, but it seems to only allow anisotropy parallel to the dimensions. Consider the case of related dimensions (for example a spatial x/y coordinate), I want the anisotropy to be able to take into account an angle by which to rotate the "ellipse" defined by the two length scales (those would become the major and minor axis).

It's not obvious to me whether this is already possible with the current kernels in Scikit-Learn, or if I should create my own kernel by subclassing one of the existing.

Consider the following example:

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, ConstantKernel

# training data
ygrid, xgrid = np.mgrid[-2:3:1, -2:3:1]
X_train = np.stack([ygrid.flat, xgrid.flat], axis=1)
y_train = np.eye(*ygrid.shape).ravel()

# define the kernel
c = ConstantKernel(constant_value=y_train.mean(), constant_value_bounds='fixed')
m = Matern(length_scale=[0.5, 0.5], nu=0.5)

kernel = c * m

# fit
gp = GaussianProcessRegressor(kernel=kernel)
gp.fit(X_train, y_train)

After fitting the kernel parameters are:

print(gp.kernel_.get_params())
{'k1': 0.447**2,
 'k2': Matern(length_scale=[0.444, 0.444], nu=0.5),
 'k1__constant_value': 0.2,
 'k1__constant_value_bounds': 'fixed',
 'k2__length_scale': array([0.4443361, 0.4443361]),
 'k2__length_scale_bounds': (1e-05, 100000.0),
 'k2__nu': 0.5}

A detailed prediction:

ygrid, xgrid = np.mgrid[-3:4:0.05, -3:4:0.05]
X_predict = np.stack([ygrid.flat, xgrid.flat], axis=1)

y_predict = gp.predict(X_predict)
y_predict = y_predict.reshape(xgrid.shape)

The left axes shows the training data (25 samples), and the right one the prediction on a finer grid:

It's already obvious from the fit that the anisotropy doesn't add anything in this case since the data correlates along the diagonal and the anisotropy from the Matern kernel doesn't capture that. The image shows that the prediction could be a lot more smooth along the diagonal if I'm able to specify this somewhere.

For my use-case, the angle can be anything (0-360 degrees), and I don't know it beforehand, so I can't transform the input coordinates.

Is there are workaround for this with the existing kernels? Or how can I add this behavior by creating a new kernel (preferably Matern-based)?

from GaussianProcessRegressor with geometric anisotropy