[ncl-talk] contour plot with 2D vertical coordinates
Laurent ROBLOU
laurent.roblou at aero.obs-mip.fr
Thu Jan 5 01:10:25 MST 2017
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170105/fd03e3ea/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: with_topo_overlaid_mesh.ps
Type: application/postscript
Size: 660911 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170105/fd03e3ea/attachment-0004.ps
-------------- next part --------------
A non-text attachment was scrubbed...
Name: with_topo_overlaid_nomesh.ps
Type: application/postscript
Size: 547191 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170105/fd03e3ea/attachment-0005.ps
-------------- next part --------------
A non-text attachment was scrubbed...
Name: without_topo_overlaid_mesh.ps
Type: application/postscript
Size: 659004 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170105/fd03e3ea/attachment-0006.ps
-------------- next part --------------
A non-text attachment was scrubbed...
Name: without_topo_overlaid_nomesh.ps
Type: application/postscript
Size: 490215 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170105/fd03e3ea/attachment-0007.ps
More information about the ncl-talk
mailing list