NCL作风场图区域出错

SG哭晕在厕所

我画了几张风场和相对湿度的叠加图，风场使用ncl自带的wrf插值函数，但是插出来以后850hpa是正常的，而到950或是925hpa时有些区域的风向标莫名消失了，各位大神能帮忙看看是什么原因吗？三张图分别是500,850和925hpa的图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

以下是代码

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/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
begin
;;read u,v,t from the data at 500hPa
files=systemfunc("ls wrfchem1/WRFV3/test/em_real/wrfout_d01_2016-11-04*")
;  f=addfile("wrfchem1/WRFV3/test/em_real/wrfout_d01_2016-11-04_06:00:00.nc","r")
  f = addfiles(files,"r")
  ListSetType(f,"cat")
  times = wrf_user_getvar(f,"times",-1)  ; get all times in the file
  fin=addfile("wrfchem1/WRFV3/test/em_real/wrfout_d01_2016-11-04_18:00:00","r")
  avg_rh = wrf_user_getvar(fin,"rh",0)        
  ntimes = dimsizes(times)         ; number of times in the file
  
  sum_u = wrf_user_getvar(f,"ua",0)        ; u averaged to mass points
  sum_v = wrf_user_getvar(f,"va",0)        ; v averaged to mass point
  p  = wrf_user_getvar(f, "pressure",0) ; pressure is our vertical coordinate
  lon2d = wrf_user_getvar(f,"XLONG",0)
  lat2d = wrf_user_getvar(f,"XLAT",0)
  avg_rh@lon2d = wrf_user_getvar(fin,"XLONG",0)
  avg_rh@lat2d = wrf_user_getvar(fin,"XLAT",0)
  do it = 1,ntimes-1,1
    u  = wrf_user_getvar(f,"ua",it)        ; u averaged to mass points
    v  = wrf_user_getvar(f,"va",it)        ; v averaged to mass point
    sum_u=sum_u+u
    sum_v=sum_v+v
    delete(u)
    delete(v)
  end do
  avg_u=sum_u/ntimes
  avg_v=sum_v/ntimes
  result_u  = wrf_user_intrp3d(avg_u,p,"h",925.,0.,False)
  result_v  = wrf_user_intrp3d(avg_v,p,"h",850.,0.,False)
  delete(avg_u)
  delete(avg_v)
  delete(sum_u)
  delete(sum_v)
  olon = fspan(104,129,100)
  olat = fspan(20,52,117)
  data1 = new((/100,117/),"float")
  olon!0          = "lon"
  olon@long_name  = "lon"
  olon@units      = "degrees-east"
  olon&lon        = olon
  olat!0          = "lat"
  olat@long_name  = "lat"
  olat@units      = "degrees_north"
  olat&lat        = olat
  grid_u = rcm2rgrid_Wrap(lat2d,lon2d,result_u,olat,olon,0)
  grid_v = rcm2rgrid_Wrap(lat2d,lon2d,result_v,olat,olon,0)
;  delete(lon2d)
;  delete(lat2d)
;  grid_u@long_name="avg of uv at 850hpa in May"
;  grid_v@long_name="avg of uv at 850hpa in May"
printVarSummary(avg_rh)

;;create plots;;
        wks = gsn_open_wks("png","overlay") ; send graphics to PNG file
;       gsn_define_colormap(wks,"BkBlAqGrYeOrReViWh200")
        res                = True
        res@gsnDraw        = False
        res@gsnFrame       = False
        res@gsnMaximize    = True
        res@tmXTOn         = False
        res@tmYROn         = False
   
        
;;set map;;
        mpres                             = res
        mpres@mpDataSetName               = "Earth..4"
        mpres@mpDataBaseVersion           = "MediumRes"
        mpres@mpOutlineOn                 = True
        mpres@mpOutlineSpecifiers         = (/"China:states","Taiwan"/)
        mpres@mpGeophysicalLineThicknessF = 2
        mpres@mpNationalLineThicknessF    = 2
        mpres@mpFillDrawOrder             = "PostDraw"
        mpres@mpFillOn                    = True
        mpres@mpFillAreaSpecifiers        = (/"water",       "land" /)
        mpres@mpSpecifiedFillColors       = (/"deepskyblue2","white"/)
