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

黄小仙儿

黄小仙儿 发表于 2014-4-17 22:23
貌似是啊：

Variable: tc (subsection)

我记得我转模式的时候设定的垂直分层就是27层，但是为什么低层没有数据啊，不懂了。不应该啊。因为一般默认就是27层。
层数设多几层不知道能不能解决这个问题？

黄小仙儿

longlivehj 发表于 2014-4-17 22:30
plot  = gsn_csm_pres_hgt_vector(wks,t_plane(23:, :),v_plane(23:, :),wscale(23:, :),res )

改成 ...

fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 38 in file /home/Huanglei/vector.ncl
line 38;         w = wrf_user_getvar(a,"W",it)      ; w wind
这个不是赋值语句吗，怎么会出现不匹配的错误。。。
坎坷的画图之路~~

longlivehj

黄小仙儿 发表于 2014-4-17 22:35
fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Ex ...

w的维数发生过改变（w := wrf_user_unstagger(w, w@stagger)），所以获取w的那句要改成：
w := wrf_user_getvar(a,"W",it)
就是等号前加个冒号。我这边只画了一个时次，所以没注意有这个情况，不好意思！

27层的问题，估计你得看看你namelist里面的设置了，是不是设置eta level的地方有问题？

黄小仙儿

longlivehj 发表于 2014-4-17 23:08
w的维数发生过改变（w := wrf_user_unstagger(w, w@stagger)），所以获取w的那句要改成：
w := wrf_user ...

哇，太感谢了。
但是还有一个问题，第一个时次之后v_plane，和t_plane又出现了维度不匹配的问题

warning:VarVarWrite: Dimension names for dimension number (0) don't match, assigning name of rhs dimension to lhs and                                                                                      overwriting coordinate variable, use "(/../)" if this change is not desired
warning:["Execute.c":8567]:Execute: Error occurred at or near line 50 in file /home/Huanglei/vector.ncl

warning:VarVarWrite: Dimension names for dimension number (0) don't match, assigning name of rhs dimension to lhs and overwriting coordinate variable, use "(/../)" if this change is not desired
warning:["Execute.c":8567]:Execute: Error occurred at or near line 51 in file /home/Huanglei/vector.ncl
50      t_plane = wrf_user_intrp3d(tc,p,"v",plane,0.,opts)
51     v_plane = wrf_user_intrp3d(v,p,"v",plane,0.,opts)
我试着把v,tc也stagger后又出现了这样的警告:
NOTE: No unstaggering required, as the input field is already on mass points.
我查了下stagger的函数说明，好像主要是用于速度变量，但是还是无法理解这个函数所起的作用，请问能不能讲解一下。

longlivehj

黄小仙儿 发表于 2014-4-18 10:47
哇，太感谢了。
但是还有一个问题，第一个时次之后v_plane，和t_plane又出现了维度不匹配的问题

这个warning没有关系，是坐标变量名称不匹配，不影响程序执行的。如果非要消除，你把下面代码addfile后的路径改改再运行试试。

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("wrfout_d01_2005-08-28_00:00:00","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 =1,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&Vertical = p_plane(:, 0)
        copy_VarCoords(t_plane, v_plane)
        copy_VarCoords(t_plane, wscale)

        ;***********************************************
        ; create plot
        ;***********************************************
        wks   = gsn_open_wks ("x11", "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(23:, :), v_plane(23:, :), wscale(23:, :), res)

        delete(t_plane&Vertical)
        delete(v_plane&Vertical)
        delete(wscale&Vertical)
    end do
end

longlivehj

黄小仙儿 发表于 2014-4-18 10:47
哇，太感谢了。
但是还有一个问题，第一个时次之后v_plane，和t_plane又出现了维度不匹配的问题

只需要对w进行stagger操作。
这方面的信息，你看下http://www.openwfm.org/wiki/How_to_interpret_WRF_variables

黄小仙儿

longlivehj 发表于 2014-4-18 11:26
只需要对w进行stagger操作。
这方面的信息，你看下http://www.openwfm.org/wiki/How_to_interpret_WRF_v ...

好的，太感谢了。自己基础知识还是不稳固。感谢家园里的热心人们~

黄小仙儿

longlivehj 发表于 2014-4-18 11:23
这个warning没有关系，是坐标变量名称不匹配，不影响程序执行的。如果非要消除，你把下面代码addfile后的 ...

你好，想再问一下，修改这个addfile路径，为什么会消除这个warning。
我如果不加+".nc"的话，会提示找不到文件
我又换了数据，warning了几个时次之后就fatal error了。不知道是什么原因

longlivehj

黄小仙儿 发表于 2014-4-23 22:27
你好，想再问一下，修改这个addfile路径，为什么会消除这个warning。
我如果不加+".nc"的话，会提示找不 ...

我在测试代码的时候，把addfile中文件换成了我电脑上的，回帖给你时忘记将文件路路径还原成你的。所以，要你把文件路径改改，并不是说warning的消失是因为这个原因。这个warning是因为那几个plane的坐标名称不一致引起的，这些warning的描述里面都已经说明了的。

至于换了数据出现错误，你给的信息量太少，无法判断是什么原因！

黄小仙儿

longlivehj 发表于 2014-4-23 23:01
我在测试代码的时候，把addfile中文件换成了我电脑上的，回帖给你时忘记将文件路路径还原成你的。所以， ...

这是我的脚本
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
   a  = addfile("/home/Huanglei/data/d032"+".nc","r")                         ; open netcdf file
;************************************************
; read needed variables from file
;************************************************
  times  = wrf_user_getvar(a,"times",-1) ; get times in the file
  ntimes = dimsizes(times)          ; number of times in the file
  mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
  nd = dimsizes(mdims)
wks   = gsn_open_wks ("pdf", "vector" )       ; open workstation
gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map

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

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

    print("Working on time: " + times(it) )
   res                 = True                     ; plot mods desired
   res@tiMainString = times(it)           ; Set Valid time to use on plots

     tc = wrf_user_getvar(a,"tc",it)     ; T in C
     v = wrf_user_getvar(a,"V",it)      ; v wind
         w := wrf_user_getvar(a,"W",it)      ; w wind
     p  = wrf_user_getvar(a, "pressure",it)     ; grid point height
    z  = wrf_user_getvar(a, "z",it)     ; grid point height
     w := wrf_user_unstagger(w, w@stagger)
    ; tc := wrf_user_unstagger(tc, tc@stagger)
    ; v := wrf_user_unstagger(v, v@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)
  
;************************************************
; Omega is significantly smaller than v, so we will
; scale it so that some vertical motion is visible
;************************************************
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)              ; copy coordinate variables
        ; printVarSummary(v_plane)
    ; printVarSummary(t_plane)
    ; printVarSummary(wscale)
        ; printVarSummary(p_plane)
         print(p_plane(:,0))
         print(p(:, 0, 0))
          print(z(:, 0, 0))
        ; print(tc(:, 0, 0))
        ; print(v(:, 0, 0))

;***********************************************
; create plot
;***********************************************


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

;*****************************************************
; draw plot from pole to pole at 170E
;*****************************************************
plot  = gsn_csm_pres_hgt_vector(wks,t_plane(29:,:),v_plane(29:,:),wscale(29:,:),res )  
end do

end