杭州锐达数字技术有限公司
查看: 2473|回复: 1

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

[复制链接]
发表于 2013-11-21 05:05 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

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

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


Python 代码,双击复制代码
# -*- 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')


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



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
     就需要重新整理数据。比如,对于下面这样的函数:

     
Python 代码,双击复制代码
def f(x,y):
    return 2*np.power(x,2)+np.power(y,2)


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

Python 代码,双击复制代码
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再画图是不行的:

Python 代码,双击复制代码
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


Python 代码,双击复制代码
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)


     结果是:
1.png


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

      
Python 代码,双击复制代码
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


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

Python 代码,双击复制代码
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)


结果是:

image.png

可以看到大致差不多了。

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

Python 代码,双击复制代码
from scipy.interpolate import griddata


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

Python 代码,双击复制代码
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)


这样结果就非常漂亮了:

3.png



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


Python 代码,双击复制代码
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


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

较密的数据点及图形:
11.png

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

22.png

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


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

QQ截图20131120230741.png

图中紫色的是压电材料:

QQ截图20131120230816.png


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

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

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



回复
分享到:

使用道具 举报

发表于 2013-12-13 16:41 | 显示全部楼层
Matplotlib的3D库很弱而且特别慢,他的作者也不建议用来做复杂的运用。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表