I have a bunch of cross plots with two sets of data and have been looking for a matploltib way of highlighting their plotted regions with smoothed polygon outlines.
At the moment i just use Adobe Illustrator and amend saved plot, but this is not ideal. Example:
I'd be grateful for any pointers/links to examples.
Cheers
Here, you have an example. I was written the main ideas, but obviously, you could do it better.
A short explanations:
1) You need to compute the convex-hull (http://en.wikipedia.org/wiki/Convex_hull)
2) With the hull, you could scale it to keep all your data inside.
3) You must to interpolate the resulting curve.
The first part was done in http://wiki.scipy.org/Cookbook/Finding_Convex_Hull. The second one is trivial. The third one is very general, and you could perform any method, there are a lot of different ways to do the same. I took the @Jaime's approach (Smooth spline representation of an arbitrary contour, f(length) --> x,y), which I think it's a very good method.
I hope it help you...
#Taken from http://wiki.scipy.org/Cookbook/Finding_Convex_Hull
import numpy as n, pylab as p, time
def _angle_to_point(point, centre):
'''calculate angle in 2-D between points and x axis'''
delta = point - centre
res = n.arctan(delta[1] / delta[0])
if delta[0] < 0:
res += n.pi
return res
def _draw_triangle(p1, p2, p3, **kwargs):
tmp = n.vstack((p1,p2,p3))
x,y = [x[0] for x in zip(tmp.transpose())]
p.fill(x,y, **kwargs)
def area_of_triangle(p1, p2, p3):
'''calculate area of any triangle given co-ordinates of the corners'''
return n.linalg.norm(n.cross((p2 - p1), (p3 - p1)))/2.
def convex_hull(points, graphic=False, smidgen=0.0075):
'''
Calculate subset of points that make a convex hull around points
Recursively eliminates points that lie inside two neighbouring points until only convex hull is remaining.
:Parameters:
points : ndarray (2 x m)
array of points for which to find hull
graphic : bool
use pylab to show progress?
smidgen : float
offset for graphic number labels - useful values depend on your data range
:Returns:
hull_points : ndarray (2 x n)
convex hull surrounding points
'''
if graphic:
p.clf()
p.plot(points[0], points[1], 'ro')
n_pts = points.shape[1]
assert(n_pts > 5)
centre = points.mean(1)
if graphic: p.plot((centre[0],),(centre[1],),'bo')
angles = n.apply_along_axis(_angle_to_point, 0, points, centre)
pts_ord = points[:,angles.argsort()]
if graphic:
for i in xrange(n_pts):
p.text(pts_ord[0,i] + smidgen, pts_ord[1,i] + smidgen, \
'%d' % i)
pts = [x[0] for x in zip(pts_ord.transpose())]
prev_pts = len(pts) + 1
k = 0
while prev_pts > n_pts:
prev_pts = n_pts
n_pts = len(pts)
if graphic: p.gca().patches = []
i = -2
while i < (n_pts - 2):
Aij = area_of_triangle(centre, pts[i], pts[(i + 1) % n_pts])
Ajk = area_of_triangle(centre, pts[(i + 1) % n_pts], \
pts[(i + 2) % n_pts])
Aik = area_of_triangle(centre, pts[i], pts[(i + 2) % n_pts])
if graphic:
_draw_triangle(centre, pts[i], pts[(i + 1) % n_pts], \
facecolor='blue', alpha = 0.2)
_draw_triangle(centre, pts[(i + 1) % n_pts], \
pts[(i + 2) % n_pts], \
facecolor='green', alpha = 0.2)
_draw_triangle(centre, pts[i], pts[(i + 2) % n_pts], \
facecolor='red', alpha = 0.2)
if Aij + Ajk < Aik:
if graphic: p.plot((pts[i + 1][0],),(pts[i + 1][1],),'go')
del pts[i+1]
i += 1
n_pts = len(pts)
k += 1
return n.asarray(pts)
if __name__ == "__main__":
import scipy.interpolate as interpolate
# fig = p.figure(figsize=(10,10))
theta = 2*n.pi*n.random.rand(1000)
r = n.random.rand(1000)**0.5
x,y = r*p.cos(theta),r*p.sin(theta)
points = n.ndarray((2,len(x)))
points[0,:],points[1,:] = x,y
scale = 1.03
hull_pts = scale*convex_hull(points)
p.plot(x,y,'ko')
x,y = [],[]
convex = scale*hull_pts
for point in convex:
x.append(point[0])
y.append(point[1])
x.append(convex[0][0])
y.append(convex[0][1])
x,y = n.array(x),n.array(y)
#Taken from https://stackoverflow.com/questions/14344099/numpy-scipy-smooth-spline-representation-of-an-arbitrary-contour-flength
nt = n.linspace(0, 1, 100)
t = n.zeros(x.shape)
t[1:] = n.sqrt((x[1:] - x[:-1])**2 + (y[1:] - y[:-1])**2)
t = n.cumsum(t)
t /= t[-1]
x2 = interpolate.spline(t, x, nt)
y2 = interpolate.spline(t, y, nt)
p.plot(x2, y2,'r--',linewidth=2)
p.show()
There are some useful papers, eg.:
http://repositorium.sdum.uminho.pt/bitstream/1822/6429/1/ConcaveHull_ACM_MYS.pdf
Also, you could try with: http://resources.arcgis.com/en/help/main/10.1/index.html#//007000000013000000
I don't know nothing about arcgis, but it looks fine.
I came across this and implemented easy to use functions as well as a couple of alternatives/improvements.
Improvements:
scipy.interpolate.spline functionAlternatives:
Hope this helps someone along the way.
import sklearn.preprocessing
import sklearn.pipeline
import scipy.spatial
import numpy as np
def calculate_hull(
X,
scale=1.1,
padding="scale",
n_interpolate=100,
interpolation="quadratic_periodic",
return_hull_points=False):
"""
Calculates a "smooth" hull around given points in `X`.
The different settings have different drawbacks but the given defaults work reasonably well.
Parameters
----------
X : np.ndarray
2d-array with 2 columns and `n` rows
scale : float, optional
padding strength, by default 1.1
padding : str, optional
padding mode, by default "scale"
n_interpolate : int, optional
number of interpolation points, by default 100
interpolation : str or callable(ix,iy,x), optional
interpolation mode, by default "quadratic_periodic"
Inspired by: https://stackoverflow.com/a/17557853/991496
"""
if padding == "scale":
# scaling based padding
scaler = sklearn.pipeline.make_pipeline(
sklearn.preprocessing.StandardScaler(with_std=False),
sklearn.preprocessing.MinMaxScaler(feature_range=(-1,1)))
points_scaled = scaler.fit_transform(X) * scale
hull_scaled = scipy.spatial.ConvexHull(points_scaled, incremental=True)
hull_points_scaled = points_scaled[hull_scaled.vertices]
hull_points = scaler.inverse_transform(hull_points_scaled)
hull_points = np.concatenate([hull_points, hull_points[:1]])
elif padding == "extend" or isinstance(padding, (float, int)):
# extension based padding
# TODO: remove?
if padding == "extend":
add = (scale - 1) * np.max([
X[:,0].max() - X[:,0].min(),
X[:,1].max() - X[:,1].min()])
else:
add = padding
points_added = np.concatenate([
X + [0,add],
X - [0,add],
X + [add, 0],
X - [add, 0]])
hull = scipy.spatial.ConvexHull(points_added)
hull_points = points_added[hull.vertices]
hull_points = np.concatenate([hull_points, hull_points[:1]])
else:
raise ValueError(f"Unknown padding mode: {padding}")
# number of interpolated points
nt = np.linspace(0, 1, n_interpolate)
x, y = hull_points[:,0], hull_points[:,1]
# ensures the same spacing of points between all hull points
t = np.zeros(x.shape)
t[1:] = np.sqrt((x[1:] - x[:-1])**2 + (y[1:] - y[:-1])**2)
t = np.cumsum(t)
t /= t[-1]
# interpolation types
if interpolation is None or interpolation == "linear":
x2 = scipy.interpolate.interp1d(t, x, kind="linear")(nt)
y2 = scipy.interpolate.interp1d(t, y, kind="linear")(nt)
elif interpolation == "quadratic":
x2 = scipy.interpolate.interp1d(t, x, kind="quadratic")(nt)
y2 = scipy.interpolate.interp1d(t, y, kind="quadratic")(nt)
elif interpolation == "quadratic_periodic":
x2 = scipy.interpolate.splev(nt, scipy.interpolate.splrep(t, x, per=True, k=4))
y2 = scipy.interpolate.splev(nt, scipy.interpolate.splrep(t, y, per=True, k=4))
elif interpolation == "cubic":
x2 = scipy.interpolate.CubicSpline(t, x, bc_type="periodic")(nt)
y2 = scipy.interpolate.CubicSpline(t, y, bc_type="periodic")(nt)
else:
x2 = interpolation(t, x, nt)
y2 = interpolation(t, y, nt)
X_hull = np.concatenate([x2.reshape(-1,1), y2.reshape(-1,1)], axis=1)
if return_hull_points:
return X_hull, hull_points
else:
return X_hull
def draw_hull(
X,
scale=1.1,
padding="scale",
n_interpolate=100,
interpolation="quadratic_periodic",
plot_kwargs=None,
ax=None):
"""Uses `calculate_hull` to draw a hull around given points.
Parameters
----------
X : np.ndarray
2d-array with 2 columns and `n` rows
scale : float, optional
padding strength, by default 1.1
padding : str, optional
padding mode, by default "scale"
n_interpolate : int, optional
number of interpolation points, by default 100
interpolation : str or callable(ix,iy,x), optional
interpolation mode, by default "quadratic_periodic"
plot_kwargs : dict, optional
`matplotlib.pyplot.plot` kwargs, by default None
ax : `matplotlib.axes.Axes`, optional
[description], by default None
"""
if plot_kwargs is None:
plot_kwargs = {}
X_hull = calculate_hull(
X, scale=scale, padding=padding, n_interpolate=n_interpolate, interpolation=interpolation)
if ax is None:
ax= plt.gca()
plt.plot(X_hull[:,0], X_hull[:,1], **plot_kwargs)
def draw_rounded_hull(X, padding=0.1, line_kwargs=None, ax=None):
"""Plots a convex hull around points with rounded corners and a given padding.
Parameters
----------
X : np.array
2d array with two columns and n rows
padding : float, optional
padding between hull and points, by default 0.1
line_kwargs : dict, optional
line kwargs (used for `matplotlib.pyplot.plot` and `matplotlib.patches.Arc`), by default None
ax : matplotlib.axes.Axes, optional
axes to plat on, by default None
"""
default_line_kwargs = dict(
color="black",
linewidth=1
)
if line_kwargs is None:
line_kwargs = default_line_kwargs
else:
line_kwargs = {**default_line_kwargs, **line_kwargs}
if ax is None:
ax = plt.gca()
hull = scipy.spatial.ConvexHull(X)
hull_points = X[hull.vertices]
hull_points = np.concatenate([hull_points[[-1]], hull_points, hull_points[[0]]])
diameter = padding * 2
for i in range(1, hull_points.shape[0] - 1):
# line
# source: https://stackoverflow.com/a/1243676/991496
norm_next = np.flip(hull_points[i] - hull_points[i + 1]) * [-1, 1]
norm_next /= np.linalg.norm(norm_next)
norm_prev = np.flip(hull_points[i - 1] - hull_points[i]) * [-1, 1]
norm_prev /= np.linalg.norm(norm_prev)
# plot line
line = hull_points[i:i+2] + norm_next * diameter / 2
ax.plot(line[:,0], line[:,1], **line_kwargs)
# arc
angle_next = np.rad2deg(np.arccos(np.dot(norm_next, [1,0])))
if norm_next[1] < 0:
angle_next = 360 - angle_next
angle_prev = np.rad2deg(np.arccos(np.dot(norm_prev, [1,0])))
if norm_prev[1] < 0:
angle_prev = 360 - angle_prev
arc = patches.Arc(
hull_points[i],
diameter, diameter,
angle=0, fill=False, theta1=angle_prev, theta2=angle_next,
**line_kwargs)
ax.add_patch(arc)
if __name__ == '__main__':
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import patches
# np.random.seed(42)
X = np.random.random((20,2))
fig, ax = plt.subplots(1,1, figsize=(10,10))
ax.scatter(X[:,0], X[:,1])
draw_rounded_hull(X, padding=0.1)
draw_hull(X)
ax.set(xlim=[-1,2], ylim= [-1,2])
fig.savefig("_out/test.png")
