<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
</head>
<body text="#000000" bgcolor="#FFFFFF">
Dear folks,
<br>
<br>
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.
<br>
<br>
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)
<br>
<br>
<br>
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:
<br>
warning:tmXBLabelFontHeightF is not a valid resource in
w.ps_contour.PlotManager at this time
<br>
warning:tmXBMode is not a valid resource in w.ps_contour.PlotManager
at this time
<br>
warning:tmXBValues is not a valid resource in
w.ps_contour.PlotManager at this time
<br>
warning:tmXBLabels is not a valid resource in
w.ps_contour.PlotManager at this time
<br>
<br>
(Figures with_topo_overlaid_mesh.ps and with_topo_overlaid_nomesh.ps
again)
<br>
<br>
<br>
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
<br>
<br>
<br>
Hereafter you will find the full script. Many thanks in advance for
any help
<br>
<br>
Best regards,
<br>
Laurent
<br>
<br>
<br>
<br>
<br>
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
<br>
; ;
<br>
; Copyright (C)
2016 ;
<br>
; Centre National de la Recherche Scientifique
CNRS ;
<br>
; ;
<br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
<br>
;
<br>
; File: croco_vslice.ncl
<br>
; Author: Laurent ROBLOU
<br>
; Centre National de la Recherche Scientifique CNRS
<br>
; Laboratoire d Aerologie, Toulouse
<br>
<br>
; this program interpolates the vertical section of a 3D variable
along <br>
; the transect line approximated by a great circle section and
overlay<br>
; the underlying sea floor topography<br>
; interpolation is made at each time step for the specified 3D var<br>
<br>
<br>
; INIT LIBS AND PLOT
<br>
; ========================
<br>
load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/gsn_code.ncl"
<br>
load "/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
<br>
load
"/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/csm/contributed.ncl"
<br>
load
"/opt/ncl/ncl-6.1.2-LA/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl"<br>
<br>
<br>
<br>
begin
<br>
<br>
<br>
;************************************
<br>
; get model variables
<br>
;************************************
<br>
<br>
f = addfile("test_noxios.nc","r")
<br>
<br>
var = f->$"w"$
<br>
if (var@units .eq. 1) then
<br>
var@units="-"
<br>
end if
<br>
<br>
if(isatt(var,"long_name"))
<br>
varname=var@long_name
<br>
else
<br>
if(isatt(var,"standard_name"))
<br>
varname=var@standard_name
<br>
else
<br>
varname=var@short_name ; last tentative
<br>
end if
<br>
end if
<br>
<br>
lat = f->lat_rho ; get lat
<br>
lon = f->lon_rho ; get lon
<br>
time = f->time(0:2) ; get time frame
<br>
<br>
vlev=f->s_rho
<br>
f = addfile ("test_noxios.nc", "r")
<br>
<br>
print("process s coordinates")
<br>
h=f->h;
<br>
h=where(h.eq.0, 1.e-14, h)
<br>
h2=h+f->hc;
<br>
h2inv=1./h2;
<br>
zeta=f->zeta
<br>
<br>
; grid depth at rho-points
<br>
sc_r=f->sc_r
<br>
Cs_r=f->Cs_r
<br>
z_r = new( dimsizes(f->rho), float, default_fillvalue("float") )
<br>
nlevels=dimsizes(f->s_rho)
<br>
do k=0,nlevels-1
<br>
z0=1.e+16*sc_r(k)+Cs_r(k)*h;
<br>
dum1=(1 + z0 * h2inv)
<br>
dum2=z0 * h * h2inv
<br>
z_r(:,k,:,:)=tofloat( conform_dims(dimsizes(zeta),
dum2, (/1,2/) ) + zeta * conform_dims(dimsizes(zeta), dum1,
(/1,2/) ) )
<br>
end do
<br>
z_r@standard_name="grid_point_depth_at_RHO-points";
<br>
z_r@long_name="grid point depth at RHO-points";
<br>
z_r@units=zeta@units;
<br>
z_r@axis="XYZT";
<br>
end if
<br>
<br>
<br>
topo=f->h * (-1) ; move to negative values
<br>
interp_topo = rcm2points (lat, lon, topo, cross_section_gc@gclat,
cross_section_gc@gclon, 2)
<br>
<br>
<br>
; PLOT
<br>
; ========================
<br>
<br>
<br>
; define cross-section
<br>
;********************************
<br>
<br>
if ( "xx" .eq. "xTruex" ) ; Turn on edges on
<br>
npts = 50
<br>
else
<br>
dum=dimsizes(lat)
<br>
d1=dum(0)
<br>
d2=dum(1)
<br>
d=sqrt(d1*d1+d2*d2)
<br>
npts = min ( (/ toint(d), 1000 /)) ; number of points in
resulting cross-section
<br>
end if
<br>
cross_section_gc = gc_latlon(35.9,-6,35.89,-5.6,npts,-2)
<br>
points = ispan(0,npts-1,1)*1.0
<br>
<br>
<br>
; create contour plot
<br>
;********************************
<br>
print("create contour plot")
<br>
wks= gsn_open_wks("x11", "figure")
<br>
gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
<br>
<br>
res = True ; Use plot options
<br>
<br>
; high-level control
<br>
res@gsnDraw = False ; dont draw
<br>
res@gsnFrame = False ; dont advance frame
<br>
res@gsnMaximize = True ; maximize plot in frame
<br>
res@gsnBoxMargin = 0.0
<br>
res@gsnLeftString = " "
<br>
res@gsnRightString = " "
<br>
res@gsnScale = False
<br>
res@trGridType = "curvilinear"; "TriangularMesh" ;
<br>
<br>
; xticks labels
<br>
nLabels=10
<br>
XBValues = toint( fspan(0,npts-1,nLabels) )
<br>
XBLabels = new(nLabels,"string")
<br>
do i=0,nLabels-1
<br>
x = cross_section_gc@gclon(XBValues(i))
<br>
y = cross_section_gc@gclat(XBValues(i))
<br>
XBLabels(i) = sprintf("%5.1f",
y)+"~S~o~N~N"+"~C~"+sprintf("%5.1f", x)+"~S~o~N~E"
<br>
end do
<br>
<br>
print(XBValues)
<br>
print(XBLabels)
<br>
<br>
<br>
res@tmXBMode = "Explicit" ; explicitly label x-axis
<br>
res@tmXBValues = XBValues
<br>
res@tmXBLabels = XBLabels
<br>
res@tmXBLabelFontHeightF = 0.01
<br>
<br>
; contours
<br>
res@cnFillOn = True
<br>
<br>
if ( "xx" .eq. "xTruex" ) ; Turn on edges on
<br>
res@cnFillMode = "CellFill"
<br>
res@cnCellFillEdgeColor= "red" ; "grey22" ; "black
<br>
res@cnLinesOn = False ; turn off countour lines
<br>
res@trGridType = "TriangularMesh" ;
<br>
end if
<br>
res@cnFillDrawOrder = "PreDraw"
<br>
res@cnLineDrawOrder = "PreDraw"
<br>
<br>
<br>
res@cnInfoLabelOn = False ; turn off contour
label
<br>
res@cnMaxLevelCount = 21
<br>
varmin=min(var)
<br>
varmax=max(var)
<br>
maxlev=50
<br>
;---Make sure contour levels are fixed for all time steps
<br>
mini=3.40282347e+38
<br>
maxi=-1.17549435e-38
<br>
i= 0
<br>
do while (i .le. 2 )
<br>
varmin=min( var(i+0,:,:,:) )
<br>
varmax=max( var(i+0,:,:,:) )
<br>
mini=min((/tofloat(mini), tofloat(varmin)/)) ; minimum value of
mini,varmin array
<br>
maxi=min((/tofloat(maxi), tofloat(varmax)/)) ; maximum value
<br>
i=i+1
<br>
end do
<br>
mnmxint = nice_mnmxintvl(varmin, varmax, maxlev, False)
<br>
<br>
res@cnLevelSelectionMode = "ManualLevels"
<br>
res@cnMinLevelValF = mnmxint(0) ; min( (<i
class="moz-txt-slash"><span class="moz-txt-tag">/</span>mnmxint(0),
mini<span class="moz-txt-tag">/</span></i>) )
<br>
res@cnMaxLevelValF = mnmxint(1)
<br>
res@cnLevelSpacingF = mnmxint(2)
<br>
<br>
<br>
; label bar
<br>
res@lbLabelAutoStride = True ; nice label bar label stride
<br>
res@gsnSpreadColors = True ; use full range of colormap
<br>
<br>
res@lbLabelBarOn = True
<br>
res@lbTitlePosition = "Bottom"
<br>
res@lbTitleDirection = "Across"
<br>
res@lbOrientation = "Vertical"
<br>
res@lbTopMarginF = 0.0
<br>
<br>
; titles
<br>
res@tiYAxisString = ("z_r" + " ["+z_r@units+"]")
<br>
res@gsnYAxisIrregular2Linear = True ; Y axis will be linearized
<br>
<br>
<br>
<br>
; another ressource for bathymetry
<br>
tres = True ; plot mods desired
<br>
tres@gsnDraw = False ; don t draw yet
<br>
tres@gsnFrame = False ; don t advance frame yet
<br>
tres@gsnMaximize = True ; maximize plot in frame
<br>
tres@xyCurveDrawOrder = "PostDraw"
<br>
tres@gsnYRefLine = min(interp_topo) ; create a
reference line
<br>
tres@gsnAboveYRefLineColor = "tan";"gray75" ; above ref
line fill color
<br>
tres@trYReverse = True ; reverse y axis
<br>
<br>
i=0
<br>
do while (i .le. 2 )
<br>
<br>
; interpolate data along a great circle
<br>
;************************************
<br>
print("interpolate data along the cross-section")
<br>
interp = new ((<i class="moz-txt-slash"><span
class="moz-txt-tag">/</span>nlevels,npts<span
class="moz-txt-tag">/</span></i>), double,
default_fillvalue("double"))
<br>
z = new ((<i class="moz-txt-slash"><span class="moz-txt-tag">/</span>nlevels,npts<span
class="moz-txt-tag">/</span></i>), double,
default_fillvalue("double"))
<br>
do k = 0,nlevels-1
<br>
interp(k,:) = rcm2points_Wrap (lat, lon, var(i+0,k,:,:),
cross_section_gc@gclat, cross_section_gc@gclon, 2)
<br>
z(k,:) = rcm2points_Wrap (lat, lon, z_r(i+0,k,:,:),
cross_section_gc@gclat, cross_section_gc@gclon, 2)
<br>
end do
<br>
<br>
interp!0 = "levels"
<br>
interp!1 = "points"
<br>
interp&levels=vlev
<br>
<br>
z!0 = "levels"
<br>
z!1 = "points"
<br>
z&levels=vlev
<br>
<br>
<br>
copy_VarAtts(var,interp) ; copy attributes
<br>
copy_VarAtts(z_r,z) ; copy attributes
<br>
<br>
; plot section
<br>
;************************************
<br>
<br>
; update title
<br>
t=time(i)
<br>
utc_date = cd_calendar(t, 0)
<br>
year = tointeger(utc_date(:,0))
<br>
month = tointeger(utc_date(:,1))
<br>
day = tointeger(utc_date(:,2))
<br>
hour = tointeger(utc_date(:,3))
<br>
minute = tointeger(utc_date(:,4))
<br>
sdate = sprinti("%0.2i/", day) + sprinti("%0.2i/",month) +
sprinti("%0.4i", year)
<br>
stime = sprinti("%0.2ih", hour) + sprinti("%0.2i", minute)
<br>
titledate = ("time = " + sdate + " - " + stime)
<br>
res@tiMainString = (varname +" ~C~" + titledate)
<br>
<br>
; print out stats as subtiltle
<br>
statsline=("min="+sprintf("%5.2f", min(interp)) + " mean="
+ sprintf("%5.2f", avg(interp)) + " max=" + sprintf("%5.2f",
max(interp)) +" ["+ interp@units +"]")
<br>
res@gsnCenterString = statsline
<br>
<br>
; plot sliced variable
<br>
<br>
res@sfXArray = tofloat(conform_dims(dimsizes(z),
points, 1) ); data for x-axis tms
<br>
res@sfYArray = z ; data for y-axis tms
<br>
<br>
; plot sliced variable
<br>
profile = gsn_csm_contour(wks,interp,res)
<br>
<br>
<br>
; add terrain topography
<br>
;************************************
<br>
topo_profile = gsn_csm_y(wks,interp_topo,tres)
<br>
overlay(profile, topo_profile)
<br>
<br>
maximize_output(wks,True)
<br>
i=i+1
<br>
<br>
; clean up variables
<br>
delete(profile)
<br>
delete(interp)
<br>
delete(z)
<br>
end do
<br>
<br>
end
</body>
</html>