NCL用wrfout数据画垂直矢量图的问题

黄小仙儿

想画垂直风矢量图，在官网上找到了例子，但是有几个变量提示说没有，比如hyam,hybm查了一下好像说hyam,hybm是含有混合系数的一维数组，实在不知道是什么东西，请各位大神帮我看看怎么才能用wrfout数据画出这种图
这是脚本：
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin

; file handling
;************************************************
  ; define filename
   in  = addfile("/home/Huanglei/data/d03"+".nc","r")                         ; open netcdf file
;************************************************
; read needed variables from file
;************************************************
   T = in->T                                    ; select variable to ave
   W = in->W
   V = in->V
   P0mb = 1000.
   hyam = in->hyam                              ; get a coefficiants
   hybm = in->hybm                              ; get b coefficiants
   PS   = in->PS                                ; get pressure
;************************************************
; define other arguments required by vinth2p
;************************************************
   interp = 2
   pnew   = (/ 900,850,800,750,700,650,600,550,500,450,400,350,300,250,200/)
   pnew@units = "mb"         
;************************************************
; interpolate to pressure levels on pressure levels
;************************************************
  t = vinth2p(T,hyam,hybm,pnew,PS,interp,P0mb,1,False)
  v = vinth2p(V,hyam,hybm,pnew,PS,interp,P0mb,1,False)
  w = vinth2p(W,hyam,hybm,pnew,PS,interp,P0mb,1,False)
;************************************************
; Omega is significantly smaller than v, so we will
; scale it so that some vertical motion is visible
;************************************************
wAve   = avg(w(0,:,:,{104}))           ; used for scaling
vAve   = avg(v(0,:,:,{104}))
scale  = fabs(vAve/wAve)
wscale = w*scale                       ; now scale

copy_VarCoords(w, wscale)              ; copy coordinate variables
;***********************************************
; create plot
;***********************************************
wks   = gsn_open_wks ("pdf", "vector" )       ; open workstation
gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map

res                 = True                     ; plot mods desired
res@tiMainString    = "Pressure/Height Vector" ; title

res@cnLineLabelsOn  = False              ; turn off line labels
res@cnFillOn        = True               ; turn on color fill
res@lbLabelStride   = 2                  ; every other color

res@gsnSpreadColors = True               ; use full range of color map

res@vcRefMagnitudeF = 3.0                ; define vector ref mag
res@vcRefLengthF    = 0.045              ; define length of vec ref
res@vcGlyphStyle    = "CurlyVector"      ; turn on curley vectors
res@vcMinDistanceF  = 0.01               ; thin out vectors
res@vcMapDirection  = False

;*****************************************************
; draw plot from pole to pole at 170E
;*****************************************************
plot  = gsn_csm_pres_hgt_vector(wks,t(0,:,:,{104}),v(0,:,:,{104}),\
                                wscale(0,:,:,{104}),res )  

end

官网上的图是这样的

longlivehj

黄小仙儿 发表于 2014-4-29 20:04
长生大神，我想把温度去掉，直接在最后的函数去掉提示有错，说必须要三个变量。然后我就把温度fillon改成 ...

取消填充，同时把等值线颜色变为背景色（白色），有必要的话，将等值线的绘图顺寻改成“PreDraw”。

longlivehj

这不是维度问题。条件表达式里面出现了missing value。你就改改addfile右面的文件路径也会出这错误啊？把代码贴上来呢！

longlivehj

你可以看一下vinth2p的帮助，它适合于CESM hybrid coordinates，和wrf的垂直坐标系统是不一样的，所以你在wrfout里面根本找不到hyam、hybm。

wrf有专属插值函数wrf_user_intrp3d。插值好后，再用gsn_csm_pres_hgt_vector做图！

黄小仙儿

longlivehj 发表于 2014-4-12 19:40
你可以看一下vinth2p的帮助，它适合于CESM hybrid coordinates，和wrf的垂直坐标系统是不一样的，所以你在w ...

好的，谢谢提醒，我再好好研究一下~

黄小仙儿

longlivehj 发表于 2014-4-12 19:40
你可以看一下vinth2p的帮助，它适合于CESM hybrid coordinates，和wrf的垂直坐标系统是不一样的，所以你在w ...

你好，我按照wrf_user_intrp3d还有gsn_csm_pres_hgt_vector的函数介绍改了一下，但是提示
    w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
  execute error。看了好久都没发现有什么错误
而且我总感觉我的这个脚本有些不对，您能不能帮我看看还有那些需要改的？万分感谢！
;----------------------------------------------------------------------
; vector_5.ncl
;
; Concepts illustrated:
;   - Drawing pressure/height vectors over filled contours
;   - Using "vinth2p" to interpolate to user specified pressure levels
;   - Drawing curly vectors
;   - Thinning vectors using a minimum distance resource
;----------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin

; file handling
;************************************************
  ; define filename
   in  = addfile("/home/Huanglei/data/d03"+".nc","r")                         ; open netcdf file
