Re: Setting cross section y axis range

From: Don Morton <Don.Morton_at_nyahnyahspammersnyahnyah>
Date: Mon Apr 25 2011 - 19:50:33 MDT

Hi Carl and NCL Colleagues,

I decided to take an interest in this because I'm just now trying to create
some better x-section plots for your area. I had to deal with this issue in
a Savoonga project last year. The script I'm pasting in is a little
different, but I think it might have the same issue. In my script (which I
stole from one of the NCL examples), zspan gets set to the top of the
"desired" scale, but then, for some reason, later on it gets reset to the
top of the original input data (or top of atmosphere, I guess). This value
ends up being used in the option opts_xy@tmYLValues.

If I comment out the second setting of zspan, then it reverts to the
pre-selected top of graph and appears to work fine.

Then, looking at your script, I see that the same thing is happening - that
it seems to be setting the scale to be over the entire span and not a
restricted span.

opts_xy@tmYLValues = ispan(0,100,10) ; over 0-100 span of Y
label

Fixing the issue "seems" to work effectively for me, though I don't feel
100% comfortable altering an official example NCL script. I recognize I
"could" be wrong in what I'm doing, but it might warrant a closer look.
Hopefully an NCL cop will jump in and tell us if I've erred here :)

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

; Include files
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

    ;;; These values are obtained from pre-set environment variables

;;;;;;;; These are used for testing - comment out when operational ;;;;;;;
; wrfoutFile =
"/wrkdir/morton/Savoonga-Study/Initial_AWIPAK/WRFRun/wrfout_d03_2009-02-27_11:30:
00.nc"
    wrfoutFile =
"/projects/ARSCWTHR/Operational/HRRR-AK/Products/wrfout/2011042518/wrfout_d01_2011-04-26_18:00:
00.nc"
    imageName = "XC" ; File name of image file (without ext)
    GRAPHICS_TYPE = "x11" ; [eps|x11|...]
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;; These are used for operational - comment out when testing ;;;;;;;
; wrfoutFile = getenv("wrfoutFile") ; Full path to the wrfout file
; imageName = getenv("imageName") ; File name of image file (without
ext)
; GRAPHICS_TYPE = "ps" ; [eps|x11|...]
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ;; Define the transect line lat/lon
    START_LAT = 62.5
    START_LON = -170.48
    END_LAT = 64.0
    END_LON = -170.48

    ;; Define top of vertical cross section (in km)
    XC_TOP_HGT = 2.0

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

    ; Need to open the input file to convert lat/lon to x/y
    f = addfile(wrfoutFile,"r") ; Opens the netCDF file

    ; Convert transect line coordinates to x/y
    xy_point = wrf_user_ll_to_ij(f, START_LON, START_LAT, False)
    STARTX = xy_point(0)
    STARTY = xy_point(1)

    xy_point = wrf_user_ll_to_ij(f, END_LON, END_LAT, False)
    ENDX = xy_point(0)
    ENDY = xy_point(1)

    ; Get the timestamp
    times = chartostring(f->Times)
    timestamp = times(0)

    ;;;;;;;;;; Set up the workstation ;;;;;;;;;;
    wks = gsn_open_wks(GRAPHICS_TYPE,imageName) ; Open workstation for
plotting

    gsn_define_colormap(wks,"wh-bl-gr-ye-re") ; Define the colormap

    resPlot = True
    resPlot@MainTitle = timestamp
    resPlot@Footer = False
    resPlot@ValidTime = False
    resPlot@InitTime = False

    tc = wrf_user_getvar(f, "tc", 0)
    u = wrf_user_getvar(f, "U", 0)
    v = wrf_user_getvar(f, "V", 0)

    z = wrf_user_getvar(f, "z", 0)
    xlat = wrf_user_getvar(f, "XLAT", 0)

    zmin = 0.0
    zmax = XC_TOP_HGT
    nz = floattoint(zmax + 1)

    ; Define the x-section vertical plane
    plane = new(4, float)
    plane = (/ STARTX, STARTY, ENDX, ENDY /) ; start x,y and end x,y
    opts = True
    ; Find the index, zmax_pos, where the top of the vertical xc is
    zz = wrf_user_intrp3d(z,z,"v",plane,0,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

;print(zmin)
;print(zmax)
;print(nz)
;print(zz(:,0))
;print(b)
;print(zmax_pos)
print(zspan)
;exit

    delete(zz)
    delete(b)

    tc_plane = wrf_user_intrp3d(tc, z, "v", plane, 0.0, opts)
    u_plane = wrf_user_intrp3d(u, z, "v", plane, 0.0, opts)
    v_plane = wrf_user_intrp3d(v, z, "v", plane, 0.0, opts)
    spd_plane = sqrt( u_plane*u_plane + v_plane*v_plane )
    spd_plane@units = "m/s"
    spd_plane@description = "Wind Speed"

    dim = dimsizes(tc_plane) ; Find the data span for use in labels

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;; COMMENTED THIS OUT - NOT SURE WHAT IT'S DOING IN HERE
;;;;;;;;; IT WAS OVER-RIDING THE PREVIOUSLY CALCULATED VALUE OF zspan
    ;zspan = dim(0)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ; set up x-axis labels (latitude in this case)
    X_plane = wrf_user_intrp2d(xlat, plane, 0.0, False)
    X_desc = "latitude"
    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 = resPlot
    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)
    opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx))
; opts_xy@tmXBLabelFontHeightF = 0.015
    opts_xy@tmYLMode = "Explicit"
    opts_xy@tmYLValues = fspan(0,zspan,nz)
    opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz))
; 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 = tc_plane@Orientation

  ; Plotting options for Temperature
    opts_tc = opts_xy
    opts_tc@cnInfoLabelOrthogonalPosF = 0.00
    opts_tc@ContourParameters = (/ 5. /)

    ; Plotting options for wind speed
    resSpd = opts_xy
    resSpd@cnFillOn = True
    resSpd@cnFillMode = "RasterFill"
    resSpd@ContourParameters = (/0., 50., 1./)

    ; Get the contour info for the rh and temp
    contour_tc = wrf_contour(f,wks,tc_plane(0:zmax_pos,:),opts_tc)
    contourSpd = wrf_contour(f, wks, spd_plane(0:zmax_pos,:), resSpd)

  ; MAKE PLOTS
    plot = wrf_overlays(f,wks,(/contourSpd,contour_tc/), resPlot)

end

-- 
Voice:  907 450 8679
Arctic Region Supercomputing Center
http://weather.arsc.edu/
http://www.arsc.edu/~morton/

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

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