;       mpres@mpSpecifiedFillColors      = (/100,0/)
        mpres@mpMaskAreaSpecifiers        = (/"China:states","Taiwan"/)

;;set area;;
   mpres@mpLimitMode             = "LatLon"
   mpres@mpMinLatF               = 26;min(lat2d)-1                        
   mpres@mpMaxLatF               = 44;max(lat2d)+1
   mpres@mpMinLonF               = 106;min(lon2d)-1
   mpres@mpMaxLonF               = 126;max(lon2d)+1            
;;set contour;;
        cnres                             = res
        cnres@cnFillDrawOrder             = "PreDraw"
        cnres@cnFillOn                    = True
        cnres@cnLinesOn                   = False
        cnres@pmLabelBarWidthF            = 0.4
        cnres@pmLabelBarHeightF           = 0.05
        cnres@pmLabelBarOrthogonalPosF    = 0.1
        cnres@lbLabelFontHeightF          = 0.006
        cnres@lbLabelAngleF               = 45
        cnres@trGridType              = "TriangularMesh"
;        cnres@cnFillMode              = "RasterFill"
; Older way to subset a color map
;       cnres@cnFillPalette               = "BkBlAqGrYeOrReViWh200"
;       cnres@gsnSpreadColorStart         = 50
;       cnres@gsnSpreadColorEnd           = 120

; Newer way to subset a color map
;        cmap = read_colormap_file("BkBlAqGrYeOrReViWh200")
;        cnres@cnFillPalette               = cmap(25:120,:)

        cnres@gsnLeftString               = "RH"
        
;;set vector;;
        res_vc                            = res
        res_vc@vcGlyphStyle               = "LineArrow"
        res_vc@vcLineArrowThicknessF      = 4
        res_vc@vcRefLengthF               = 0.04
        res_vc@vcRefMagnitudeF         = 10.                ; make vectors larger
        res_vc@vcMinDistanceF        = 0.03                ; thin out vectors  
;       res_vc@vcWindBarbColor     = "Blue"
        res_vc@vcRefAnnoString1On = True
        res_vc@vcRefAnnoString1   = "10m/s"
;        res_vc@vcMonoLineArrowColor = False
;        res_vc@vcLevelColors = True
;        res_vc@vpXF = 0   
;        res_vc@vpYF = 1.0              ; height and width of plot
;        res_vc@vpHeightF = 0.8


;;wind barb resources don't apply
;;      res_vc@vcGlyphStyle               = "WindBarb"
;;      res_vc@vcWindBarbLineThicknessF   = 5
;;      res_vc@vcWindBarbColor            = "Gray40"

      
        res_vc@vcRefAnnoSide             = "Top"
        res_vc@vcRefAnnoString2On        = False
        res_vc@vcRefAnnoPerimOn          = False
        res_vc@vcRefAnnoOrthogonalPosF   = -0.12
        res_vc@vcRefAnnoParallelPosF     = 0.999
        res_vc@vcRefAnnoBackgroundColor  = "Purple"
        res_vc@vcVectorDrawOrder         = "PostDraw"
        res_vc@gsnRightString            = "Wind"
        
;;plot;;
        map     = gsn_csm_map(wks,mpres)
        contour = gsn_csm_contour(wks,avg_rh(0,:,:),cnres)
        vector  = gsn_csm_vector(wks,grid_u,grid_v,res_vc)

;;overlay filled contours and vectors on the map;;
        overlay(map,contour)
        overlay(map,vector)



;;drawing "map" will draw everything: map, contours, vectors, and text;;
        draw(map)
        frame(wks)
end

随缘

说明那块风速很小，或者地形已经高过925hpa了

SG哭晕在厕所

随缘 发表于 2017-4-26 14:30
说明那块风速很小，或者地形已经高过925hpa了

原来如此。。。谢谢

xiaoguoguo327

请教楼主，grid_u = rcm2rgrid_Wrap(lat2d,lon2d,result_u,olat,olon,0)
  grid_v = rcm2rgrid_Wrap(lat2d,lon2d,result_v,olat,olon,0)这两句用这个函数插值是要达到什么目的呢？我看到好多帖子都提到了这个函数，但是还是不太明白它是做什么的？