[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