;************************************************
; read needed variables from file
;************************************************
  times  = wrf_user_getvar(in,"times",-1) ; get times in the file
  ntimes = dimsizes(times)          ; number of times in the file
  FirstTime = True

  mdims = getfilevardimsizes(in,"P") ; get some dimension sizes for the file
  nd = dimsizes(mdims)

;---------------------------------------------------------------

  do it =3,ntimes-1,2                  ; TIME LOOP

    print("Working on time: " + times(it) )
   

    tc  = wrf_user_getvar(in,"tc",it)     ; T in C
    v = wrf_user_getvar(in,"V",it)      ; v wind
        w = wrf_user_getvar(in,"W",it)      ; w wind
    z   = wrf_user_getvar(in, "z",it)     ; grid point height

    if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = max(z)/1000.
      nz   = floattoint(zmax/2 + 1)
      FirstTime = False
    end if

       plane = new(4,float)
       plane = (/ 2,2, mdims(nd-1)-2, mdims(nd-2)-2 /)   
     opts = True  
     t_plane = wrf_user_intrp3d(tc,z,"v",plane,0.,opts)
     v_plane = wrf_user_intrp3d(v,z,"v",plane,0.,opts)
     w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)

     dim = dimsizes(t_plane)                      ; Find the data span - for use in labels
      zspan = dim(0)
                ; Options for XY Plots
        opts_xy                         = res
        opts_xy@tiYAxisString           = "Height (km)"
        opts_xy@cnMissingValPerimOn     = True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
        opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
        opts_xy@tiXAxisFontHeightF      = 0.020
        opts_xy@tiYAxisFontHeightF      = 0.020
        opts_xy@tmXBMajorLengthF        = 0.02
        opts_xy@tmYLMajorLengthF        = 0.02
        opts_xy@tmYLLabelFontHeightF    = 0.015
        opts_xy@PlotOrientation         = tc_plane@Orientation
;************************************************
; Omega is significantly smaller than v, so we will
; scale it so that some vertical motion is visible
;************************************************
wAve   = avg(w_plane(it,:,:,{104}))           ; used for scaling
vAve   = avg(v_plane(it,:,:,{104}))
scale  = fabs(vAve/wAve)
wscale = w*scale                       ; now scale

copy_VarCoords(w, wscale)              ; copy coordinate variables
;***********************************************
; create plot
;***********************************************
wks   = gsn_open_wks ("pdf", "vector" )       ; open workstation
gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map

res                 = True                     ; plot mods desired
res@tiMainString    = "Pressure/Height Vector" ; title
res@TimeLabel = times(it)           ; Set Valid time to use on plots
res@cnLineLabelsOn  = False              ; turn off line labels
res@cnFillOn        = True               ; turn on color fill
res@lbLabelStride   = 2                  ; every other color

res@gsnSpreadColors = True               ; use full range of color map

res@vcRefMagnitudeF = 3.0                ; define vector ref mag
res@vcRefLengthF    = 0.045              ; define length of vec ref
res@vcGlyphStyle    = "CurlyVector"      ; turn on curley vectors
res@vcMinDistanceF  = 0.01               ; thin out vectors
res@vcMapDirection  = False

;*****************************************************
; draw plot from pole to pole at 170E
;*****************************************************
plot  = gsn_csm_pres_hgt_vector(wks,t(it,:,:,{104}),v(it,:,:,{104}),\
                                wscale(it,:,:,{104}),res )  
end do
end

longlivehj

黄小仙儿 发表于 2014-4-14 11:59
你好，我按照wrf_user_intrp3d还有gsn_csm_pres_hgt_vector的函数介绍改了一下，但是提示
    w_plane = ...

有点小乱哦，我仔细看看哈，需要点时间。
目前发现：
1. w在剖面之前要进行处理，w := wrf_user_unstagger(w, w@stagger)
2. 代码是在示例基础上改的，没有完全更新，到后面还有还有错，你好好看看。

黄小仙儿

这个程序改了好久，还是走不通，后来给ncltalk发了邮件，回复说gsn_csm_pres_hgt_vector函数必须要求lev的单位是hpa,p啊，mb之类的，所以wrfout数据不能用这个函数画
n order to use any of the gsn_csm_pres_hgt routines, you must have data arrays that are rectilinear and contain a 1D level coordinate array.

The wrf interpolation functions you are using do not return the data on a rectilinear grid, and hence you can’t use gsn_csm_pres_hgt_vector to plot it.

longlivehj

黄小仙儿 发表于 2014-4-17 17:31
这个程序改了好久，还是走不通，后来给ncltalk发了邮件，回复说gsn_csm_pres_hgt_vector函数必须要求lev的 ...

o(∩∩)o...哈哈，差点忘了这事儿了，刚把代码修改一下！图片是用官网katrina例子的输出画的，看来是可以实现的哦！注意标红的部分。



