[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