Python和科学计算认证群组  - 讨论区

标题:请高手解答科学3D画图技巧

2011年03月28日 星期一 12:01

能否将一个风速数据画成如下三维图:

我想用如下的数据画成上面这个3D图,带箭头的,然后将每个箭头末端连线(黄色的)映射到下面的一个极坐标中,极坐标画图我已经在上一个帖子中说了,请 张若愚 高手能给个提示或者有个实例吗,谢谢了。

一个数据样本如下:

--------------------
  高度  风向  风速
    m   deg   knot
--------------------
   317   215     20
   522   198     27
   610   190     30
   797   195     30
   914   195     28
  1219   195     22
  1540   210     28
  1571   210     28
  1632   210     28
  1829   210     28
  2134   215     30
  2438   215     29
  2743   230     29
  2998   234     27
  3048   235     26
  3080   235     26
  3139   235     26
  3199   235     26
  3658   230     24
  4165   218     24
  4267   215     24
  4271   215     24
  4572   210     24
  4610   211     24
  4877   220     23
  4947   223     23
  5486   250     22
  5940   255     19
  5987   255     19

2011年03月28日 星期一 16:25

可以用 风向  风速 计算每一个XY平面上X Y的值,然后用mlab.plot3d画出三维连线图,

箭头的起始点都知道,相当于画一个矢量,然后用mlab.quiver3d画的箭头

 

呵呵,都是从用Python做科学计算试读版例子凑出来的

但是Z(高度)变化比较大,画出的效果不太好,取对数之后也不好,可能要单独设置Z轴的间距吧

 

还是等能人吧~

2011年03月28日 星期一 18:10

急呀,还没人吗,在线等

2011年03月28日 星期一 19:11

查了一下,不打算用mayavi,好像是收费的吧

2011年03月28日 星期一 20:20

mayavi不是收费的。如王振所说,高度变化很大,效果可能不太好,不知道你准备怎么处理。

这个图还可以用matplotlib的三维功能,或者VPython做。

2011年03月28日 星期一 20:50

下面是VPython的程序:

 

import numpy as np

from visual import *

data = np.loadtxt("data.txt", delimiter=",")

data[:,0] /= 60

c = curve(color = color.yellow, radius=1)

for m, deg, knot in data:

    rad = np.deg2rad(deg)

    x, y = knot*np.cos(rad), knot*np.sin(rad)

    arrow(pos=(0,0,m), axis=(x,y,0), shaftwidth=1)

    arrow(pos=(0,0,0), axis=(0,0,np.max(data[:,0])+10), shaftwidth=0.5)

    c.append(pos=(x,y,m))    

scene.autocenter = True    

下面是显示效果

 

如果你觉得曲线不平滑,可以用插值:

import numpy as np

from scipy.interpolate import interp1d

from visual import *

data = np.loadtxt("data.txt", delimiter=",")

data[:,0] /= 60

rad = np.deg2rad(data[:,1])

data[:,1],data[:,2] = data[:,2]*np.cos(rad), data[:,2]*np.sin(rad)

for m, x, y in data:

    arrow(pos=(0,0,m), axis=(x,y,0), shaftwidth=1)

    arrow(pos=(0,0,0), axis=(0,0,np.max(data[:,0])+10), shaftwidth=0.5)

 

zs = np.linspace(np.min(data[:,0]), np.max(data[:,0]), 200)

xs = interp1d(data[:, 0], data[:, 1], 2)(zs)

ys = interp1d(data[:, 0], data[:, 2], 2)(zs)    

c = curve(color = color.yellow, radius=1, pos=np.c_[xs,ys,zs])

scene.autocenter = True

2011年03月30日 星期三 14:23

 

花了点时间,用matplotlib画了个图,但是没有箭头,主要是我不会,表示遗憾,

