<div dir="ltr"><div class="gmail_default" style="font-size:small">Hi Laurent,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thanks for providing the data file.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">This email was internally sent to several people, so unfortunately everybody assumed somebody else was working on it.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I decided to take a stab at it, but I was unable to run your script because 1) there&#39;s a stray &quot;end if&quot; statement, and 2) cross_section_gc is undefined in the script you sent.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Instead of copying-and-pasting a script in an email, can you send a clean script that works as an attachment?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Thanks,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">--Mary</div><div class="gmail_default" style="font-size:small"><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Jan 5, 2017 at 1:10 AM, Laurent ROBLOU <span dir="ltr">&lt;<a href="mailto:laurent.roblou@aero.obs-mip.fr" target="_blank">laurent.roblou@aero.obs-mip.fr</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
  

    
  
  <div 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&#39;t match the last sigma layer
    and I would bet on an error overlaying the topography. (Figures
    <a href="http://with_topo_overlaid_mesh.ps" target="_blank">with_topo_overlaid_mesh.ps</a> and <a href="http://with_topo_overlaid_nomesh.ps" target="_blank">with_topo_overlaid_nomesh.ps</a> with
    respect to <a href="http://without_topo_overlaid.ps" target="_blank">without_topo_overlaid.ps</a>)
    <br>
    <br>
    <br>
    2) Another error is that I &quot;loose&quot; 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 <a href="http://with_topo_overlaid_mesh.ps" target="_blank">with_topo_overlaid_mesh.ps</a> and <a href="http://with_topo_overlaid_nomesh.ps" target="_blank">with_topo_overlaid_nomesh.ps</a>
    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
    <a href="http://without_topo_overlaid_nomesh.ps" target="_blank">without_topo_overlaid_nomesh.<wbr>ps</a>
    <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>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<wbr>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<wbr>;;;;;;;;;;;;
    <br>
    ; ;
    <br>
    ;                       Copyright (C) 
    2016                          <wbr>  ;
    <br>
    ;        Centre National de la Recherche Scientifique
    CNRS             ;
    <br>
    ; ;
    <br>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<wbr>;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;<wbr>;;;;;;;;;;;;
    <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 &quot;/opt/ncl/ncl-6.1.2-LA/lib/<wbr>ncarg/nclscripts/csm/gsn_code.<wbr>ncl&quot;
    <br>
    load &quot;/opt/ncl/ncl-6.1.2-LA/lib/<wbr>ncarg/nclscripts/csm/gsn_csm.<wbr>ncl&quot;
    <br>
    load
    &quot;/opt/ncl/ncl-6.1.2-LA/lib/<wbr>ncarg/nclscripts/csm/<wbr>contributed.ncl&quot;
    <br>
    load
