沙尘模式地球坐标系三维动图

今年沙尘暴目前发生比较频繁，很多沙尘暴模式能够比较好地模拟预报沙尘浓度的时空变化，这里利用中国气象局的沙尘暴模式 CUACE/Dust 预报结果绘制一个地球坐标系中的三维动图。在MeteoInfoLab中，axes3d(projection='earth') 创建一个三维地球坐标系，缺省绘制一个偏暗色的地形纹理图像，可以加上 image 参数指定地球表面的纹理图像（在MeteoInfo的map文件夹中有几个jpg文件可用）。模式预报结果是一个netCDF文件，读取沙尘浓度和风场U/V分量的三维数组。沙尘浓度可以用等值面（isosurface）来表示，这里用了两个浓度：100微克每立方米和500微克每立方米，高浓度用更深的颜色和更小的透明度可以展示出浓度的层次。风场用streamslice函数绘制850 hPa风场的流线图，geoshow函数绘制地图行政界线。

单个时次的示例脚本程序如下：

#Set date
sdate = datetime.datetime(2023, 4, 9, 0)

#Set directory
datadir = r'V:\Asian-dust\model\CMA'
datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))
#datadir = r'D:\Temp\cuace\20221212'

#Read data
fn = os.path.join(datadir, 'CUACE-DUST_CMA_{}.nc'.format(sdate.strftime('%Y%m%d%H')))
f = addfile(fn)
st = f.gettime(0)
t = 3
dust = f['CONC_DUST'][t]
u = f['U'][t]
v = f['V'][t]
speed = sqrt(u*u + v*v)
levels = dust.dimvalue(0)
#dust[dust<5] = 0
height = meteolib.pressure_to_height_std(levels)
height = height / 10
lat = dust.dimvalue(1)
lon = dust.dimvalue(2)

#Plot
ax = axes3d(projection='earth', image='shadedrelief.jpg')
lighting(True)
ax.set_rotation(160)
ax.set_elevation(-50)
ax.set_head(0)
ax.set_pitch(-61.5)

geoshow('cn_province', edgecolor='gray', lighting=False)
geoshow('country', edgecolor='gray', lighting=False)
#Beijing location
plot3([116.39,116.39], [39.91,39.91], [0,1200], color='w', lighting=False)
#ax.add_zaxis(116.39, 39.91)
#streamslice
ss = slice(None, None, 4)
levs = arange(2, 20, 2)
hidx = 4    # 850hPa
#hidx = 11   # 500hPa
gss = streamslice(lon[ss], lat[ss], height, u[:,ss,ss], v[:,ss,ss], u[:,ss,ss],
    speed[:,ss,ss], levs=levs, zslice=[height[hidx]], interval=10)
colorbar(gss[0], aspect=30, tickcolor='w', xshift=150, label='m/s')
#isosurface
isosurface(lon, lat, height, dust, 500, facecolor='r', edgecolor=None, 
    alpha=0.6, nthread=4)
isosurface(lon, lat, height, dust, 100, facecolor=[255,180,0], 
    edgecolor=None, alpha=0.4, nthread=4)

v = 1200
axis([-v,v,-v,v,-v,v])

tt = st + datetime.timedelta(hours=t*3)
title('Dust concentration higher than 100/500 ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')),
    color='w')
#savefig('D:/Temp/test/dust_3d_{}.png'.format(tt.strftime('%Y%m%d%H')), dpi=150)

读取多个时次的数据可以绘制时间动画 gif 图，主要用到 gifanimation, gifaddframe 等函数。

 

#Set date
sdate = datetime.datetime(2023, 4, 9, 0)

#Set directory
datadir = r'V:\Asian-dust\model\CMA'
datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))
outdir = 'V:/Asian-dust/figures'
outdir = os.path.join(outdir, 'CMA', sdate.strftime('%Y%m%d'))
if not os.path.exists(outdir):
    os.makedirs(outdir)

#Read data
fn = os.path.join(datadir, 'CUACE-DUST_CMA_{}.nc'.format(sdate.strftime('%Y%m%d%H')))
f = addfile(fn)
st = f.gettime(0)
t = 10
dust = f['CONC_DUST'][t]
levels = dust.dimvalue(0)
#dust[dust<5] = 0
height = meteolib.pressure_to_height_std(levels)
height = height / 10
lat = dust.dimvalue(1)
lon = dust.dimvalue(2)

#Plot
ax = axes3d(projection='earth', image='shadedrelief.jpg')
lighting(True)
ax.set_rotation(160)
ax.set_elevation(-50)
ax.set_head(0)
ax.set_pitch(-61.5)

geoshow('cn_province', edgecolor='gray', lighting=False)
geoshow('country', edgecolor='gray', lighting=False)
#Beijing location
plot3([116.39,116.39], [39.91,39.91], [0,1200], color='w', lighting=False)
#streamslice
ss = slice(None, None, 4)
levs = arange(2, 20, 2)

#Loop
giffn = os.path.join(outdir, 'Dust_3d_isosurface_relief_{}--loop_1.gif'.format(st.strftime('%Y%m%d%H')))
print giffn
animation = gifanimation(giffn)
tn = f.timenum()
tn = 24
for t in range(1, tn):    
    print(t)
    if t > 1:
        cll()
        cll()
        cll()
    
    dust = f['CONC_DUST'][t]
    u = f['U'][t]
    v = f['V'][t]
    speed = sqrt(u*u + v*v)
    
    gss = streamslice(lon[ss], lat[ss], height, u[:,ss,ss], v[:,ss,ss], u[:,ss,ss],
        speed[:,ss,ss], levs=levs, zslice=[height[4]], interval=10)
    colorbar(gss[0], aspect=30, tickcolor='w', xshift=135, label='m/s')
    #isosurface
    isosurface(lon, lat, height, dust, 500, facecolor='r', edgecolor=None, 
        alpha=0.6, nthread=4)
    isosurface(lon, lat, height, dust, 100, facecolor=[255,180,0], 
        edgecolor=None, alpha=0.4, nthread=4)

    vv = 1200
    axis([-vv,vv,-vv,vv,-vv,vv])
    
    tt = f.gettime(t)
    title('Dust concentration higher than 100/500 ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')),
        color='w')
    #Add frame to gif animation
    gifaddframe(animation, width=1095, height=647)

#Finish gif animation
animation.finish()
print 'Finished...'