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

MeteoInfo

无尽的兰博 发表于 2016-10-12 16:10
老师您好，我用您的脚本用grib2的数据画了涡度图，但是我想选择200hPa那层的，我的设置如下：

f = addfi ...

你目前的用法是取第200层的数据（应该没有那么多层），数字是序号，如果要取200hPa，试试用[200]。

无尽的兰博

MeteoInfo 发表于 2016-10-12 16:28
你目前的用法是取第200层的数据（应该没有那么多层），数字是序号，如果要取200hPa，试试用[200]。

谢谢老师，但是我没明白是哪里改成[200]呢？
是改成
u = f['u-component_of_wind_pressure_difference_layer'][0,[200],:,:]
v = f['v-component_of_wind_pressure_difference_layer'][0,[200],:,:]
这样吗？

无尽的兰博

老师我还有一个问题，我想在图里加入一条0值的白线，需要怎么设置啊？找了好久都没找到

MeteoInfo

无尽的兰博 发表于 2016-10-12 16:58
老师我还有一个问题，我想在图里加入一条0值的白线，需要怎么设置啊？找了好久都没找到

参考此贴：http://bbs.06climate.com/forum.p ... &extra=page%3D1

自定义等值区范围和颜色。

无尽的兰博

MeteoInfo 发表于 2016-10-12 17:48
参考此贴：http://bbs.06climate.com/forum.php?mod=viewthread&tid=36422&extra=page%3D1

自定义等值 ...

多谢老师，我去学习学习

半颗麦芽糖

王老师，我对这个脚本进行了一点修改想画我模式的涡度图，但是画出来不仅量级有问题，而且图中只显示了一种颜色，但是我用grads使用相同的数据画图却是正确的，所以我不太明白问题在哪里，请王老师指导，附上我的ctl、脚本和图。
模式ctl：
dset E:\gongzuo\3.8.1\summer0801.dat
options  byteswapped
undef 1.e30
title  OUTPUT FROM WRF V3.8.1 MODEL
pdef  174 120 lcc  37.383   90.551   87.500   60.500  60.00000  30.00000   90.55100  30000.000  30000.000
xdef  576 linear   51.60618   0.13513513
ydef  284 linear   16.84722   0.13513513
zdef   19 levels  
1000.00000
950.00000
900.00000
850.00000
800.00000
750.00000
700.00000
650.00000
600.00000
550.00000
500.00000
450.00000
400.00000
350.00000
300.00000
250.00000
200.00000
150.00000
100.00000
tdef  737 linear 00Z01JUN2008     180MN      
VARS    4
U             19  0  x-wind component (m s-1)
V             19  0  y-wind component (m s-1)
pressure      19  0  Model pressure (hPa)
height        19  0  Model height (km)
ENDVARS
@ global String comment TITLE =  OUTPUT FROM WRF V3.8.1 MODEL
@ global String comment START_DATE = 2008-06-01_00:00:00
@ global String comment SIMULATION_START_DATE = 2008-06-01_00:00:00
@ global String comment WEST-EAST_GRID_DIMENSION =   175
@ global String comment SOUTH-NORTH_GRID_DIMENSION =   121
@ global String comment BOTTOM-TOP_GRID_DIMENSION =    30
@ global String comment DX =     30000.00
@ global String comment DY =     30000.00
@ global String comment SKEBS_ON =     0
@ global String comment SPEC_BDY_FINAL_MU =     1
@ global String comment USE_Q_DIABATIC =     0
@ global String comment GRIDTYPE = C
@ global String comment DIFF_OPT =     1
@ global String comment KM_OPT =     4
@ global String comment DAMP_OPT =     0
@ global String comment DAMPCOEF =         0.01
@ global String comment KHDIF =         0.00
@ global String comment KVDIF =         0.00
@ global String comment MP_PHYSICS =     6
@ global String comment RA_LW_PHYSICS =     1
@ global String comment RA_SW_PHYSICS =     1
@ global String comment SF_SFCLAY_PHYSICS =     2
@ global String comment SF_SURFACE_PHYSICS =     2
@ global String comment BL_PBL_PHYSICS =     2
@ global String comment CU_PHYSICS =     1
@ global String comment SF_LAKE_PHYSICS =     0
@ global String comment SURFACE_INPUT_SOURCE =     1
@ global String comment SST_UPDATE =     0
@ global String comment GRID_FDDA =     0
@ global String comment GFDDA_INTERVAL_M =     0
@ global String comment GFDDA_END_H =     0
@ global String comment GRID_SFDDA =     0
@ global String comment SGFDDA_INTERVAL_M =     0
@ global String comment SGFDDA_END_H =     0
@ global String comment HYPSOMETRIC_OPT =     2
@ global String comment USE_THETA_M =     0
@ global String comment SF_URBAN_PHYSICS =     0
@ global String comment SHCU_PHYSICS =     0
@ global String comment MFSHCONV =     0
@ global String comment FEEDBACK =     1
@ global String comment SMOOTH_OPTION =     0
@ global String comment SWRAD_SCAT =         1.00
@ global String comment W_DAMPING =     0
@ global String comment DT =       180.00
@ global String comment RADT =        30.00
@ global String comment BLDT =         0.00
@ global String comment CUDT =         5.00
@ global String comment AER_OPT =     0
@ global String comment SWINT_OPT =     0
@ global String comment AER_TYPE =     1
@ global String comment AER_AOD550_OPT =     1
@ global String comment AER_ANGEXP_OPT =     1
@ global String comment AER_SSA_OPT =     1
@ global String comment AER_ASY_OPT =     1
@ global String comment AER_AOD550_VAL =         0.12
@ global String comment AER_ANGEXP_VAL =         1.30
@ global String comment AER_SSA_VAL =         0.00
@ global String comment AER_ASY_VAL =         0.00
@ global String comment MOIST_ADV_OPT =     1
@ global String comment SCALAR_ADV_OPT =     1
@ global String comment TKE_ADV_OPT =     1
@ global String comment DIFF_6TH_OPT =     0
@ global String comment DIFF_6TH_FACTOR =         0.12
@ global String comment OBS_NUDGE_OPT =     0
@ global String comment BUCKET_MM =        -1.00
@ global String comment BUCKET_J =        -1.00
@ global String comment PREC_ACC_DT =         0.00
@ global String comment SF_OCEAN_PHYSICS =     0
@ global String comment ISFTCFLX =     0
@ global String comment ISHALLOW =     0
@ global String comment ISFFLX =     1
@ global String comment ICLOUD =     1
@ global String comment ICLOUD_CU =     0
@ global String comment TRACER_PBLMIX =     1
@ global String comment SCALAR_PBLMIX =     0
@ global String comment YSU_TOPDOWN_PBLMIX =     0
@ global String comment GRAV_SETTLING =     0
@ global String comment DFI_OPT =     0
@ global String comment SIMULATION_INITIALIZATION_TYPE = REAL-DATA CASE
@ global String comment WEST-EAST_PATCH_START_UNSTAG =     1
@ global String comment WEST-EAST_PATCH_END_UNSTAG =   174
@ global String comment WEST-EAST_PATCH_START_STAG =     1
@ global String comment WEST-EAST_PATCH_END_STAG =   175
@ global String comment SOUTH-NORTH_PATCH_START_UNSTAG =     1
@ global String comment SOUTH-NORTH_PATCH_END_UNSTAG =   120
@ global String comment SOUTH-NORTH_PATCH_START_STAG =     1
@ global String comment SOUTH-NORTH_PATCH_END_STAG =   121
@ global String comment BOTTOM-TOP_PATCH_START_UNSTAG =     1
@ global String comment BOTTOM-TOP_PATCH_END_UNSTAG =    29
@ global String comment BOTTOM-TOP_PATCH_START_STAG =     1
@ global String comment BOTTOM-TOP_PATCH_END_STAG =    30
@ global String comment GRID_ID =     1
@ global String comment PARENT_ID =     1
@ global String comment I_PARENT_START =     1
@ global String comment J_PARENT_START =     1
@ global String comment PARENT_GRID_RATIO =     1
@ global String comment CEN_LAT =        37.38
@ global String comment CEN_LON =        90.55
@ global String comment TRUELAT1 =        30.00
@ global String comment TRUELAT2 =        60.00
@ global String comment MOAD_CEN_LAT =        37.38
@ global String comment STAND_LON =        90.55
@ global String comment POLE_LAT =        90.00
@ global String comment POLE_LON =         0.00
@ global String comment GMT =         0.00
@ global String comment JULYR =  2008
@ global String comment JULDAY =   153
@ global String comment MAP_PROJ =     1
@ global String comment MAP_PROJ_CHAR = Lambert Conformal
@ global String comment MMINLU = MODIFIED_IGBP_MODIS_NOAH
@ global String comment NUM_LAND_CAT =    21
@ global String comment ISWATER =    17
@ global String comment ISLAKE =    21
@ global String comment ISICE =    15
@ global String comment ISURBAN =    13
@ global String comment ISOILWATER =    14

