;------------------------------------------------------------------------
;
; Temporary changes to add_shapefile_poly functions,
; to compensate for swapped X/Y bug in some builds of NCL 6.6.2.
;
; Original function versions are in ncarg/nclscripts/csm/gsn_code.ncl.
; These versions replace the originals via the "load" statement.
;
; 2022 April 6, by Dave Allured, NOAA/PSL/CIRES.
;
;------------------------------------------------------------------------

;***********************************************************************;
; Procedure : gsn_shapefile_polylines                                   ;
;                   wks: workstation object                             ;
;                plotid: plot object                                    ;
;                 fname: Name of shapefile ("xxxx.shp")                 ;
;               resources: optional resources                           ;
;                                                                       ;
; This code is deprecated. gsn_add_shapefile_polylines should be used.  ;
; ANY CHANGES MADE TO THIS CODE SHOULD POTENTIALLY BE MADE TO           ;
; gsn_add_shapefile_polylines BELOW!! Note that in V6.2.0, the "add"    ;
; version of this code was sped up significantly, so this procedure is  ;
; not needed so much anymore.                                           ;
;                                                                       ;
; This function draws shapefile polylines on the plot "plotid".         ;
; See gsn_add_shapefile_polylines if you want to add them instead (this ;
; can be much slower).                                                  ;
;                                                                       ;
; In version 6.1.0, some code was added to add checks if the lat/lon    ;
; segments are within the range of the map. This works best for a C.E.  ;
; map. You have to set the special min/max/lat/lon attributes for this  ;
; to work.  I won't advertise this yet, because the interface could     ;
; probably be made better. See note about V6.2.0 above.                 ;
;***********************************************************************;
undef("gsn_shapefile_polylines")
procedure gsn_shapefile_polylines(wks,plot,fname:string,lnres)
local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, i, lat, lon, \
startSegment, numSegments, seg, startPT, endPT
begin
;---Open the shapefile
  f = addfile(fname,"r")

;---Error checking
  if(ismissing(f)) then
    print("Error: gsn_shapefile_polylines: Can't open shapefile '" + fname + "'")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  elseif(.not.isatt(f,"geometry_type")) then
    print("Error: gsn_shapefile_polylines: this shapefile has no 'geometry_type' attribute.")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  elseif(f@geometry_type.eq."point") then
    print("Error: gsn_shapefile_polylines: geometry_type attribute is equal to 'point'.")
    print("       Use gsn_add_shapefile_polymarkers to plot this data.")
    return(new(1,graphic))
  elseif(.not.any(f@geometry_type.eq.(/"polygon","polyline"/))) then
    print("Warning: gsn_shapefile_polylines: geometry_type attribute is not equal to 'polygon'")
    print("         or 'polyline'. Will attempt to add shapefile outlines anyway.")
  end if

;---Read data off the shapefile
  geomDims    = getfilevardimsizes(f,"geometry")
  numFeatures = geomDims(0)
  if(numFeatures.eq.0) then
    print("Error: gsn_shapefile_polylines: the number of features in this file is 0.")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  end if

  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

;---Section to attach polylines to plot.

;;  lon = f->x		; ORIGINAL in NCL 6.6.2
;;  lat = f->y

  lon = f->y		; REVERSED TO FIX SWAPPED X/Y, 2022 April 6
  lat = f->x

; Special check for minlat/maxlat/minlon/maxlon attributes.
;
; If set, then each lat/lon segment will be checked if it's
; in the range.  This can speed up plotting, but I need to
; verify this!
; 
  if(isatt(lnres,"minlon").and.isatt(lnres,"maxlon").and.\
     isatt(lnres,"minlat").and.isatt(lnres,"maxlat")) then
    do i=0, numFeatures-1  
       startSegment = geometry(i, geom_segIndex)
       numSegments  = geometry(i, geom_numSegs)
       do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
          lat_sub = lat(startPT:endPT)
          lon_sub = lon(startPT:endPT) 
          if(.not.(all(lon_sub.lt.lnres@minlon).or. \
                   all(lon_sub.gt.lnres@maxlon).or. \
                   all(lat_sub.lt.lnres@minlat).or. \
                   all(lat_sub.gt.lnres@maxlat))) then
            gsn_polyline(wks, plot, lon_sub, lat_sub, lnres)
          end if
          delete([/lat_sub,lon_sub/])
       end do
    end do
  else       ; Don't do any range checking. 
    do i=0, numFeatures-1  
       startSegment = geometry(i, geom_segIndex)
       numSegments  = geometry(i, geom_numSegs)
       do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
          gsn_polyline(wks, plot, lon(startPT:endPT),  \
                                  lat(startPT:endPT), lnres)
       end do
    end do
  end if
