RE: [ncl-talk] NCL not plotting first 200 m in the vertical cross section

From: Basit Khan <basit.khan_at_nyahnyahspammersnyahnyah>
Date: Thu Mar 17 2011 - 23:31:43 MDT

Hi Mary,
In fact the problem was not resolved and it still exists. When you spotted the error in my code i thought it was the cause of not plotting first 200 m but even after fixing the bug, the problem remained there. With further investigation i found that function wrf_user_intrp3d was playing up. Vertical levels in my wrf output starts from 6 m agl (terrain following sigma level)but after interpolation the first level i got was at 200 m. At that time I had to complete some analysis so i used GrADS instead for vertical cross sections. But it is good that this issue popped-up again and i got a chance to mention this problem.

Cheers,

Basit A. Khan
Centre for Atmospheric Research
Department of Geography
University of Canterbury, Private Bag 4800
Christchurch 8140, New Zealand.
Tel: +643 3642987 ext. 3088
Email. basit.khan@canterbury.ac.nz

-----Original Message-----
From: ncl-talk-bounces@ucar.edu on behalf of Mary Haley
Sent: Thu 3/17/2011 3:38 AM
To: Song Guiting
Cc: ncl forum
Subject: Re: NCL not plotting first 200 m in the vertical cross section
 
Guiting,

The problem in Basit's case, if I remember correctly, is that all the data were missing in that location, so the plot was being drawn correctly.

If you think NCL is not plotting your data correctly, then you need to provide us with your script and any necessary data files. Dennis Shea has already asked for these offline.

Thanks,

--Mary

On Mar 14, 2011, at 2:00 AM, Song Guiting (Dr) wrote:

