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')

Do the artifacts also appear for a Cartesian (non-polar) plot? Contourf also interpolates. A [mcve] would help debug this. — Andras Deak -- Слава Україні, Apr 22 '21 at 08:14
No, the artifacts do not appear for the Cartesian plot. None of contourf's methods seems to do the trick. I think *did* make a minimial reproductible example! :) — ngc1535, Apr 23 '21 at 04:27
Very good MCVE! (There are still styling things you could remove but this is concise enough to be useful.) I'm not sure what's going on but others will be able to help you now. — Andras Deak -- Слава Україні, Apr 23 '21 at 13:33
Thank you Andras. To create the cartesian result- the ax.set polar stylings are turned off and cax = ax.contourf(xgrid, ygrid, nsb_grid, levels=levels, cmap=colors). This has NO artifacts regardless of contour levels. This is why I have focused on the polar projection as being an issue. — ngc1535, Apr 23 '21 at 19:14
Oh yeah, I get that. I meant that I have no idea where the problem is coming from. I tried removing the nans from the interpolated data (didn't seem to help), I tried interpolating the polar values directly (I couldn't get it right within the ten minutes I could spend on all this). — Andras Deak -- Слава Україні, Apr 23 '21 at 19:15
I have given up... all I want to figure out now is how to overlay an empty polar axes over the Cartesian contours... and get on with life! — ngc1535, Apr 25 '21 at 07:16
Note that it's possible that the rendering bug is a matplotlib bug (or shortcoming). You might also try opening an issue on their repo using the above example. It's not perfect because there's a scipy.interpolate step in generating the data, but it's clearly reproducible. — Andras Deak -- Слава Україні, Apr 26 '21 at 20:27

ngc1535 · Answer 1 · 2021-04-28T19:15:40.367

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)

I never meant for you to delete this question; leaving this answer is the ideal solution. Good job getting it figured out. I'd also post a link to the discussion on github, assuming that's where you got your answer. — Andras Deak -- Слава Україні, Apr 27 '21 at 18:36
Oh..I thought about that- but thought that "external" links were frowned on.I will include int he above... — ngc1535, Apr 28 '21 at 19:14
They are frowned upon when they are the meat of a post. "My problem is there", "the solution to your problem is there" etc. But in this case this just adds some context, and the github link is as stable as it gets :) — Andras Deak -- Слава Україні, Apr 28 '21 at 19:16

Griddata and Contourf produce artifacts with increasing steps/levels

1 Answers1