MeteoInfoLab脚本示例：绘制污染物三维粒子图

MeteoInfo · 发表于 2020-5-28 08:33:17

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

Axes3DGL类中包含了 plot_particles 方法来绘制三维粒子图，适合用来绘制污染物浓度三维分布。

绘制沙尘浓度三维分布：

#Set date
sdate = datetime.datetime(2020, 5, 11, 0)
#Set directory
datadir = r'W:\Asia-dust\CMA'
datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))
#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 = 4
dust = f['CONC_DUST'][t]
levels = dust.dimvalue(0)
#dust[dust<5] = 0
height = meteolib.pressure_to_height_std(levels)
lat = dust.dimvalue(1)
lon = dust.dimvalue(2)
#Relief data
rfn = 'D:/Temp/nc/elev.0.25-deg.nc'
rf = addfile(rfn)
elev = rf['data'][0,'15:65','65:155']
elev[elev<0] = -1
lon1 = elev.dimvalue(1)
lat1 = elev.dimvalue(0)
lon1, lat1 = meshgrid(lon1, lat1)
#Map
lchina = shaperead('cn_province')
clon = lchina.x_coord
clat = lchina.y_coord
calt = zeros(len(clon))
h = interp2d(elev, clon, clat)
calt = calt + h
lworld = shaperead('country')
wlon = lworld.x_coord
wlat = lworld.y_coord
walt = zeros(len(wlon))
h = interp2d(elev, wlon, wlat)
walt = walt + h
#Plot
ax = axes3dgl(bbox=True)
ax.set_elevation(-20)
ax.set_rotation(335)
rlevs = arange(0, 6000, 200)
cols = makecolors(len(rlevs) + 1, cmap='MPL_gist_yarg', alpha=1)
cols[0] = [51,153,255]
ax.plot_surface(lon1, lat1, elev, rlevs, colors=cols, edge=False)
ax.plot(clon, clat, calt, color=[255,153,255])
ax.plot(wlon, wlat, walt, color='b')
#Beijing location
ax.plot([116.39,116.39], [39.91,39.91], [0,12000])
#ax.set_lighting(True, position=[1,1,1,1], mat_specular=[0.5,0.5,0.5,1])
levs = [50,100,200,300,400,500]
#levs = [100,200,300,400,500]
cmap='WhiteBlueGreenYellowRed'
#cmap = 'MPL_Oranges'
pp = ax.plot_particles(lon, lat, height, dust, levs, vmin=20, s=2, \
cmap=cmap, alpha_min=0.1, alpha_max=0.7, density=1)
colorbar(pp, aspect=30)
xlim(65, 155)
xlabel('Longitude')
ylim(15, 65)
ylabel('Latitude')
zlim(0, 12000)
zlabel('Height (m)')
tt = st + datetime.timedelta(hours=t*3)
title('Dust concentration ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')))
#savefig('D:/Temp/test/dust_3d_{}.png'.format(tt.strftime('%Y%m%d%H')), dpi=150)

#Set date 
sdate = datetime.datetime(2020, 5, 11, 0) 
 
#Set directory 
datadir = r'W:\Asia-dust\CMA' 
datadir = os.path.join(datadir, sdate.strftime('%Y%m%d')) 
 
#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 = 4 
dust = f['CONC_DUST'][t] 
levels = dust.dimvalue(0) 
#dust[dust<5] = 0 
height = meteolib.pressure_to_height_std(levels) 
lat = dust.dimvalue(1) 
lon = dust.dimvalue(2) 
 
#Relief data 
rfn = 'D:/Temp/nc/elev.0.25-deg.nc' 
rf = addfile(rfn) 
elev = rf['data'][0,'15:65','65:155'] 
elev[elev<0] = -1 
lon1 = elev.dimvalue(1) 
lat1 = elev.dimvalue(0) 
lon1, lat1 = meshgrid(lon1, lat1) 
 
#Map 
lchina = shaperead('cn_province') 
clon = lchina.x_coord 
clat = lchina.y_coord 
calt = zeros(len(clon)) 
h = interp2d(elev, clon, clat) 
calt = calt + h 
lworld = shaperead('country') 
wlon = lworld.x_coord 
wlat = lworld.y_coord 
walt = zeros(len(wlon)) 
h = interp2d(elev, wlon, wlat) 
walt = walt + h 
 
