[ncl-talk] plot longitude in a cross section
Jesús Garcia Rosales
jesus21gr at gmail.com
Wed Mar 30 07:37:54 MDT 2016
Hi ncl users,
I have a problem when I want to plot longitude in axis x in my cross
section. Please help me.
The error that appears is:
fatal:fspan: number of elements parameter is less-than-or-equal-to one,
can't continue
fatal:["Execute.c":8575]:Execute: Error occurred at or near line 281 in
file prueba.ncl
This is my script:
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
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/usuario/TOSHIBA/WRF_alan/9km/exp_re/
wrfout_2002F.cu_betts.3km.nc","r")
; We generate plots, but what kind do we prefer?
type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
; type = "png"
wks = gsn_open_wks(type,"cross_section_2")
;gsn_define_colormap(wks,"BlueWhiteOrangeRed")
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
; wks2 = gsn_open_wks(type,"plan_map")
;;;;;;;;;; Dominio
; si se quiere disminuir el dominio
bordh=8 ;8 ;h=haut(hight)
bordb= 20 ;20 ;b=bas(bajo)
bordd= 15 ;15 ;d=droite(derecha)
bordg= 17;17 ;g=gauche (izquierda)
tmin= 72 ; la seleccion empieza el 01/01/2011 a las 00hs (UTC)
tmax= 744 ; la seleccion se termina el 31/01/2011 a las 21hs (UTC)
; Set some basic resources
res = True
res at MainTitle = "Hurricane Irene 27 km WRF (P1)"
res at Footer = False
pltres = True
; Set the path the cross section will take
lat_start = -12.
lon_start = -76.00
lat_end = -12.
lon_end = -74.5
; Limited vertical extent? (True or False)
LimitedExtent = False
z_top = 10. ; top in km (ignored if above is False)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times = wrf_user_getvar(a,"times",-1) ; get times in the file
ntimes = dimsizes(times) ; number of times in the file
FirstTime = True
mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
nd = dimsizes(mdims)
;---------------------------------------------------------------
FirstTimeMap = True
; do it = 0,ntimes-1,72 ; TIME LOOP
; print("Working on time: " + times(it) )
; res at TimeLabel = times(it) ; Set Valid time to use on plots
; tc = wrf_user_getvar(a,"tc",it) ; temperature (deg C)
; rh = wrf_user_getvar(a,"rh",it) ; relative humidity (%)
z = wrf_user_getvar(a, "z",0) ; grid point height (m)
p = wrf_user_getvar(a,"pressure",0) ; full model pressure (hPa)
u = wrf_user_getvar(a,"ua",72) ; x-component wind (m/s)
v = wrf_user_getvar(a,"va",72) ; y-component wind (m/s)
lat = wrf_user_getvar(a,"lat",0) ; latitude grid
lon = wrf_user_getvar(a,"lon",0) ; longitude grid
ter = wrf_user_getvar(a,"HGT",0) ; model terrain height (m)
w = wrf_user_getvar(a,"wa",72)
q= wrf_user_getvar(a,"QVAPOR",72)*1000
; Compute the relative vorticity from the u and v wind fields
printVarSummary(p)
printVarSummary(u)
printVarSummary(lat)
printVarSummary(w)
;tim_u=u(tmin:tmax,:,:,:)
;tim_v=v(tmin:tmax,:,:,:)
;printVarSummary(tim_u)
;vm=dim_avg_n(v,0)
;um=dim_avg_n(u,0)
;wm=dim_avg_n(w,0)
; Sacando los promedios de las varaibles u,v,w
;um=dim_avg_n(u(tmin:tmax,:,:,:),0)
;vm=dim_avg_n(v(tmin:tmax,:,:,:),0)
;wm=dim_avg_n(w(tmin:tmax,:,:,:),0)
;qm=dim_avg_n(q(tmin:tmax,:,:,:),0)
vort = uv2vr_cfd(u,v,lat(:,0),lon(0,:),2)
vort at description = "Relative Vorticity"
vort at units = "1/s"
if ( FirstTime ) then ; get height info for labels
zmin = 0. ; bottom of plot
if ( LimitedExtent ) then
angle = 0.
zmax = z_top
nz = floattoint(zmax + 1)
else
angle = 0.
zmax = max(z) / 1000.
nz = floattoint(zmax / 2 + 1)
FirstTime = False
end if
end if
;---------------------------------------------------------------
; Plot a cross session that run from point A to point B
; starting point for cross-section
llres = True
llres at ReturnInt = False
locij_a = wrf_user_ll_to_ij(a, lon_start, lat_start, llres)
locij_a = locij_a - 1
; array pointers in NCL space
locX_a = locij_a(0)
locY_a = locij_a(1)
; ending point for cross-section
locij_b = wrf_user_ll_to_ij(a, lon_end, lat_end, llres)
locij_b = locij_b - 1
; array pointers in NCL space
locX_b = locij_b(0)
locY_b = locij_b(1)
plane = new(4,float)
plane = (/ locX_a,locX_b , locY_a,locY_b /) ; start x;y & end
x;y point
opts = True ; start and end
points specified
;rh_plane = wrf_user_intrp3d(rh,z,"v",plane,0.,opts) ; relative
humidity plane
;tc_plane = wrf_user_intrp3d(tc,z,"v",plane,0.,opts) ; temperature
plane
p_plane = wrf_user_intrp3d(p,z,"v",plane,0.,opts) ; pressure plane
vort_plane = wrf_user_intrp3d(vort,z,"v",plane,0.,opts) ; relative
vorticity plane
u_plane = wrf_user_intrp3d(u,z,"v",plane,0.,opts)
w_plane = wrf_user_intrp3d(w,z,"v",plane,0.,opts)
q_plane = wrf_user_intrp3d(q,z,"v",plane,0.,opts)
printVarSummary(u_plane)
;************************************************
; Omega is significantly smaller than v, so we will
; scale it so that some vertical motion is visible
;************************************************
wAve = avg(w_plane(:,:)) ; used for scaling
uAve = avg(u_plane(:,:))
scale = fabs(uAve/wAve)
w_plane = w_plane*scale ; now scale
; w_plane=w_plane*100
vort_plane = vort_plane * 100000. ; normalize vorticity to 10^-5
s^-1
vort_plane at units = "10^5 s^-1"
dim = dimsizes(vort_plane) ; Find the data
span - for use in labels
zspan = dim(0)
if ( FirstTime ) then
angle = 0.
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
printVarSummary(zmax_pos)
else
zspan = b(0)
end if
delete(zz)
delete(b)
FirstTime = False
end if
x_plane = wrf_user_intrp2d(lon,plane,angle,opts)
x_plane at description = "Longitude"
dimsX = dimsizes(x_plane)
xmin = x_plane(0)
xmax = x_plane(dimsX(0)-1)
xspan = dimsX(0)-1
nx = floattoint((xmax-xmin)/2+1)
; printVarSummary(zmax_pos)
; Options for XY Plots
opts_xy = res
opts_xy at tiYAxisString = "Height (km)"
opts_xy at cnMissingValPerimOn = True
opts_xy at cnMissingValFillColor = 0
opts_xy at cnMissingValFillPattern = 11
opts_xy at tmYLMode = "Explicit"
opts_xy at tmYLValues = fspan(0,zspan,nz)
; Create tick marks
opts_xy at tmYLLabels =
sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels
opts_xy at tiXAxisFontHeightF = 0.020
opts_xy at tiYAxisFontHeightF = 0.020
opts_xy at tmXBMajorLengthF = 0.02
opts_xy at tmYLMajorLengthF = 0.02
opts_xy at tmYLLabelFontHeightF = 0.015
opts_xy at PlotOrientation = p_plane at Orientation
printVarSummary(xspan)
printVarSummary(nx)
opts_xy at tmXBMode = "Explicit"
opts_xy at tmXBValues = fspan(0,xspan,nx)
opts_xy at tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx))
; Plotting options for vorticity
opts_vort = opts_xy
opts_vort at ContourParameters = (/ -20., 20., 0.5 /)
opts_vort at pmLabelBarOrthogonalPosF = -0.07
opts_vort at cnFillOn = True
opts_vort at cnLinesOn = False
; Plotting options for Pressure
opts_p = opts_xy
opts_p at ContourParameters = (/ 0.,1000.,50. /)
opts_p at cnInfoLabelOn = False
vres = True
vres at gsnMaximize = True ; Maximize plot in frame
;vres at pmLabelBarDisplayMode = "Always" ; Turn on a label bar.
;vres at pmLabelBarWidthF = 0.075
vres at lbPerimOn = False
vres at tiMainString = "No resources set"
;vres at NumVectors = 60
vres at Footer = False
vres at InitTime = False
vres at vcGlyphStyle = "CurlyVector";"LineArrow";"CurlyVector"
vres at vcMonoLineArrowColor = True ; color arrows based on magnitude
;vres at vcLevelSelectionMode = "ExplicitLevels"
;vres at vcLevels =
(/5,10,15,20,25,30,35,40,45,50,55,60,65/)
;vres at vcLevelColors = (//)
vres at vcRefLengthF = 0.1
vres at vcFillArrowWidthF = 50
vres at vcLineArrowThicknessF = 2
vres at vcMinFracLengthF = 0.5
vres at vcMinAnnoOn = True
vres at vcMinAnnoPerimOn = False
vres at vcMinAnnoString1On = False
;vres at vcMinAnnoString2 = "0.2"
vres at vcMinAnnoFontThicknessF = 10
vres at vcMinAnnoArrowLineColor = "Black"
;vres at vcMinAnnoArrowUseVecColor = "False "
;vres at MainTitle = "Vertically-integrated humidity flux"
;vres at MainTitlePos = "Center"
vres at vcMinAnnoArrowSpaceF = 4.0
opts_q = res
;opts_q at ContourParameters = (/ 10., 90., 10. /)
opts_q at pmLabelBarOrthogonalPosF = -0.1
opts_q at cnFillOn = True
opts_q at ContourParameters = (/ .2, 10., 0.5 /)
; Get the contour info for the rh and temp
if ( LimitedExtent ) then
contour_p = wrf_contour(a,wks,p_plane(0:zmax_pos,:),opts_p)
contour_vort = wrf_contour(a,wks,vort_plane(0:zmax_pos,:),opts_vort)
vector = wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:),vres)
contour_q = wrf_contour(a,wks,q_plane(0:zmax_pos,:),opts_q)
else
contour_p = wrf_contour(a,wks,p_plane,opts_p)
contour_vort = wrf_contour(a,wks,vort_plane,opts_vort)
vector = wrf_vector(a,wks,u_plane,w_plane,vres)
contour_q = wrf_contour(a,wks,q_plane,opts_q)
end if
; vector =
wrf_vector(a,wks,u_plane(0:zmax_pos,:),w_plane(0:zmax_pos,:),vres)
; MAKE PLOTS
plot = wrf_overlays(a,wks,(/contour_q,contour_p,vector/),pltres)
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160330/4f4c10ad/attachment.html
More information about the ncl-talk
mailing list