> Dear Mary,
>
> I also encounter the same problem. My NCL version is 6.00 beta.
>
> Best regards,
> Guiting
>
>
> From: ncl-talk-bounces@ucar.edu [mailto:ncl-talk-bounces@ucar.edu] On Behalf Of Mary Haley
> Sent: Monday, February 14, 2011 11:18 PM
> To: Basit Khan
> Cc: ncl-talk@ucar.edu
> Subject: Re: NCL not plotting first 200 m in the vertical cross section
>
> This might be due to a bug we fixed awhile back. What version of NCL (ncl -V) are you running?
> If you can provide me with the file, then I can try your script and see if the bug is fixed.
>
> --Mary
>
>
>
> On Feb 13, 2011, at 3:15 PM, Basit Khan wrote:
>
>
> Hi,
> I am trying to plot vertical cross section (height x latitude) from a WRF output file. The problem is that NCL not plotting anything in the first 200 m while I have at least 11 vertical levels in the first 200 m asl. I am using terrain following coordinates system and most of the cross section area is either sea or terrain less than 50 m high asl. I will greatly appreciate if somebody could help me fix this problem. The ncl code is copied below and the output of the ncl code is attached.
>
> Regards,
>
> Basit
> Centre for Atmospheric Research
> Department of Geography,
> University of Canterbury, Private Bag 4800,
> Christchurch 8140, New Zealand.
> Tel: +643 364 2987 ext. 3088
> email: basit.khan@canterbury.ac.nz
>
> -------------------------------------------------------------------------------
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
>
> begin
> ;
> a = addfile("/hpc/bluefern/WRF/basit/basit-archive/b-haze/bh30/wrfout_d03_2006-05-30_12: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,"plt_CrossSection4")
>
>
> ; Set some basic resources
> res = True
> res@MainTitle = "REAL-TIME WRF"
> res@Footer = False
>
> pltres = True
>
> ;;;;;;;;;;;;;;;; HIGHRES RANGS/GSHHS COASTLINE PLOT RESOURCES ;;;;;;;;;;;;;;;;;;;;;;;;
> mpres =True
> mpres@gsnMaximize = True
> mpres@mpOutlineOn = True
> mpres@mpFillOn = False
> mpres@mpProjection = "Mercator"
> mpres@mpLimitMode = "LatLon"
> mpres@mpMinLonF = 174.4061
> mpres@mpMaxLonF = 175.2502
> mpres@mpMinLatF = -37.27655
> mpres@mpMaxLatF = -36.60195
> mpres@mpOutlineOn = True
> mpres@mpGridLineColor = -1
> mpres@gsnDraw = False ; don't draw the plots
> mpres@gsnFrame = False ; don't advance the frame
> mpres@tmYROn = False
> mpres@tmXTOn = False
> mpres@mpGeophysicalLineColor= "Red"
> mpres@mpLimbLineColor = "Blue"
> mpres@mpGeophysicalLineThicknessF = 1.5
> mpres@mpDataBaseVersion = "HighRes" ; Use the high-res database
> mpres@pmTickMarkDisplayMode = "Always" ; Turn on map tickmarks.
> ; mpres@tiMainString = "Using the RANGS-GSHHS coastal database"
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> FirstTime = True
> FirstTimeMap = True
> times = wrf_user_list_times(a) ; 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; print(mdims)=97x50x150x150
> nd = dimsizes(mdims) ; incase of 'P' the print(nd)=4 cause P has four dims
>
> xlat = wrf_user_getvar(a, "XLAT",0)
> xlon = wrf_user_getvar(a, "XLONG",0)
> ter = wrf_user_getvar(a, "HGT",0)
>
> ;---------------------------------------------------------------
>
> do it = 58,74,1 ; TIME LOOP
>
> print("Working on time: " + times(it) )
> res@TimeLabel = times(it) ; Set Valid time to use on plots
>
> tc = wrf_user_getvar(a,"tc",it) ; T in C
> rh = wrf_user_getvar(a,"rh",it) ; relative humidity
> z = wrf_user_getvar(a, "height",it) ; grid point height
>
> if ( FirstTime ) then ; get height info for labels
> zmin = 0
> zmax = 1. ; We are only interested in the first 1000 m
> nz = floattoint(zmax + 10)
> end if
>
>
> ;---------------------------------------------------------------
>
> plane = new(2,float)
> plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y)
> opts = False
>
>
> angle = 0.
> X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
> X_desc = "latitude"
>
> rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts) ; "v" is for vertical and "h" is for horinzontal interporlation
> tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) ; Array of vertical coordninates to interpolate
>
>
> ; 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) ; printVarSummary(zz) = dimsizesV x Dimsize H = 100 x 149
>
> b = ind(zz(:,0) .gt. zmax*1000. ) ; (zmax *1000 --> conver km into metres
>
> 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 + 10)
> ; nx = (xmax-xmin)/2
> ;---------------------------------------------------------------
> ; print("the xmin value is "+xmin)
> ; print("the xmax value is "+xmax)
>
> ; 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@tmXBValues = fspan(0,xspan,nx) ; Create tick marks
> opts_xy@tmXBLabels = sprintf("%6.2f",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("%3.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 = tc_plane@Orientation
>
>
> ; Plotting options for RH
> opts_rh = opts_xy
> opts_rh@ContourParameters = (/ 10., 90., 10. /)
> 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 Temperature
> opts_tc = opts_xy
> opts_tc@cnInfoLabelZone = 1
> opts_tc@cnInfoLabelSide = "Top"
> opts_tc@cnInfoLabelPerimOn = True
> opts_tc@cnInfoLabelOrthogonalPosF = -0.00005
> opts_tc@ContourParameters = (/ 1. /)
>
>
> ; Get the contour info for the rh and temp
> contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc)
> contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh)
>
> ;---------------------------------------------------------------
>
> ; 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_rh,contour_tc/),pltres) ; plot x-section
>
> ; Delete options and fields, so we don't have carry over
> delete(opts_xy)
> delete(opts_tc)
> delete(opts_rh)
> delete(tc_plane)
> delete(rh_plane)
> delete(X_plane)
>
> end do ; make next cross section
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> FirstTimeMap = False
> ; end do ; END OF TIME LOOP
>
> end
>
> This email may be confidential and subject to legal privilege, it may
> not reflect the views of the University of Canterbury, and it is not
> guaranteed to be virus free. If you are not an intended recipient,
> please notify the sender immediately and erase all copies of the message
> and any attachments.
>
> Please refer to http://www.canterbury.ac.nz/emaildisclaimer for more
> information.
>
> <vcs_ncl.eps>_______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
> CONFIDENTIALITY: This email is intended solely for the person(s) named and may be confidential and/or privileged. If you are not the intended recipient, please delete it, notify us and do not copy, use, or disclose its content. Thank you.
>
> Towards A Sustainable Earth: Print Only When Necessary
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

This email may be confidential and subject to legal privilege, it may
not reflect the views of the University of Canterbury, and it is not
guaranteed to be virus free. If you are not an intended recipient,
please notify the sender immediately and erase all copies of the message
and any attachments.

Please refer to http://www.canterbury.ac.nz/emaildisclaimer for more
information.
Received on Thu Mar 17 23:31:58 2011

This archive was generated by hypermail 2.1.8 : Wed Mar 23 2011 - 16:15:59 MDT