#Plot 
ax = axes3dgl(bbox=True) 
ax.set_elevation(-20) 
ax.set_rotation(335) 
rlevs = arange(0, 6000, 200) 
cols = makecolors(len(rlevs) + 1, cmap='MPL_gist_yarg', alpha=1) 
cols[0] = [51,153,255] 
ax.plot_surface(lon1, lat1, elev, rlevs, colors=cols, edge=False) 
ax.plot(clon, clat, calt, color=[255,153,255]) 
ax.plot(wlon, wlat, walt, color='b') 
#Beijing location 
ax.plot([116.39,116.39], [39.91,39.91], [0,12000]) 
#ax.set_lighting(True, position=[1,1,1,1], mat_specular=[0.5,0.5,0.5,1]) 
levs = [50,100,200,300,400,500] 
#levs = [100,200,300,400,500] 
cmap='WhiteBlueGreenYellowRed' 
#cmap = 'MPL_Oranges' 
pp = ax.plot_particles(lon, lat, height, dust, levs, vmin=20, s=2, \ 
    cmap=cmap, alpha_min=0.1, alpha_max=0.7, density=1) 
colorbar(pp, aspect=30) 
xlim(65, 155) 
xlabel('Longitude') 
ylim(15, 65) 
ylabel('Latitude') 
zlim(0, 12000) 
zlabel('Height (m)') 
tt = st + datetime.timedelta(hours=t*3) 
title('Dust concentration ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00'))) 
#savefig('D:/Temp/test/dust_3d_{}.png'.format(tt.strftime('%Y%m%d%H')), dpi=150)

输出沙尘三维浓度随时间变化的gif动画图：

import time
#Set date
sdate = datetime.datetime(2020, 5, 10, 0)
#Set directory
datadir = r'W:\Asia-dust\CMA'
datadir = os.path.join(datadir, sdate.strftime('%Y%m%d'))
#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 = 0
dust = f['CONC_DUST'][t]
levels = dust.dimvalue(0)
#dust[dust<5] = 0
height = meteolib.pressure_to_height_std(levels)
lat = dust.dimvalue(1)
lon = dust.dimvalue(2)
#Relief data
rfn = 'D:/Temp/nc/elev.0.25-deg.nc'
rf = addfile(rfn)
elev = rf['data'][0,'15:65','65:155']
elev[elev<0] = -1
lon1 = elev.dimvalue(1)
lat1 = elev.dimvalue(0)
lon1, lat1 = meshgrid(lon1, lat1)
#Map
lchina = shaperead('cn_province')
clon = lchina.x_coord
clat = lchina.y_coord
calt = zeros(len(clon))
h = interp2d(elev, clon, clat)
calt = calt + h
lworld = shaperead('country')
wlon = lworld.x_coord
wlat = lworld.y_coord
walt = zeros(len(wlon))
h = interp2d(elev, wlon, wlat)
walt = walt + h
#Plot
ax = axes3dgl(bbox=True)
ax.set_elevation(-20)
ax.set_rotation(335)
rlevs = arange(0, 6000, 200)
cols = makecolors(len(rlevs) + 1, cmap='MPL_gist_yarg', alpha=1)
cols[0] = [51,153,255]
ax.plot_surface(lon1, lat1, elev, rlevs, colors=cols, edge=False)
ax.plot(clon, clat, calt, color=[255,153,255])
ax.plot(wlon, wlat, walt, color='b')
#Beijing location
ax.plot([116.39,116.39], [39.91,39.91], [0,12000])
#ax.set_lighting(True, position=[1,1,1,1], mat_specular=[0.5,0.5,0.5,1])
levs = [50,100,150,200,300,400,500,1000]
cmap='WhiteBlueGreenYellowRed'
#cmap = 'MPL_Oranges'
#Loop
giffn = os.path.join(datadir, 'Dust_3d_particles_relief_{}--loop-1.gif'.format(st.strftime('%Y%m%d')))
print giffn
animation = gifanimation(giffn)
tn = f.timenum()
tn = 24
for t in range(1, tn):
print t
if t > 1:
cll()
dust = f['CONC_DUST'][t]
pp = ax.plot_particles(lon, lat, height, dust, levs, vmin=20, s=2, \
cmap=cmap, alpha_min=0.1, alpha_max=0.7, density=1)
colorbar(pp, newlegend=False, aspect=30)
xlim(65, 155)
xlabel('Longitude')
ylim(15, 65)
ylabel('Latitude')
zlim(0, 12000)
zlabel('Height (m)')
#zticks(arange(len(levels))[1:], levels[1:])
tt = f.gettime(t)
title('Dust concentration ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00')))
#Add frame to gif animation
gifaddframe(animation, width=1024, height=516)
#Finish gif animation
animation.finish()
print 'Finished...'

import time 
 
#Set date 
sdate = datetime.datetime(2020, 5, 10, 0) 
 
#Set directory 
datadir = r'W:\Asia-dust\CMA' 
datadir = os.path.join(datadir, sdate.strftime('%Y%m%d')) 
 
#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 = 0 
dust = f['CONC_DUST'][t] 
levels = dust.dimvalue(0) 
#dust[dust<5] = 0 
height = meteolib.pressure_to_height_std(levels) 
lat = dust.dimvalue(1) 
lon = dust.dimvalue(2) 
 
