Two-dimensional interpolation/smoothing of unevenly sampled angular values

I have some data that consist of unevenly sampled 2D spatial locations, where each x, y coordinate has an associated phase value theta between 0 and 2pi. I'd like to be able to interpolate the theta values onto a regular x, y grid. The data is degenerate in the sense that the same (or very nearby) x, y locations may be associated with multiple phase values, and vice versa for values of theta, so this is strictly speaking a smoothing problem rather than straight interpolation.

I've briefly experimented with scipy's radial basis functions, but these give nasty edge effects because of the discontinuity in the theta values from 2pi --> 0.

Here's a toy example (the real spatial distribution of phases is a lot messier):

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import colorbar
from matplotlib.colors import Normalize
from scipy import interpolate

# randomly sampled spatial locations
x, y = np.random.uniform(-1, 1, size=(2, 1000))
# theta varies smoothly with location apart from the singularity at 0, 0
z = np.arctan2(x, y) % (2 * np.pi)

# smooth with a simple linear RBF
rbf = interpolate.Rbf(x, y, z, function='linear', smooth=0.1)

# resample on a finer grid
xi, yi = np.mgrid[-1:1:100j, -1:1:100j].reshape(2, -1)
zi = rbf(xi, yi) % (2 * np.pi)

# plotting
fig, ax = plt.subplots(1, 1, subplot_kw={'aspect': 'equal'})
ax.hold(True)
norm = Normalize(0, 2 * np.pi)
im = ax.imshow(zi.reshape(100, 100).T, extent=(-1, 1, -1, 1),
               origin='lower', cmap=plt.cm.hsv, norm=norm)
sc = ax.scatter(x, y, s=30, c=z, cmap=im.cmap, norm=norm)
cax, kw = colorbar.make_axes_gridspec(ax)
cb = plt.colorbar(im, cax=cax, **kw)
ax.set_xlabel(r'$X_0$', fontsize='x-large')
ax.set_ylabel(r'$Y_0$', fontsize='x-large')
cb.set_ticks(np.arange(0, 2.1*np.pi, np.pi/2.))
cb.set_ticklabels([r'$0$', r'$\frac{\pi}{2}$', r'$\pi$',
                   r'$\frac{3\pi}{2}$', r'$2\pi$'])
cb.set_label(r'$\theta$', fontsize='x-large')
cb.ax.tick_params(labelsize='x-large')
plt.show()

enter image description here

What would be a good way to go about interpolating angular quantities like this? Does scipy have any built in interpolation method that will deal with angles nicely, or will I have to write my own?

Solution

I feel pretty stupid now!

The answer was very simple - this answer on MathOverflow clued me in. There is no problem with discontinuity provided that I convert from a polar coordinate space to a Cartesian one, then interpolate the x and y components of the vector independently:

x, y = np.random.uniform(-1, 1, size=(2, 1000))
z = np.arctan2(y, x) % (2*np.pi)

# convert from polar --> cartesian
u, v = np.cos(z), np.sin(z)

# interpolate x and y components separately
rbf_u = interpolate.Rbf(x, y, u, function='linear', smooth=0.1)
rbf_v = interpolate.Rbf(x, y, v, function='linear', smooth=0.1)
xi, yi = np.mgrid[-1:1:100j, -1:1:100j].reshape(2, -1)
ui = rbf_u(xi, yi)
vi = rbf_v(xi, yi)

# convert from cartesian --> polar
zi = np.arctan2(ui, vi) % (2*np.pi)

enter image description here

It would be nice performance-wise if there was a way to avoid performing two separate interpolations on the x and y components, but I don't really see a way around this.