;----------------------------------------------------------------------
;
; 2017-mar-31	shapefiles_15.ncl: Original example from NCL website.
; 2017-may-30	Modified by Dave Allured.
;		Overlay NCL state outlines and shapefile on single plot.
;		Add point marker.
;
; Concepts illustrated:
;   - Comparing NCL state and county outlines to shapefile map outlines
;   - Zooming in on a cylindrical equidistant map
;   - Using functions for cleaner code
;
; The USA_adm files were downloaded from gadm.org/country.
;
;----------------------------------------------------------------------

; 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"

;---------------------------------------------------------------------------
; This function creates a cylindrical equidistant map in a specified area.
;---------------------------------------------------------------------------

function create_colorado_map(wks,res,draw_ncl_outlines)
local mpres, map
begin
  mpres                       = res

  mpres@gsnMaximize           = True
  mpres@gsnPaperOrientation   = "portrait"

  mpres@gsnDraw               = False
  mpres@gsnFrame              = False

  mpres@mpFillOn              = False

  mpres@tmXBMode              = "Automatic"
  mpres@tmYLMode              = "Automatic"
  mpres@tmXBLabelFontHeightF  = 0.012      ; smaller tickmark labels

;---Zoom in on area of interest

  mpres@mpMinLatF             = 36.4
  mpres@mpMaxLatF             = 36.8
  mpres@mpMinLonF             = -79.75
  mpres@mpMaxLonF             = -79.15

  mpres@mpLimitMode           = "LatLon"
  mpres@mpFillOn              = False

; Add NCL built-in state and county outlines.

  mpres@mpOutlineOn           = True
  mpres@mpOutlineBoundarySets = "AllBoundaries"
  mpres@mpDataBaseVersion     = "MediumRes"
  mpres@mpDataSetName         = "Earth..4"      ; U.S. counties

;---Create map.

  map = gsn_csm_map(wks,mpres)

  return(map)
end

;--------------------------------------------------
; Main code
;--------------------------------------------------

begin
  wtype = "pdf"
  shapefile = "$DD/shapefiles/us/gadm/USA_adm2.shp"

; Set reference point.  Use double precision for accurate location.
; (For multiple points, could also be arrays of lats and lons.)

  lat =  36.5459089d		; southwest of Danville, Virginia
  lon = -79.4585847d

  wks = gsn_open_wks (wtype, "shapefiles_15.danville")

  res                   = True
  res@tiMainFontHeightF = 0.015

  res@tiMainString    = "State and county outlines; GADM shapefile + NCL " \
                         + get_ncl_version()
  res@gsnCenterString =  "Near Danville, Virginia"

  map = create_colorado_map (wks, res, True)

  lnres               = True
  lnres@gsLineColor   = "black"

  id = gsn_add_shapefile_polylines (wks, map, shapefile, lnres)

; Add location dot(s) at the specified point location(s).

  mkres               = True
  mkres@gsMarkerIndex = 17     ; Filled circle
  mkres@gsMarkerSizeF = 0.03

  map@markers = gsn_add_polymarker (wks, map, lon, lat, mkres)

  draw (map)
  frame (wks)
end
exit