&quot;/opt/ncl/ncl-6.1.2-LA/lib/<wbr>ncarg/nclscripts/contrib/<wbr>calendar_decode2.ncl&quot;<br>
    <br>
    <br>
    <br>
    begin
    <br>
    <br>
    <br>
    ;*****************************<wbr>*******
    <br>
    ; get model variables
    <br>
    ;*****************************<wbr>*******
    <br>
    <br>
    f = addfile(&quot;<a href="http://test_noxios.nc" target="_blank">test_noxios.nc</a>&quot;,&quot;r&quot;)
    <br>
    <br>
    var  = f-&gt;$&quot;w&quot;$
    <br>
    if (var@units .eq. 1) then
    <br>
        var@units=&quot;-&quot;
    <br>
    end if
    <br>
    <br>
    if(isatt(var,&quot;long_name&quot;))
    <br>
        varname=var@long_name
    <br>
    else
    <br>
        if(isatt(var,&quot;standard_name&quot;))
    <br>
        varname=var@standard_name
    <br>
        else
    <br>
        varname=var@short_name ; last tentative
    <br>
        end if
    <br>
    end if
    <br>
    <br>
    lat = f-&gt;lat_rho  ; get lat
    <br>
    lon = f-&gt;lon_rho  ; get lon
    <br>
    time = f-&gt;time(0:2)  ; get time frame
    <br>
    <br>
    vlev=f-&gt;s_rho
    <br>
    f = addfile (&quot;<a href="http://test_noxios.nc" target="_blank">test_noxios.nc</a>&quot;, &quot;r&quot;)
    <br>
    <br>
     print(&quot;process s coordinates&quot;)
    <br>
     h=f-&gt;h;
    <br>
     h=where(h.eq.0, 1.e-14, h)
    <br>
     h2=h+f-&gt;hc;
    <br>
     h2inv=1./h2;
    <br>
     zeta=f-&gt;zeta
    <br>
    <br>
     ; grid depth at rho-points
    <br>
     sc_r=f-&gt;sc_r
    <br>
     Cs_r=f-&gt;Cs_r
    <br>
     z_r = new( dimsizes(f-&gt;rho), float, default_fillvalue(&quot;float&quot;) )
    <br>
     nlevels=dimsizes(f-&gt;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=&quot;grid_<wbr>point_depth_at_RHO-points&quot;;
    <br>
     z_r@long_name=&quot;grid point depth at RHO-points&quot;;
    <br>
     z_r@units=zeta@units;
    <br>
     z_r@axis=&quot;XYZT&quot;;
    <br>
    end if
    <br>
    <br>
    <br>
    topo=f-&gt;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>
    ;*****************************<wbr>***
    <br>
    <br>
    if ( &quot;xx&quot; .eq. &quot;xTruex&quot; ) ; 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,<wbr>npts,-2)
    <br>
    points   = ispan(0,npts-1,1)*1.0
    <br>
    <br>
    <br>
    ; create contour plot
    <br>
    ;*****************************<wbr>***
    <br>
    print(&quot;create contour plot&quot;)
    <br>
    wks= gsn_open_wks(&quot;x11&quot;, &quot;figure&quot;)
    <br>
    gsn_define_colormap(wks,&quot;<wbr>WhiteBlueGreenYellowRed&quot;)
    <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 = &quot; &quot;
    <br>
    res@gsnRightString = &quot; &quot;
    <br>
    res@gsnScale        = False
    <br>
    res@trGridType            = &quot;curvilinear&quot;; &quot;TriangularMesh&quot; ;
    <br>
    <br>
    ; xticks labels
    <br>
    nLabels=10
    <br>
    XBValues    = toint( fspan(0,npts-1,nLabels) )
    <br>
    XBLabels    = new(nLabels,&quot;string&quot;)
    <br>
    do i=0,nLabels-1
    <br>
        x = cross_section_gc@gclon(<wbr>XBValues(i))
    <br>
        y = cross_section_gc@gclat(<wbr>XBValues(i))
    <br>
        XBLabels(i) = sprintf(&quot;%5.1f&quot;,
    y)+&quot;~S~o~N~N&quot;+&quot;~C~&quot;+sprintf(&quot;%<wbr>5.1f&quot;, x)+&quot;~S~o~N~E&quot;
    <br>
    end do
    <br>
    <br>
    print(XBValues)
    <br>
    print(XBLabels)
    <br>
    <br>
    <br>
    res@tmXBMode            = &quot;Explicit&quot;    ; 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 ( &quot;xx&quot; .eq. &quot;xTruex&quot; ) ; Turn on edges on
    <br>
       res@cnFillMode            = &quot;CellFill&quot;
    <br>
       res@cnCellFillEdgeColor= &quot;red&quot; ; &quot;grey22&quot;     ; &quot;black
    <br>
       res@cnLinesOn           = False        ; turn off countour lines
    <br>
       res@trGridType            = &quot;TriangularMesh&quot; ;
    <br>
    end if
    <br>
    res@cnFillDrawOrder = &quot;PreDraw&quot;
    <br>
    res@cnLineDrawOrder = &quot;PreDraw&quot;
    <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 = &quot;ManualLevels&quot;
    <br>
    res@cnMinLevelValF       = mnmxint(0) ; min( (<i class="m_8118306918613670801moz-txt-slash"><span class="m_8118306918613670801moz-txt-tag">/</span>mnmxint(0),
      mini<span class="m_8118306918613670801moz-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 = &quot;Bottom&quot;
    <br>
    res@lbTitleDirection = &quot;Across&quot;
    <br>
    res@lbOrientation = &quot;Vertical&quot;
    <br>
    res@lbTopMarginF = 0.0
    <br>
    <br>
    ; titles
    <br>
    res@tiYAxisString       = (&quot;z_r&quot; + &quot; [&quot;+z_r@units+&quot;]&quot;)
    <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 = &quot;PostDraw&quot;
    <br>
    tres@gsnYRefLine           = min(interp_topo)             ; create a
    reference line
    <br>
    tres@gsnAboveYRefLineColor = &quot;tan&quot;;&quot;gray75&quot;              ; 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>
        ;*****************************<wbr>*******
    <br>
        print(&quot;interpolate data along the cross-section&quot;)
    <br>
        interp = new ((<i class="m_8118306918613670801moz-txt-slash"><span class="m_8118306918613670801moz-txt-tag">/</span>nlevels,npts<span class="m_8118306918613670801moz-txt-tag">/</span></i>), double,
    default_fillvalue(&quot;double&quot;))
    <br>
        z = new ((<i class="m_8118306918613670801moz-txt-slash"><span class="m_8118306918613670801moz-txt-tag">/</span>nlevels,npts<span class="m_8118306918613670801moz-txt-tag">/</span></i>), double,
    default_fillvalue(&quot;double&quot;))
    <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 = &quot;levels&quot;
    <br>
        interp!1 = &quot;points&quot;
    <br>
        interp&amp;levels=vlev
    <br>
    <br>
        z!0 = &quot;levels&quot;
    <br>
        z!1 = &quot;points&quot;
    <br>
        z&amp;levels=vlev
    <br>
    <br>
    <br>
        copy_VarAtts(var,interp)      <wbr>    ; copy attributes
    <br>
        copy_VarAtts(z_r,z)          ; copy attributes
    <br>
    <br>
        ; plot section
    <br>
        ;*****************************<wbr>*******
    <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(&quot;%0.2i/&quot;, day) + sprinti(&quot;%0.2i/&quot;,month)  +
    sprinti(&quot;%0.4i&quot;, year)
    <br>
       stime = sprinti(&quot;%0.2ih&quot;, hour)  + sprinti(&quot;%0.2i&quot;, minute)
    <br>
       titledate = (&quot;time = &quot; + sdate + &quot; - &quot; + stime)
    <br>
           res@tiMainString      = (varname +&quot; ~C~&quot; + titledate)
    <br>
    <br>
            ; print out stats as subtiltle
    <br>
           statsline=(&quot;min=&quot;+sprintf(&quot;%5.<wbr>2f&quot;, min(interp)) + &quot;    mean=&quot;
    + sprintf(&quot;%5.2f&quot;, avg(interp)) + &quot;    max=&quot; + sprintf(&quot;%5.2f&quot;,
    max(interp)) +&quot;    [&quot;+ interp@units +&quot;]&quot;)
    <br>
           res@gsnCenterString = statsline
    <br>
    <br>
            ; plot sliced variable
    <br>
    <br>
        res@sfXArray              = tofloat(conform_dims(dimsizes(<wbr>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,<wbr>res)
    <br>
    <br>
    <br>
    ; add terrain topography
    <br>
    ;*****************************<wbr>*******
    <br>
    topo_profile = gsn_csm_y(wks,interp_topo,<wbr>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
  </div>

<br>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/<wbr>mailman/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>