Eulerian grid hydro sims output their fluid variables using tables of grid cells. In the case of 2D LDI simulations, these grids are in terms of radius and **angle**. This presents an annoying problem for plotting, since typical plotting libraries work in Cartesian.

I spent an hour coding up a hackish and slow projection function, but then it occurred to me that Python must do its own projections when generating display images. Our use case is weird though, because we don’t want a full polar plot — just a range of radii and angles. Matplotlib can’t do this by default, but this stackoverflow post gives a handy function using axisartist. Put this at the top of your code:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from matplotlib.transforms import Affine2D
import mpl_toolkits.axisartist.floating_axes as floating_axes
import mpl_toolkits.axisartist.angle_helper as angle_helper
from matplotlib.projections import PolarAxes
from mpl_toolkits.axisartist.grid_finder import MaxNLocator
# define how your plots look:
def setup_axes(fig, rect, theta, radius):
# PolarAxes.PolarTransform takes radian. However, we want our coordinate
# system in degree
tr = Affine2D() + PolarAxes.PolarTransform()
# Find grid values appropriate for the coordinate (degree).
# The argument is an approximate number of grids.
grid_locator1 = angle_helper.LocatorD(2)
# And also use an appropriate formatter:
tick_formatter1 = angle_helper.FormatterDMS()
# set up number of ticks for the r-axis
grid_locator2 = MaxNLocator(4)
# the extremes are passed to the function
grid_helper = floating_axes.GridHelperCurveLinear(tr,
extremes=(theta[0], theta[1], radius[0], radius[1]),
grid_locator1=grid_locator1,
grid_locator2=grid_locator2,
tick_formatter1=tick_formatter1,
tick_formatter2=None,
)
ax1 = floating_axes.FloatingSubplot(fig, rect, grid_helper=grid_helper)
fig.add_subplot(ax1)
# adjust axis
# the axis artist lets you call axis with
# "bottom", "top", "left", "right"
ax1.axis["left"].set_axis_direction("bottom")
ax1.axis["right"].set_axis_direction("top")
ax1.axis["bottom"].set_visible(False)
ax1.axis["top"].set_axis_direction("bottom")
ax1.axis["top"].toggle(ticklabels=True, label=True)
ax1.axis["top"].major_ticklabels.set_axis_direction("top")
ax1.axis["top"].label.set_axis_direction("top")
ax1.axis["left"].label.set_text("R")
ax1.axis["top"].label.set_text(ur"$\theta$")
# create a parasite axes
aux_ax = ax1.get_aux_axes(tr)
aux_ax.patch = ax1.patch # for aux_ax to have a clip path as in ax
ax1.patch.zorder=0.9 # but this has a side effect that the patch is
# drawn twice, and possibly over some other
# artists. So, we decrease the zorder a bit to
# prevent this.
return ax1, aux_ax

Here’s some example code (append this to the above to test it) which generates an example density function and then distorts it from -0.1 to 0.1 radians, and from 1 to 10 in radius.

# input file grid dimensions
xmax = 2000
ymax = 50
thetamin = -0.1
thetamax = 0.1
rmin = 1.0
rmax = 10.0
# generate a fake density variable on a 2000x50 grid
img = np.zeros( (xmax,ymax) )
for x in range(0,xmax):
for y in range(0,ymax):
img[x,y] = np.sin(0.05*x) * np.sin(0.05*y) + 10
theta,rad = np.meshgrid(np.linspace(thetamin,thetamax,ymax),
np.linspace(rmin,rmax,xmax))
X = theta
Y = rad
# let's show the starting materials
plt.imshow(img.T)
plt.show()
plt.clf()
# now let's plot the finished product!
fig = plt.figure(1, figsize=(8, 4))
fig.subplots_adjust(wspace=0.3, left=0.05, right=0.95)
ax1, aux_ax1 = setup_axes(fig, 111, theta=[thetamin,thetamax], radius=[rmin, rmax])
aux_ax1.pcolormesh( X,Y, img )
plt.show()

You first get a Cartesian plot.

After our changes, we should get the desired result.