end

;***********************************************************************;
; Function : gsn_add_shapefile_polylines                                ;
;                   wks: workstation object                             ;
;                plotid: plot object                                    ;
;                 fname: Name of shapefile ("xxxx.shp")                 ;
;               resources: optional resources                           ;
;                                                                       ;
; ANY CHANGES MADE TO THIS CODE SHOULD POTENTIALLY BE MADE TO           ;
; gsn_shapefile_polylines ABOVE!!                                       ;
;                                                                       ;
; This function attaches shapefile polylines to the plot "plotid".      ;
; See gsn_shapefile_polylines if you want to just draw them.            ; 
;                                                                       ;
; In version 6.2.0 this function was significantly sped-up via the use  ;
; of a new resource called "gsSegments". Only one primitive object id   ;
; is now created. In some cases, we saw speed-ups of 80x (for example,  ;
; with the France_adm5 shapefile). Because of this speed-up, the code   ;
; added in V6.1.0 for checking for special minlat/maxlat/minlon/maxlon  ;
; is no longer needed. It is commented out for now. It may eventually   ;
; be removed once we are certain it doesn't serve any benefit.          ;
;                                                                       ;
; In version 6.1.0, some code was added to add checks if the lat/lon    ;
; segments are within the range of the map. This works best for a C.E.  ;
; map. You have to set the special min/max/lat/lon attributes for this  ;
; to work.  I won't advertise this yet, because the interface could     ;
; probably be made better. THIS CODE HAS BEEN COMMENTED OUT in V6.2.0.  ;
;                                                                       ;
; In version 6.2.1, this function was modified to allow multiple plots. ;
;***********************************************************************;
undef("gsn_add_shapefile_polylines")
function gsn_add_shapefile_polylines(wks,plots[*]:graphic,fname:string,lnres)
local f, geomDims, numFeatures, lnres2, valid_shapefile
begin
;---Open the shapefile
  f = addfile(fname,"r")

  valid_shapefile = is_valid_shapefile(f,"polyline")
  if(.not.valid_shapefile) then
    return(new(1,graphic))
  end if

  if(lnres) then
    lnres2 = lnres      ; Make a copy so that we don't keep gsSegments
  else
    lnres2 = True
  end if

;;  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
;;
;;;---Create array to hold all polylines
;;  npoly = sum(geometry(:,geom_numSegs)) 
;;  poly  = new(npoly,graphic)
;;
;;;---Section to attach polylines to plot.
;;  lon = f->x
;;  lat = f->y
;;  npl = 0     ; polyline counter
;;;
;;; Special check for minlat/maxlat/minlon/maxlon attributes.
;;;
;;; If set, then each lat/lon segment will be checked if it's
;;; in the range.  This can speed up plotting, but I need to
;;; verify this!
;;; 
;;  if(isatt(lnres,"minlon").and.isatt(lnres,"maxlon").and.\
;;     isatt(lnres,"minlat").and.isatt(lnres,"maxlat")) then
;;    do i=0, numFeatures-1  
;;       startSegment = geometry(i, geom_segIndex)
;;       numSegments  = geometry(i, geom_numSegs)
;;       do seg=startSegment, startSegment+numSegments-1
;;          startPT = segments(seg, segs_xyzIndex)
;;          endPT   = startPT + segments(seg, segs_numPnts) - 1
;;          lat_sub = lat(startPT:endPT)
;;          lon_sub = lon(startPT:endPT) 
;;          if(.not.(all(lon_sub.lt.lnres@minlon).or. \
;;                   all(lon_sub.gt.lnres@maxlon).or. \
;;                   all(lat_sub.lt.lnres@minlat).or. \
;;                   all(lat_sub.gt.lnres@maxlat))) then
;;            poly(npl) = gsn_add_polyline(wks, plot, lon_sub, lat_sub, lnres)
;;            npl = npl + 1
;;          end if
;;          delete([/lat_sub,lon_sub/])
;;       end do
;;    end do
;;  else       ; Don't do any range checking. 
;;    do i=0, numFeatures-1  
;;       startSegment = geometry(i, geom_segIndex)
;;       numSegments  = geometry(i, geom_numSegs)
;;       do seg=startSegment, startSegment+numSegments-1
;;          startPT = segments(seg, segs_xyzIndex)
;;          endPT   = startPT + segments(seg, segs_numPnts) - 1
;;          poly(npl) = gsn_add_polyline(wks, plot, lon(startPT:endPT),  \
;;                                       lat(startPT:endPT), lnres)
;;          npl = npl + 1
;;       end do
;;    end do
;;  end if
;;  return(poly(0:npl-1))