load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin

    ; file handling
    ;************************************************
    ; define filename
    in  = addfile("/home/Huanglei/data/d03"+".nc","r")                         ; open netcdf file
    ;************************************************
    ; read needed variables from file
    ;************************************************
    times  = wrf_user_getvar(in,"times",-1) ; get times in the file
    ntimes = dimsizes(times)          ; number of times in the file
    mdims = getfilevardimsizes(in,"P") ; get some dimension sizes for the file
    nd = dimsizes(mdims)

    do it =3,ntimes-1,2                  ; TIME LOOP

        print("Working on time: " + times(it) )
      
        tc  = wrf_user_getvar(in,"tc",it)     ; T in C
        v = wrf_user_getvar(in,"V",it)      ; v wind
        w = wrf_user_getvar(in,"W",it)      ; w wind
        p = wrf_user_getvar(in, "pressure",it)

        w := wrf_user_unstagger(w, w@stagger)

        plane = new(4,float)
        plane = (/ 2,2, mdims(nd-1)-2, mdims(nd-2)-2 /)   

        opts = True
        t_plane = wrf_user_intrp3d(tc,p,"v",plane,0.,opts)
        v_plane = wrf_user_intrp3d(v, p,"v",plane, 0.,opts)
        w_plane = wrf_user_intrp3d(w, p,"v",plane, 0.,opts)
        p_plane = wrf_user_intrp3d(p, p,"v",plane, 0.,opts)

        ;************************************************
        wAve   = avg(w_plane(:,104))           ; used for scaling
        vAve   = avg(v_plane(:,104))
        scale  = fabs(vAve/wAve)
        wscale = w_plane*scale                       ; now scale

        t_plane!0="level"
        t_plane&level = p_plane(:, 0)
        copy_VarCoords(t_plane, v_plane)
        copy_VarCoords(t_plane, wscale)

        ;***********************************************
        ; create plot
        ;***********************************************
        wks   = gsn_open_wks ("pdf", "vector" )       ; open workstation
        gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map

        res                 = True                    ; plot mods desired
        res@tiMainString    = "Pressure/Height Vector" ; title
        res@tiYAxisString           = "Pressure (mb)"
        res@cnLineLabelsOn  = False              ; turn off line labels
        res@cnFillOn        = True               ; turn on color fill
        res@lbLabelStride   = 2                  ; every other color

        res@gsnSpreadColors = True               ; use full range of color map

        res@vcRefMagnitudeF = 3.0                ; define vector ref mag
        res@vcRefLengthF    = 0.045              ; define length of vec ref
        res@vcGlyphStyle    = "CurlyVector"      ; turn on curley vectors
        res@vcMinDistanceF  = 0.03               ; thin out vectors
        res@vcMapDirection  = False

        plot  = gsn_csm_pres_hgt_vector(wks,t_plane,v_plane,wscale,res )  
    end do
end

黄小仙儿

longlivehj 发表于 2014-4-17 19:25
o(∩∩)o...哈哈，差点忘了这事儿了，刚把代码修改一下！图片是用官网katrina例子的输出画的，看来是可以 ...

哇，太感谢了。
但是我刚刚试着运行了下，还是提示维度不匹配

fatal:The result of the conditional expression yields a missing value. NCL can not determine branch, see ismissing function
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 12154 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl

fatal:["Execute.c":8567]:Execute: Error occurred at or near line 12566 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl

fatal:["Execute.c":8567]:Execute: Error occurred at or near line 90 in file /home/Huanglei/vector.ncl
第90行就是plot  = gsn_csm_pres_hgt_vector(wks,t_plane,v_plane,wscale,res )  

黄小仙儿

longlivehj 发表于 2014-4-17 19:25
o(∩∩)o...哈哈，差点忘了这事儿了，刚把代码修改一下！图片是用官网katrina例子的输出画的，看来是可以 ...

我刚查看了数据维度，好像没有错啊

Variable: v_plane
Type: float
Total Size: 99736 bytes
            24934 values
Number of Dimensions: 2
Dimensions and sizes:   [level | 91] x [Horizontal | 274]
Coordinates:
            level: [9.96921e+36..70]
Number Of Attributes: 4
  _FillValue :  9.96921e+36
  units :       m s-1
  description : y-wind component
  Orientation : Cross-Section: (2,2) to (195,195)

Variable: t_plane
Type: float
Total Size: 99736 bytes
            24934 values
Number of Dimensions: 2
Dimensions and sizes:   [level | 91] x [Horizontal | 274]
Coordinates:
            level: [9.96921e+36..70]
Number Of Attributes: 4
  _FillValue :  9.96921e+36
  units :       C
  description : Temperature
  Orientation : Cross-Section: (2,2) to (195,195)

Variable: wscale
Type: float
Total Size: 99736 bytes
            24934 values
Number of Dimensions: 2
Dimensions and sizes:   [level | 91] x [Horizontal | 274]
Coordinates:
            level: [9.96921e+36..70]
Number Of Attributes: 1
  _FillValue :  9.96921e+36
fatal:The result of the conditional expression yields a missing value. NCL can not determine branch, see ismissing function
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 12154 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl