请教用MI做时间剖面的方法

lhy277969385

今天看了一天示例文件，
但剖面好像都是经度为横坐标
编程试了好久，都没办法实现时间为横坐标做数据的剖面……
王老师让我参考：http://www.meteothink.org/docs/m ... etag=s_pctim_aiomsg
但看了好久，这里边的图也都是经度为横坐标的图
用contour()函数的时候，只要经度固定住，contour()函数就报错啊……
唉……不知道该怎么用contour函数了……
然后又发现MeteoInfoMap可以直接画出时间剖面图！
但是……时间不可选择……
然后又去试着用MI截取nc文件的时间段……感觉跑偏了……作为一个新手我表示很头大
求大佬帮忙解决下吧……跪谢……
想处理的文件就是NCEP的每日四次的这种数据
f
File Name: D:\Meteoinfolab_Project\readnc\omega.2019.nc


File type: Hierarchical Data Format, version 5 (NetCDF-4)


Dimensions: 4


        level = 12;


        lat = 73;


        lon = 144;


        time = 232;


X Dimension: Xmin = 0.0; Xmax = 357.5; Xsize = 144; Xdelta = 2.5


Y Dimension: Ymin = -90.0; Ymax = 90.0; Ysize = 73; Ydelta = 2.5


Global Attributes:


        : Conventions = "COARDS"


        : title = "4x daily NMC reanalysis (2014)"


        : history = "created 2013/12 by Hoop (netCDF2.3)"


        : description = "Data is from NMC initialized reanalysis\n(4x/day).  It consists of most variables interpolated to\npressure surfaces from model (sigma) surfaces."


        : platform = "Model"


        : References = "http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html"


        : dataset_title = "NCEP-NCAR Reanalysis 1"


Variations: 5


        float level(level);


                level: units = "millibar"


                level: actual_range = 1000.0f, 100.0f


                level: long_name = "Level"


                level: positive = "down"


                level: GRIB_id = 100S


                level: GRIB_name = "hPa"


                level: axis = "Z"


        float lat(lat);


                lat: units = "degrees_north"


                lat: actual_range = 90.0f, -90.0f


                lat: long_name = "Latitude"


                lat: standard_name = "latitude"


                lat: axis = "Y"


        float lon(lon);


                lon: units = "degrees_east"


                lon: long_name = "Longitude"


                lon: actual_range = 0.0f, 357.5f


                lon: standard_name = "longitude"


                lon: axis = "X"


        double time(time);


                time: long_name = "Time"


                time: delta_t = "0000-00-00 06:00:00"


                time: standard_name = "time"


                time: axis = "T"


                time: units = "hours since 1800-01-01 00:00:0.0"


                time: actual_range = 1919712.0, 1921098.0


                time: _ChunkSizes = 1


        float omega(time,level,lat,lon);


                omega: long_name = "4xDaily Omega"


                omega: units = "Pascal/s"


                omega: precision = 3S


                omega: least_significant_digit = 3S


                omega: GRIB_id = 39S


                omega: GRIB_name = "VVEL"


                omega: var_desc = "Omega (dp/dt)"


                omega: level_desc = "Multiple levels"


                omega: statistic = "Individual Obs"


                omega: parent_stat = "Other"


                omega: missing_value = -9.96921E36f


                omega: valid_range = -3.0f, 7.0f


                omega: dataset = "NCEP Reanalysis"


                omega: actual_range = -1.7110002f, 1.173f


                omega: _ChunkSizes = 1, 1, 73, 144


直接就能下载
ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/pressure/omega.2019.nc
再次跪谢各位大佬！求指教！

lhy277969385

贴上代码……

1. clf()
   
2. figure(figsize=[700,550], newfig=False)
   
3. f = addfile(r'D:\Meteoinfolab_Project\readnc\omega.2019.nc')
   
4. omega = f['omega']#[:,0:5,'50','120']
   
5. time = f['time'][0:8]
   
6. lon = f['lon'][:]
   
7. lat = f['lat'][:]
   
8. level = f['level']
   
9. 
   
10. layer = contour(omega)
    
11. clabel(layer)

复制代码

到这里就一直报错不知道该怎么改下去了……
报错信息：
>>> run script...
Traceback (most recent call last):
  File "D:\Meteoinfolab_Project\readnc\????.py", line 14, in <module>
    layer = contour(omega)
  File "D:\MeteoInfo\pylib\mipylib\plotlib\miplot.py", line 2063, in contour
    r = gca.contour(*args, **kwargs)
  File "D:\MeteoInfo\pylib\mipylib\plotlib\axes.py", line 1841, in contour
    ls = LegendManage.createLegendScheme(gdata.min(), gdata.max(), cmap)
AttributeError: 'NoneType' object has no attribute 'min'

MeteoInfo

lhy277969385 发表于 2019-3-2 19:37
贴上代码……

代码没有贴全

lhy277969385

不好意思王老师~我补全~还请大佬指教！

1. clf()
   
2. figure(figsize=[700,550], newfig=False)
   
3. f = addfile(r'D:\Meteoinfolab_Project\readnc\omega.2019.nc')
   
4. data = f['omega']#[0:5,:,'50','120']
   
5. lev1 = tdata.dimvalue(0)
   
6. lev1=lev1[::-1]
   