#Relief data 
rfn = 'D:/Temp/nc/elev.0.25-deg.nc' 
rf = addfile(rfn) 
elev = rf['data'][0,'15:65','65:155'] 
elev[elev<0] = -1 
lon1 = elev.dimvalue(1) 
lat1 = elev.dimvalue(0) 
lon1, lat1 = meshgrid(lon1, lat1) 
 
#Map 
lchina = shaperead('cn_province') 
clon = lchina.x_coord 
clat = lchina.y_coord 
calt = zeros(len(clon)) 
h = interp2d(elev, clon, clat) 
calt = calt + h 
lworld = shaperead('country') 
wlon = lworld.x_coord 
wlat = lworld.y_coord 
walt = zeros(len(wlon)) 
h = interp2d(elev, wlon, wlat) 
walt = walt + h 
 
#Plot 
ax = axes3dgl(bbox=True) 
ax.set_elevation(-20) 
ax.set_rotation(335) 
rlevs = arange(0, 6000, 200) 
cols = makecolors(len(rlevs) + 1, cmap='MPL_gist_yarg', alpha=1) 
cols[0] = [51,153,255] 
ax.plot_surface(lon1, lat1, elev, rlevs, colors=cols, edge=False) 
ax.plot(clon, clat, calt, color=[255,153,255]) 
ax.plot(wlon, wlat, walt, color='b') 
#Beijing location 
ax.plot([116.39,116.39], [39.91,39.91], [0,12000]) 
#ax.set_lighting(True, position=[1,1,1,1], mat_specular=[0.5,0.5,0.5,1]) 
levs = [50,100,150,200,300,400,500,1000] 
cmap='WhiteBlueGreenYellowRed' 
#cmap = 'MPL_Oranges' 
 
#Loop 
giffn = os.path.join(datadir, 'Dust_3d_particles_relief_{}--loop-1.gif'.format(st.strftime('%Y%m%d'))) 
print giffn 
animation = gifanimation(giffn) 
tn = f.timenum() 
tn = 24 
for t in range(1, tn): 
    print t 
    if t > 1: 
        cll() 
    dust = f['CONC_DUST'][t] 
    pp = ax.plot_particles(lon, lat, height, dust, levs, vmin=20, s=2, \ 
        cmap=cmap, alpha_min=0.1, alpha_max=0.7, density=1) 
    colorbar(pp, newlegend=False, aspect=30) 
    xlim(65, 155) 
    xlabel('Longitude') 
    ylim(15, 65) 
    ylabel('Latitude') 
    zlim(0, 12000) 
    zlabel('Height (m)') 
    #zticks(arange(len(levels))[1:], levels[1:]) 
    tt = f.gettime(t) 
    title('Dust concentration ug/m3 ({}UTC)'.format(tt.strftime('%Y-%m-%d %H:00'))) 
    #Add frame to gif animation 
    gifaddframe(animation, width=1024, height=516) 
 
#Finish gif animation 
animation.finish() 
print 'Finished...'

开心的玉米超人 · 发表于 2020-6-2 09:16:18

王老师，您的数据（CUACE-DUST_CMA_2020051000.nc和elev.0.25-deg.nc）在哪里下载的，谢谢

随缘 · 发表于 2020-5-28 08:58:29

{:eb502:}{:eb502:}

午后阳光 · 发表于 2020-5-29 10:03:01

厉害，好漂亮的图，喜欢了

kc121 · 发表于 2020-6-6 19:29:13

这图画的太好了！

guo066600 · 发表于 2020-7-18 19:35:00

开心的玉米超人发表于 2020-6-2 09:16
王老师，您的数据（CUACE-DUST_CMA_2020051000.nc和elev.0.25-deg.nc）在哪里下载的，谢谢

请问您找到CUACE-DUST_CMA_2020051000.nc和elev.0.25-deg.nc的下载地址了吗？如果找到的话，能分享一下吗？

handuoduo0215 · 发表于 2020-11-4 20:41:47

图画的好炫酷啊，需要什么数据，才能画出来呢？

森南 · 发表于 2020-11-9 16:09:17

所以数据大家有找到吗

突破的渴望 · 发表于 2020-12-5 11:07:39

王老师，没有数据啊，怎么画？提供下数据下载链接啊，要不这个帖子会沉下去的

陈浩 · 发表于 2020-12-17 10:37:47

王老师，您的数据（CUACE-DUST_CMA_2020051000.nc和elev.0.25-deg.nc）在哪里下载的，谢谢

		自动登录	找回密码
密码			立即注册

MeteoInfoLab脚本示例：绘制污染物三维粒子图

登录后查看更多精彩内容~

评分