[数据展示]绘制三维曲面，二维投影以及相关问题
在这个贴子里，我主要想跟大家讨论如下几个问题：1，一个我自己写的绘制三维曲面和二维投影（等高线）等的函数，早在帖子：用matplotlib及mplot3d绘的图 的时候就开始使用，自己用得挺顺手，分享给大家：
# -*- coding: cp936 -*-
import numpy as math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
def Surf_and_Z_Cont(X,Y,Z,
lables = ('X','Y','Z'),
title = '',
Is3Dsurf = True,
IsContour = True,
IsFillContour = True,
IsCheck = True,
IsSave = False
):
XMax = math.max(X)
XMin = math.min(X)
Xost = XMin - (XMax - XMin)*0.1
YMax = math.max(Y)
YMin = math.min(Y)
Yost = YMax + (YMax - YMin)*0.1
ZMax = math.max(Z)
ZMin = math.min(Z)
Zost = ZMin - (ZMax - ZMin)*0.1
if IsCheck:
print 'Surf_and_Z_Cont Check information'
print 'X.shape = ', X.shape
print 'Y.shape = ', Y.shape
print 'Z.shape = ', Z.shape
print 'Surf_and_Z_Cont Check information'
if Is3Dsurf:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_title(title)
ax.plot_surface(X,Y,Z,rstride=1, cstride=1,alpha=0.2, lw = 0.2)
cset = ax.contour( X,Y, Z, zdir='x', offset=Xost,levels = math.linspace(XMin,XMax,20))
cset = ax.contour( X,Y, Z, zdir='y', offset=Yost,levels = math.linspace(YMin,YMax,20))
cset = ax.contour( X,Y, Z, zdir='z', offset=Zost,levels = math.linspace(ZMin,ZMax,20))
ax.set_xlabel(lables)
ax.set_xlim3d(Xost,XMax)
ax.set_ylabel(lables)
ax.set_ylim3d(YMin,Yost)
ax.set_zlabel(lables)
ax.set_zlim3d(Zost,ZMax)
if IsSave:
fig.savefig('3D_'+ title, dpi=600)
if IsFillContour:
fig = plt.figure()
ax = fig.add_subplot(111)
CS = ax.contourf(X,Y,Z,25,cmap=cm.jet)
ax.set_xlabel(lables)
ax.set_ylabel(lables)
CBI = fig.colorbar(CS, orientation='vertical', shrink=0.8,drawedges=True)
if IsContour:
fig = plt.figure()
ax = fig.add_subplot(111)
CS = ax.contour(X,Y,Z,25,cmap=cm.jet)
ax.set_xlabel(lables)
ax.set_ylabel(lables)
CB = fig.colorbar(CS, shrink=0.8, extend='both')
这些代码不难写出，但是如果每次遇到都要重写一遍的话，就没必要了。
2，绘制此类图形之前的数据组织问题，因为一般来讲，最终传递给绘图库函数的数据（我上面给出的函数也是这样）总是要组织成三个矩阵（二维序列）的形式：X[:,:],Y[:,:],Z[:,:]，其对应关系是
F(X,Y)=Z, i=0,1,2,...,M-1, j=0,1,2,...,N-1
但在一些情形下，并不能直接获得这样的数据，比如数据是从ANSYS中导出的，储存在文件中，读入程序中的初始状态是三个一维序列 X[:], Y[:], Z[:]，对应关系是
F(X,Y)=Z, i = 0, 1, 2, .....,M*N-1
就需要重新整理数据。比如，对于下面这样的函数：
def f(x,y):
return 2*np.power(x,2)+np.power(y,2)
如果我们的取值点是随机生成的，比如这样：
def Get_XYZ_RandXY(Cxlines,Cylines):
Xlist = []
Ylist = []
Zlist = []
for j in range(0,Cxlines):
for i in range(0,Cylines):
x = np.random.random()*2-1
y = np.random.random()*2-1
Xlist.append(x)
Ylist.append(y)
Zlist.append(f(x,y))
return Xlist,Ylist,Zlist
那么直接reshape再画图是不行的：
def ReshapeData(xlist,ylist,zlist,Cx,Cy):
x_t = np.array(xlist)
y_t = np.array(ylist)
z_t = np.array(zlist)
X = np.reshape(x_t,(Cx,Cy))
Y = np.reshape(y_t,(Cx,Cy))
Z = np.reshape(z_t,(Cx,Cy))
return X,Y,Z
Cxlines = 50
Cylines = 50
Xlist, Ylist, Zlist = Get_XYZ_RandXY(Cxlines,Cylines)
X,Y,Z = ReshapeData(Xlist, Ylist, Zlist,Cxlines,Cylines)
Surf_and_Z_Cont(X,Y,Z,
lables = ('X','Y','Z'),
Is3Dsurf = False,
IsContour = False,
IsFillContour = True,
IsCheck = True)
plt.plot(Xlist,Ylist, 'k.', ms=1)
结果是：
如果不想修改原始数据，可以将上述数据进行两次排序，主要是让空间中相近的点储存在X,Y,Z矩阵中的临近位置中。
看代码：
def ReorgnizeDataForContourPlot(xlist,ylist,zlist,Cx,Cy):
x_t = np.array(xlist)
y_t = np.array(ylist)
z_t = np.array(zlist)
ag = x_t.argsort()
x_t = x_t
y_t = y_t
z_t = z_t
X = np.reshape(x_t,(Cx,Cy))
Y = np.reshape(y_t,(Cx,Cy))
Z = np.reshape(z_t,(Cx,Cy))
t = list(np.ogrid[])
t = X.argsort(0)
X = X
Y = Y
Z = Z
t = list(np.ogrid[])
t[-1] = Y.argsort(-1)
X = X
Y = Y
Z = Z
return X,Y,Z
在绘图之前调用这个函数即可：
Cxlines = 50
Cylines = 50
Xlist, Ylist, Zlist = Get_XYZ_RandXY(Cxlines,Cylines)
X,Y,Z = ReorgnizeDataForContourPlot(Xlist, Ylist, Zlist,Cxlines,Cylines)
Surf_and_Z_Cont(X,Y,Z,
lables = ('X','Y','Z'),
Is3Dsurf = False,
IsContour = False,
IsFillContour = True,
IsCheck = True)
plt.plot(Xlist,Ylist, 'k.', ms=1)
结果是：
可以看到大致差不多了。
另一种方法实际上是对原始数据进行了插值，用均匀的网格对原始数据的差值函数进行采样，实际上已经有库函数干了这个活儿，它的名字在PYTHON.SCIPY中和MATLAB中是一样的——griddata：
from scipy.interpolate import griddata
因此我们直接使用就可以了，注意在用之前需要先生成均匀的二维序列：
Cxlines = 100
Cylines = 100
Xlist, Ylist, Zlist = Get_XYZ_RandXY(Cxlines,Cylines)
grid_x, grid_y = np.mgrid[-1.0:1.0:50j, -1.0:1.0:50j]
grid_z0 = griddata((Xlist,Ylist),Zlist,(grid_x, grid_y),method='cubic')
print np.max(grid_z0)
Surf_and_Z_Cont(grid_x,grid_y,grid_z0,
lables = ('X','Y','Z'),
Is3Dsurf = False,
IsContour = False,
IsFillContour = True,
IsCheck = True)
plt.plot(grid_x,grid_y, 'k.', ms=1)
这样结果就非常漂亮了：
如果这些储存在1维序列中的数据点在空间的分布本身是均匀的，比如用如下点进行采样：
def Get_XYZ(Cxlines,Cylines):
Xlist = []
Ylist = []
Zlist = []
for x in np.linspace(-1,1,Cxlines):
for y in np.linspace(-1,1,Cylines):
Xlist.append(x)
Ylist.append(y)
Zlist.append(f(x,y))
return Xlist,Ylist,Zlist
那么上述两种方法所得的结果就是一样的了。由于在第二种方法中重新采样的点数不一定是原数据点数，可以自己设定，这样如果原始数据点太多，也不失为一种减少数据点数的方法？
较密的数据点及图形：
差值后重新取的数据点及图形：
最后再贴一个上面的代码与ANSYS进行后处理的例子。
在ANSYS中，计算这样的一个压电悬臂梁在端部受力时的静变形：
图中紫色的是压电材料：
我想确定压电梁的横截面正应变S1是否还满足经典欧拉梁的中性面假设，即S1的分布是否还呈一个过Y轴的斜面，这样最好是将某个横截面的正应力提取出来，再画成三维图形。虽然网格是均匀的，但提取出来的数据储存在文件中，却是一维序列的形式，这就涉及上面所说的问题了。
由于横截面的网格本身就是均匀的，只是顺序乱了，所以不需要进行差值，通过两次排序就可以了，结果是：
仍然满足欧拉伯努利梁的基本假设。
Matplotlib的3D库很弱而且特别慢，他的作者也不建议用来做复杂的运用。
页:
[1]