ncl画wrfout数据风场图

风之牧语

我想用wrfout数据画850hpa风场图，我按照官网上的程序改的，画几天的数据没有问题！但是数据大之后就会报内存不够的错误，经纬度范围又不能缩小了，请教大神们帮忙看一下程序，该如何优化啊！
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

  files=systemfunc("ls -1 wrfout_d01_2013-05-01*")
  f=addfiles(files,"r")
  ListSetType(f,"cat")

  wks = gsn_open_wks("x11","wind_May_850hpa")

; Set some Basic Plot options
   res=True
   res@gsnAddCyclic = False
   res@gsnMaximize             = True

   res@mpMinLatF               = 16.                        
   res@mpMaxLatF               = 54.
   res@mpMinLonF               = 73.
   res@mpMaxLonF               = 135.
   res@mpDataBaseVersion       = "MediumRes"
   res@mpDataSetName           = "Earth..4"
   res@mpOutlineOn           = True         ; Turn on map outlines
   res@mpOutlineSpecifiers   = (/"China","Taiwan"/)       ;China:states
   res@mpAreaMaskingOn = True   ;cover it
   res@mpMaskAreaSpecifiers = (/"China","Taiwan"/)   ;China:states
   res@mpLandFillColor         = "white"
   res@mpInlandWaterFillColor  = "white"
   res@mpOceanFillColor        = "white"
   res@mpOutlineBoundarySets   = "NoBoundaries"
   res@mpFillDrawOrder = "PostDraw"
   res@vcRefMagnitudeF         = 4.                ; make vectors larger
   res@vcRefLengthF            = 0.032              ; ref vec length
   ;res@vcGlyphStyle            = "WindBarb"         ; select wind barbs
   res@vcMinDistanceF        = 0.015                ; thin out vectors
   res@vcGlyphStyle          = "CurlyVector"   
   res@vcWindBarbColor     = "Blue"
   res@vcRefAnnoString1On = True
   res@vcRefAnnoString1   = "4m/s"
   ;res@vcLevelColors = True
   res@vpXF = 0     ;左边距
   res@vpYF = 0      ;上边距
   res@vpWidthF  = 1.0              ; height and width of plot
   res@vpHeightF = 0.8
   ;res@vcVectorDrawOrder = "PostDraw"

; The specific pressure levels that we want the data interpolated to.
  pressure_levels = (/ 850./)   ; pressure levels to plot
  nlevels         = dimsizes(pressure_levels)     ; number of pressure levels
  u  = wrf_user_getvar(f[:],"ua",-1)        ; u averaged to mass points
  v  = wrf_user_getvar(f[:],"va",-1)        ; v averaged to mass points
  p  = wrf_user_getvar(f[:], "pressure",-1) ; pressure is our vertical coordinate
  lon2d = wrf_user_getvar(f[:],"XLONG",0)
  lat2d = wrf_user_getvar(f[:],"XLAT",0)
  avg_u=dim_avg_n_Wrap(u,0)
  avg_v=dim_avg_n_Wrap(v,0)
  avg_p=dim_avg_n_Wrap(p,0)

