Re: multiple overlays

From: Dave Allured <dave.allured_at_nyahnyahspammersnyahnyah>
Date: Wed, 11 Jul 2007 13:12:27 -0600

Joe,

We had great results recently with almost the identical plot that
you requested.

We took advantage of the high level gsn routine
gsn_csm_vector_scalar_map_ce, which draws shaded contours and curly
vectors together in one function. This simplifies the code. The
only difference from yours is that the top overlay is the line
contours, rather than the curly vectors.

We thought this was a reasonable compromise, since those are both
line plots. However, you could instead draw the vectors as a
separate overlay and get the exact layering that you want.

The very simplified code sequence is:

     res_at_gsnDraw = False
     res_at_gsnFrame = False
     res2_at_gsnDraw = False
     res2_at_gsnFrame = False

     plot1 = gsn_csm_vector_scalar_map_ce (wks, u, v, grid, res)
     plot2 = gsn_csm_contour (wks, sfunc, res2)
     overlay (plot1, plot2)
     draw (plot1)
     frame (wks)

Attached is the full working script. This one also adds a fourth
"overlay" on top, in the form of a circle symbol using wmlabs.
Apologies for the clutter, this version has been through a lot.

--Dave A.
CU/CIRES Climate Diagnostics Center (CDC)
NOAA/ESRL/PSD, Climate Analysis Branch (CAB)

Joe Metzger wrote:
> I've searched the website and archives for this capability and may have
> missed it, but I'm wondering if NCL can plot *three* overlays? I want a
> color-filled plot on the bottom, then line contours on top of that, and
> finally curly vectors on the very top. Is this possible? Thanks for any
> advice or sample scripts.
>
> Joe
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

;----------------------------------------------------------------------------
;
; plotcomp.ncl -- Make offset plot series of composite fields.
;
; 2006-jun-22 bl, convert to south america. Adapted from sstcor.ncl.
; 2006-jun-29 dra, add both binary and Netcdf readers for mulitple offsets.
; 2007-jan-23 bl, add vector capability (program copied from plotcomp.ncl)
; 2007-jun-05 dra, add stream function from composite U and V winds.
; 2007-jun-25 dra, add contour line overlay for stream function, on
; same plot as shaded contour data field and wind vectors.
; Add command output selectors "ps=1" and "pdf=1", default = X11.
; Fix overlay dots.
;
;----------------------------------------------------------------------------

print ("Load ncl scripts")
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
print ("Load complete")

begin

; Parameters.

   olrdata = False ; select data field, set only one to True
   raindata = True

   if(olrdata)then
   infile = "/mac_home/bl/transient/composite/composite.olr.nc"
   varname = "compolr"
   end if

   if(raindata)then
   infile = "/mac_home/bl/transient/composite/composite.rain.lanczos.nc"
   varname = "comprain"
   end if

  idir = "/mac_home/bl/transient/composite/"
  uwindfile = idir + "composite.uwnd.lanczos.nc"
  varnameu = "compuwnd"

  vwindfile = idir + "composite.vwnd.lanczos.nc"
  varnamev = "compvwnd"

  sf_scale = "1e6" ; scale factor for stream function
;; sf_units = fs_scale + " x unity" ; units label for stream function
  
;; output = "x" ; x=window, ps=postscript, pdf=PDF
;; ; not used any more; use "ps=1" or "pdf=1"
; ; in command line; default = X11 window

  outfile = "kindexrain850-2.5N95W-3day10w"
; outfile = "kindexstream200-largest"
                                  ; name of output file; leave off .nc extension
  
  title = "Composite from index"

; Plot borders. Independent from data grid borders.

  lonw = -120 ;
  lone = -5
  lats = -50
  latn = 30
  
; lat_values = (/ -40, -30, -20, -10, 0, 10 /) ; lat labels with "EQ"
; lat_labels = (/ "40S", "30S", "20S", "10S", " EQ", "10N" /)

; contour_levels = (/ -7,-6,-5,-4,-3,-2, 2, 3, 4, 5, 6, 7 /)