;---This is all that's needed in V6.2.0.
  lnres2@gsSegments = f->segments(:,0)
  nplots = dimsizes(plots)
  poly   = new(nplots,graphic)
  do i=0,nplots-1

;;  poly(i) = gsn_add_polyline(wks, plots(i), f->x, f->y, lnres2)
				; ORIGINAL in NCL 6.6.2

    poly(i) = gsn_add_polyline(wks, plots(i), f->y, f->x, lnres2)
				; REVERSED TO FIX SWAPPED X/Y, 2022 April 6
  end do
  return(poly)
end

;***********************************************************************;
; Function : gsn_add_shapefile_polygons                                 ;
;                   wks: workstation object                             ;
;                plotid: plot object                                    ;
;                 fname: Name of shapefile ("xxxx.shp")                 ;
;               resources: optional resources                           ;
;                                                                       ;
; This function attaches shapefile polygons to the plot "plotid".       ;
;                                                                       ;
; In version 6.2.0 this function was significantly sped-up via the use  ;
; of a new resource called "gsSegments". Only one primitive object id   ;
; is now created. In some cases, we saw speed-ups of 80x (for example,  ;
; with the France_adm5 shapefile). Because of this speed-up, the code   ;
; added in V6.1.0 for checking for special minlat/maxlat/minlon/maxlon  ;
; is no longer needed. It is commented out for now. It may eventually   ;
; be removed once we are certain it doesn't serve any benefit.          ;
;                                                                       ;
; In version 6.1.0, some code was added to add checks if the lat/lon    ;
; segments are within the range of the map. This works best for a C.E.  ;
; map. You have to set the special min/max/lat/lon attributes for this  ;
; to work.  I won't advertise this yet, because the interface could     ;
; probably be made better. THIS CODE HAS BEEN COMMENTED OUT in V6.2.0.  ;
;                                                                       ;
; In version 6.2.1, this function was modified to allow multiple plots. ;
;***********************************************************************;
undef("gsn_add_shapefile_polygons")
function gsn_add_shapefile_polygons(wks,plots[*]:graphic,fname:string,gnres)
local f, geomDims, numFeatures, gnres2
begin
;---Open the shapefile
  f = addfile(fname,"r")

  valid_shapefile = is_valid_shapefile(f,"polygon")
  if(.not.valid_shapefile) then
    return(new(1,graphic))
  end if

;---Read data off the shapefile
  geomDims    = getfilevardimsizes(f,"geometry")
  numFeatures = geomDims(0)
  if(numFeatures.eq.0) then
    print("Error: gsn_add_shapefile_polygons: the number of features in this file is 0.")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  end if

  if(gnres) then
    gnres2 = gnres      ; Make a copy so that we don't keep gsSegments
  else
    gnres2 = True
  end if

;---Get the number of colors
  if(.not.gnres2.or.(.not.isatt(gnres2,"gsFillColor").and.\
                     .not.isatt(gnres2,"gsColors"))) then
    getvalues wks
      "wkColorMapLen"  : cmap_len
    end getvalues
    gnres2          = True   ; Make sure this is True
    gnres2@gsColors = toint(random_uniform(2,cmap_len-2,numFeatures))
  end if

