我想做一个 matplotlib contour 或 contourf 非网格 3D 数据 (x, y, z) 的图，它在 x 和 y 中呈某种 C 形（见草图）——因此部分围绕数据的封闭外壳在 x 和 y 上是凹的。
通常我会先用插值来绘制非网格 3D 数据的图
from matplotlib.mlab import griddata
griddata...
但这会在数据的凹入部分产生伪影，从而使凹入部分被插值填充。
是否可以进行插值或轮廓/等高线图，使得 尊重数据的凹入部分？
【问题讨论】：
标签： python matplotlib plot
下面是一个示例，说明如何使用 tricontourf 和遮罩来获得没有数据外插值部分的凹形。它依赖于根据条件屏蔽数据的能力。
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
# create some data
rawx = np.random.rand(500)
rawy = np.random.rand(len(rawx))
cond01 = (rawx-1)**2 + rawy**2 <=1
cond02 = (rawx-0.7)**2 + rawy**2 >0.3
x = rawx[cond01 & cond02]
y = rawy[cond01 & cond02]
f = lambda x,y: np.sin(x*4)+np.cos(y)
z = f(x,y)
# now, x,y are points within a partially concave shape
triang0 = tri.Triangulation(x, y)
triang = tri.Triangulation(x, y)
x2 = x[triang.triangles].mean(axis=1)
y2 = y[triang.triangles].mean(axis=1)
#note the very obscure mean command, which, if not present causes an error.
#now we need some masking condition.
# this is easy in this case where we generated the data according to the same condition
cond1 = (x2-1)**2 + y2**2 <=1
cond2 = (x2-0.7)**2 + (y2)**2 >0.3
mask = np.where(cond1 & cond2,0,1)
# apply masking
triang.set_mask(mask)
fig, (ax, ax2) = plt.subplots(ncols=2, figsize=(6,3))
ax.set_aspect("equal")
ax2.set_aspect("equal")
ax.tricontourf(triang0, z, cmap="Oranges")
ax.scatter(x,y, s=3, color="k")
ax2.tricontourf(triang, z, cmap="Oranges")
ax2.scatter(x,y, s=3, color="k")
ax.set_title("tricontourf without mask")
ax2.set_title("tricontourf with mask")
ax.set_xlim(0,1)
ax.set_ylim(0,1)
ax2.set_xlim(0,1)
ax2.set_ylim(0,1)
plt.show()
如果您无法获得确切的条件，但具有点之间的最大边长（距离），则以下将是一个解决方案。它将掩盖至少一侧长于某个最大距离的所有三角形。如果点密度相当高，这可以很好地应用。
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
# create some data
rawx = np.random.rand(500)
rawy = np.random.rand(len(rawx))
cond01 = (rawx-1)**2 + rawy**2 <=1
cond02 = (rawx-0.7)**2 + rawy**2 >0.3
x = rawx[cond01 & cond02]
y = rawy[cond01 & cond02]
f = lambda x,y: np.sin(x*4)+np.cos(y)
z = f(x,y)
# now, x,y are points within a partially concave shape
triang1 = tri.Triangulation(x, y)
triang2 = tri.Triangulation(x, y)
triang3 = tri.Triangulation(x, y)
def apply_mask(triang, alpha=0.4):
# Mask triangles with sidelength bigger some alpha
triangles = triang.triangles
# Mask off unwanted triangles.
xtri = x[triangles] - np.roll(x[triangles], 1, axis=1)
ytri = y[triangles] - np.roll(y[triangles], 1, axis=1)
maxi = np.max(np.sqrt(xtri**2 + ytri**2), axis=1)
# apply masking
triang.set_mask(maxi > alpha)
apply_mask(triang2, alpha=0.1)
apply_mask(triang3, alpha=0.3)
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(9,3))
ax1.tricontourf(triang1, z, cmap="Oranges")
ax1.scatter(x,y, s=3, color="k")
ax2.tricontourf(triang2, z, cmap="Oranges")
ax2.scatter(x,y, s=3, color="k")
ax3.tricontourf(triang3, z, cmap="Oranges")
ax3.scatter(x,y, s=3, color="k")
ax1.set_title("tricontourf without mask")
ax2.set_title("with mask (alpha=0.1)")
ax3.set_title("with mask (alpha=0.3)")
for ax in (ax1, ax2, ax3):
ax.set(xlim=(0,1), ylim=(0,1), aspect="equal")
plt.show()
可以看出，在这里找到正确的参数 (alpha) 可能需要一些调整。
【讨论】：
TheImportanceOfBeingErnest 提供的答案给了我一个完美的起点，我操纵了上面的代码，为凹壳/阿尔法形状的等高线图提供了一个通用解决方案。我使用 Python 包形状地创建了凹壳的多边形。这不是一个需要添加的内置函数，我取自this post。一旦我们有了多边形，我们只需检查三角点的平均值是否在多边形内，并将其用作我们形成掩码的条件。
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
# Start Using SHAPELY
import shapely.geometry as geometry
from shapely.geometry import Polygon, MultiPoint, Point
from shapely.ops import triangulate
from shapely.ops import cascaded_union, polygonize
from scipy.spatial import Delaunay
from descartes.patch import PolygonPatch
import math
# create some data
rawx = np.random.rand(500)
rawy = np.random.rand(len(rawx))
cond01 = (rawx-1)**2 + rawy**2 <=1
cond02 = (rawx-0.7)**2 + rawy**2 >0.3
x = rawx[cond01 & cond02]
y = rawy[cond01 & cond02]
f = lambda x,y: np.sin(x*4)+np.cos(y)
z = f(x,y)
# now, x,y are points within a partially concave shape
triang0 = tri.Triangulation(x, y)
triang = tri.Triangulation(x, y)
# Function for finding an alpha function
def alpha_shape(points, alpha):
"""
Compute the alpha shape (concave hull) of a set
of points.
@param points: Iterable container of points.
@param alpha: alpha value to influence the
gooeyness of the border. Smaller numbers
don't fall inward as much as larger numbers.
Too large, and you lose everything!
"""
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull
def add_edge(edges, edge_points, coords, i, j):
"""
Add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(coords[ [i, j] ])
coords = np.array([point.coords[0]
for point in points])
tri = Delaunay(coords)
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the
# triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = math.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
#print circum_r
if circum_r < 1.0/alpha:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
#import ipdb; ipdb.set_trace()
return cascaded_union(triangles), edge_points
# create array of points from reduced exp data to convert to Polygon
crds=np.array( [x,y]).transpose()
# Adjust the length of acceptable sides by adjusting the alpha parameter
concave_hull, edge_points = alpha_shape(MultiPoint(crds), alpha=2.3)
# View the polygon and adjust alpha if needed
def plot_polygon(polygon):
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
margin = .3
x_min, y_min, x_max, y_max = polygon.bounds
ax.set_xlim([x_min-margin, x_max+margin])
ax.set_ylim([y_min-margin, y_max+margin])
patch = PolygonPatch(polygon, fc='#999999',
ec='#000000', fill=True,
zorder=-1)
ax.add_patch(patch)
return fig
plot_polygon(concave_hull); plt.plot(x,y,'.'); #plt.show()
# Use the mean distance between the triangulated x & y poitns
x2 = x[triang.triangles].mean(axis=1)
y2 = y[triang.triangles].mean(axis=1)
##note the very obscure mean command, which, if not present causes an error.
##now we need some masking condition.
# Create an empty set to fill with zeros and ones
cond = np.empty(len(x2))
# iterate through points checking if the point lies within the polygon
for i in range(len(x2)):
cond[i] = concave_hull.contains(Point(x2[i],y2[i]))
mask = np.where(cond,0,1)
# apply masking
triang.set_mask(mask)
fig, (ax, ax2) = plt.subplots(ncols=2, figsize=(6,3))
ax.set_aspect("equal")
ax2.set_aspect("equal")
ax.tricontourf(triang0, z, cmap="Oranges")
ax.scatter(x,y, s=3, color="k")
ax2.tricontourf(triang, z, cmap="Oranges")
ax2.scatter(x,y, s=3, color="k")
ax.set_title("tricontourf without mask")
ax2.set_title("tricontourf with mask")
ax.set_xlim(0,1)
ax.set_ylim(0,1)
ax2.set_xlim(0,1)
ax2.set_ylim(0,1)
plt.show()
【讨论】：