MeteoInfoLab脚本示例:计算垂直螺旋度

尝试编写MeteoInfoLab脚本计算垂直螺旋度,结果未经验证。

脚本程序:

print ‘Open data files...‘
f_uwnd = addfile(‘D:/Temp/nc/uwnd.2011.nc‘)
f_vwnd = addfile(‘D:/Temp/nc/vwnd.2011.nc‘)
f_omega = addfile(‘D:/Temp/nc/omega.2011.nc‘)

print ‘Calculate vertical helicity...‘
tidx = 173    # Jun 23, 2011
t = f_uwnd.gettime(tidx)
level = [1000, 100]
lat = [15,55]
lon = [70,135]
uwnd = f_uwnd[‘uwnd‘][tidx,level,lat,lon][:,::-1,:]
vwnd = f_vwnd[‘vwnd‘][tidx,level,lat,lon][:,::-1,:]
omega = f_omega[‘omega‘][tidx,level,lat,lon][:,::-1,:]
wd = hcurl(uwnd, vwnd)
lx = -(wd*omega*10.)/12.64*1e6
lx1 = lx[:,[40],:]
lev1 = lx1.dimvalue(0)
#lev2 = 1000 - lev1
lev2 = p2h(lev1)
levels = []
for i in range(0, len(lev1)):
    levels.append(‘%i‘ % lev1[i])
lx1.setdimvalue(0, lev2)

print ‘Plot...‘
layer = contourf(lx1, 20)
title(‘Vertical helicity (‘ + t.strftime(‘%Y-%m-%d‘) + ‘)‘)
yticks(lev2, levels)
xlabel(‘Longitude‘)
ylabel(‘Pressure (hPa)‘)
colorbar(layer)

时间: 2024-10-07 04:11:01

MeteoInfoLab脚本示例:计算垂直螺旋度的相关文章

MeteoInfoLab脚本示例:计算水平螺旋度

尝试了用MeteoInfoLab编写计算水平螺旋度的脚本,结果未经验证. 脚本程序: print 'Open data files...' f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc') f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc') print 'Calculate average wind field from 850 to 600 hpa...' tidx = 173 # Jun 23, 2011 t = f_

MeteoInfoLab脚本示例:FY-3C全球火点HDF数据

FY-3C全球火点HDF数据包含一个FIRES二维变量,第一维是火点数,第二维是一些属性,其中第3.4列分别是火点的纬度和经度.下面的脚本示例读出所有火点经纬度并绘图. 脚本程序: #Add data file fn = 'D:/Temp/hdf/FY3C_VIRRX_GBAL_L2_GFR_MLT_GLL_20150811_POAD_1000M_MS.HDF' f = addfile(fn) #Get data variable v = f['FIRES'] #Get data array d

MeteoInfoLab脚本示例:计算温度平流

需要温度和风场U/V分量格点数据,计算中主要用到cdiff函数,结果正确与否未经验证(希望有人用GrADS验证一下). 脚本程序: print 'Open data files...' f_air = addfile('D:/Temp/nc/air.2011.nc') f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc') f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc') print 'Read data array...'

MeteoInfoLab脚本示例:利用比湿、温度计算相对湿度

利用比湿和温度计算相对湿度的函数是qair2rh(qair, temp, press=1013.25),三个参数分别是比湿.温度和气压,气压有一个缺省值1013.25,因此计算地面相对湿度的时候也可以不给气压参数.下面的例子打开CFS模式的地面结果,并从比湿和温度数组计算出相对湿度数组并绘图. 脚本程序: #Open and read data fn = 'U:/data/climate/CFS/20150712/flxf2015081500.01.2015071200.grb2' f = ad

MeteoInfoLab脚本示例:计算涡度、散度

用U/V分量数据计算涡度和散度,计算涡度的函数是hcurl,计算散度的函数是hdivg,参数都是U, V. 脚本程序: f = addfile('D:/Temp/GrADS/model.ctl') u = f['U'][0,0,:,:] v = f['V'][0,0,:,:] vort = hcurl(u, v) divg = hdivg(u, v) axesm() mlayer = shaperead('D:/Temp/map/country1.shp') geoshow(mlayer) #l

MeteoInfoLab脚本示例:创建netCDF文件(合并文件)

在MeteoInfoLab中增加了创建netCDF文件并写入数据的功能,这里利用合并多个netCDF文件为一个新的netCDF文件为例. 1.创建一个可写入的netCDF文件对象(下面用ncfile表示),用addfile函数,第一个参数是文件名,第二次参数'c'表示创建新的netCDF文件.ncfile = addfile(outfn, 'c') 2.添加维(Dimensions),用ncfile的adddim函数,两个参数分别是维名称和维长度.stn = 26564stdim = ncfil

MeteoInfoLab脚本示例:合并数组

对于全球数据来说,经度要么是-180 - 180,要么是0 - 360,都会存在边界数据不连续的问题.比如0 - 360的数据,怎么得到 -20 - 30度的连续格点数据就是个问题(跨越了数据的经度边界),在MeteoInfoLab中可以用DimArray或MIArray的join函数来将两个数组合并为一个,参数分别是另一个数组和合并的维的序号,比如下面例子中的二维数组(y, x),经度维是第二维(序号为1),也就是说两个数组沿着经度维合并.还需要给合并后的数组的经度维设置新的经度值. 脚本程序

MeteoInfoLab脚本示例:AMSR-E卫星数据投影

AMSR-E(http://nsidc.org/data/amsre/index.html)数据中的Land3数据是HDF-EOS4格式,投影是Cylindrical_Equal_Area.这里示例读取数据并投影至等经纬度投影. 脚本程序: #Add data file fn = 'AMSR_E_L3_DailyLand_V06_20091231.hdf' f = addfile(os.path.join('D:/Temp/hdf', fn)) #vname = 'D_Soil_Moisture

MeteoInfoLab脚本示例:多Y轴图

数据范围相差比较大的数据序列进行对比的时候多Y轴图就很重要了.MeteoInfoLab中提供了一个twinx函数来根据已有的坐标系(Axes)生成一个新的Axes,这个命令会使得已有的Axes不绘制右边的Y轴,而新生成的Axes只绘制右边的Y轴.yaxis函数可以对某个Axes的Y轴进行相应的设置,第一个参数是某个Axes的对象,可以设置Y轴的颜色(color)和偏移(shift),对于超过2个Y轴的图形来说第3个Y轴就必须设置偏移以避免压盖. 双Y轴图: 多Y轴图: 脚本程序: ax1 = a