还有个急于需要解决的问题,就是如何让z轴变为semilog方式的显示方式,就是让z轴变为对数的显示方式。

  #打开图片
  fig = plt.figure(figsize=(12.6,9.2))  # (宽度,高度)
  axes = Axes3D(fig)

  zaxis_min=-2000
  zaxis_max=max(ltEt_Height1)

  #Z轴坐标 标签(标准层17层)
  ltyticks = [1013.,925.,850.,700.,600.,500.,400.,300.,\
              250.,200.,150.,100.,70.,50.,30.,20.,10.]
  #最大半径
  fMax_Radius=60.

  #X-y轴坐标
  ltxtick = np.arange(-fMax_Radius,fMax_Radius,5.)  #-60->60
  ltytick = np.arange(-fMax_Radius,fMax_Radius,5.)  #-60->60

  #画中心垂直线
  axes.plot([0,0],[0,0],[zaxis_min,zaxis_max],linestyle ='-',color='black',linewidth=1.0 )
  axes.scatter3D([0],[0],[zaxis_max],c=tp_cOrange,s=40) #顶部帽子

  #画风向杆
  for u,v,h in zip(ltEt_Uw1,ltEt_Vw1,ltEt_Height1):
    axes.plot([0,u], [0,v],[h,h],linestyle ='-',color='blue',linewidth=1.5 )
  #最大风向杆
  axes.plot([0,ltEt_Uw1[iWSmax_Index]], [0,ltEt_Vw1[iWSmax_Index]],
            [ltEt_Height1[iWSmax_Index],ltEt_Height1[iWSmax_Index]],
            linestyle ='-',color='red',linewidth=3.0 )

  #连接风向矢端
  axes.plot(ltEt_Uw1,ltEt_Vw1,ltEt_Height1,linestyle ='-',color='yellow',linewidth=2.0 )

  #底面投影风矢量(x=u,y=风速)
  axes.plot(ltEt_Uw1,ltEt_WindS1,[0]*len(ltEt_WindS1),
            linestyle ='-',color='green',linewidth=0.5 )

  #底部十字线(x-y面)
  axes.plot([-fMax_Radius,fMax_Radius],[0,0],[zaxis_min,zaxis_min],linestyle ='--',color='black',linewidth=0.5 )
  axes.plot([0,0],[-fMax_Radius,fMax_Radius],[zaxis_min,zaxis_min],linestyle ='--',color='black',linewidth=0.5 )
  #底部十字斜线(x-y面)
  x=y=[-fMax_Radius/np.sqrt(2.),fMax_Radius/np.sqrt(2.)]
  axes.plot(x,y,[zaxis_min,zaxis_min],linestyle ='--',color='black',linewidth=0.5 )
  x=[-fMax_Radius/np.sqrt(2.),fMax_Radius/np.sqrt(2.)]
  y=[fMax_Radius/np.sqrt(2.),-fMax_Radius/np.sqrt(2.)]
  axes.plot(x,y,[zaxis_min,zaxis_min],linestyle ='--',color='black',linewidth=0.5 )

  #画里圈
  for rd in np.arange(10.,fMax_Radius,10.):
    x=np.arange(-rd, rd+1, 1.)
    y=np.sqrt(rd*rd-np.square(x))
    z=[zaxis_min]*len(x)
    axes.plot(x,y,z,linestyle ='--',color='black',linewidth=0.3 )
    axes.plot(x,-y,z,linestyle ='--',color='black',linewidth=0.3 )
  #画最外圈
  x=np.arange(-fMax_Radius, fMax_Radius+1, 1.)
  y=np.sqrt(fMax_Radius*fMax_Radius-np.square(x))
  z=[zaxis_min]*len(x)
  axes.plot(x,y,z,linestyle ='-',color='black',linewidth=0.5 )
  axes.plot(x,-y,z,linestyle ='-',color='black',linewidth=0.5 )


  drift=0.3  #漂移值
  #标签坐标
  labelx = [0,
            (fMax_Radius+drift)/np.sqrt(2.),   #45
            fMax_Radius+drift,                 #90
            (fMax_Radius+drift)/np.sqrt(2.),   #135
            0,                                 #180
            (-fMax_Radius-drift)/np.sqrt(2.),  #225
            -fMax_Radius-drift,                #270
            (-fMax_Radius-drift)/np.sqrt(2.)]  #315
  labely = [fMax_Radius+1,
            (fMax_Radius+drift)/np.sqrt(2.),   #45
            0,                                 #90
            (-fMax_Radius-drift)/np.sqrt(2.),  #135
            -fMax_Radius-drift,                #180
            (-fMax_Radius-drift)/np.sqrt(2.),  #225
            0,                                 #270
            (fMax_Radius+drift)/np.sqrt(2.)]   #315
  #标签
  labeldeg = ['0',
              '45',
              '90',
              '135',
              '180',
              '225',
              '270',
              '315']
  labeldir = ['N',
              'NE',
              'E',
              'SE',
              'S',
              'SW',
              'W',
              'NW']
  labelstr=[]
  for deg,dirn in zip(labeldeg,labeldir):
    labelstr.append(deg+r'$^\circ$'+dirn)

  #标签对齐方式
  holalt=['center',
          'left',      #45
          'left',      #90
          'left',      #135
          'center',    #180
          'right',     #225
          'right',     #270
          'right']     #315
  velalt=['bottom',
          'bottom',    #45
          'center',    #90
          'top',       #135
          'top',       #180
          'top',       #225
          'center',    #270
          'bottom']    #315
  #画标签
  for x,y,lb,ha,va in zip(labelx,labely,labelstr,holalt,velalt):
    axes.text(x,y,zaxis_min,lb,
              fontsize='small',
              fontweight='bold',
              horizontalalignment=ha,
              verticalalignment=va)


  axes.set_xlim3d(-fMax_Radius,fMax_Radius)
  axes.set_ylim3d(-fMax_Radius,fMax_Radius)
  axes.set_zlim3d(zaxis_min,zaxis_max)

  plt.show()