; Plotting options for Wind Speed  
  do level = 0,nlevels-1                 ; LOOP OVER LEVELS
      pressure = pressure_levels(level)
      result_u  = wrf_user_intrp3d( avg_u,avg_p,"h",pressure,0.,False)
      result_v  = wrf_user_intrp3d( avg_v,avg_p,"h",pressure,0.,False)

      ;chazhi

      olon = fspan(73,135,248)
      olat = fspan(16,54,152)
      data1 = new((/248,152/),"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)
      grid_u@long_name="avg of uv at 850hpa in May"
      grid_v@long_name="avg of uv at 850hpa in May"

        if ( pressure .eq. 850 ) then   
          ;plot = wrf_map_overlays(f[0],wks,(/vector/),pltres,mpres)
           plot = gsn_csm_vector_map(wks,grid_u,grid_v,res)
        end if


    end do      ; END OF LEVEL LOOP
end

风之牧语

不好意思，确实弄错！

风之牧语

这个问题解决了，传一个程序很大家分享一下，可能大家觉得没有什么高端的地方的，就当是感谢论坛的朋友的帮助吧！

longlivehj

tmp  = wrf_user_getvar(f[:], "ua", -1)
avg_u = dim_avg_n_Wrap(tmp, 0)
tmp  = wrf_user_getvar(f[:], "va", -1)
avg_v = dim_avg_n_Wrap(tmp, 0)

类似于这样吧！结构相同的数组尽量重复使用，大块的数据不用了赶紧delete掉。

又是那隻貓

优化你的数组对内存的占据

anlizheng

学习学习{:5_275:}{:5_275:}{:5_275:}

风之牧语

longlivehj 发表于 2014-8-21 16:45
tmp  = wrf_user_getvar(f[:], "ua", -1)
avg_u = dim_avg_n_Wrap(tmp, 0)
tmp  = wrf_user_getvar(f[:], ...

还是不行啊！
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  files=systemfunc("ls -1 wrfout_d01_2013-05*")
  f=addfiles(files,"r")
  ListSetType(f,"cat")
  ;printVarSummary(f)
  wks = gsn_open_wks("png","wind_May_850hpa")

; Set some Basic Plot options
   res=True
   res@gsnAddCyclic = False
   res@gsnMaximize             = True
   ;res@gsnDraw                 = False
   ;res@gsnFrame                = False
   res@mpMinLatF               = 16.                        
   res@mpMaxLatF               = 54.
   res@mpMinLonF               = 73.
   res@mpMaxLonF               = 135.
   res@mpDataBaseVersion       = "MediumRes"
   res@mpDataSetName           = "Earth..4"
   res@mpOutlineOn           = True         ; Turn on map outlines
   res@mpOutlineSpecifiers   = (/"China","Taiwan"/)       ;China:states
   res@mpAreaMaskingOn = True   ;cover it
   res@mpMaskAreaSpecifiers = (/"China","Taiwan"/)   ;China:states
   res@mpLandFillColor         = "white"
   res@mpInlandWaterFillColor  = "white"
   res@mpOceanFillColor        = "white"
   res@mpOutlineBoundarySets   = "NoBoundaries"
   res@mpFillDrawOrder = "PostDraw"
   res@vcRefMagnitudeF         = 4.                ; make vectors larger
   res@vcRefLengthF            = 0.032              ; ref vec length
   ;res@vcGlyphStyle            = "WindBarb"         ; select wind barbs
   res@vcMinDistanceF        = 0.015                ; thin out vectors
   res@vcGlyphStyle          = "CurlyVector"   
   res@vcWindBarbColor     = "Blue"
   res@vcRefAnnoString1On = True
   res@vcRefAnnoString1   = "4m/s"
   ;res@vcLevelColors = True
   res@vpXF = 0     ;左边距
   res@vpYF = 0      ;上边距
   res@vpWidthF  = 1.0              ; height and width of plot
   res@vpHeightF = 0.8
   ;res@vcVectorDrawOrder = "PostDraw"

; The specific pressure levels that we want the data interpolated to.
   olon = fspan(73,135,248)
   olat = fspan(16,54,152)
   data1 = new((/248,152/),"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
   lon2d = wrf_user_getvar(f[:],"XLONG",0)
   lat2d = wrf_user_getvar(f[:],"XLAT",0)
   temp  = wrf_user_getvar(f[:],"ua",-1)        ; u averaged to mass points
   avg_u=dim_avg_n_Wrap(temp(:,0:8,:,:),0)
   temp  = wrf_user_getvar(f[:], "pressure",-1) ; pressure is our vertical coordinate
   avg_p=dim_avg_n_Wrap(temp(:,0:8,:,:),0)
   result_u  = wrf_user_intrp3d( avg_u,avg_p,"h",850.,0.,False)
   delete(avg_u)
   grid_u = rcm2rgrid_Wrap(lat2d,lon2d,result_u,olat,olon,0)
   delete(result_u)  
   temp  = wrf_user_getvar(f[:],"va",-1)        ; v averaged to mass points
   avg_v=dim_avg_n_Wrap(temp(:,0:8,:,:),0)
   delete(temp)
   result_v  = wrf_user_intrp3d( avg_v,avg_p,"h",850.,0.,False)
   delete(avg_v)
   delete(avg_p)
   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"
   ;spd     = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
   ;spd@description = "Wind Speed"
   ;spd@units = "m/s"
   ;u_plane = u_plane*1.94386     ; kts
   ;v_plane = v_plane*1.94386     ; kts
   ;u_plane@units = "kts"
   ;v_plane@units = "kts"              
   plot = gsn_csm_vector_map(wks,grid_u,grid_v,res)

end

风之牧语

我发现在 u  = wrf_user_getvar(f[:],"ua",-1)        ; u averaged to mass points读取一个变量的这个时候就会出现内存不足的错误，请问大婶们有其他的方式来读取吗？（主要是44*200*200）的格点

longlivehj

风之牧语 发表于 2014-8-25 10:26
这个问题解决了，传一个程序很大家分享一下，可能大家觉得没有什么高端的地方的，就当是感谢论坛的朋友的帮 ...

附件是不是贴错了，怎么是一堆xml的东西？！

lideqin

没往下看，果断的丢了积分！