Cross Section with Wind Vectos

The NCL graphics package.

Cross Section with Wind Vectos

Postby Behrouz » Fri Apr 20, 2012 9:30 am

Hi

I've tried to study on mountainous terrain circulations and I need to plot Cross sections with Wind Vector and W as a contour. has any one worked with similar script? is there any NCL script to recommend?
thanks.
Behrouz
 
Posts: 2
Joined: Thu Jul 21, 2011 10:50 am

Re: Cross Section with Wind Vectos

Postby Orome » Mon Apr 23, 2012 6:08 am

I was about to ask a very similar question - so instead of making another thread I'll just ask here.

Would it be possible for someone to check my code below for plotting a cross section including vertical velocity. The graph made at the end looks ok, but am totally unsure whether it is right or not, and whether it is actually plotting what I want it too!

Thanks

Code: Select all
;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   Plot data on a cross section
;   This script will plot data at a set angle through a specified point
;   This script adds lon/lat info along X-axis

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

begin
;
; The WRF ARW input file. 
; This needs to have a ".nc" appended, so just do it.
  a = addfile("../WRFV3/test/em_real/wrfout_d03_2011-02-09_00:00:00.nc","r")


; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"CrossSection_d03_RH")
  gsn_define_colormap(wks,"gui_default")

; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"
  res@Footer = False
 
  pltres = True


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  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 = 0,ntimes-1,2             ; TIME LOOP

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

    th   = wrf_user_getvar(a,"th",it)      ; potentialT in C
    rh   = wrf_user_getvar(a,"rh",it)      ; relative humidity
    z    = wrf_user_getvar(a, "z",it)      ; grid point height
    pvo = wrf_user_getvar(a,"pvo",it)                     
    u   = wrf_user_getvar(a,"ua",it)                         
    v   = wrf_user_getvar(a,"va",it)                       
    p   = wrf_user_getvar(a, "pressure",it)                 
    w   = wrf_user_getvar(a,"wa",it)
    qv  = wrf_user_getvar(a,"QVAPOR",it)   ; Qv

    if ( FirstTime ) then                ; get height info for labels
      zmin = 0.
      zmax = 16.                          ; We are only interested in the first 6km
      nz   = floattoint(zmax + 1)
    end if


;---------------------------------------------------------------

    do ip = 1, 3         ; we are doing 3 plots
         ; all with the pivot point (plane) in the center of the domain
         ; at angles 0, 45 and 90
 ; 
 ;                   |
 ;       angle=0 is  |
 ;                   |
 ;

        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

        rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)
        th_plane = wrf_user_intrp3d(th,p,"v",plane,angle,opts)
   qv_plane = wrf_user_intrp3d(qv,z,"v",plane,angle,opts)
   u_plane  = wrf_user_intrp3d(u,p,"v",plane,0.,opts)
        v_plane  = wrf_user_intrp3d(v,p,"v",plane,0.,opts)
        w_plane  = wrf_user_intrp3d(w,p,"v",plane,0.,opts)

   qv = qv*1000
   qv@units = "g/kg"

         u_plane = u_plane*1.94386                  ; kts conversion
         v_plane = v_plane*1.94386                   ; kts conversion
      ;  w_plane = w_plane*1.94386

         u_plane@units = "kts"
         v_plane@units = "kts"
    w_plane@units = "kts"

      ; 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,nx)                    ; Create tick marks
        opts_xy@tmXBLabels              = sprintf("%.1f",fspan(xmin,xmax,nx))  ; Create labels
        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         = th_plane@Orientation


      ; Plotting options for RH
        opts_rh = opts_xy
        opts_rh@ContourParameters       = (/ 10., 90., 5. /)
        opts_rh@pmLabelBarOrthogonalPosF = -0.1
        opts_rh@cnFillOn                = True
       ; opts_rh@cnFillColors            = (/"White","White","White", \
       ;                                     "White","Chartreuse","Green", \
       ;                                     "Green3","Green4", \
        ;                                    "ForestGreen","PaleGreen4"/)

      ; Plotting options for Qv
   opts_qv = opts_xy
   opts_qv@ContorParameters   = (/1., 100., 10. /)
   opts_qv@pmLabelBarOrthogonalPosF = -0.1
   opts_qv@cnFillOn      = True


      ; Plotting options for Potential Temperature
        opts_th = opts_xy
        opts_th@cnInfoLabelZone = 1
        opts_th@cnInfoLabelSide = "Top"
        opts_th@cnInfoLabelPerimOn = True
        opts_th@cnInfoLabelOrthogonalPosF = -0.00005
        opts_th@ContourParameters  = (/ 1. /)


     ; Vertical Velocity
       w_plane = 100.*w_plane
       opts_w   = res
       ;opts_w@FieldTitle            = w@description
       opts_w@UnitLabel             = "m/3"
       ;opts_w@PlotLevelID           = 0.001*height + " km"
       ;opts_w@cnFillOn              = True
       opts_w@gsnSpreadColorEnd     = -3
       contour_w = wrf_contour(a,wks, w_plane,opts_w)


      ;Plotting options for Wind Vectors

        opts_w = res         
       opts_w@FieldTitle = "Wind"        ; Overwrite Field Title
       opts_w@NumVectors = 27     ; Wind barb density
         

      ; Get the contour info for the rh and temp
        contour_th = wrf_contour(a,wks,th_plane(0:zmax_pos,:),opts_th)
        contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh)
   contour_qv = wrf_contour(a,wks,qv_plane(0:zmax_pos,:),opts_qv)
        vector = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:),opts_w)

      ;---------------------------------------------------------------

  ; 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
          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"
          do ii = 0,dimsX(0)-2
            gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres)
          end do
          frame(wks)
          delete(lon_plane)
          delete(lat_plane)
          pltres@FramePlot = True
       end if

     plot = wrf_overlays(a,wks,(/contour_qv,contour_th,vector/),pltres)    ; plot x-section

;plot = wrf_overlays(a,wks,(/contour_th,contour_rh/),pltres)   

  ; Delete options and fields, so we don't have carry over
        delete(opts_xy)
        delete(opts_th)
        delete(opts_rh)
   delete(opts_qv)
   delete(qv_plane)
        delete(th_plane)
        delete(rh_plane)
   delete(opts_w)
   delete(u_plane)
   delete(v_plane)
   delete(w_plane)
        delete(X_plane)

    end do  ; make next cross section

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

    FirstTimeMap = False
  end do        ; END OF TIME LOOP

end
Orome
 
Posts: 35
Joined: Mon Dec 05, 2011 6:09 am

Re: Cross Section with Wind Vectos

Postby bingxuewu » Tue Aug 04, 2020 1:32 am

I think when the angle is 45,u_plane(0:zmax_pos,:)need to change, like u_plane=u*cos()+vu*cos()
:?
vector = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:),opts_w)
bingxuewu
 
Posts: 1
Joined: Tue Aug 04, 2020 1:08 am


Return to NCL

Who is online

Users browsing this forum: No registered users and 6 guests