Griddata and Contourf produce artifacts with increasing steps/levels

Question

I am using SciPy Griddata to interpolate data in its Cartesian form and then plot these data using contourf with a polar projection. When the Cartesian interpolated data is plotted with contourf there are no artifacts. However, when the projection is polar, artifacts develop with increasing "levels".

The artifacts are polygons or rays that form near regions of steep gradients. The code below plots the brightness of the sky with the moon. With graphlevels of "12" there isn't an issue. Artifacts develop with graphlevel of "25." My desired level is 80 or more - which shows terrible artifacts. The below is example real data from one night. These artifacts always occur. See images with Levels = 12 and Levels = 80

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata


gridsize =150
graphlevels =12

plt.figure(figsize=(12,10))
ax = plt.subplot(111,projection='polar')


x = [72.90,68.00,59.14,44.38,29.63,63.94,59.68,51.92,38.98,26.03,47.34,44.20,38.46,28.89,19.31,23.40,20.40,15.34,10.28,-0.18,-0.14,-0.09,-0.04,0.02,-25.39,-23.66,-20.57,-15.40,-10.23,-47.56,-44.34,-38.54,-28.89,-19.22,-64.01,-59.68,-51.89,-38.90,-25.90,-72.77,-67.84,-58.98,-44.21,-29.44,-72.75,-67.83,-58.96,-44.18,-29.41,-59.63,-51.82,-38.83,-25.84,-47.42,-44.20,-38.40,-28.76,-19.12,-23.40,-20.32,-15.19,-10.08,0.27,0.25,0.23,0.20,23.92,20.80,15.63,10.46,47.93,44.67,38.86,29.17,19.48,64.40,60.03,52.20,39.18,26.15,73.08,68.12,59.26,44.47,29.68,-4.81]
y = [12.93,12.01,10.38,7.67,4.99,37.03,34.49,29.93,22.33,14.77,56.60,52.75,45.82,34.26,22.72,64.60,56.14,42.02,27.90,73.66,68.67,59.68,44.68,29.68,69.12,64.45,56.00,41.92,27.84,56.26,52.45,45.56,34.08,22.61,36.59,34.11,29.61,22.11,14.62,12.48,11.62,10.04,7.43,4.83,-13.33,-12.31,-10.78,-8.21,-5.58,-34.84,-30.36,-22.87,-15.36,-57.04,-53.20,-46.31,-34.83,-23.34,-65.20,-56.72,-42.62,-28.53,-69.33,-60.31,-45.31,-30.31,-65.09,-56.63,-42.55,-28.47,-56.81,-52.99,-46.13,-34.69,-23.23,-36.99,-34.53,-30.08,-22.66,-15.22,-12.73,-11.93,-10.44,-7.94,-5.40,-1.22,]
skybrightness = [19.26,19.31,19.21,19.65,19.40,19.26,19.23,19.43,19.57,19.52,19.19,19.31,19.33,19.68,19.50,19.29,19.45,19.50,19.23,18.98,19.28,19.46,19.54,19.22,19.03,19.18,19.35,19.37,19.08,18.99,18.98,19.26,19.36,19.08,18.79,18.85,19.13,19.17,19.05,18.51,18.64,18.88,18.92,18.93,18.12,18.34,18.72,18.82,18.74,18.22,18.46,18.76,18.26,18.13,18.24,18.46,18.58,17.30,18.38,18.08,18.24,17.68,18.34,18.46,18.65,18.23,18.70,18.52,18.79,18.83,18.18,18.51,19.01,19.08,19.08,18.99,19.02,19.07,19.20,19.27,19.06,19.01,19.28,19.46,19.30,18.94]

xgrid = np.linspace(min(x), max(x),gridsize)
ygrid = np.linspace(min(y), max(y),gridsize)

xgrid, ygrid = np.meshgrid(xgrid, ygrid, indexing='ij')

nsb_grid = griddata((x,y),skybrightness,(xgrid, ygrid), method='linear')

r = np.sqrt(xgrid**2 + ygrid**2)
theta = np.arctan2(ygrid, xgrid)

plt.rc('ytick', labelsize=16)
ax.set_facecolor('#eeddcc')

colors = plt.cm.get_cmap('RdYlBu')
levels,steps = np.linspace(min(skybrightness), max(skybrightness)+0.3,graphlevels, retstep=True)
ticks = np.linspace(min(skybrightness), max(skybrightness)+0.3,12)

cax = ax.contourf(theta, r, nsb_grid, levels=levels, cmap=colors)

cbar = plt.colorbar(cax, fraction=0.046, pad=0.04, ticks=ticks)
cbar.set_label(r'mag/arcsec$^2$')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.set_rmax(75)
ax.set_yticks(range(10, 80, 20))
ax.set_xticklabels([r'N', r'NE', r'E', r'SE', r'S', r'SW', r'W', r'NW'])
ax.grid(alpha=0.3)
plt.savefig('StackOverflowHELP.png')

Answer 1

