Setting cross section y axis range

From: Carl F Dierking <carl.dierking_at_nyahnyahspammersnyahnyah>
Date: Mon Apr 25 2011 - 18:22:25 MDT

  I have modified one of the cross section example scripts to display a
WRF cross section of theta along a line between two lat/lon points. The
vertical coordinate is height in km. Currently the script displays the
entire vertical depth of all WRF levels, but I'd like to display only
the lowest 4 km or so. What changes do I need to make to do that?

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Script to plot WRF forecast theta on a cross section
; This script will plot data from a a given point A to point B
; Vertical coordinate is height

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("/media/disk-1/taku_00011612/wrfout_d03_2000-01-16_21: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,"CrossSection2")
   gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

; Set some Basic Plot options
   ARWres = True
   ARWres@MainTitle = "REAL-TIME WRF"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
   FirstTime = True
   times = wrf_user_list_times(a) ; get times in the file
   ntimes = dimsizes(times) ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   do it = 0,ntimes-1,2 ; TIME LOOP
     print("Working on time: " + times(it) )
     ARWres@TimeLabel = times(it) ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
     w = wrf_user_getvar(a,"wa",it) ; T in C
     th = wrf_user_getvar(a,"th",it) ; relative humidity
     z = wrf_user_getvar(a, "z",it) ; grid point height
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     if (FirstTime) then
       ; THIS IS NEEDED FOR Y LABLES - ALWAYS DO
       zz = z/100.
       dims = dimsizes(zz)
       zmax = -999.
       do imax = 0,dims(0)-1
         if ( .not.ismissing(zz(imax,1,1)) .and. zz(imax,1,1) .gt. zmax
) then
           zmax = zz(imax,1,1)
         end if
       end do
       imax = floattoint(zmax)
       zmax = imax/10.
       print("Y labels set to: 0 - " + zmax + " km")
       FirstTime = False
       ; END OF ALWAYS DO
     end if

     do ip = 1, 1
         opts = True
         plane = new(4,float)

         if(ip .eq. 1) then ; Specify the start/end points of the
cross section in lat/lon
           lats = (/58.20, 58.38/) ;
           lons = (/-134.46, -134.18/) ;
           loc = wrf_user_ll_to_ij(a, lons, lats, True)
           ; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y)
           ; subtract one since we want to use it as an index in NCL
           x_start = loc(0,0) - 1
           x_end = loc(0,1) - 1
           y_start = loc(1,0) - 1
           y_end = loc(1,1) - 1
           plane = (/ x_start,y_start , x_end,y_end /) ; approx. start
x;y and end x;y point
         end if
         if(ip .eq. 2) then
           plane = (/ 130,1 , 130,162 /) ; approx. start x;y and end
x;y point
         end if
         if(ip .eq. 3) then
           plane = (/ 49,1 , 210,162 /) ; approx. start x;y and end
x;y point
         end if

         ; Interpolate data vertically (in z)
         th_plane = wrf_user_intrp3d(th,z,"v",plane,0.,opts)
         w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)

   ; Options for XY Plots
         opts_xy = ARWres
         opts_xy@tiXAxisString = "Number of Grid Points"
         opts_xy@tiYAxisString = "Height (km)"
         opts_xy@AspectRatio = .75
         opts_xy@cnMissingValPerimOn = True
         opts_xy@cnMissingValFillColor = 0
         opts_xy@cnMissingValFillPattern = 11
         opts_xy@tmYLMode = "Explicit"
         opts_xy@tmYLValues = ispan(0,100,10) ; over 0-100
span of Y label
         opts_xy@tmYLLabels = fspan(0,zmax,11);
Corresponding 11 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 theta
         opts_th = opts_xy
         opts_th@cnLevelSelectionMode = "ManualLevels" ;set manual
contour lvls
         opts_th@cnMinLevelValF = 260. ; set min
contour level
         opts_th@cnMaxLevelValF = 290. ; set max
contour level
         opts_th@cnLevelSpacingF = 1 ; set contour spacing
         opts_th@cnFillOn = True
         opts_th@gsnSpreadColors = True

   ; Plotting options for vertical vel
         opts_w = opts_xy
         opts_w@ContourParameters = (/ 1. /)
   ; Get the contour info for the th and temp
         contour_th = wrf_contour(a,wks,th_plane,opts_th)
         contour_w = wrf_contour(a,wks,w_plane,opts_w)
   ; MAKE PLOTS
         wrf_overlay(wks,(/contour_th,contour_w/),True)
   ; Delete options and fields, so we don't have carry over
         delete(opts_w)
         delete(opts_th)
         delete(w_plane)
         delete(th_plane)
     end do ; make next cross section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   end do ; END OF TIME LOOP
end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Apr 25 18:22:41 2011

This archive was generated by hypermail 2.1.8 : Tue May 03 2011 - 14:47:35 MDT