;----------------------------------------------------------------------
; This procedure fills each feature in a given shapefile in a color
; based on a data average. A shapefile is used to get the 
; lat/lon outlines for each feature. 
;
; This code currently only works on rectilinear, curvilinear, or
; unstructured grids.
;
;----------------------------------------------------------------------
undef("plot_averages_using_shapefile_features")
procedure plot_averages_using_shapefile_features(wks,plot,data,fname:string,\
                                                 levels,fill_colors)
local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, i, lat, lon, \
startSegment, numSegments, seg, startPT, endPT, dims, 
minlat, maxlat, minlon, maxlon
begin
;
; Get the lat/lon area of interest so we can cull the lat/lon
; pairs that we have to check later.
;
  getvalues plot
    "mpMinLatF" : minlat
    "mpMaxLatF" : maxlat
    "mpMinLonF" : minlon
    "mpMaxLonF" : maxlon
  end getvalues

  dims = dimsizes(data)
  rank = dimsizes(dims)

;
; Determine the grid type.
;
; Get every lat/lon pair of the original data array, and
; turn them into 1D arrays.  We will use these 1D arrays
; when checking which lat/lon pairs fall in the given
; feature in the shapefile.
;
  grid_type = ""
  if(rank.eq.2.and.\
     isdimnamed(data,0).and.iscoord(data,data!0).and.\
     isdimnamed(data,1).and.iscoord(data,data!1)) then
    lat1d = ndtooned(conform_dims(dims,data&$data!0$,0))
    lon1d = ndtooned(conform_dims(dims,data&$data!1$,1))
    grid_type = "rectilinear"
  else if(rank.eq.2.and.all(isatt(data,(/"lat2d","lon2d"/)))) then
;---Curvilinear
    lat1d = ndtooned(data@lat2d)
    lon1d = ndtooned(data@lon2d)
    if(product(dims).eq.dimsizes(lat1d).and.\
       product(dims).eq.dimsizes(lon1d)) then
      grid_type = "curvilinear"
    end if
  else if(rank.eq.1.and.all(isatt(data,(/"lat1d","lon1d"/)))) then
;---Unstructured
    lat1d = data@lat1d
    lon1d = data@lon1d
    if(dims.eq.dimsizes(lat1d).and.\
       product(dims).eq.dimsizes(lon1d)) then
      grid_type = "unstructured"
    end if
  end if
  end if
  end if

  if(grid_type.eq."") then
    print("plot_averages_using_shapefile_features: Error")
    print("       Not a valid rectilinear, curvilinear, or unstructured grid")
    exit
  end if

  data1d  = ndtooned(data)
  nlatlon = dimsizes(lat1d)

;
; Get the approximate index values that contain the area of interest.
;
; This will make our gc_inout loop below go faster, because we're
; not checking every single lat/lon point to see if it's within
; the area of interest.
;
  ii_latlon = ind(lat1d.ge.minlat.and.lat1d.le.maxlat.and.\
                  lon1d.ge.minlon.and.lon1d.le.maxlon)
  nii_latlon = dimsizes(ii_latlon)

;---Open the shapefile
  f = addfile(fname,"r")

;---Read data off the shapefile
  geomDims    = getfilevardimsizes(f,"geometry")
  numFeatures = geomDims(0)
  segments    = f->segments
  geometry    = f->geometry
  segsDims    = dimsizes(segments)

;---Read global attributes  
  geom_segIndex = f@geom_segIndex
  geom_numSegs  = f@geom_numSegs
  segs_xyzIndex = f@segs_xyzIndex
  segs_numPnts  = f@segs_numPnts

;---Read the lat/lon data off the shapefile
  lat = f->y
  lon = f->x

  if( (min(lon).lt.0.and.max(lon1d).gt.180).or.\
      (min(lon1d).lt.0.and.max(lon).gt.180)) then
    print("plot_averages_using_shapefile_features: Error: your data longitude array and")
    print("shapefile longitude values are not in the same range. You may need to use 'lonFlip'")
  end if

;---Create array to hold new data mask, and set all values to 0 initially.
  data_avg_1d = new(nlatlon,integer)
  skip_check  = new(nlatlon,logical)
  skip_check  = False

;---Make sure "data" has a missing value.
  if(.not.isatt(data,"_FillValue")) then
    data@_FillValue = default_fillvalue(typeof(data))
  end if

  data_sum = new(1,typeof(data))    ; Scalar to hold data sums for average.

;
; This is the loop that loops across every shapefile feature
; and sums all the data whose lat/lon points fall inside the 
; given feature.
;
; It then takes an average of the data values for this feature
; and attaches it as a filled polygon to the given plot.
;
  gnres = True                       ; polygon resource list
  nfill = dimsizes(fill_colors)

  do i=0, numFeatures-1
;---Some features have multiple segments
    startSegment = geometry(i, geom_segIndex)
    numSegments  = geometry(i, geom_numSegs)

;
; Loop through each segment of a climate division, and for each point
; inside a segment, add it to the running sum.
;
    data_sum  = 0
    ndata_sum = 0
    do seg=startSegment, startSegment+numSegments-1
      startPT = segments(seg, segs_xyzIndex)
      endPT   = startPT + segments(seg, segs_numPnts) - 1

;---Loop through each point on data grid and see if it's in this feature's segment
      do n=0,nii_latlon-1
        nn = ii_latlon(n)     ; Get index of point we're checking
;
; This "if" statement *significantly speeds up this code, by making sure
; we don't needlessly check a lat/lon point if:
;
;   - We've already found it in another feature.
;   - It doesn't fall within the general lat/lon box that covers this feature.
; 
        if(skip_check(nn).or.\
           lat1d(nn).lt.min(lat(startPT:endPT)).or.\
           lat1d(nn).gt.max(lat(startPT:endPT)).or.\
           lon1d(nn).lt.min(lon(startPT:endPT)).or.\
           lon1d(nn).gt.max(lon(startPT:endPT))) 
           continue
        end if

;---Here's the check if the point is in this lat/lon segment.
        if(gc_inout(lat1d(nn),lon1d(nn),\
                    lat(startPT:endPT),lon(startPT:endPT))) then
          ndata_sum = ndata_sum + 1
          data_sum  = data_sum + data1d(nn)
          skip_check(nn)  = True   ; Don't check this point again
        end if
      end do
    end do    ; End of looping across segment to collect points for this feature.

;
; If there are > 0 values found in this feature, then calculate
; the average, set a fill color for this feature, and do the
; fill for each segment.
;
    if(ndata_sum.gt.0) then
      data_avg = data_sum / ndata_sum
;; Debug prints
;;     print("Average over this feature = " + data_avg)

;---Set the fill color for this feature to the appropriate color
      iiavg := ind(data_avg.lt.levels)
      if(.not.any(ismissing(iiavg))) then
        gnres@gsFillColor = fill_colors(iiavg(0))
      else
        gnres@gsFillColor = fill_colors(nfill-1)
      end if
      do seg=startSegment, startSegment+numSegments-1
        startPT = segments(seg, segs_xyzIndex)
        endPT   = startPT + segments(seg, segs_numPnts) - 1

        str = unique_string("polygon")
        plot@$str$ = gsn_add_polygon(wks,plot,lon(startPT:endPT),\
                                              lat(startPT:endPT),gnres)
      end do
    end if
  end do
end