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

longlivehj

黄小仙儿 发表于 2014-4-17 21:38
我刚查看了数据维度，好像没有错啊

Variable: v_plane

level坐标变量有问题。
还是要看你的代码，估计pressure插值的地方错了。

黄小仙儿

longlivehj 发表于 2014-4-17 21:40
level坐标变量有问题。
还是要看你的代码，估计pressure插值的地方错了。

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/d03"+".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)

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

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

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


     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
     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)
  
;************************************************
; 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)
;***********************************************
; create plot
;***********************************************
wks   = gsn_open_wks ("pdf", "vector" )       ; open workstation
gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map

res                 = True                     ; plot mods desired
res@TimeLabel = times(it)           ; Set Valid time to use on plots
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,v_plane,wscale,res )  
end do

end

黄小仙儿

longlivehj 发表于 2014-4-17 21:40
level坐标变量有问题。
还是要看你的代码，估计pressure插值的地方错了。

我发现wscale的attribute只有一个fillvalue，其他都有4个。
但是已经 copy_VarCoords(t_plane, wscale)了啊，是不是跟t的单位这些不匹配造成的？

longlivehj

黄小仙儿 发表于 2014-4-17 22:02
我发现wscale的attribute只有一个fillvalue，其他都有4个。
但是已经 copy_VarCoords(t_plane, wscale) ...

level: [9.96921e+36..70]第一个值是有问题的，也就是missing value。所以纵坐标会出错。
你把p_plane的信息打印一个贴出来看看，另外输出p_plane(:,0)所有值贴出来。

黄小仙儿

longlivehj 发表于 2014-4-17 22:05
level: [9.96921e+36..70]第一个值是有问题的，也就是missing value。所以纵坐标会出错。
你把p_plane的 ...

Variable: p_plane (subsection)
Type: float
Total Size: 364 bytes
            91 values
Number of Dimensions: 1
Dimensions and sizes:   [Vertical | 91]
Coordinates:
Number Of Attributes: 4
  Orientation : Cross-Section: (2,2) to (195,195)
  description : Pressure
  units :       hPa
  _FillValue :  9.96921e+36
(0)     9.96921e+36
(1)     9.96921e+36
(2)     9.96921e+36
(3)     9.96921e+36
(4)     9.96921e+36
(5)     9.96921e+36
(6)     9.96921e+36
(7)     9.96921e+36
(8)     9.96921e+36
(9)     9.96921e+36
(10)    9.96921e+36
(11)    9.96921e+36
(12)    9.96921e+36
(13)    9.96921e+36
(14)    9.96921e+36
(15)    9.96921e+36
(16)    9.96921e+36
(17)    9.96921e+36
(18)    9.96921e+36
(19)    9.96921e+36
(20)    9.96921e+36
(21)    9.96921e+36
(22)    9.96921e+36
(23)    740
(24)    730
(25)    720
(26)    710
(27)    700
(28)    690
(29)    680
(30)    670
(31)    660
(32)    650
(33)    640
(34)    630
(35)    620
(36)    610
(37)    600
(38)    590
(39)    580
(40)    570
(41)    560
(42)    550
(43)    540
(44)    530
(45)    520
(46)    510
(47)    500
(48)    490
(49)    480
(50)    470
(51)    460
(52)    450
(53)    440
(54)    430
(55)    420
(56)    410
(57)    400
(58)    390
(59)    380
(60)    370
(61)    360
(62)    350
(63)    340
(64)    330
看来低层pressure是错的，这是为什么？有什么修改办法？

longlivehj

黄小仙儿 发表于 2014-4-17 22:10
Variable: p_plane (subsection)
Type: float
Total Size: 364 bytes

print(p(:, 0, 0))看看输出结果，怀疑文件里面的pressure不对。

黄小仙儿

longlivehj 发表于 2014-4-17 22:12
print(p(:, 0, 0))看看输出结果，怀疑文件里面的pressure不对。

Variable: p (subsection)
Type: float
Total Size: 108 bytes
            27 values
Number of Dimensions: 1
Dimensions and sizes:   [bottom_top | 27]
Coordinates:
Number Of Attributes: 6
  FieldType :   104
  MemoryOrder : XYZ
  description : Pressure
  units :       hPa
  stagger :
  coordinates : XLONG XLAT
(0)     734.5763
(1)     728.6539
(2)     720.5775
(3)     710.4689
(4)     698.0468
(5)     682.7861
(6)     663.9152
(7)     636.6399
(8)     601.5903
(9)     566.8007
(10)    532.2227
(11)    484.9306
(12)    428.3757
(13)    377.5354
(14)    331.8683
(15)    290.9242
(16)    254.2611
(17)    221.5859
(18)    192.4622
(19)    166.6089
(20)    143.7232
(21)    123.4688
(22)    105.766
(23)    90.2112
(24)    76.62739
(25)    64.83072
(26)    54.69307

longlivehj

哇，真的是少低层数据哦！
看看t、v的信息，是不是垂直都只有27层？

黄小仙儿

longlivehj 发表于 2014-4-17 22:17
哇，真的是少低层数据哦！
看看t、v的信息，是不是垂直都只有27层？

貌似是啊：

Variable: tc (subsection)
Type: float
Total Size: 108 bytes
            27 values
Number of Dimensions: 1
Dimensions and sizes:   [bottom_top | 27]
Coordinates:
Number Of Attributes: 6
  coordinates : XLONG XLAT
  stagger :
  units :       C
  MemoryOrder : XYZ
  FieldType :   104
  description : Temperature
(0)     15.73135
(1)     15.39087
(2)     14.76349
(3)     14.14319
(4)     13.31049
(5)     12.22168
(6)     10.72394
(7)     8.699066
(8)     6.098785
(9)     3.867065
(10)    1.590088
(11)    -2.299255
(12)    -8.191071
(13)    -14.4848
(14)    -21.32133
(15)    -28.58134
(16)    -35.94986
(17)    -43.39331
(18)    -51.00345
(19)    -58.81908
(20)    -66.00589
(21)    -70.42078
(22)    -73.3201
(23)    -73.6741
(24)    -72.00471
(25)    -69.03786
(26)    -64.97437


Variable: v (subsection)
Type: float
Total Size: 108 bytes
            27 values
Number of Dimensions: 1
Dimensions and sizes:   [bottom_top | 27]
Coordinates:
Number Of Attributes: 6
  FieldType :   104
  MemoryOrder : XYZ
  description : y-wind component
  units :       m s-1
  stagger :     Y
  coordinates : XLONG_V XLAT_V
(0)     0.8362358
(1)     1.365672
(2)     1.593711
(3)     2.053607
(4)     2.742851
(5)     3.284138
(6)     3.329784
(7)     3.668312
(8)     2.785434
(9)     0.3634361
(10)    -2.049346
(11)    -4.532928
(12)    -5.971051
(13)    -6.456445
(14)    -8.979906
(15)    -11.7275
(16)    -11.22292
(17)    -8.77671
(18)    -6.193141
(19)    -5.025019
(20)    -3.98056
(21)    -4.298276
(22)    -7.213253
(23)    -8.591394
(24)    -7.177768
(25)    -4.694066
(26)    -2.382471

longlivehj

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

Variable: tc (subsection)

plot  = gsn_csm_pres_hgt_vector(wks,t_plane(23:, :),v_plane(23:, :),wscale(23:, :),res )

改成这个，应该可以出图了。