; colormap = "WhViBlGrYeOrRe"

   colormap = "testcmap"
   black = 1 ; colors are specific to each color map
   blue = 35
   green = 80
   red = 200

;----------------------------------------------------------------------------

; Get command line parameters.

  if (isdefined ("ps")) then ; get output selector
    output = "ps" ; default = X11 window
  else
    if (isdefined ("pdf")) then
      output = "pdf"
    else
      output = "x11"
    end if
  end if

;----------------------------------------------------------------------------

; Derived parameters.

  chars = stringtochar (infile) ; determine file type
  len = dimsizes (chars) - 1 ; shortcut: assume 3 char ".xx" exten.
  e1 = max ( (/ 0, len-3 /) )
  e2 = max ( (/ e1, len-1 /) )
  suffix = chartostring (chars(e1:e2))

  netcdf = (suffix .eq. ".nc")
  binary = .not. netcdf

  if (binary) then
    nt = offset2 - offset_base + 1
  end if

;; Replaced by command line output selector.
;; if (output .ne. "ps" .and. output .ne. "pdf") then
;; output = "x11"
;; end if

  sf_scale_F = stringtofloat (sf_scale)

;----------------------------------------------------------------------------

; Read composite main input file.

; if (netcdf) then
    print ("Read main input file.")
    ifile = addfile (infile, "r")
    in = ifile->$varname$
; end if
  
    if (any (isnan_ieee (in))) then ; check input grids for NaN's
      print ("*** Abort, input contains NaN's.")
      exit
    end if
    
    if (.not. isdefined ("_FillValue")) then
      in@_FillValue = in_at_missing_value
    end if

; printVarSummary (in)
  
  print ("Min data = " + sprintf ("%0.3f", min (in)))
  print ("Max data = " + sprintf ("%0.3f", max (in)))
  print ("")

;----------------------------------------------------------------------------

  print ("Read composite wind files.")

  wfu = addfile (uwindfile, "r")
  inu = wfu->$varnameu$

  wfv = addfile (vwindfile, "r")
  inv = wfv->$varnamev$

;----------------------------------------------------------------------
; Common setup for all plots.
;----------------------------------------------------------------------

; print ("Start output file or window.")

   wks = gsn_open_wks (output, outfile) ; open output "device"

   gsn_define_colormap (wks, colormap)

  res = True ; plot mods desired
  res_at_gsnDraw = False ; hold draw and auto frame for
  res_at_gsnFrame = False ; overlays (contours and dots)

  res_at_cnFillOn = True ; turn on color fill
  res_at_cnLinesOn = True ; turn on contour lines
  
  res_at_gsnSpreadColors = True ; choose colors across entire range
; res_at_gsnSpreadColorStart = -1 ; reverse colors
; res_at_gsnSpreadColorEnd = 2
; res_at_gsnSpreadColorEnd = 10
; res_at_gsnSpreadColorEnd = 30

   res_at_lbOrientation = "horizontal" ; horizontal label bars
; res_at_lbOrientation = "vertical" ; vertical label bars

  res_at_pmLabelBarOrthogonalPosF = 0.3 ; label bar farther from plot
  res_at_tiMainFontHeightF = 0.018 ; adjust main title font size

  if(olrdata)then ; set contours for OLR data

    res_at_cnLevelSelectionMode = "ExplicitLevels"
    contour_levels = (/ -30,-25.,-20.,-15.,-10.,-5.,5.,10.,15.,20.,25.,30. /)
    res_at_cnLevels = contour_levels

; res_at_cnLevelSelectionMode = "ManualLevels" ; manual levels
; res_at_cnMinLevelValF = -10 ; min level
; res_at_cnMaxLevelValF = 10 ; max level
; res_at_cnLevelSpacingF = 2 ; interval

; res_at_cnMinLevelValF = -7 ; min level ; good for rain
; res_at_cnMaxLevelValF = 7 ; max level
; res_at_cnLevelSpacingF = 2 ; interval

; res_at_cnMinLevelValF = -30 ; min level
; res_at_cnMaxLevelValF = 30 ; max level
; res_at_cnLevelSpacingF = 5 ; interval

  end if
  
  if(raindata)then ; set contours for rain data

    res_at_cnLevelSelectionMode = "ExplicitLevels"
    contour_levels = (/ -10.,-8.,-6.,-4.,-2.,2.,4.,6.,8.,10. /)
    res_at_cnLevels = contour_levels