;;  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
;;
;;;---Create array to hold all polylines
;;  npoly = sum(geometry(:,geom_numSegs)) 
;;  poly  = new(npoly,graphic)
;;
;;;---Section to attach polygons to plot.
;;  lon = f->x
;;  lat = f->y
;;  npl = 0     ; polyline counter
;;;
;;; Special check for minlat/maxlat/minlon/maxlon attributes.
;;;
;;; If set, then each lat/lon segment will be checked if it's
;;; in the range.  This can speed up plotting, but I need to
;;; verify this!
;;; 
;;  if(isatt(gnres,"minlon").and.isatt(gnres,"maxlon").and.\
;;     isatt(gnres,"minlat").and.isatt(gnres,"maxlat")) then
;;    do i=0, numFeatures-1  
;;       startSegment = geometry(i, geom_segIndex)
;;       numSegments  = geometry(i, geom_numSegs)
;;       do seg=startSegment, startSegment+numSegments-1
;;          startPT = segments(seg, segs_xyzIndex)
;;          endPT   = startPT + segments(seg, segs_numPnts) - 1
;;          lat_sub = lat(startPT:endPT)
;;          lon_sub = lon(startPT:endPT) 
;;          if(set_fill_color) then
;;;---Pick a random color 
;;            gnres@gsFillColor = toint(random_uniform(2,cmap_len-2,1))
;;          end if
;;          if(.not.(all(lon_sub.lt.gnres@minlon).or. \
;;                   all(lon_sub.gt.gnres@maxlon).or. \
;;                   all(lat_sub.lt.gnres@minlat).or. \
;;                   all(lat_sub.gt.gnres@maxlat))) then
;;;---Attach the line segment
;;            poly(npl) = gsn_add_polygon(wks, plot, lon_sub, lat_sub, gnres)
;;            npl = npl + 1
;;          end if
;;          delete([/lat_sub,lon_sub/])
;;       end do
;;    end do
;;  else
;;    do i=0, numFeatures-1  
;;       startSegment = geometry(i, geom_segIndex)
;;       numSegments  = geometry(i, geom_numSegs)
;;       do seg=startSegment, startSegment+numSegments-1
;;          startPT = segments(seg, segs_xyzIndex)
;;          endPT   = startPT + segments(seg, segs_numPnts) - 1
;;          if(set_fill_color) then
;;;---Pick a random color 
;;            gnres@gsFillColor = toint(random_uniform(2,cmap_len-2,1))
;;          end if
;;;---Attach the line segment
;;          poly(npl) = gsn_add_polygon(wks, plot, lon(startPT:endPT),  \
;;                                      lat(startPT:endPT), gnres)
;;          npl = npl + 1
;;       end do
;;    end do
;;  end if
;;  return(poly)

;---This is all that's needed in V6.2.0.
  gnres2@gsSegments = f->segments(:,0)
  nplots = dimsizes(plots)
  poly   = new(nplots,graphic)
  do i=0,nplots-1
;;  poly(i) = gsn_add_polygon(wks, plots(i), f->x, f->y, gnres2)
				; ORIGINAL in NCL 6.6.2

    poly(i) = gsn_add_polygon(wks, plots(i), f->y, f->x, gnres2)
				; REVERSED TO FIX SWAPPED X/Y, 2022 April 6
  end do
  return(poly)
end

;***********************************************************************;
; Function : gsn_add_shapefile_polymarkers                              ;
;                   wks: workstation object                             ;
;                plotid: plot object                                    ;
;                 fname: Name of shapefile ("xxxx.shp")                 ;
;               resources: optional resources                           ;
;                                                                       ;
; This function attaches shapefile point data to the plot "plotid".     ;
;                                                                       ;
; In version 6.1.0, some code was added to add checks if the lat/lon    ;
; segments are within the range of the map. This works best for a C.E.  ;
; map. You have to set the special min/max/lat/lon attributes for this  ;
; to work.  I won't advertise this yet, because the interface could     ;
; probably be made better.                                              ;
;                                                                       ;
; In version 6.2.1, this function was modified to allow multiple plots. ;
;***********************************************************************;
undef("gsn_add_shapefile_polymarkers")
function gsn_add_shapefile_polymarkers(wks,plots[*]:graphic,fname:string,mkres)
local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, i, lat, lon, \
startSegment, numSegments, seg, startPT, endPT, npoly, npl
begin
;---Open the shapefile
  f = addfile(fname,"r")

  valid_shapefile = is_valid_shapefile(f,"point")
  if(.not.valid_shapefile) then
    return(new(1,graphic))
  end if

