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 © 2024
京ICP备05028076号