雷达剖面图和水平图对应不上

我不喜欢吃生姜

沿着图1所示剖线做dbz的剖面图，结果发现和图1对应不上（例如31-31.5N在剖面图上应该有dbz的大值）。看了官网的脚本，也不知道是哪地方出错了~请各位大神指教！

图1

我不喜欢吃生姜

脚本如下

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

  a = addfile("/wrfout_24_00:00:00.nc","r")
; We generate plots, but what kind do we prefer?
  type = "png"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"dbz115.95-12km")
  gsn_define_colormap(wks,"radar_1")
  
      ; Cross section
  
       minlat= 30.5
       maxlat=  32.0
       minlon = 115.95
       maxlon = 115.95

  
   loc = wrf_user_ll_to_ij(a,(/minlon,maxlon/),(/minlat,maxlat/),True)

FirstTime = True
opts = True
angle= 0.
plane = new(4,float)
plane = (/ loc(0,0)-1, loc(0,1)-1, loc(1,0)-1,loc(1,1)-1 /)


; Set some basic resources
  res = True
  res@MainTitle = ""
  pltres = True


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  FirstTime = True
  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)

;---------------------------------------------------------------
xlon = wrf_user_getvar(a, "XLONG",0)
xlat = wrf_user_getvar(a,"XLAT",0)

do it = 31,31              ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots

   ; u = wrf_user_getvar(a,"ua",it)        ; theta
   ; w = wrf_user_getvar(a,"wa",it)        ; relative humidity
    z  = wrf_user_getvar(a, "z",it)        ; grid point height
    p  = wrf_user_getvar(a, "pressure",it) ; grid point height
    dbz= wrf_user_getvar(a,"dbz",it)
    ter = wrf_user_getvar(a, "HGT",0)
     
      ;  u_plane = wrf_user_intrp3d(u,z,"v",plane,0.,opts)
      ;  w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
        dbz_plane= wrf_user_intrp3d(dbz,z,"v",plane,0.,opts)
      

      ; Let's create nice labels - only have to do this once
       if ( FirstTime ) then
          zmin = 0.
          zmax = 12.
          nz = floattoint(zmax+1)
        end if

     if ( FirstTime ) then
          zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
          b = ind(zz(:,0) .gt. zmax*1000. )
          zmax_pos = b(0) - 1
          if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
            zspan = b(0) - 1
          else
            zspan = b(0)
          end if
          delete(zz)
          delete(b)
          FirstTime = False
        end if

         ; X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
         ; X_desc = "longitude"
          X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
          X_desc = "latitude"
          dimsX = dimsizes(X_plane)
          printVarSummary(X_plane)
          xmin  = X_plane(0)
         xmax  = X_plane(dimsX-1)
          xspan = dimsX(0)-1
          nx    = floattoint( (xmax-xmin)/0.5 + 1)
      
      ; Options for XY Plots
        opts_xy                         = res
        opts_xy@tiYAxisString           = "Height (km)"
        opts_xy@AspectRatio             = 0.75
        opts_xy@cnMissingValPerimOn     = True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmXBMode                = "Explicit"
        opts_xy@tmXBValues              = fspan(0,xspan,nx)                    ; Create tick marks
        opts_xy@tmXBLabels              = sprintf("%.02f",fspan(xmin,xmax,nx))  ; Create labels
       ; opts_xy@trXReverse               = True
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
        opts_xy@tmYLLabels              = sprintf("%.02f",fspan(zmin,zmax,nz))  ; Create labels
        opts_xy@tiXAxisFontHeightF      = 0.050
        opts_xy@tiYAxisFontHeightF      = 0.050
        opts_xy@tmXBMajorLengthF        = 0.05
        opts_xy@tmYLMajorLengthF        = 0.05
        opts_xy@tmYLLabelFontHeightF    = 0.05
        opts_xy@PlotOrientation         = dbz_plane@Orientation

        opts_dbz = opts_xy
        opts_dbz@pmLabelBarOrthogonalPosF = -0.06
        opts_dbz@cnFillOn                = True

        opts_dbz@cnLevelSelectionMode = "ExplicitLevels"
        opts_dbz@cnLevels    = (/10,15,20,25,30,35,40,45,50,55,60,65,70/)
        opts_dbz@cnFillColors = (/ 0,11,10,13,16,17,18,19,20,21,22,23,4,24/)
       ; opts_dbz@lbLabelFontHeightF        =0.035
       ; opts_dbz@lbTitleFontHeightF        =0.035

        contour_dbz = wrf_contour(a,wks,dbz_plane(0:zmax_pos,:),opts_dbz)
       ; MAKE PLOTS         
        plot = wrf_overlays(a,wks,(/contour_dbz/),pltres)


       ; Delete options and fields, so we don't have carry over
        delete(opts_dbz)
        delete(dbz_plane)

      

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

end

1649518749

你确定你没剖错？我看你第一张图好像是在120度的长江入海口。

我不喜欢吃生姜

我剖的是115.95~

我不喜欢吃生姜

1649518749 发表于 2016-11-17 21:40
你确定你没剖错？我看你第一张图好像是在120度的长江入海口。

刚才的剖面图放的范围太小了，我又更新了下，在楼下~

尽头的尽头

仔细看你的红线，并没有从最大值中心穿过

葫芦娃小星星

求问楼主 你这个问题中用到的nc资料是由雷达基数据处理得到的么 如果是 麻烦能教一下我么

我不喜欢吃生姜

葫芦娃小星星 发表于 2017-4-22 09:27
求问楼主 你这个问题中用到的nc资料是由雷达基数据处理得到的么 如果是 麻烦能教一下我么

不是的，是wrfout

lsucky

色标非常漂亮呢~楼主的grb文件可以看一下吗?

lgz

楼主的grb文件可以看一下吗?