;---Read data off the shapefile
  geomDims    = getfilevardimsizes(f,"geometry")
  numFeatures = geomDims(0)
  if(numFeatures.eq.0) then
    print("Error: gsn_add_shapefile_polymarkers: the number of features in this file is 0.")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  end if

  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

;---Create array to hold all polymarkers
  npoly  = sum(geometry(:,geom_numSegs)) 
  nplots = dimsizes(plots)
  poly   = new(npoly*nplots,graphic)

;---Section to attach polymarkers to plots.

;;  lon = f->x		; ORIGINAL in NCL 6.6.2
;;  lat = f->y

  lon = f->y		; REVERSED TO FIX SWAPPED X/Y, 2022 April 6
  lat = f->x

  npl = 0     ; polyline counter
;
; Special check for minlat/maxlat/minlon/maxlon attributes.
;
; If set, then each lat/lon segment will be checked if it's
; in the range.  This can speed up plotting, but I need to
; verify this!
; 
  if(isatt(mkres,"minlon").and.isatt(mkres,"maxlon").and.\
     isatt(mkres,"minlat").and.isatt(mkres,"maxlat")) then
    do i=0, numFeatures-1  
       startSegment = geometry(i, geom_segIndex)
       numSegments  = geometry(i, geom_numSegs)
       do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
          lat_sub = lat(startPT:endPT)
          lon_sub = lon(startPT:endPT) 
          if(.not.(all(lon_sub.lt.mkres@minlon).or. \
                   all(lon_sub.gt.mkres@maxlon).or. \
                   all(lat_sub.lt.mkres@minlat).or. \
                   all(lat_sub.gt.mkres@maxlat))) then
;---Attach the markers
            do n=0,nplots-1
              poly(npl) = gsn_add_polymarker(wks, plots(n), lon_sub, lat_sub, mkres)
              npl = npl + 1
            end do
           end if
           delete([/lat_sub,lon_sub/])
         end do
    end do
  else       ; Don't do any range checking. 
    do i=0, numFeatures-1  
       startSegment = geometry(i, geom_segIndex)
       numSegments  = geometry(i, geom_numSegs)
       do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
;---Attach the markers
          do n=0,nplots-1
            poly(npl) = gsn_add_polymarker(wks, plots(n), lon(startPT:endPT),  \
                                           lat(startPT:endPT), mkres)
            npl = npl + 1
          end do
       end do
    end do
  end if
  return(poly(0:npl-1))
end

;***********************************************************************;
; Function : gsn_add_shapefile_text                                     ;
;                   wks: workstation object                             ;
;                plotid: plot object                                    ;
;                 fname: Name of shapefile ("xxxx.shp")                 ;
;                 vname: Name of string variable containing text strings;
;               resources: optional resources                           ;
;                                                                       ;
; This function attaches shapefile text strings to the plot "plotid".   ;
; The assumption is that there are "num_features" text strings, and this;
; routine gets the approximate mid lat/lon area for each text string.   ;
;***********************************************************************;
undef("gsn_add_shapefile_text")
function gsn_add_shapefile_text(wks,plot,fname[1]:string,vname[1]:string,txres)
local f, segments, geometry, segsDims, geomDims, geom_segIndex, \
geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures, i, lat, lon, \
startSegment, numSegments, seg, startPT, endPT, ntxt, \
minlat, maxlat, minlon, maxlon
begin
;---Open the shapefile
  f = addfile(fname,"r")
  if(.not.isfilevar(f,vname))
    print("Error: gsn_add_shapefile_text: '" + vname + "' is not a variable")
    print("       in file '" + fname + "'.")
    return(new(1,graphic))
  end if
;---Error checking
  if(ismissing(f)) then
    print("Error: gsn_add_shapefile_text: Can't open shapefile '" + \
           fname + "'")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  end if

