Workaround for spherical polar coordinate Jacobian division

I am having data visualization issue. Short summup, I'm working on a project involving a polar spherical coordinate mesh, and trying to solved coupled system of ODE (chemical reactions) for each cell. For a specific reason I need my state vector to be of the form (r^2*sin(theta)*n_i), i={1,2,3...}.

I rewrote a quick example of my issue below, you can run it as is. Why is it that cst2, which I would have assumed equal to a np.ones(a[0].shape), doesn't show a uniform pcolormesh. And even more surprising to me, why does adding the colorbar make this issue vanish ?

My best guess is that dividing by r^2*sin(theta) causes numerical issues, but how can I workaround this ? (I need to see my data without the curvature term in order to interpret it -> the jacobian division seems mandatory to me at first thought).

import numpy as np
import matplotlib.pylab as plt
fig, ax = plt.subplots()

### Edges
r = np.logspace(np.log10(1), np.log10(4.6), num=14) #cell edges
theta = np.linspace(0+0.001,np.pi-0.001,num=10) 
b = np.meshgrid(r,theta)

### Center
r_c = r[0:-1] + np.ediff1d(r)/2 #get the cell center
theta_c = theta[0:-1] + np.ediff1d(theta)/2
a = np.meshgrid(r_c,theta_c)

### The jacobian division
cst = pow(a[0],2)*np.sin(a[1])
cst2 = np.copy(cst)/pow(a[0],2)/np.sin(a[1])
pcm = ax.pcolormesh(b[0]*np.cos(b[1]),\
                    b[0]*np.sin(b[1]), \
                    cst2,cmap='seismic',edgecolor='black')
# clb = fig.colorbar(pcm, ax=ax, orientation='horizontal') 

Tweaking the formula slightly seems to work. Also used a lighter color for cmap.

cst2 = np.copy(cst) / (pow(a[0],2)*np.sin(a[1]))
pcm = ax.pcolormesh(b[0]*np.cos(b[1]), \
                    b[0]*np.sin(b[1]), \
                    cst2,cmap='coolwarm',edgecolor='black')