7. lev2 = meteo.pressure_to_height_std(lev1)
   
8. lev2 = lev2[:]/1000
   
9. tdata.setdimvalue(0, lev2)
   
10. 
    
11. layer = contour(data,20)
    
12. clabel(layer)
    
13. colorbar(layer)
    
14. yticks(lev2, levels)
    
15. xlabel('Time')
    
16. ylabel('Pressure (hPa)')
    
17. title('Time Profile')

复制代码

报错提示：
>>> run script...
Traceback (most recent call last):
  File "D:\Meteoinfolab_Project\readnc\????.py", line 11, in <module>
    layer = contour(data,20)
  File "D:\MeteoInfo\pylib\mipylib\plotlib\miplot.py", line 2063, in contour
    r = gca.contour(*args, **kwargs)
  File "D:\MeteoInfo\pylib\mipylib\plotlib\axes.py", line 1835, in contour
    ls = LegendManage.createLegendScheme(gdata.min(), gdata.max(), cn, cmap)
AttributeError: 'NoneType' object has no attribute 'min'

MeteoInfo

lhy277969385 发表于 2019-3-4 16:09
不好意思王老师~我补全~还请大佬指教！

参考此脚本：

1. fn = 'D:/Temp/nc/omega.2011.nc'
   
2. f = addfile(fn)
   
3. data = f['omega'][0:5,:,'50','120']
   
4. data = data.T
   
5. lev1 = data.dimvalue(0)
   
6. lev2 = meteo.pressure_to_height_std(lev1)
   
7. lev2 = lev2[:]/1000
   
8. data.setdimvalue(0, lev2)
   
9. 
   
10. layer = contour(data,20)
    
11. clabel(layer)
    
12. colorbar(layer)
    
13. yticks(lev2, lev1)
    
14. xlabel('Time')
    
15. ylabel('Pressure (hPa)')
    
16. title('Time Profile')

fn = 'D:/Temp/nc/omega.2011.nc'<br /> f = addfile(fn)<br /> data = f['omega'][0:5,:,'50','120']<br /> data = data.T<br /> lev1 = data.dimvalue(0)<br /> lev2 = meteo.pressure_to_height_std(lev1)<br /> lev2 = lev2[:]/1000<br /> data.setdimvalue(0, lev2) <br /> <br /> layer = contour(data,20)<br /> clabel(layer)<br /> colorbar(layer)<br /> yticks(lev2, lev1)<br /> xlabel('Time')<br /> ylabel('Pressure (hPa)')<br /> title('Time Profile')

lhy277969385

MeteoInfo 发表于 2019-3-4 22:55
参考此脚本：

太感谢王老师了！！！

920daibingbing

老师您好，我想请教下上面这个程序，如果想画水汽通量应该怎么改，才学习不久，请老师指教，谢谢

MeteoInfo

920daibingbing 发表于 2020-8-2 22:31
老师您好，我想请教下上面这个程序，如果想画水汽通量应该怎么改，才学习不久，请老师指教，谢谢

可以看看脚本汇总帖：http://bbs.06climate.com/forum.p ... &extra=page%3D1

920daibingbing

MeteoInfo 发表于 2020-8-3 08:57
可以看看脚本汇总帖：http://bbs.06climate.com/forum.php?mod=viewthread&tid=36151&extra=page%3D1

老师你好，我参考您的脚本，会画某一层的水汽通量散度，也会画单一变量的垂直剖面，就是不知道怎么把两个合起来，画水汽通量散度的垂直剖面，因为公式里有一个levs，您有时间可以回复下我吗，谢谢！

920daibingbing

老师您好，我已经画出图了，但是想请教时间坐标轴怎么控制，我设的时间范围并没有作用，还有就是坐标轴标记怎么改，谢谢。
fn = 'D:/nc/4/uwnd.2019.nc'
fm = 'D:/nc/4/vwnd.2019.nc'
fo = 'D:/nc/4/air.2019.nc'
fp = 'D:/nc/4/rhum.2019.nc'
f1 = addfile(fn)
f2 = addfile(fm)
f3 = addfile(fo)
f4 = addfile(fp)
data1 = f1['uwnd'][872:886,2:8,'46','100']
data2 = f2['vwnd'][872:886,2:8,'46','100']
data3 = f3['air'][872:886,2:8,'46','100']
data4 = f4['rhum'][872:886,2:8,'46','100']
data1 = data1.T
data2 = data2.T
data3 = data3.T
data4 = data4.T
lev1 = data1.dimvalue(0)
lev2 = meteo.pressure_to_height_std(lev1)
lev2 = lev2[:]/1000
data1.setdimvalue(0, lev2)
data2.setdimvalue(0, lev2)
data3.setdimvalue(0, lev2)
data4.setdimvalue(0, lev2)
es = 6.112*exp(17.67*(data3-273.16)/(data3-29.65))
qv = data4*es/100
uv = magnitude(data1, data2)
data = uv*qv/(9.8*1000)
levs=arange(0,100,5)
layer = contour(data*10000,levs,color='black')
clabel(layer)
#xticks([07/00,07/06,07/12,07/18,08/00,08/06,08/12,08/18,09/00,09/06,09/12,09/18,10/00,10/06,10/12,10/18])
yticks(lev2, lev1)
xlabel('Time')
ylabel('Pressure (hPa)')