[ncl-talk] Adding a Legend to Shapefile Plot

[Carlos J. Valle-Diaz]

Hello,

I've been running WRF simulations to observe 24hr precipitation in Puerto
Rico from (July 12, 2011 12UTC to July 13, 2011 12UTC). I want to compare
the model runs with the observed from the National Weather Service. I have
their hourly shapefiles, see link to access them:
http://www.srh.noaa.gov/ridge2/Precip/qpehourlyshape/2011/201107/20110712/).
In the shapefiles, Globvalue refers to precipitation.

I've been trying all day to add a legend to my plot, but haven't been
successful. I've been using example #7 in the legends section of NCL as a
reference. See leg_7.ncl at http://ncl.ucar.edu/Applications/legend.shtml

I have attached the script I have at the moment. It plots fine when you
take out the lines that codes for the legend and text. What am I missing
and/or doing wrong?

Thanks beforehand for the help,

P.S. I set up the array(arr) data for 16 values because I'm using a color
map that displays 17 colors.
-- 
Carlos J. Valle Diaz
Ph.D. Chemistry Student
University of Puerto Rico Rio Piedras Campus
Department of Graduate Chemistry and Environmental Sciences
cj.vallediaz at gmail.com
Tel. 787-764-0000 x88192
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151218/ab1663a5/attachment.html 
-------------- next part --------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
  wks  = gsn_open_wks("X11","NWS-Precip-2011071215")

  f = addfile("nws_precip_2011071215.shp","r")

  minLat = 17.30
  maxLat = 19
  latRange = maxLat - minLat
  minLon = 292
  maxLon = 295
  lonRange = maxLon - minLon

  globMin = min(f->Globvalue)
  globMax = max(f->Globvalue)
  globDiv = globMax / 16
  globRange = globMax - globMin

  arr = (/globMin, globMin + globDiv, globMin + 2 * globDiv, globMin + 3 * globDiv, globMin + 4 * globDiv, globMin + 5 * globDiv, globMin + 6 * globDiv, globMin + 7 * globDiv, globMin + 8 * globDiv, globMin + 9 * globDiv, globMin + 10 * globDiv, globMin + 11 * globDiv, globMin + 12 * globDiv, globMin + 13 * globDiv, globMin + 14 * globDiv, globMin + 15 * globDiv/)
  labels = new(dimsizes(arr)+1,string)
  num_distinct_markers = dimsizes(arr)+1      ; number of distinct markers

 ; Group the points according to which range they fall in. At the
 ; same time, create the label that will be used later in the legend.

  do i = 0, num_distinct_markers-1
    if (i.eq.0) then
      indexes = ind(R.lt.arr(0))
      labels(i) = "x < " + arr(0)
    end if
    if (i.eq.num_distinct_markers-1) then
      indexes = ind(R.ge.max(arr))
      labels(i) = "x >= " + max(arr)
    end if
    if (i.gt.0.and.i.lt.num_distinct_markers-1) then
      indexes = ind(R.ge.arr(i-1).and.R.lt.arr(i))
      labels(i) = arr(i-1) + " <= x < " + arr(i)
    end if

  mres                       = True
  mres at mpDataBaseVersion     = "MediumRes"
  mres at mpLimitMode           = "Corners"     ; corner method of zoom
  mres at mpLeftCornerLatF      = minLat - (latRange * .01)
  mres at mpLeftCornerLonF      = minLon - (lonRange * .01)
  mres at mpRightCornerLatF     = maxLat + (latRange * .01)
  mres at mpRightCornerLonF     = maxLon + (lonRange * .01)

  mres at gsnDraw               = True          ; don't draw yet
  mres at gsnFrame              = False         ; don't advance frame yet
  mres at gsnMaximize           = False
  mres at pmTickMarkDisplayMode = "Always"      ; turn on tickmarks
  mres at tiMainString          = "NWS Precipitation July 12 2011 at 15 UTC"

  plot = gsn_csm_map(wks,mres)

  cmap = read_colormap_file("prcp_1")
  numColors = dimsizes(cmap(:,0)) - 1        ; skip the first color in the table

  print(globMin + " " + globRange)

  txres = True
  txres at txFontHeightF = 0.015

  ; Loop through each grouping of markers, and draw them one set at
  ; a time, assigning the proper color and size with gsn_marker.
  ; At the same time, draw a legend showing the meaning of the markers.

  gsres = True

  xleg = (/0.07,0.07,0.32,0.32,0.57,0.57,0.82,0.82/)   ; Location of
  xtxt = (/0.16,0.16,0.41,0.41,0.66,0.66,0.91,0.91/)   ; legend markers
  yleg = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/)   ; and text
  ytxt = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/)   ; strings.

  do i = 0, num_distinct_markers-1
    if (.not.ismissing(f->Lat(i)))
      gsres at gsMarkerColor      = numColors
      gsres at gsMarkerThicknessF = 0.7*(i+1)
      gsn_polymarker(wks, plot, f->Lon(i), f->Lat(i), gsres)

; Add marker and text for the legend.

      gsn_polymarker_ndc(wks,          xleg(i),yleg(i),gsres)
      gsn_text_ndc      (wks,labels(i),xtxt(i),ytxt(i),txres)
    end if
  end do

   do i=0, dimsizes(f->Globvalue) - 1

      indx = doubletoint( (f->Globvalue(i)-globMin) / globRange * numColors) + 1  ; skip first color in table
      if (indx.ge.numColors) then
          indx = numColors - 1
      end if

       gsres at gsMarkerIndex = 16
       gsres at gsMarkerSizeF = 0.005
       gsres at gsMarkerColor = cmap(indx, :)
       gsn_polymarker(wks, plot, f->Lon(i), f->Lat(i), gsres)

   end do

 draw(plot)
 frame(wks)

end