I am going to leave my question and this answer on StackOverflow... because I did get an answer from the developers of Matploblib. The problem is Contourf . In its attempt to project data in polar dimensions there are overlaps and extensions of polygons at the cyclic boundaries that cause problems. The only way to avoid this is to add points at the boundary. To quote the developer:

The workaround is a lot of effort and has to be tuned to each particular problem, so is a very long way from being ideal. We (Matplotlib) should do better in these situations. Inserting extra points into the triangulation isn't the right approach, we should instead correct the lines/polygons that traverse the discontinuity to provide a general solution.

See https://github.com/matplotlib/matplotlib/issues/20060 for the full discussion

The answer I settled on is to interpolate and render the result in Cartesian space. Then I format an empty polar plot with axes and labels to overlay on the top... and get on with my life!

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata


gridsize =150
graphlevels = 200

fig = plt.figure(figsize=(12,10))

ax = fig.add_subplot(111, aspect='equal')
pax = fig.add_subplot(111,projection='polar')
pax.set_facecolor('none')
ax.set_axis_off()
ax.set_xlim([-75,75])
ax.set_ylim([-75,75])

x = [72.90,68.00,59.14,44.38,29.63,63.94,59.68,51.92,38.98,26.03,47.34,44.20,38.46,28.89,19.31,23.40,20.40,15.34,10.28,-0.18,-0.14,-0.09,-0.04,0.02,-25.39,-23.66,-20.57,-15.40,-10.23,-47.56,-44.34,-38.54,-28.89,-19.22,-64.01,-59.68,-51.89,-38.90,-25.90,-72.77,-67.84,-58.98,-44.21,-29.44,-72.75,-67.83,-58.96,-44.18,-29.41,-59.63,-51.82,-38.83,-25.84,-47.42,-44.20,-38.40,-28.76,-19.12,-23.40,-20.32,-15.19,-10.08,0.27,0.25,0.23,0.20,23.92,20.80,15.63,10.46,47.93,44.67,38.86,29.17,19.48,64.40,60.03,52.20,39.18,26.15,73.08,68.12,59.26,44.47,29.68,-4.81]
y = [12.93,12.01,10.38,7.67,4.99,37.03,34.49,29.93,22.33,14.77,56.60,52.75,45.82,34.26,22.72,64.60,56.14,42.02,27.90,73.66,68.67,59.68,44.68,29.68,69.12,64.45,56.00,41.92,27.84,56.26,52.45,45.56,34.08,22.61,36.59,34.11,29.61,22.11,14.62,12.48,11.62,10.04,7.43,4.83,-13.33,-12.31,-10.78,-8.21,-5.58,-34.84,-30.36,-22.87,-15.36,-57.04,-53.20,-46.31,-34.83,-23.34,-65.20,-56.72,-42.62,-28.53,-69.33,-60.31,-45.31,-30.31,-65.09,-56.63,-42.55,-28.47,-56.81,-52.99,-46.13,-34.69,-23.23,-36.99,-34.53,-30.08,-22.66,-15.22,-12.73,-11.93,-10.44,-7.94,-5.40,-1.22,]
skybrightness = [19.26,19.31,19.21,19.65,19.40,19.26,19.23,19.43,19.57,19.52,19.19,19.31,19.33,19.68,19.50,19.29,19.45,19.50,19.23,18.98,19.28,19.46,19.54,19.22,19.03,19.18,19.35,19.37,19.08,18.99,18.98,19.26,19.36,19.08,18.79,18.85,19.13,19.17,19.05,18.51,18.64,18.88,18.92,18.93,18.12,18.34,18.72,18.82,18.74,18.22,18.46,18.76,18.26,18.13,18.24,18.46,18.58,17.30,18.38,18.08,18.24,17.68,18.34,18.46,18.65,18.23,18.70,18.52,18.79,18.83,18.18,18.51,19.01,19.08,19.08,18.99,19.02,19.07,19.20,19.27,19.06,19.01,19.28,19.46,19.30,18.94]

xgrid = np.linspace(min(x), max(x),gridsize)
ygrid = np.linspace(min(y), max(y),gridsize)

xgrid, ygrid = np.meshgrid(xgrid, ygrid, indexing='ij')

nsb_grid = griddata((x,y),skybrightness,(xgrid, ygrid), method='linear')

plt.rc('ytick', labelsize=16) #colorbar font

colors = plt.cm.get_cmap('RdYlBu')
levels,steps = np.linspace(min(skybrightness), max(skybrightness)+0.3,graphlevels, retstep=True)
ticks = np.linspace(min(skybrightness), max(skybrightness)+0.3,12)

cax = ax.contourf(xgrid, ygrid, nsb_grid, levels=levels, cmap=colors)

cbar = plt.colorbar(cax, fraction=0.046, pad=0.04, ticks=ticks)
cbar.set_label(r'mag/arcsec$^2$')
pax.set_theta_zero_location('N')
pax.set_theta_direction(-1)
pax.set_rmax(75)
pax.set_yticks(range(10, 80, 20))
pax.set_xticklabels([r'N', r'NE', r'E', r'SE', r'S', r'SW', r'W', r'NW'])
pax.grid(alpha=0.3)