; res_at_cnLevelSelectionMode = "ManualLevels" ; manual levels
; res_at_cnMinLevelValF = -8 ; min level
; res_at_cnMaxLevelValF = 8 ; max level
; res_at_cnLevelSpacingF = 1 ; interval

  end if

; Display color bar without end labels.

; res_at_cnLevels = contour_levels

; Display color bar with forced end labels.
; Warning: Data beyond end values will be plotted to appear as if inside
; the end values. No warning is printed unless code is added.

; This version: Forced bottom end, open top end.
  
; nclevs = dimsizes (contour_levels)
;; res_at_cnLevels = contour_levels(1:nclevs-2) ; hide end vals from NCL
; res_at_cnLevels = contour_levels(1:) ; hide only bottom value
; res_at_cnExplicitLabelBarLabelsOn = True
; res_at_cnLabelBarEndLabelsOn = True
; label_strings = new (nclevs+1, string, "")
; label_strings(0:nclevs-1) = sprintf ("%1.0f", contour_levels)
; res_at_lbLabelStrings = label_strings

; end of forced bottom and top labels
 
   res_at_cnMissingValPerimOn = True
  res_at_cnMissingValPerimDashPattern = 1

   res_at_mpFillOn = False ; turn off gray continents
  
  res_at_gsnMaximize = True ; makes figure fit in frame properly

  res_at_gsnPaperOrientation = "portrait" ; always portrait

; Plot boundaries.

  res_at_mpMinLonF = lonw
  res_at_mpMaxLonF = lone
  res_at_mpMinLatF = lats
  res_at_mpMaxLatF = latn

  res_at_mpCenterLonF = 0.5 * (lonw + lone) ; get accurate center longitude

  res_at_gsnAddCyclic = False ; must be false if not wrapping
; res_at_gsnAddCyclic = True ; must be false if not wrapping
                                                  ; data at Greenwich meridian

; res_at_mpOutlineBoundarySets = "AllBoundaries"
; res_at_mpOutlineBoundarySets = "USStates"

; res_at_mpOutlineBoundarySets = "geophysicalandusstates"; turn on states
  res_at_mpOutlineBoundarySets = "national" ; countries, continents, islands

  res_at_mpDataBaseVersion = "mediumres" ; select database
  res_at_mpDataSetName = "Earth..2"

;------------------------------------------------------------------------------

; Common setup for stream function overlay plot.

  res2 = True ; make resource object for plot2
  res2_at_gsnDraw = False ; also hold draw and auto frame
  res2_at_gsnFrame = False ; for overlay plot

  res2_at_cnLevelSelectionMode = "ManualLevels" ; manual contour levels
  res2_at_lbLabelAutoStride = True
  res2_at_cnMinLevelValF = -6 ; min level
  res2_at_cnMaxLevelValF = 6 ; max level
  res2_at_cnLevelSpacingF = 0.5 ; interval

  res2_at_gsnLeftString = "" ; no auto labels
  res2_at_gsnRightString = ""
  res2_at_gsnCenterString = ""
  res2_at_cnInfoLabelOrthogonalPosF = 0.15 ; move info label down
                                                  ; (to be fixed later)

;----------------------------------------------------------------------
; Loop on offset, for each frame.
;----------------------------------------------------------------------

  offset1= in_at_first_offset
  offset2= in_at_last_offset
  dims = dimsizes(in)
  noffsets = dims(0)
  
  do ti = 0, noffsets -1
    offset = offset1+ti
    print("offset in loop = " + offset)

    grid = in(ti,:,:) ; get composite grid for this offfset
    
    u = inu(ti,::2,::2) ; get U and V grids for this offset
    v = inv(ti,::2,::2)
; v = 0. ;NOTE SETTING V COMPONENT to 0

    res_at_gsnCenterString = sprinti ("offset %0i", offset)
    res_at_gsnStringFontHeightF = 0.018

    printVarSummary(in)
    degrees_north = in_at_key_lat
    degrees_west = 360 - in_at_key_lon
    comp_days = in_at_comp_count
    mblevel = in_at_level
    