2011年03月30日 星期三 20:18

平行于X-Y平面的箭头可以这样画,对数坐标系可以自己先取对数然后在一般坐标系中绘图。只是网格之类的不是对数坐标。

import matplotlib.pyplot as plt

from matplotlib.patches import Arrow

from mpl_toolkits.mplot3d import Axes3D

import mpl_toolkits.mplot3d.art3d as art3d

import numpy as np

fig = plt.figure()

ax = fig.add_subplot(111, projection="3d")

 

for z in np.linspace(-2*np.pi, 2*np.pi, 40):

    a = Arrow(0, 0, 5*np.sin(z), 5*np.cos(z))

    ax.add_patch(a)

    art3d.pathpatch_2d_to_3d(a, z=z, zdir="z")

 

ax.set_xlim3d(-10,10)

ax.set_ylim3d(-10,10)

ax.set_zlim3d(-10,10)

ax.set_xlabel("x")

ax.set_ylabel("y")

ax.set_zlabel("z")

plt.show()

2011年03月30日 星期三 23:02

我上面可能没有说清楚,这个对数坐标是下面大,上面小,如下图

 

  这个是一个二维的,在二维中只用调用axes_11.semilogy(x, y) 画图后,再将坐标限制翻转,就可以

  二维坐标中设置:ax.set_ylim(1013,50) 

下面是 需要画图的z坐标刻度(从底(大)==》顶(小))

 ltzticks = [1013.,925.,850.,700.,600.,500.,400.,300.,\
                 250.,200.,150.,100.,70.,50.,30.,20.,10.]

我试验了设置,但是会出错,画不出图

ax.set_zlim3d(10,-10)

 请再帮忙解释一下,看看有解决方案吗?

2011年03月30日 星期三 23:59

三维空间是可以旋转的,换一个位置看就可以让Z轴上面小,下面大。

关于不等距网格,我想可以通过设置Z轴的Locator和Formatter手工做。你试试看。z轴通过下面的属性获得:

ax.w_zaxis

2011年03月31日 星期四 08:52

没有听懂 换一个位置看就可以让Z轴上面小,下面大。

不是旋转的问题,是数据的问题,  z一开始就是大的,到了高层才变小,就是气压。气压岁高度是变小的

2011年03月31日 星期四 10:31

不是要下图的效果吗?

2011年03月31日 星期四 15:05

对,就是这个怎么弄的来着?

2011年03月31日 星期四 15:11

数据是这种的,气压随高度增加而减小,但是又想把气压标注在z轴上,而有时候没有高度项,所以只能以气压为z,风速风向的分解量为x,y。我不知道你上面的图是不是这样做的

-----------------------------
   气压   高度   风向  风速
    hPa     m    deg   knot
-----------------------------
  976.0    317    215     20
  954.0    522    198     27
  944.6    610    190     30
  925.0    797    195     30
  912.8    914    195     28
  881.6   1219    195     22
  850.0   1540    210     28
  847.0   1571    210     28
  841.0   1632    210     28
  821.9   1829    210     28
  793.1   2134    215     30
  765.5   2438    215     29
  738.6   2743    230     29
  717.0   2998    234     27
  712.7   3048    235     26
  710.0   3080    235     26
  705.0   3139    235     26
  700.0   3199    235     26
  662.3   3658    230     24
  623.0   4165    218     24
  615.3   4267    215     24
  615.0   4271    215     24
  592.7   4572    210     24
  590.0   4610    211     24
  570.9   4877    220     23
  566.0   4947    223     23
  529.2   5486    250     22
  500.0   5940    255     19

2011年03月31日 星期四 16:32

那个图只是旋转了一下观察角度,让它上下颠倒。

我也不知道如何同时标注高度和气压。

如下红色区域有误,请重新填写。

    你的回复:

    请 登录 后回复。还没有在Zeuux哲思注册吗?现在 注册 !

    Zeuux © 2024

    京ICP备05028076号