How to make four-way logarithmic plot in Matplotlib?

Question

Four-way logarithmic plot is a very often used graph for vibration control and earthquake protection. I am quite interesting in how this plot can be plotted in Matplotlib instead of adding axes in Inkscape. A sample of Four-way logarithmic plot is here.

A quick and dirty Python code can generate main part of the figure, but I cannot add the two axes onto the figure. http://matplotlib.org/examples/axes_grid/demo_curvelinear_grid.html provides an example of adding axes, but I fails to make it working. Anyone has similar experience on adding axes to Matplotlib figure?

from pylab import *
from  mpl_toolkits.axisartist.grid_helper_curvelinear import GridHelperCurveLinear
from mpl_toolkits.axisartist import Subplot
beta=logspace(-1,1,500)
Rd={}
for zeta in [0.01,0.1,0.2,0.7,1]:
    Rd[zeta]=beta/sqrt((1-beta*beta)**2+(2*beta*zeta)**2)
    loglog(beta,Rd[zeta])
ylim([0.1,10])
xlim([0.1,10])
grid('on',which='minor')

Update: Thank you all! I use Inkscape to modify the figure above. I think the result is just fine. However, I am still looking for methods to draw this figure in Matplotlib.

Answer 1

Here is a partial solution. I am still working on how to do all of this in a natural loglog() plot rather than scaling the data. (To complete this example you would have to define custom tick-lables so that they display 10**x rather than x.)

%matplotlib inline                   # I am doing this in an IPython notebook.
from matplotlib import pyplot as plt
import numpy as np
from numpy import log10

# Generate the data
beta = np.logspace(-1, 1, 500)[:, None]
zeta = np.array([0.01,0.1,0.2,0.7,1])[None, :]
Rd = beta/np.sqrt((1 - beta*beta)**2 + (2*beta*zeta)**2)

def draw(beta=beta, Rd=Rd):
    plt.plot(log10(beta), log10(Rd))
    plt.ylim([log10(0.1), log10(10)])
    plt.xlim([log10(0.1), log10(10)])
    plt.grid('on',which='minor')
    ax = plt.gca()
    ax.set_aspect(1)

from mpl_toolkits.axisartist import GridHelperCurveLinear
from matplotlib.transforms import Affine2D
from mpl_toolkits.axisartist import SubplotHost
from mpl_toolkits.axisartist import Subplot

#tr = Affine2D().rotate(-np.pi/2)
#inv_tr = Affine2D().rotate(np.pi/2)

class Transform(object):
    """Provides transforms to go to and from rotated grid.

    Parameters
    ----------
    ilim : (xmin, xmax, ymin, ymax)
       The limits of the displayed axes (in physical units)
    olim : (xmin, xmax, ymin, ymax)
       The limits of the rotated axes (in physical units)
    """
    def __init__(self, ilim, olim):
        # Convert each to a 3x3 matrix and compute the transform
        # [x1, y1, 1] = A*[x0, y0, 1]
        x0, x1, y0, y1 = np.log10(ilim)
        I = np.array([[x0, x0, x1],
                      [y0, y1, y1],
                      [ 1,  1,  1]])

        x0, x1, y0, y1 = np.log10(olim)
        x_mid = (x0 + x1)/2
        y_mid = (y0 + y1)/2
        O = np.array([[   x0, x_mid, x1],
                      [y_mid,    y1, y_mid],
                      [    1,     1,     1]])
        self.A = np.dot(O, np.linalg.inv(I))
        self.Ainv = np.linalg.inv(self.A)

    def tr(self, x, y):
        """From "curved" (rotated) coords to rectlinear coords"""
        x, y = map(np.asarray, (x, y))
        return np.dot(self.A, np.asarray([x, y, 1]))[:2]

    def inv_tr(self, x, y):
        """From rectlinear coords to "curved" (rotated) coords"""
        x, y = map(np.asarray, (x, y))
        return np.dot(self.Ainv, np.asarray([x, y, 1]))[:2]

ilim = (0.1, 10)
olim = (0.01, 100)
tr = Transform(ilim + ilim, olim + olim)

grid_helper = GridHelperCurveLinear((tr.tr, tr.inv_tr))

fig = plt.gcf()
ax0 = Subplot(fig, 1, 1, 1)
ax1 = Subplot(fig, 1, 1, 1, grid_helper=grid_helper, frameon=False)
ax1.set_xlim(*np.log10(olim))
ax1.set_ylim(*np.log10(olim))
ax1.axis["left"] = ax1.new_floating_axis(0, 0.)
ax1.axis["bottom"] = ax1.new_floating_axis(1, 0.0)
fig.add_subplot(ax0)
fig.add_subplot(ax1)
ax0.grid('on', which='both')
ax1.grid('on', which='both')

plt.plot(log10(beta), log10(Rd))
plt.ylim(np.log10(ilim))
plt.xlim(np.log10(ilim))

Answer 2

This seems to be a bit tricker than it should. There are ways to center the spines (axis lines), and ways to rotate them, but those do not work together. Adding a normal axis on a line (a la mpl demos) results in a curved axis (because it is logarithmic). Here is a [poor] example of how to draw -- as in, like you would with Inkscape something to look like an additional pair of axis spines with the example data.

import matplotlib.pyplot as plt
import numpy as np

#data
b = np.logspace(-1, 1, 500)
Rd = {}
for zeta in [0.01, 0.1, 0.2, 0.7, 1]:
    Rd[zeta] = b / np.sqrt((1 - b * b) ** 2 + (2 * b * zeta) ** 2)

#plot
fig = plt.figure()   
ax1 = fig.add_subplot(111)

for z in Rd:  
    ax1.loglog(b, Rd[z])

ax1.set_xlim([0.1, 10])
ax1.set_ylim([0.1, 10])
ax1.set_aspect(1.)

#draw lines to look like diagonal spines (axes)
xmin, xmax = ax1.get_xlim()  # xlim == ylim

a = np.log10(xmin)
b = np.log10(xmax)
span = b - a
period_points = 3  # number of points/ticks per decade
npts = (span * period_points) + 1  # +1 for even powers of 10
x1 = np.logspace(a, b, num=npts)
x2 = np.logspace(b, a, num=npts)

ax1.plot(x1, x1, color='k', marker='x', ms='9')
ax1.plot(x1, x2, color='k', marker='x', ms='9')
#NOTE: v1.2.1 lacks 'TICKUP' and similar - these may be
#   a better choice in v1.3x and beyond

ax1.text(0.97, 0.9,
         "axis label: A",
         size='large',
         horizontalalignment='right',
         verticalalignment='top',
         rotation=45,
         transform=ax1.transAxes,
         #bbox={'facecolor': 'white', 'alpha': 0.5, 'pad': 10},
         )

ax1.text(0.03, 0.9,
         "axis label: B",
         size='large',
         horizontalalignment='left',
         verticalalignment='top',
         rotation=-45,
         transform=ax1.transAxes,
         #bbox={'facecolor': 'white', 'alpha': 0.5, 'pad': 10},
         )

plt.savefig("example.pdf")