MeteoInfoLab脚本示例：计算温度平流(附全球七大洲边界地图shp)

本文源码作者：王亚强 研究员 中国气象科学研究院

---

温度和风场U/V分量格点数据，

计算中主要用到cdiff函数（结果已用GrADS验证一致）。

print 'Open data files...'
f_air = addfile('F:/RMeteoInfo/data/air.2011.nc')
f_uwnd = addfile('F:/RMeteoInfo/data/uwnd.2011.nc')
f_vwnd = addfile('F:/RMeteoInfo/data/vwnd.2011.nc')
print 'Read data array...'
tidx = 173    # Jun 23, 2011
t = f_air.gettime(tidx)
lidx = 0    # 1000 hPa
air = f_air['air'][tidx,lidx,:,:]
uwnd = f_uwnd['uwnd'][tidx,lidx,:,:]
vwnd = f_vwnd['vwnd'][tidx,lidx,:,:]
lon = f_air['lon'][:]
lat = f_air['lat'][:]
lon, lat = meshgrid(lon, lat)
# Calculate
print 'Calculate...'
dtx = cdiff(air, 1)
dty = cdiff(air, 0)
dx = cdiff(lon, 1) * pi / 180
dy = cdiff(lat, 0) * pi / 180
tadv = -1*((uwnd*dtx)/(cos(lat*pi/180)*dx)+vwnd*dty/dy)/6.37e6
figure(figsize=[700,400], newfig=False)
ax=axesm(tickfontsize=11,axis=True)
# axesm()
world = shaperead('F:/RMeteoInfo/data/map/country1.shp') 
lchina = shaperead('F:/RMeteoInfo/data/map/bou2_4p.shp')
lchina2 = shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp')
#geoshow(world)
geoshow(lchina,edgecolor='k',size=0.3)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.3)
china_layer = geoshow('china', visible=False)
levs = [-1*10e-5,-0.8*10e-5,-0.6*10e-5,-0.4*10e-5,-0.2*10e-5,0,0.2*10e-5,0.4*10e-5,0.6*10e-5,0.8*10e-5,1*10e-5]
layer = contourfm(tadv, levs, cmap='grads_rainbow')
masklayer(lchina, [layer])
title(u'Temperature advection (' + t.strftime('%Y-%m-%d') + ')',bold=True,fontsize=11)
colorbar(layer,orientation='horizontal',ticklen=0,extendrect=False,shrink=1,pad=0.01,aspect=50,label='mm',bold=True,fontsize=11)
yaxis(tickvisible=True,location='left',tickwidth=2,linewidth=2,ticklength=3) #ticklength刻度线长度，tickwidth刻度线宽度，linewidth边框宽度
yaxis(tickvisible=False,location='right',tickwidth=2,linewidth=2,ticklength=4) #分别调试tick的宽度，边框线宽和tick的长度
xaxis(tickvisible=False,location='top',tickwidth=2,linewidth=2,ticklength=4)
xaxis(tickvisible=True,location='bottom',tickwidth=2,linewidth=2,ticklength=3)
xlim(70, 140)
ylim(15, 55)
xticks(arange(70, 140.1, 10),bold=True,fontsize=11)
yticks(arange(15, 55.1, 10),bold=True,fontsize=11)
xlabel('Longitude',bold=True,fontsize=11)
ylabel('Latitude',bold=True,fontsize=11)
ax.scale_bar(0.2,0.35,width=165,linewidth=1,bold=True,fontsize=10,bartype='alternating_bar')
ax.north_arrow(0.2,0.9,width=40,height=40,linewidth=1,bold=True,fontsize=10)
#Add south China Sea
axesm(position=[0.7305,0.238,0.15,0.2], axison=False, frameon=True)
geoshow(lchina,edgecolor='k',size=0.2)
bou1_layer=geoshow(lchina2,edgecolor='k',size=0.1)
layer = contourfm(tadv, levs, cmap='grads_rainbow')
masklayer(lchina, [layer])
xlim(106, 123)
ylim(2, 23)
savefig('F:/RMeteoInfo/test/plot128.png',dpi=800)

---

获取本文全球七大洲边界地图shp数据的途径：

气象水文科研猫公众号后台回复：

“全球七大洲边界地图shp”，

获取百度云免费下载链接。

---