;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; track_map.ncl / ibtracs_2.ncl
;; Carl Schreck (carl@cicsnc.org) / Lily Kailyn (lkailyn@ucar.edu)
;; February 2012 / October 2016
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Description: Draw a blank map and put genesis locations on it
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; These files are loaded by default in NCL V6.2.0 and newer
; 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/contributed.ncl"
;
; These files still have to be loaded manually
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

;;****************************************************************
;; These are some parameters that could be useful to have up top
;;****************************************************************
  plotType = "pdf"                 ; send graphics to PNG file
  plotName = "Tracks"
  plotDpi  = 150

;;****************************************************************
;;****************************************************************
  ibtDir  = "./"
  ibtPath = "Allstorms_ibtracs_wmo_v03r04.nc" 

;;****************************************************************
;; region
;;****************************************************************
  minLon = 100
  maxLon = 180
  minLat = 0
  maxLat = 30
  centerLon = 180

  basin = 2 ;; -1=All, 0=North Atlantic, 1=South Atlantic, 2=West Pacific,
             ;; 3=East Pacific, 4=South Pacific, 5=North Indian, 6=South Indian


;;****************************************************************
;; Read the data
;;****************************************************************
  inFile = addfile( ibtPath, "r" )
  Lat = short2flt( inFile->lat_wmo )
  Lon = short2flt( inFile->lon_wmo )
  Lon = where( Lon.lt.0, Lon+360, Lon )
  wind = inFile->wind_wmo
  type = inFile->nature_wmo

  printMinMax(Lat, True)
  printMinMax(Lon, True)
  

;;****************************************************************
;; Customize base plot
;;****************************************************************
  res = True
  res@gsnMaximize            = True
  res@gsnDraw  = False
  res@gsnFrame = False
  res@gsnPaperOrientation = "Portrait"
  res@tiMainString = "Tropical Cyclones over the Oceans (1848-2012)"
  res@pmLegendDisplayMode = "NoCreate"
 ; res@tiXAxisString         = "latitude"           ; x-axis label
 ; res@tiYAxisString         = "longitude"         ; y-axis label 

  ;; ...set the spacing of the lat/lon labels
  res@gsnMajorLonSpacing = 20
  res@gsnMinorLonSpacing = 10
  res@gsnMajorLatSpacing = 20
  res@gsnMinorLatSpacing = 10

  ;; ...map gridlines
  res@mpGridAndLimbOn = False
  res@mpGridLatSpacingF = 0.0
  res@mpGridLonSpacingF = 0.0
  ;res@mpGridMaskMode    = "MaskLand"
  res@mpGridLineColor   = "Gray30"

  ;; ...set the bounds of a map plot
  res@mpMinLatF    = minLat
  res@mpMaxLatF    = maxLat
  res@mpMinLonF    = minLon
  res@mpMaxLonF    = maxLon
  res@mpCenterLonF = centerLon

  ;; ...set map resources
  res@mpFillOn = False
  res@mpGeophysicalLineThicknessF = 5
  res@mpGeophysicalLineColor      = "black"
  res@mpNationalLineColor         = "blue"
  res@mpUSStateLineColor          = "blue"
  ;;res@mpOutlineBoundarySets       = "AllBoundaries"
; res@mpDataBaseVersion           = "MediumRes"
; res@mpDataSetName               = "Earth..1"


  ;; Set up an x-y line graph
  xyRes = True
  xyRes@gsnDraw  = False
  xyRes@gsnFrame = False
  ;xyRes@xyMarkLineMode     = "MarkLines"
  ;xyRes@xyMonoDashPattern  = True
  xyRes@xyLineColor        = red
  xyRes@xyLineThicknessF   = 1
  ;xyRes@xyMarker           = 1
  ;xyRes@xyMarkerColor      = "black"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
  
     res_poly = True ; polyline mods desired 
    res_poly@gsEdgesOn = True ; draw border around polygons
  res_poly@gsEdgeColor = "black" 

   TypeOn                    = True  ; If you just want the legend to have the basic (H, TS, and TD)
  PlotSubTropicalSegments   = True   ; If you wish to exclude all subtropical line segments, set to False 
   
; Note: if TypeOn is set to False, the following (W, L, E) cases won't be plotted regardless of their value   
   PlotWaveSegments          = True   ; If you wish to exclude all wave segments, set to False
   PlotLowSegments	     = True   ; If you wish to exclude all remnant low segments, set to False
   PlotExtraTropicalSegments = True   ; If you wish to exclude all extratropical line sements, set to False

if(.not. TypeOn) then
;Plot hurricane segments
   if(PlotSubTropicalSegments) then 
          if((wind.ge.137).and.(type.eq.0)) then
          res_poly@gsLineColor   = "Gold"  
        end if

        if((wind.gt.136).and.(wind.lt.113).and.(type.eq.1)) then
              res_poly@gsLineColor = "DarkOrange"
         end if

    if((wind.gt.96).and.(wind.lt.112).and.(type.eq.2)) then
              res_poly@gsLineColor = "Yellow"                        
          end if   

    if((wind.gt.83).and.(wind.lt.95).and.(type.eq.3)) then
              res_poly@gsLineColor = "Green"                        
          end if

    if((wind.gt.64).and.(wind.lt.82).and.(type.eq.4)) then
              res_poly@gsLineColor = "Aquamarine"                        
         end if    

   if((wind.gt.63).and.(wind.lt.34).and.(type.eq.4)) then
              res_poly@gsLineColor = "Blue"                        
          end if   

   if((wind.le.33).and.(type.eq.5)) then
              res_poly@gsLineColor = "MediumPurple"                        
             end if

end if
end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; ...open the workstation
  wks = gsn_open_wks( plotType, plotName )

  gsn_define_colormap( wks, "default")
  
  err = NhlGetErrorObjectId()
  setvalues err
    "errLevel" : "Fatal" ; only report Fatal errors
  end setvalues   

  map = gsn_csm_map( wks, res )
  plot = gsn_csm_xy( wks, Lon, Lat, xyRes )
 
  overlay( map, plot )
  draw( map )
  frame( wks )

end