[FIXED] Ellipse and Polygon Patches in log-log scale

Issue

I want to do some routines to automate the generation of Ashby charts. This type of plots are, essentially, scatter plots. It is useful for material selection, and it is common to add some envelope shape to groups of materials. It is also common to plot the axis in logarithmic scale due to the changes in several orders of magnitude.

To generate the envelope I have done two different things:

• Find the convex hull of the group of points; and
• Find the ellipsoid that represent 2 standar deviations in each dataset.

When I use these approaches in linear-scaled plots everything is fine, they also work sometimes in the log-log-scale. Here there are two examples:

In other cases, my approaches don't work, i.e., the shaded regions are calculated fine but their plots in log-log-scale is wrong (The plots in linear-scale is always fine). Here some examples of what happens:

Questions:

• Is there a problem with the plot of Paths/Patches in matplotlib when using logarithmic axis?

• Is there a way to do the intended plots?

Edit

I added a working code as suggested by @farenorth. I noticed that negative values are problematic. In the case of polygons (convex-hull) this is not a problem, since the material properties tend to be positive, but there are some exceptions. For this cases, I think that the use of symlog scaling is enough.

In the case of ellipses, I need to think more about my implementation. It is based in the use of ellipses centered in the arithmetic mean with 2*inc standard deviations as semi-axes. This can lead to negative values for some region (since data is not necessarily symmetric around the mean). The data files can be downloaded here.

import numpy as np
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
from matplotlib import rcParams

rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16


def poly_enclose(points, color, inc=1.2, rad=0.3, lw=2):
    """
    Plot the convex hull around a set of points as a 
    shaded polygon.
    """
    hull = ConvexHull(points)


    cent = np.mean(points, 0)
    pts = []
    for pt in points[hull.simplices]:
        pts.append(pt[0].tolist())
        pts.append(pt[1].tolist())

    pts.sort(key=lambda p: np.arctan2(p[1] - cent[1],
                                    p[0] - cent[0]))
    pts = pts[0::2]  # Deleting duplicates
    pts.insert(len(pts), pts[0])


    verts = inc*(np.array(pts)- cent) + cent
    verts2 = np.zeros((3*verts.shape[0]-2,2))
    verts2[0::3] = verts
    verts2[1::3,:] = (1-rad)*verts[0:-1,:] + rad*verts[1:,:]
    verts2[2::3,:] = rad*verts[0:-1,:] + (1-rad)*verts[1:,:]
    verts2[0:-1] = verts2[1:]
    verts2[-1] = verts2[0]

    codes = [Path.MOVETO, Path.LINETO, Path.CURVE3,]
    for j in range(len(pts)-2):
        codes.extend([Path.CURVE3, Path.LINETO, Path.CURVE3,])
    codes.append(Path.CURVE3)


    path = Path(verts2, codes)
    patch = patches.PathPatch(path, facecolor=color, lw=0, alpha=0.2)
    edge = patches.PathPatch(path, edgecolor=color, facecolor='none', lw=lw)
    plt.gca().add_patch(patch)
    plt.gca().add_patch(edge)


def ellip_enclose(points, color, inc=1.2, lw=2, nst=2):
    """
    Plot the minimum ellipse around a set of points.

    Based on: 
    https://github.com/joferkington/oost_paper_code/blob/master/error_ellipse.py
    """

    def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:,order]

    x = points[:,0]
    y = points[:,1]
    cov = np.cov(x, y)
    vals, vecs = eigsorted(cov)
    theta = np.degrees(np.arctan2(*vecs[:,0][::-1]))
    w, h = 2 * nst * np.sqrt(vals)        
    center = np.mean(points, 0)
    ell = patches.Ellipse(center, width=inc*w, height=inc*h, angle=theta,
                          facecolor=color, alpha=0.2, lw=0)
    edge = patches.Ellipse(center, width=inc*w, height=inc*h, angle=theta,
                          facecolor='none', edgecolor=color, lw=lw)
    plt.gca().add_artist(ell)
    plt.gca().add_artist(edge)


inc = 1.2
rad = 0.3
lw = 2
colors = ['blue', 'green', 'red', 'orange']

## Example 1
plt.figure()
for k in range(4):
    points = 1.5*(np.random.rand(20, 2) - 0.5) + k + 3
    plt.plot(points[:,0], points[:,1], 'o', ms=8, color=colors[k],
             mfc="white", mec=colors[k])
    poly_enclose(points, colors[k], inc=inc, rad=rad, lw=lw)
##    ellip_enclose(points, colors[k], inc=inc, lw=lw)

#plt.xscale('symlog')
#plt.yscale('symlog')
plt.grid(True)

##  Example 2
E = {}
E["poly"] = np.loadtxt('young_poly.txt')
E["metals"] = np.loadtxt('young_metals.txt')
E["comp"] = np.loadtxt('young_comp.txt')
E["ceramic"] = np.loadtxt('young_ceramic.txt')

rho = {}
rho["poly"] = np.loadtxt('dens_poly.txt')
rho["metals"] = np.loadtxt('dens_metals.txt')
rho["comp"] = np.loadtxt('dens_comp.txt')
rho["ceramic"] = np.loadtxt('dens_ceramic.txt')

plt.figure()
for k, key  in enumerate(E.keys()):
    x = rho[key][:,0] * 1000
    y = E[key][:,0] * 1e9
    points = np.vstack([x,y]).T
    poly_enclose(points, colors[k], inc=inc, rad=0.3, lw=lw)
##    ellip_enclose(points, colors[k], inc=1, lw=lw)
    plt.plot(x, y, 'o', ms=8, color=colors[k], mfc="white", mec=colors[k])

##plt.xscale('symlog')
##plt.yscale('symlog')
plt.grid(True)

plt.show()

Solution

There are two places in your code where you are assuming a linear plot: calculating the center of the point cloud as a linear average (numpy.mean), and the constant coefficient (inc) multiplied to the difference vector between each point and the calculated center.

If you want visual proof, try the simple exercise of plotting the mean center of a point cloud on a logarithmic plot, it should appear well off-center, too far towards the top and the right. Similarly, if you multiply your point cloud by a constant coefficient and graph both versions, you will not see the expected linear dilation of a "bigger" cloud, but rather actually something that appears more like the points going down a drain near the top right corner of the original cloud.

If you convert your points to a logarithmic form before taking the mean and scaling the difference vectors and then change the points back to linear form before plotting, then you would have equivalent results to what you see on a linear plot. This could be accomplished either by adding an extra flag to your functions or by having your original functions return the shape coordinates rather than plotting them directly. The latter case seems cleaner to me.

For example:

def poly_enclose(...):
    ....
    return patch, edge

log_points = np.log(points)
log_patch, log_edge = poly_enclose(log_points, colors[k], inc=inc, rad=0.3, lw=lw)
lin_patch = np.exp(log_patch)
lin_edge = np.exp(log_edge)
plt.gca().add_patch(lin_patch)
plt.gca().add_patch(lin_edge)

Answered By - cheepychappy