[ncl-talk] contour plot with 2D vertical coordinates

Mary Haley haley at ucar.edu
Thu Jan 12 11:21:47 MST 2017


Hi Laurent,

Thanks for providing the data file.

This email was internally sent to several people, so unfortunately
everybody assumed somebody else was working on it.

I decided to take a stab at it, but I was unable to run your script because
1) there's a stray "end if" statement, and 2) cross_section_gc is undefined
in the script you sent.

Instead of copying-and-pasting a script in an email, can you send a clean
script that works as an attachment?

Thanks,

--Mary


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/20170112/43d3efc5/attachment.html 


More information about the ncl-talk mailing list