Python不规则x，y数据到原始域上的等高线图

Question

我在“x-cord”、“y-cord”、“value”列下有包含点的文件。这些是不规则间隔的。我正在尝试制作“价值”的等高线图并将其覆盖在原始域上。我放弃了在 pgfplots 和 matlab 中尝试这样做，并认为我会给 python 一个机会。任何这些脚本中的答案都可以。 python脚本如下

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import numpy.ma as ma
from numpy.random import uniform, seed
from scipy.spatial import ConvexHull
#
# Loading data
filename = "strain.dat"
coordinates = []
x_c = []
y_c = []
z_c = []
xyz = open(filename)
title = xyz.readline()
for line in xyz:
    x,y,z = line.split()
    coordinates.append([float(x), float(y), float(z)])
    x_c.append([float(x)])
    y_c.append([float(y)])
    z_c.append([float(z)])
xyz.close()
#
# Rehaping and translating data
x_c=np.ravel(np.array(x_c))
y_c=np.ravel(np.array(y_c))
z_c=np.ravel(np.array(z_c))
x_c = x_c-100.0
y_c = y_c-100.0
#
# Checking the convex hull
points=np.column_stack((x_c,y_c))
hull = ConvexHull(points);
plt.plot(points[hull.vertices,0], points[hull.vertices,1], 'r--', lw=2)
plt.scatter(x_c, y_c, marker='o', s=5, zorder=10)
#
# Mapping the irregular data onto a regular grid and plotting
xic = np.linspace(min(x_c), max(x_c), 1000)
yic = np.linspace(min(y_c), max(y_c), 1000)
zic = griddata((x_c, y_c), z_c, (xic[None,:], yic[:,None]))
CS = plt.contour(xic,yic,zic,15,linewidths=0.5,colors='k')
CS = plt.contourf(xic,yic,zic,15,cmap=plt.cm.summer)
plt.colorbar() # draw colorbar
#
#plt.scatter(x_c, y_c, marker='o', s=5, zorder=10)
plt.axis('equal')
plt.savefig('foo.pdf', bbox_inches='tight')
plt.show()

输出看起来像

问题是griddata使用了一个凸包，而这个凸包超出了不规则数据的边缘。有什么办法可以将原始点边界边缘之外的griddata点的值设置为零？

编辑

最后我认输了，回到了 Matlab。我必须将数据导出到 pgfplots 以获得漂亮的图。我想出的代码是

x = strain.x;
y = strain.y;
z = strain.eps;
% Get the alpha shape (couldn't do this in python easily)
shp = alphaShape(x,y,.001);
% Get the boundary nodes
[bi, xy] = boundaryFacets(shp);
no_grid = 500;
xb=xy(:,1);
yb=xy(:,2);
[X,Y] = ndgrid(linspace(min(x),max(x),no_grid),linspace(min(y),max(y),no_grid));
Z = griddata(x,y,z,X,Y,'v4');

% Got through the regular grid and set the values which are outside the boundary of the original domain to Nans
for j = 1:no_grid
    [in,on] = inpolygon(X(:,j),Y(:,j),xb,yb);
    Z(~in,j) = NaN;
end

contourf(X,Y,Z,10),axis equal
colorbar
hold on
plot(xb,yb)
axis equal
hold off

这是生成的图像。

如果有人可以在 Python 中做类似的事情，我会很乐意接受答案。

Answer 1

我必须在复杂几何图形上绘制插值数据（参见图中的蓝点）P(x,z)（z 是水平坐标）。我使用了掩码操作，效果很好。没有遮罩，整个正方形 (x=0..1 ; z=0..17.28) 被轮廓覆盖。

## limiting values for geometry
xmax1=0.408
zmin1=6.
xmax2=0.064
zmin2=13.12
xmin=0.
xmax=1.
zmin=0.
zmax=17.28

# Grid for points
x1 = np.arange(xmin,xmax+dx,dx)
z1 = np.arange(zmin,zmax+dz,dz)
zi2,xi2 = np.meshgrid(z1,x1)
mask = (((zi2 > zmin2) & (xi2 > xmax2)) | ((zi2 > zmin1) & (zi2 <= zmin2) & (xi2 > xmax1)))
zim=np.ma.masked_array(zi2,mask)
xim=np.ma.masked_array(xi2,mask)

# Grid for P values
# npz=z coordinates of data, npx is the x coordinates and npp is P values
grid_p = scipy.interpolate.griddata((npz, npx), npp, (zim,xim),method='nearest')
pm=np.ma.masked_array(grid_p,mask)

# plot 
plt.contour(zim, xim, pm, 25, linewidths=0.5, colors='k',corner_mask=False)
plt.contourf(zim, xim, pm, 25,vmax=grid_p.max(), vmin=grid_p.min(),corner_mask=False)
plt.colorbar()

# Scatter plot to check 
plt.scatter(npz,npr, marker='x', s=2)

plt.show()

enter image description here