;***************************************************************************;
; Load libraries.
;***************************************************************************;
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

;-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
; This function creates a cylindrical equidistant map of the area of
; interest.
;-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
function create_SD_map(wks,res,draw_ncl_outlines)
local a, mpres
begin
mpres = res
mpres@gsnMaximize = True
mpres@gsnPaperOrientation = "portrait"
mpres@gsnDraw = False
mpres@gsnFrame = False
mpres@mpFillOn = False
;-- -Turn on fancier tickmark labels.
mpres@pmTickMarkDisplayMode = "Always"
mpres@tmXBLabelFontHeightF = 0.008 ; smaller tickmark labels
;-- -Zoom in on area of interest
mpres@mpLimitMode = "LatLon"
mpres@mpMaxLatF = 43.518145
mpres@mpMinLonF = -97.592181
mpres@mpMinLatF = 43.365088
mpres@mpMaxLonF = -97.067584
mpres@mpFillOn = False
if(draw_ncl_outlines) then
mpres@mpOutlineOn = True
mpres@mpOutlineBoundarySets = "AllBoundaries"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpDataSetName = "Earth..4" ; U.S. counties
else
mpres@mpOutlineOn = False
end if
;-- -Create map.
map = gsn_csm_map(wks,mpres)
return(map)
end
;-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -;
; MAIN CODE
;
; Basemap for tornado location plot on July 06 2015 in southeastern SD.
;-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -;
begin
;-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- ;
; Plot!
;-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- ;
type = "png"
type@wkWidth = 2000
type@wkHeight = 2000
wks = gsn_open_wks(type,"tornado")
res = True
res@mpFillOn = False
map_sd = create_SD_map(wks,res,False)

;-- Add USA county shapefile
cores = True
cores@gsLineColor = "NavyBlue"
cores@gsLineThicknessF = 2.0
co = gsn_add_shapefile_polylines(wks,map_sd,"USA_adm2.shp",cores)

;-- Add roads
cores@gsLineColor = "Red"
cores@gsLineThicknessF = 1.0
ro = gsn_add_shapefile_polylines(wks,map_sd,"roads.shp",cores)

draw(map_sd)
frame(wks)

end