; res_at_gsnLeftString = in_at_title3 + " (" + comp_days + " days)"
    title3_plus = in_at_title3 + " (" + comp_days + " days)"
    res_at_gsnLeftString = ""
    
    if (in_at_key_lat .eq. 0) then
       print_lat = "Equator"
    else
       if (in_at_key_lat .gt. 0) then
          print_lat = degrees_north + " deg N"
       else
          print_lat = -degrees_north + " deg S"
       end if
    end if
    
    res_at_tiMainString = in_at_title1 + "~C~ " \
      + in_at_title2 + " " + degrees_west + " deg W, " + print_lat \
      + "~C~ " + in_at_title4 + " " + in_at_title5 \
      + "~C~" + title3_plus
    
    res_at_tiMainOffsetYF = 0.005

    print (sprinti("offset %3i",offset) + sprintf (" min %5.2f", min(grid)) \
      + sprintf (" max %5.2f", max(grid)))

    res_at_gsnScalarContour = True ; contours desired

; res_at_vcRefMagnitudeF = 8.0 ; define vector ref mag ;200 mb
; res_at_vcRefMagnitudeF = 2.0 ; define vector ref mag 1000 mb
    res_at_vcRefMagnitudeF = 4.0 ; define vector ref mag 1000 mb
    res_at_vcRefLengthF = 0.05 ; define length of vec ref
    res_at_vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
    res_at_vcRefAnnoArrowLineColor = "black" ; change ref vector color
    res_at_vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref

    res_at_vcGlyphStyle = "CurlyVector" ; turn on curley vectors
    res_at_vcLineArrowColor = "black" ; change vector color
    res_at_vcLineArrowThicknessF = 2.0 ; change vector thickness
    res_at_vcVectorDrawOrder = "PostDraw" ; draw vectors last

    res_at_vcMinDistanceF = 0.005 ; Thin the vectors

    plot1 = gsn_csm_vector_scalar_map_ce (wks, u, v, grid, res)
; plot1 = gsn_csm_vector_scalar_map_ce (wks, u, v, u, res)
                                ; make contour plot with wind vector overlay.

; plot1 = gsn_csm_contour_map_ce (wks, grid, res) ; make next plot frame

;----------------------------------------------------------------------

; Add overlay contour lines for stream function.

    sf_vp = uv2sfvpF (inu(ti,:,:), inv(ti,:,:)) ; compute stream function
                                                    ; for current offset
    sfunc = sf_vp(0,:,:) / sf_scale_F ; get stream function alone
    ugrid = inu(ti,:,:) ; copy over lat and lon coords
    copy_VarCoords (ugrid, sfunc)
;; grid_at_units = sf_units ; SF not labeled in this version

    plot2 = gsn_csm_contour (wks, sfunc, res2) ; add stream func overlay
    overlay (plot1, plot2) ; on top of main plot

    draw (plot1)

;----------------------------------------------------------------------

; Draw dot at specified location.

; Multiple dots may be drawn with multiple wmlabs commands.
; Color and size may be changed between dots.

    dot_lat = degrees_north
    dot_lon = -degrees_west ; must be -180 to +180 ???

    wmsetp ("EZF", 1)
    wmsetp ("DBC", blue) ; set color of dot outer ring
    wmsetp ("DTC", red) ; set color of dot center
    wmsetp ("DTS", 0.012) ; set dot size

    wmlabs (wks, dot_lat, dot_lon, "DT") ; draw dot at given lat, lon

    print("Dot lat = " + dot_lat + ", lon = " + dot_lon)
    
;----------------------------------------------------------------------

    frame (wks) ; send final plot to output

    delete (grid) ; must delete working grid
    delete (plot1) ; and working plot objects
    delete (plot2) ; between loops
  
 end do

  print ("Done.")
  
end

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jul 11 2007 - 13:12:27 MDT

This archive was generated by hypermail 2.2.0 : Tue Jul 17 2007 - 06:52:00 MDT