[数据展示]绘制三维曲面，二维投影以及相关问题

Rainyboy

在这个贴子里，我主要想跟大家讨论如下几个问题：

1，一个我自己写的绘制三维曲面和二维投影（等高线）等的函数，早在帖子：用matplotlib及mplot3d绘的图 的时候就开始使用，自己用得挺顺手，分享给大家：


[code=Python width=600px]# -*- 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[0])
        ax.set_xlim3d(Xost,XMax)
        ax.set_ylabel(lables[1])
        ax.set_ylim3d(YMin,Yost)
        ax.set_zlabel(lables[2])
        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[0])
        ax.set_ylabel(lables[1])
        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[0])
        ax.set_ylabel(lables[1])
        CB = fig.colorbar(CS, shrink=0.8, extend='both')[/code]

这些代码不难写出，但是如果每次遇到都要重写一遍的话，就没必要了。



2，绘制此类图形之前的数据组织问题，因为一般来讲，最终传递给绘图库函数的数据（我上面给出的函数也是这样）总是要组织成三个矩阵（二维序列）的形式：X[:,:],Y[:,:],Z[:,:]，其对应关系是


      F(X[i,j],Y[i,j])=Z[i,j],   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
     就需要重新整理数据。比如，对于下面这样的函数：

     [code=Python width=600px]def f(x,y):
    return 2*np.power(x,2)+np.power(y,2)[/code]

    如果我们的取值点是随机生成的，比如这样：

[code=Python width=600px]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[/code]

    那么直接reshape再画图是不行的：

[code=Python width=600px]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[/code]

[code=Python width=600px]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)[/code]

     结果是：



如果不想修改原始数据，可以将上述数据进行两次排序，主要是让空间中相近的点储存在X,Y,Z矩阵中的临近位置中。
      看代码：

      [code=Python width=600px]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[ag]
    y_t = y_t[ag]
    z_t = z_t[ag]
   
    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[[slice(xt) for xt in X.shape]])
    t[0] = X.argsort(0)
    X = X[t]
    Y = Y[t]
    Z = Z[t]
      
    t = list(np.ogrid[[slice(xt) for xt in X.shape]])
    t[-1] = Y.argsort(-1)
    X = X[t]
    Y = Y[t]
    Z = Z[t]
   
    return X,Y,Z[/code]

      在绘图之前调用这个函数即可：

[code=Python width=600px]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)[/code]

结果是：



可以看到大致差不多了。

    另一种方法实际上是对原始数据进行了插值，用均匀的网格对原始数据的差值函数进行采样，实际上已经有库函数干了这个活儿，它的名字在PYTHON.SCIPY中和MATLAB中是一样的——griddata：

[code=Python width=600px]from scipy.interpolate import griddata
[/code]

因此我们直接使用就可以了，注意在用之前需要先生成均匀的二维序列：

[code=Python width=600px]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)[/code]

这样结果就非常漂亮了：





    如果这些储存在1维序列中的数据点在空间的分布本身是均匀的，比如用如下点进行采样：


[code=Python width=600px]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[/code]

    那么上述两种方法所得的结果就是一样的了。由于在第二种方法中重新采样的点数不一定是原数据点数，可以自己设定，这样如果原始数据点太多，也不失为一种减少数据点数的方法？

较密的数据点及图形：


差值后重新取的数据点及图形：



最后再贴一个上面的代码与ANSYS进行后处理的例子。


在ANSYS中，计算这样的一个压电悬臂梁在端部受力时的静变形：



图中紫色的是压电材料：




我想确定压电梁的横截面正应变S1是否还满足经典欧拉梁的中性面假设，即S1的分布是否还呈一个过Y轴的斜面，这样最好是将某个横截面的正应力提取出来，再画成三维图形。虽然网格是均匀的，但提取出来的数据储存在文件中，却是一维序列的形式，这就涉及上面所说的问题了。

由于横截面的网格本身就是均匀的，只是顺序乱了，所以不需要进行差值，通过两次排序就可以了，结果是：


仍然满足欧拉伯努利梁的基本假设。

mayaview

Matplotlib的3D库很弱而且特别慢，他的作者也不建议用来做复杂的运用。