1. 
   
2. <P>
   
3. <P>#脚本</P>
   
4. <P>f = addfile('E:/gongzuo/3.8.1/summer0801.ctl')
   
5. t = datetime.datetime(2008,06,01,00)</P>
   
6. <P>u = f['U'][[t],[500],:,:]
   
7. v = f['V'][[t],[500],:,:]
   
8. vort = hcurl(u, v)
   
9. axesm()
   
10. lchina = shaperead('F:/11/MeteoInfo/map/bou2_4p.shp')
    
11. ltibet = shaperead('F:/11/MeteoInfo/map/DBATP/DBATP_Polygon.shp')
    
12. geoshow(lchina, edgecolor='k')
    
13. geoshow(ltibet, edgecolor='r',size=2)
    
14. layer = contourfm(vort, 20)
    
15. colorbar(layer,fontsize=17)
    
16. xlim(70,110)
    
17. ylim(25,45)</P></P>

复制代码

图：

MeteoInfo

半颗麦芽糖 发表于 2017-1-1 15:13
王老师，我对这个脚本进行了一点修改想画我模式的涡度图，但是画出来不仅量级有问题，而且图中只显示了一种 ...

这个数据是有投影的（Lambert），hcurl, hdivg等函数是针对经纬度投影设计的，因此需要将Lambert投影的数据投影到经纬度投影坐标系中然后再计算涡度散度等。GrADS其实在读取数据时就对数据自动进行了投影转换，所以在代码里不要再进行投影处理。

投影的代码：
toproj = projinfo()
u = u.project(toproj=toproj)
u.dims[0].setDimType(DimensionType.Y)
u.dims[1].setDimType(DimensionType.X)
v = v.project(toproj=toproj)

第3和4行语句显得多余，后续的版本会在project函数中设置好，就可以去掉这2行语句了。

半颗麦芽糖

MeteoInfo 发表于 2017-1-1 17:27
这个数据是有投影的（Lambert），hcurl, hdivg等函数是针对经纬度投影设计的，因此需要将Lambert投影的数 ...

如王老师所说，修改以后能够正确出图了

马里奥大叔

很有用，谢谢分享

guo066600

无尽的兰博 发表于 2016-10-12 16:42
谢谢老师，但是我没明白是哪里改成[200]呢？
是改成
u = f['u-component_of_wind_pressure_difference ...

请问200hPa的图怎么实现的？