Rainyboy 发表于 2013-11-21 05:05

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

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

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轴的斜面,这样最好是将某个横截面的正应力提取出来,再画成三维图形。虽然网格是均匀的,但提取出来的数据储存在文件中,却是一维序列的形式,这就涉及上面所说的问题了。

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


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



mayaview 发表于 2013-12-13 16:41

Matplotlib的3D库很弱而且特别慢,他的作者也不建议用来做复杂的运用。
页: [1]
查看完整版本: [数据展示]绘制三维曲面,二维投影以及相关问题