ncl可以进行曲线插值吗？

denglin19901030

请问想对wrfout数据，如图进行曲线的垂直剖面，有可用的命令吗？

又是那隻貓

没有现成命令，不过可以自己写啊

denglin19901030

又是那隻貓 发表于 2015-4-10 10:05
没有现成命令，不过可以自己写啊

大神！我后来是想着用两个点之间做一个循环，一段一段的插值，但是画出来的图就真的变成一段一段的了，这怎么办？？？脚本如下：
;ncon_1.ncl
;
; Concepts illustrated:
;   - Drawing pressure/height contours on top of another set of contours
;   - Drawing negative contour lines as dashed lines
;   - Drawing the zero contour line thicker
;   - Changing the color of a contour line
;   - Overlaying dashed contours on solid line contours
;*************************************************
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/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;*************************************************

begin





f = addfile("wrfout_d02_2014-07-17_18:00:00","r")
hdf4_file = addfile("2A12.20140717.94948.7.HDF","r")
times = chartostring(f->Times)
lon1=hdf4_file->Longitude
lat1=hdf4_file->Latitude
qc_1=new((/100,5000/),float)
qr_1=new((/100,5000/),float)
qi_1=new((/100,5000/),float)
qs_1=new((/100,5000/),float)
qg_1=new((/100,5000/),float)
ta_1=new((/100,5000/),float)
j=0

do i=88,204
loc =wrf_user_latlon_to_ij(f,lat1(1000,i),lon1(1000,i))
loc1=wrf_user_latlon_to_ij(f,lat1(1000,i+1),lon1(1000,i+1))


pic_name = "satellite"

mdims = getfilevardimsizes(f,"P")
nd = dimsizes(mdims)
;print(mdims(0))
  qc     = f->QCLOUD(0,:,:,:)*1000                             ; get qc January
  qr     = f->QRAIN(0,:,:,:)*1000
  qi     = f->QICE(0,:,:,:)*1000
  qs     = f->QSNOW(0,:,:,:)*1000
  qg     = f->QGRAUP(0,:,:,:)*1000
  ta      = wrf_user_getvar(f,"tk",0)
  t      = ta(:,:,:)
  z      = wrf_user_getvar(f, "z",0)     ; grid point height
  printVarSummary(qc)
  Times  = chartostring(f->Times)
  xlon = wrf_user_getvar(f,"lon",0)
FirstTime = True
if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = floattoint(max(z)/1000.)
      nz   = floattoint(zmax/2 + 1)
      FirstTime = False
  end if
   angle=0.0
   plane = (/loc(1),loc(0),loc1(1),loc1(0)/)
  ; plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /)
;   X_plane = wrf_user_intrp2d(xlon,plane,angle,False)
;  X_desc = "longitude"
   opts = True         
   qc_plane = wrf_user_intrp3d(qc,z,"v",plane,angle,opts)
   qr_plane = wrf_user_intrp3d(qr,z,"v",plane,angle,opts)
   qi_plane = wrf_user_intrp3d(qi,z,"v",plane,angle,opts)
   qs_plane = wrf_user_intrp3d(qs,z,"v",plane,angle,opts)
   qg_plane = wrf_user_intrp3d(qg,z,"v",plane,angle,opts)
   ta_plane = wrf_user_intrp3d(t,z,"v",plane,angle,opts)

ndd = dimsizes(qc_plane(0,:))
   do m=0,99
   do n=0,ndd-1
   qc_1(m,j+n)=qc_plane(m,n)
   qr_1(m,j+n)=qr_plane(m,n)
   qi_1(m,j+n)=qi_plane(m,n)
   qs_1(m,j+n)=qs_plane(m,n)
   qg_1(m,j+n)=qg_plane(m,n)
   ta_1(m,j+n)=ta_plane(m,n)
   end do
   end do
j=j+ndd+1

delete(qc_plane)
delete(qr_plane)
delete(qi_plane)
delete(qs_plane)
delete(qg_plane)
delete(ta_plane)


end do

print(j)
; printVarSummary(qc_1)
   dim = dimsizes(qc_1)                      ; Find the data span - for use in labels
  zspan = dim(0)
  ; dimsX=dimsizes(X_plane)
  ; xmin = X_plane(0)
  ; xmax = X_plane(dimsX(0)-1)
  ; xspan = dimsX(0)-1
  ; nx = floattoint((xmax-xmin)+1)
  
  res                         = True
  ;      res@tiXAxisString           = X_desc
  ;      res@tiYAxisString           = "Height (km)"
;       res@cnMissingValPerimOn     = True
;        res@cnMissingValFillColor   = 0
;       res@cnMissingValFillPattern = 11
   ;     res@tmXTOn =False
    ;    res@tmYROn = False
     ;   res@tmXBMode = "Explicit"
      ;  res@tmXBValues = fspan(0,xspan,nx)
       ; res@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx))
       ; res@tmXBLabelFontHeightF = 0.015

        res@tmYLMode                = "Explicit"
        res@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
        res@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
        res@tmYLLabelStride         = 3
        res@tiXAxisFontHeightF      = 0.020
        res@tiYAxisFontHeightF      = 0.020
        res@tmXBMajorLengthF        = 0.02
        res@tmYLMajorLengthF        = 0.02
        res@tmYLLabelFontHeightF    = 0.015
;       res@PlotOrientation         = qc_plane@Orientation
  
  
  
  
  
  
  wks   = gsn_open_wks ("x11", pic_name )        ; open workstation
  

  res@gsnDraw      = False                        ; don't draw yet
  res@gsnFrame     = False                        ; don't advance frame yet
  res@tiMainString = Times(0)   ; add title
; res@gsnContourZeroLineThicknessF = 2.

freekiller

你试试把插值后的数据放到一个数组，分段画就可能不连续

denglin19901030

freekiller 发表于 2015-4-11 18:29
你试试把插值后的数据放到一个数组，分段画就可能不连续

恩恩，已经成功啦~~~~

freekiller

denglin19901030 发表于 2015-4-11 20:35  恩恩，已经成功啦~~~~

恭喜，期待分享

zzttll

请教，这个图是用什么软件画的？

denglin19901030

zzttll 发表于 2015-4-14 17:14
请教，这个图是用什么软件画的？

用ncl画的

zzttll

denglin19901030 发表于 2015-4-14 20:14
用ncl画的

哦，看过别人画过风格类似的图，没用过ncl，只会matlab。

denglin19901030

zzttll 发表于 2015-4-15 09:01
哦，看过别人画过风格类似的图，没用过ncl，只会matlab。

matlab应该都可以实现的，都是画图语言嘛