NCL绘制剖面图出错

刀客白刃

我第三重区域设置为90km*45km垂直为35层
在绘制剖面图过程中只绘制两张就出现了问题，红色为报错行
(0)        Working on time: 2014-01-10_16:00:00
(0)        Max terrain height in plot 2460.46
warning:cnMissingValPerimOn is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillColor is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillPattern is not a valid resource in uw_theta_vector at this time
warning:cnMissingValPerimOn is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillColor is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillPattern is not a valid resource in uw_theta_vector at this time
fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8640]:Execute: Error occurred at or near line 71 in file poumiantu.ncl

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"
begin
  a = addfile("/home/public2/data/18qinc/wrfout_d03_2014-01-10_00:00:00.nc","r")
  wks = gsn_open_wks("png","uw_theta")

; Set some basic resources
  res = True
  res@Footer = False

  pltres = True
  ter_res = True
  opts_ter = ter_res
  opts_ter@gsnYRefLine = 0.0
  opts_ter@gsnAboveYRefLineColor = "black"
  opts_ter@gsnDraw = False
  opts_ter@gsnFrame = False
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  FirstTime = True
  FirstTimeMap = 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)
  xlat = wrf_user_getvar(a, "XLAT",0)
  xlon = wrf_user_getvar(a, "XLONG",0)
  ter = wrf_user_getvar(a, "HGT",0)
;---------------------------------------------------------------
  do it = 16,39,1              ; TIME LOOP

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

    u = wrf_user_getvar(a,"ua",it)
;   v = wrf_user_getvar(a,"va",it)
    w = wrf_user_getvar(a,"wa",it)
    theta = wrf_user_getvar(a,"theta",it)
    z = wrf_user_getvar(a, "z",it)
    if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = 6.                          ; We are only interested in the first 6km
      nz   = floattoint(zmax + 1)
    end if
        opts_ter@trYMaxF = zmax*1000

;---------------------------------------------------------------
   do ip = 1, 3       ; we are doing 3 plots

        plane = new(2,float)
        plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /)    ; pivot point is center of domain (x,y)

        opts = False
        if(ip .eq. 1) then
          angle = 90.
          X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
          X_desc = "longitude"
        end if
        if(ip .eq. 2) then
          angle = 0.
          X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
          X_desc = "latitude"
        end if
        if(ip .eq. 3) then
          angle = 45.
          X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
          X_desc = "longitude"
        end if
        theta_plane = wrf_user_intrp3d(theta,z,"v",plane,angle,opts)
        u_plane = wrf_user_intrp3d(u,z,"v",plane,angle,opts)
;        v_plane = wrf_user_intrp3d(v,z,"v",plane,angle,opts)
        w_plane = wrf_user_intrp3d(w,z,"v",plane,angle,opts)
        ter_plane = wrf_user_intrp2d(ter,plane,angle,opts)
        print("Max terrain height in plot " + max(ter_plane))

        theta_plane2 = theta_plane
        w_plane2 = w_plane
        cross_dims = dimsizes(theta_plane2)
        rank = dimsizes(cross_dims)
        iz_do = 25
        do iz = 0,24
          iz_do = iz_do-1
          do ix = 0,cross_dims(rank-1)-1
            if ( ismissing(theta_plane2(iz_do,ix)) ) then
              theta_plane2(iz_do,ix) = theta_plane2(iz_do+1,ix)
            end if
            if ( ismissing(w_plane2(iz_do,ix)) ) then
              w_plane2(iz_do,ix) = w_plane2(iz_do+1,ix)
            end if
          end do
        end do

      ; Find the index where 6km is - only need to do this once
        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-axis lables
      dimsX = dimsizes(X_plane)
      xmin  = X_plane(0)
      xmax  = X_plane(dimsX(0)-1)
      xspan = dimsX(0)-1
      nx    = floattoint( (xmax-xmin)/2 + 1)
      ;---------------------------------------------------------------

      ; Options for XY Plots
        opts_xy                         = res
        opts_xy@tiXAxisString           = X_desc
        opts_xy@tiYAxisString           = "Height (km)"
        opts_xy@cnMissingValPerimOn     = True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmXTOn                  = False
        opts_xy@tmYROn                  = False
        opts_xy@tmXBMode                = "Explicit"
        opts_xy@tmXBValues              = fspan(0,xspan,6)                    ; Create tick marks
        opts_xy@tmXBLabels              = sprintf("%.1f",fspan(xmin,xmax,6))
        opts_xy@tmXBLabelFontHeightF    = 0.015
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,zspan,nz)                    ; Create tick marks
        opts_xy@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz))  ; Create labels
        opts_xy@tiXAxisFontHeightF      = 0.020
        opts_xy@tiYAxisFontHeightF      = 0.020
        opts_xy@tmXBMajorLengthF        = 0.02
        opts_xy@tmYLMajorLengthF        = 0.02
        opts_xy@tmYLLabelFontHeightF    = 0.015
        opts_xy@PlotOrientation         = w_plane@Orientation

      ; Plotting options for RH
        opts_theta = opts_xy
        opts_theta@pmLabelBarOrthogonalPosF = -0.18
        opts_theta@cnFillOn                = True
        opts_theta@cnLevelSelectionMode = "ExplicitLevels"
        opts_theta@cnLevels = (/286,287,288,289,290,292,294,296,300,304,308,312/)

      ; Plotting options for Temperature
        ;opts_vec = opts_xy
        vcres = opts_xy
        vcres@vcGlyphStyle       = "LineArrow"
        vcres@vcPositionMode     = "ArrowTail"
        vcres@vcRefAnnoOn        = True
        vcres@vcMinDistanceF     = 0.03
        vcres@vcRefMagnitudeF    = 8.
        vcres@vcLineArrowThicknessF = 3.0
        vcres@vcMapDirection        = False
        vcres@vcFillArrowsOn     = True
        vcres@vcRefLengthF       = 0.02
        vcres@vcFillArrowHeadXF  = 0.2
        vcres@vcFillArrowHeadYF  = 0.2
        vcres@vcFillArrowHeadInteriorXF = 0.15
        ;opts_vec@cnInfoLabelZone = 1
        ;opts_vec@cnInfoLabelSide = "Top"
        ;opts_vec@cnInfoLabelPerimOn = True
        ;opts_vec@cnInfoLabelOrthogonalPosF = -0.00005
        ;opts_vec@ContourParameters  = (/ 5. /)
