[ncl-talk] plot longitude in a cross section
Jesús Garcia Rosales
jesus21gr at gmail.com
Tue Apr 5 12:53:38 MDT 2016
Thanks for your answer Mary.
Greetings,
Alan
2016-04-02 16:01 GMT-05:00 Mary Haley <haley at ucar.edu>:
> 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/20160405/67df9c8e/attachment.html
More information about the ncl-talk
mailing list