[ncl-talk] plot longitude in a cross section
Mary Haley
haley at ucar.edu
Sat Apr 2 15:01:16 MDT 2016
Jesús,
It's a bit hard to debug your script because you are setting *lots* of
resources and it's hard to follow which plots they are applying to.
The error says the third argument to "fspan" is equal to 0 or 1, which is
not allowed. I believe this is occurring on this line:
opts_xy at tmXBValues = fspan(0,xspan,nx)
so you need to look at how "nx" is being calculated:
nx = floattoint((xmax-xmin)/2+1)
There's probably an issue with xmin and/or xmax. Maybe they are both equal?
--Mary
On Wed, Mar 30, 2016 at 7:37 AM, Jesús Garcia Rosales <jesus21gr at gmail.com>
wrote:
> 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
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160402/94c3cbb9/attachment.html
More information about the ncl-talk
mailing list