[ncl-talk] contour plot with 2D vertical coordinates
Dennis Shea
shea at ucar.edu
Thu Jan 5 09:24:16 MST 2017
Could you please ftp the file? ... or place it somewhere we can get access
to it?
ftp ftp.cgd.ucar.edu
anonymous
your_email
cd incoming
put test_noxios.nc
quit
Merci
On Thu, Jan 5, 2017 at 1:10 AM, Laurent ROBLOU <
laurent.roblou at aero.obs-mip.fr> wrote:
> Dear folks,
>
> Although I have carefully read around the corresponding pages on the topic
> (.../2Dvertcoord.shtml), I am getting into troubles for plotting a vertical
> contour section of the 3D variable w from ROMS model and overlay the model
> topography.
>
> 1) It looks that the topography doesn't match the last sigma layer and I
> would bet on an error overlaying the topography. (Figures
> with_topo_overlaid_mesh.ps and with_topo_overlaid_nomesh.ps with respect
> to without_topo_overlaid.ps)
>
>
> 2) Another error is that I "loose" my xticks labels when switching
> trGridType from triangularMesh (to turn the edges on as in 2dvertcoords_2)
> to curvilinear grid (to get a nicer plot) with that error message:
> warning:tmXBLabelFontHeightF is not a valid resource in
> w.ps_contour.PlotManager at this time
> warning:tmXBMode is not a valid resource in w.ps_contour.PlotManager at
> this time
> warning:tmXBValues is not a valid resource in w.ps_contour.PlotManager at
> this time
> warning:tmXBLabels is not a valid resource in w.ps_contour.PlotManager at
> this time
>
> (Figures with_topo_overlaid_mesh.ps and with_topo_overlaid_nomesh.ps
> again)
>
>
> 3) Last error of mine, there is no frame to my contour plot when no grid
> nor tickmarks when no topography overlaid. Do see figure
> without_topo_overlaid_nomesh.ps
>
>
> Hereafter you will find the full script. Many thanks in advance for any
> help
>
> Best regards,
> Laurent
>
>
>
>
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; ;
> ; Copyright (C) 2016 ;
> ; Centre National de la Recherche Scientifique CNRS ;
> ; ;
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;
> ; File: croco_vslice.ncl
> ; Author: Laurent ROBLOU
> ; Centre National de la Recherche Scientifique CNRS
> ; Laboratoire d Aerologie, Toulouse
>
> ; this program interpolates the vertical section of a 3D variable along
> ; the transect line approximated by a great circle section and overlay
> ; the underlying sea floor topography
> ; interpolation is made at each time step for the specified 3D var
>
>
> ; INIT LIBS AND PLOT
> ; ========================
> load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/contrib/
> calendar_decode2.ncl"
>
>
>
> begin
>
>
> ;************************************
> ; get model variables
> ;************************************
>
> f = addfile("test_noxios.nc","r")
>
> var = f->$"w"$
> if (var at units .eq. 1) then
> var at units="-"
> end if
>
> if(isatt(var,"long_name"))
> varname=var at long_name
> else
> if(isatt(var,"standard_name"))
> varname=var at standard_name
> else
> varname=var at short_name ; last tentative
> end if
> end if
>
> lat = f->lat_rho ; get lat
> lon = f->lon_rho ; get lon
> time = f->time(0:2) ; get time frame
>
> vlev=f->s_rho
> f = addfile ("test_noxios.nc", "r")
>
> print("process s coordinates")
> h=f->h;
> h=where(h.eq.0, 1.e-14, h)
> h2=h+f->hc;
> h2inv=1./h2;
> zeta=f->zeta
>
> ; grid depth at rho-points
> sc_r=f->sc_r
> Cs_r=f->Cs_r
> z_r = new( dimsizes(f->rho), float, default_fillvalue("float") )
> nlevels=dimsizes(f->s_rho)
> do k=0,nlevels-1
> z0=1.e+16*sc_r(k)+Cs_r(k)*h;
> dum1=(1 + z0 * h2inv)
> dum2=z0 * h * h2inv
> z_r(:,k,:,:)=tofloat( conform_dims(dimsizes(zeta), dum2,
> (/1,2/) ) + zeta * conform_dims(dimsizes(zeta), dum1, (/1,2/) ) )
> end do
> z_r at standard_name="grid_point_depth_at_RHO-points";
> z_r at long_name="grid point depth at RHO-points";
> z_r at units=zeta at units;
> z_r at axis="XYZT";
> end if
>
>
> topo=f->h * (-1) ; move to negative values
> interp_topo = rcm2points (lat, lon, topo, cross_section_gc at gclat,
> cross_section_gc at gclon, 2)
>
>
> ; PLOT
> ; ========================
>
>
> ; define cross-section
> ;********************************
>
> if ( "xx" .eq. "xTruex" ) ; Turn on edges on
> npts = 50
> else
> dum=dimsizes(lat)
> d1=dum(0)
> d2=dum(1)
> d=sqrt(d1*d1+d2*d2)
> npts = min ( (/ toint(d), 1000 /)) ; number of points in resulting
> cross-section
> end if
> cross_section_gc = gc_latlon(35.9,-6,35.89,-5.6,npts,-2)
> points = ispan(0,npts-1,1)*1.0
>
>
> ; create contour plot
> ;********************************
> print("create contour plot")
> wks= gsn_open_wks("x11", "figure")
> gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
>
> res = True ; Use plot options
>
> ; high-level control
> res at gsnDraw = False ; dont draw
> res at gsnFrame = False ; dont advance frame
> res at gsnMaximize = True ; maximize plot in frame
> res at gsnBoxMargin = 0.0
> res at gsnLeftString = " "
> res at gsnRightString = " "
> res at gsnScale = False
> res at trGridType = "curvilinear"; "TriangularMesh" ;
>
> ; xticks labels
> nLabels=10
> XBValues = toint( fspan(0,npts-1,nLabels) )
> XBLabels = new(nLabels,"string")
> do i=0,nLabels-1
> x = cross_section_gc at gclon(XBValues(i))
> y = cross_section_gc at gclat(XBValues(i))
> XBLabels(i) = sprintf("%5.1f", y)+"~S~o~N~N"+"~C~"+sprintf("%5.1f",
> x)+"~S~o~N~E"
> end do
>
> print(XBValues)
> print(XBLabels)
>
>
> res at tmXBMode = "Explicit" ; explicitly label x-axis
> res at tmXBValues = XBValues
> res at tmXBLabels = XBLabels
> res at tmXBLabelFontHeightF = 0.01
>
> ; contours
> res at cnFillOn = True
>
> if ( "xx" .eq. "xTruex" ) ; Turn on edges on
> res at cnFillMode = "CellFill"
> res at cnCellFillEdgeColor= "red" ; "grey22" ; "black
> res at cnLinesOn = False ; turn off countour lines
> res at trGridType = "TriangularMesh" ;
> end if
> res at cnFillDrawOrder = "PreDraw"
> res at cnLineDrawOrder = "PreDraw"
>
>
> res at cnInfoLabelOn = False ; turn off contour label
> res at cnMaxLevelCount = 21
> varmin=min(var)
> varmax=max(var)
> maxlev=50
> ;---Make sure contour levels are fixed for all time steps
> mini=3.40282347e+38
> maxi=-1.17549435e-38
> i= 0
> do while (i .le. 2 )
> varmin=min( var(i+0,:,:,:) )
> varmax=max( var(i+0,:,:,:) )
> mini=min((/tofloat(mini), tofloat(varmin)/)) ; minimum value of
> mini,varmin array
> maxi=min((/tofloat(maxi), tofloat(varmax)/)) ; maximum value
> i=i+1
> end do
> mnmxint = nice_mnmxintvl(varmin, varmax, maxlev, False)
>
> res at cnLevelSelectionMode = "ManualLevels"
> res at cnMinLevelValF = mnmxint(0) ; min( (*/mnmxint(0), mini/*) )
> res at cnMaxLevelValF = mnmxint(1)
> res at cnLevelSpacingF = mnmxint(2)
>
>
> ; label bar
> res at lbLabelAutoStride = True ; nice label bar label stride
> res at gsnSpreadColors = True ; use full range of colormap
>
> res at lbLabelBarOn = True
> res at lbTitlePosition = "Bottom"
> res at lbTitleDirection = "Across"
> res at lbOrientation = "Vertical"
> res at lbTopMarginF = 0.0
>
> ; titles
> res at tiYAxisString = ("z_r" + " ["+z_r at units+"]")
> res at gsnYAxisIrregular2Linear = True ; Y axis will be linearized
>
>
>
> ; another ressource for bathymetry
> tres = True ; plot mods desired
> tres at gsnDraw = False ; don t draw yet
> tres at gsnFrame = False ; don t advance frame yet
> tres at gsnMaximize = True ; maximize plot in frame
> tres at xyCurveDrawOrder = "PostDraw"
> tres at gsnYRefLine = min(interp_topo) ; create a
> reference line
> tres at gsnAboveYRefLineColor = "tan";"gray75" ; above ref line
> fill color
> tres at trYReverse = True ; reverse y axis
>
> i=0
> do while (i .le. 2 )
>
> ; interpolate data along a great circle
> ;************************************
> print("interpolate data along the cross-section")
> interp = new ((*/nlevels,npts/*), double,
> default_fillvalue("double"))
> z = new ((*/nlevels,npts/*), double, default_fillvalue("double"))
> do k = 0,nlevels-1
> interp(k,:) = rcm2points_Wrap (lat, lon, var(i+0,k,:,:),
> cross_section_gc at gclat, cross_section_gc at gclon, 2)
> z(k,:) = rcm2points_Wrap (lat, lon, z_r(i+0,k,:,:),
> cross_section_gc at gclat, cross_section_gc at gclon, 2)
> end do
>
> interp!0 = "levels"
> interp!1 = "points"
> interp&levels=vlev
>
> z!0 = "levels"
> z!1 = "points"
> z&levels=vlev
>
>
> copy_VarAtts(var,interp) ; copy attributes
> copy_VarAtts(z_r,z) ; copy attributes
>
> ; plot section
> ;************************************
>
> ; update title
> t=time(i)
> utc_date = cd_calendar(t, 0)
> year = tointeger(utc_date(:,0))
> month = tointeger(utc_date(:,1))
> day = tointeger(utc_date(:,2))
> hour = tointeger(utc_date(:,3))
> minute = tointeger(utc_date(:,4))
> sdate = sprinti("%0.2i/", day) + sprinti("%0.2i/",month) +
> sprinti("%0.4i", year)
> stime = sprinti("%0.2ih", hour) + sprinti("%0.2i", minute)
> titledate = ("time = " + sdate + " - " + stime)
> res at tiMainString = (varname +" ~C~" + titledate)
>
> ; print out stats as subtiltle
> statsline=("min="+sprintf("%5.2f", min(interp)) + " mean=" +
> sprintf("%5.2f", avg(interp)) + " max=" + sprintf("%5.2f", max(interp))
> +" ["+ interp at units +"]")
> res at gsnCenterString = statsline
>
> ; plot sliced variable
>
> res at sfXArray = tofloat(conform_dims(dimsizes(z), points,
> 1) ); data for x-axis tms
> res at sfYArray = z ; data for y-axis tms
>
> ; plot sliced variable
> profile = gsn_csm_contour(wks,interp,res)
>
>
> ; add terrain topography
> ;************************************
> topo_profile = gsn_csm_y(wks,interp_topo,tres)
> overlay(profile, topo_profile)
>
> maximize_output(wks,True)
> i=i+1
>
> ; clean up variables
> delete(profile)
> delete(interp)
> delete(z)
> end do
>
> 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/20170105/3853eacd/attachment.html
More information about the ncl-talk
mailing list