;---Read data off the shapefile
  geomDims    = getfilevardimsizes(f,"geometry")
  numFeatures = geomDims(0)
  if(numFeatures.eq.0) then
    print("Error: gsn_add_shapefile_text: the number of features in this file is 0.")
    print("       No shapefile information will be added.")
    return(new(1,graphic))
  end if

  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
  numFeatures   = geomDims(0)

;---Create array to hold all text
  text  = new(numFeatures,graphic)

;---Section to attach text to plot.
;;  lon = f->x		; ORIGINAL in NCL 6.6.2
;;  lat = f->y

  lon = f->y		; REVERSED TO FIX SWAPPED X/Y, 2022 April 6
  lat = f->x

  ntxt = 0     ; text counter
;
; Special check for minlat/maxlat/minlon/maxlon attributes.
;
; If set, then each lat/lon segment will be checked if it's
; in the range.  This can speed up plotting, but I need to
; verify this!
;
  if(isatt(txres,"minlon").and.isatt(txres,"maxlon").and.\
     isatt(txres,"minlat").and.isatt(txres,"maxlat")) then
    do i=0, numFeatures-1
       startSegment = geometry(i, geom_segIndex)
       numSegments  = geometry(i, geom_numSegs)
       minlat = new(1,typeof(lat))
       maxlat = new(1,typeof(lat))
       minlon = new(1,typeof(lon))
       maxlon = new(1,typeof(lon))
       do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
          lat_sub = lat(startPT:endPT)
          lon_sub = lon(startPT:endPT)
          if(.not.(all(lon_sub.lt.txres@minlon).or. \
                   all(lon_sub.gt.txres@maxlon).or. \
                   all(lat_sub.lt.txres@minlat).or. \
                   all(lat_sub.gt.txres@maxlat))) then
            if(any((/ismissing(minlat),ismissing(maxlat), \
                     ismissing(minlon),ismissing(maxlon)/))) then
              minlat = min(lat_sub)
              maxlat = max(lat_sub)
              minlon = min(lon_sub)
              maxlon = max(lon_sub)
            else
              minlat = min((/minlat,min(lat_sub)/))
              maxlat = max((/maxlat,max(lat_sub)/))
              minlon = min((/minlon,min(lon_sub)/))
              maxlon = max((/maxlon,max(lon_sub)/))
            end if
          end if
          delete([/lat_sub,lon_sub/])
       end do
;---Attach the text string
       if(.not.any((/ismissing(minlat),ismissing(maxlat), \
                     ismissing(minlon),ismissing(maxlon)/))) then
         avglat = (minlat+maxlat)/2.
         avglon = (minlon+maxlon)/2.
         print("Text = '" + f->$vname$(i) + "'")
         print("Location = " + avglat + "/" + avglon)
         text(ntxt) = gsn_add_text(wks, plot,f->$vname$(i),avglon,avglat,txres)
         ntxt = ntxt + 1
       end if
    end do
  else
    do i=0, numFeatures-1  
       startSegment = geometry(i, geom_segIndex)
       numSegments  = geometry(i, geom_numSegs)
       minlat = new(1,typeof(lat))
       maxlat = new(1,typeof(lat))
       minlon = new(1,typeof(lon))
       maxlon = new(1,typeof(lon))
       do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
          if(any((/ismissing(minlat),ismissing(maxlat), \
                   ismissing(minlon),ismissing(maxlon)/))) then
            minlat = min(lat(startPT:endPT))
            maxlat = max(lat(startPT:endPT))
            minlon = min(lon(startPT:endPT))
            maxlon = max(lon(startPT:endPT))
          else
            minlat = min((/minlat,min(lat(startPT:endPT))/))
            maxlat = max((/maxlat,max(lat(startPT:endPT))/))
            minlon = min((/minlon,min(lon(startPT:endPT))/))
            maxlon = max((/maxlon,max(lon(startPT:endPT))/))
          end if
       end do
;---Attach the text string
       if(.not.any((/ismissing(minlat),ismissing(maxlat), \
                     ismissing(minlon),ismissing(maxlon)/))) then
         avglat = (minlat+maxlat)/2.
         avglon = (minlon+maxlon)/2.
         text(ntxt) = gsn_add_text(wks,plot,f->$vname$(i),avglon,avglat,txres)
         ntxt = ntxt + 1
       end if
    end do
  end if
  return(text)
end