How do I create curved bins in matplotlib polar 2D histogram?

I am plotting a polar 2d histogram in Python 3.7 using matplotlib and the following code (adapted from this answer to another question):

import numpy as np
import matplotlib.pyplot as plt

# input data
azimut = np.random.rand(3000)*2*np.pi
radius = np.random.rayleigh(9, size=3000)

# binning
rbins = np.linspace(0, radius.max(), 10)
abins = np.linspace(0, 2*np.pi, 10)

# histogram
hist, _, _ = np.histogram2d(azimut, radius, bins=(abins, rbins))
A, R = np.meshgrid(abins, rbins)

# plot
fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))

pc = ax.pcolormesh(A, R, hist.T, cmap='inferno')
fig.colorbar(pc)

plt.show()

To produce the following plot:

Due to the larger bin sizes, the polar projection is appearing more like a polygon rather than a circle.

Is there any way to plot this so that the bins appear curved rather than straight? I.E. so that the plot is always circular, regardless of the bin size and doesn't become polygon-like when bins are larger?

A matplotlib solution would be preferable, but others are welcome.

Thanks very much for any help.

To get a rounded look, the mesh can be subdivided into more angles. Note that np.linspace(0, 2 * np.pi, 10) creates 9 bins (and 10 boundaries). For the subdivided mesh you need e.g. 90 bins, so 91 boundaries. The histogram values need to be repeated by the same factor.

The code below uses a different colormap for debugging purposes. An optional grid highlights the original boundaries.

import numpy as np
import matplotlib.pyplot as plt

# input data
azimut = np.random.rand(3000) * 2 * np.pi
radius = np.random.rayleigh(9, size=3000)

# binning
rbins = np.linspace(0, radius.max(), 7)
abins = np.linspace(0, 2 * np.pi, 10)
subdivs = 10
abins2 = np.linspace(0, 2 * np.pi, (len(abins) - 1) * subdivs + 1)

# histogram
hist, _, _ = np.histogram2d(azimut, radius, bins=(abins, rbins))
A1, R1 = np.meshgrid(abins, rbins)
A2, R2 = np.meshgrid(abins2, rbins)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10, 4), subplot_kw=dict(projection="polar"))

# plot with original mesh
pc1 = ax1.pcolormesh(A1, R1, hist.T, cmap='hsv')
ax1.tick_params(axis='y', labelcolor='white')
ax1.set_xticks(abins[:-1])
fig.colorbar(pc1, ax=ax1)

# plot with subdivided mesh
pc2 = ax2.pcolormesh(A2, R2, np.repeat(hist.T, subdivs, axis=1), cmap='hsv')
ax2.tick_params(axis='y', labelcolor='white')
ax2.set_xticks(abins[:-1])
ax2.set_yticks(rbins, minor=True)
ax2.grid(axis='x', color='white')
ax2.grid(axis='y', which='minor', color='white')
fig.colorbar(pc2, ax=ax2)

plt.tight_layout()
plt.show()