;Contour terrain cross section
        contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter)

      ; Get the contour info for the rh and temp
        vector = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:)*20.,vcres)
        contour_theta = wrf_contour(a,wks,theta_plane(0:zmax_pos,:),opts_theta)
        vector2 = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane2(0:zmax_pos,:)*20.,vcres)
        contour_theta2 = wrf_contour(a,wks,theta_plane2(0:zmax_pos,:),opts_theta)
      ;---------------------------------------------------------------
  ; MAKE PLOTS
        if (FirstTimeMap) then
          lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
          lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
          mpres = True
          pltres = True
          pltres@FramePlot = False
          optsM = res
          optsM@NoHeaderFooter = True
          optsM@cnFillOn = True
          optsM@lbTitleOn = False
          optsM@cnLevelSelectionMode = "ExplicitLevels"
          optsM@cnLevels = (/1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600/)
          contour  = wrf_contour(a,wks,ter,optsM)
          plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
          lnres = True
          lnres@gsLineThicknessF = 3.0
          lnres@gsLineColor = "Red"
          gsn_polyline(wks,plot,(/lon_plane(0),lon_plane(45)/),(/30.58,36.58/),lnres)
;          gsn_polyline(wks,plot,(/103.67,103.67/),(/lat_plane(0),lat_plane(90)/),lnres)

          frame(wks)
          delete(lon_plane)
          delete(lat_plane)
          pltres@FramePlot = True
       end if
       ;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)
       ;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc,contour_ter/),pltres)
       plot = wrf_overlays(a,wks,(/contour_theta2,contour_ter,vector/),pltres)
  ; Delete options and fields, so we don't have carry over
        delete(opts_xy)
        delete(vcres)
        delete(opts_theta)
        delete(w_plane)
        delete(theta_plane)
        delete(w_plane2)
        delete(theta_plane2)
        delete(X_plane)
        delete(ter_plane)
    end do  ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    FirstTimeMap = False
  end do        ; END OF TIME LOOP
end

Adam-Lee

u_plane = wrf_user_intrp3d(v,z,"v",plane,angle,opts)
错误行改成上面这样，似乎变量写错了???

刀客白刃

Adam-Lee 发表于 2018-5-8 14:42
u_plane = wrf_user_intrp3d(v,z,"v",plane,angle,opts)
错误行改成上面这样，似乎变量写错了???

改完之后还是同样的错误？

Adam-Lee

刀客白刃 发表于 2018-5-8 14:50
改完之后还是同样的错误？

71行是哪一行？是你标红的地方吗

刀客白刃

Adam-Lee 发表于 2018-5-8 14:51
71行是哪一行？是你标红的地方吗

是的
(0)        Working on time: 2014-01-10_06:00:00
(0)        Max terrain height in plot 2338.19
warning:cnMissingValPerimOn is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillColor is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillPattern is not a valid resource in uw_theta_vector at this time
warning:cnMissingValPerimOn is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillColor is not a valid resource in uw_theta_vector at this time
warning:cnMissingValFillPattern is not a valid resource in uw_theta_vector at this time
fatal:Dimension sizes of left hand side and right hand side of assignment do not match
fatal:["Execute.c":8640]:Execute: Error occurred at or near line 71 in file poumian.ncl
报错为左右变量维数不一致

Jmy·Mavericks

学习学习，学习学习