NCL, problem to plot lat/long as x-axis in crosssections

The NCL graphics package.

NCL, problem to plot lat/long as x-axis in crosssections

Postby anisetty » Wed Apr 27, 2016 3:54 am

Hello,
I am new user of NCL. I could able to plot vertical cross sections. At the x-axix, it is printing grids instead of latitude or longitude. Can some one correct my script or suggest it. I tried wrf_user_ij_to_ll, but did not succeed. Please any of you help me. my script is attached below
Thank you very much in advance

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

wrffiles = systemfunc("ls wrfout_d01_2010-10-15_1*")
numFiles = dimsizes(wrffiles)
do i = 0, numFiles -1
wrffiles(i) = wrffiles(i)+".nc"
end do
inpFiles = addfiles(wrffiles,"r")

; We generate plots, but what kind do we prefer?
type = "x11"
type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"plt_CrossSectionW-U")


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

pltres = True


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

FirstTime = True
do ifile = 0, numFiles-1
a = inpFiles[ifile]
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
nd = dimsizes(mdims)

;---------------------------------------------------------------
do it = 0,ntimes-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,"ua",it) ; zonal wind at mass point
rh = wrf_user_getvar(a,"wa",it) ; vertical velocity at mass point
z = wrf_user_getvar(a, "z",it) ; grid point height

if ( FirstTime ) then ; get height info for labels
zmin = 0.
zmax = max(z)/1000.
nz = floattoint(zmax/2 + 1)
FirstTime = False
end if

; ajout pour moi zmax et plot birmanie
zmax=10.

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

ip = 1 ; Just do the one (constant y coord) plot
opts = True ; setting start and end times
plane = new(4,float)

if(ip .eq. 1) then
plane = (/ 0,84, 200,84 /) ; start x;y & end x;y point
; The jangmi data is on a grid from 0 to 200
end if
if(ip .eq. 2) then
plane = (/ 130,1, 130,162 /) ; start x;y & end x;y point
end if
if(ip .eq. 3) then
plane = (/ 49,1, 210,162 /) ; start x;y & end x;y point
end if

; plan N-S entre y=30 et y=80 pour x fixe a 20
plane = (/ 100,100, 250,250 /) ; start x;y & end x;y point
; plan W-E entre x=10 et x=40 pour y fixe a 65
; plane = (/ 10,65, 40,65/) ; start x;y & end x;y point


rh_plane = wrf_user_intrp3d(rh,z,"v",plane,0.,opts)
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,0.,opts)

dim = dimsizes(rh_plane) ; Find the data span - for use in labels
zspan = dim(0)

; Options for XY Plots
opts_xy = res
opts_xy@tiYAxisString = "Height (km)"
opts_xy@AspectRatio = 0.75
opts_xy@cnMissingValPerimOn = True
opts_xy@cnMissingValFillColor = 0
opts_xy@cnMissingValFillPattern = 11
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 = tc_plane@Orientation


; Plotting options for RH
opts_rh = opts_xy
opts_rh@pmLabelBarOrthogonalPosF = -0.07
opts_rh@ContourParameters = (/ -4.5, 4.5, .1 /)
opts_rh@cnFillOn = True
; opts_rh@cnFillColors = (/"White","White","White", \
; "White","Chartreuse","Green", \
; "Green3","Green4", \
; "ForestGreen","PaleGreen4","Red"/)

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


; Get the contour info for the rh and temp
contour_tc = wrf_contour(a,wks,tc_plane,opts_tc)
contour_rh = wrf_contour(a,wks,rh_plane,opts_rh)

; MAKE PLOTS
plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)

; Delete options and fields, so we don't have carry over
delete(opts_tc)
delete(opts_rh)
delete(tc_plane)
delete(rh_plane)

; end do ; make next cross section (loop over ip was removed)

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

end do ; END OF TIME LOOP

end do ; END OF FILE LOOP

end
anisetty
 
Posts: 6
Joined: Mon Oct 05, 2015 11:25 am

Return to NCL

Who is online

Users browsing this forum: No registered users and 2 guests