;----------------------------------------------------------------------
; shapefiles_8.ncl
;
; Concepts illustrated:
;   - Reading shapefiles
;   - Plotting data from shapefiles
;   - Using shapefile data to draw areas of interest in India
;   - Zooming in on India on a cylindrical equidistant map
;----------------------------------------------------------------------
; This example shows how to read geographic data from a shapefile
; and plot it on a map created by NCL.
;
; This particular example plots data for India
;----------------------------------------------------------------------
; The shapefiles for this example were obtained from the
; "Global Administratie Areas" website:
;
;      http://www.gadm.org/country
;----------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

;----------------------------------------------------------------------
; This procedure attaches polyline data from a shapefile
;----------------------------------------------------------------------
undef("attach_shapefile_polyline")
function attach_shapefile_polyline(f:file,wks,plot,line_color)
local lnres
begin
  gsn_define_colormap(wks,(/"white","black","tan","LightBlue",\
                            "brown","yellow","navyblue","green"/))

  lnres                  = True      ; resources for polylines
  lnres@gsLineThicknessF = 2.0       ; 2x as thick
  lnres@gsLineColor      = line_color
;
; Loop through files that we want to read geographic information from.
;
; If this loop is extremely slow, consider using gsn_polyline instead
; of gsn_add_polyline. This can have a significant effect. Remember
; that gsn_polyline is a procedure, not a function, and it draws the
; lines right when you call it, so you need to make sure your plot is
; already drawn to the frame.
;
;---Read data off the shapefile
  segments = f->segments
  geometry = f->geometry
  segsDims = dimsizes(segments)
  geomDims = dimsizes(geometry)

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

;---Section to draw polylines on plot.
  lon = f->x
  lat = f->y
  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
;
; This call adds the line segment.
;
; Can use gsn_polyline here to make it faster.
;
        dumstr = unique_string("primitive")
        plot@$dumstr$ = gsn_add_polyline(wks, plot, lon(startPT:endPT),  \
                                        lat(startPT:endPT), lnres)
     end do
  end do

;---Clean up before we read in same variables again.
  delete(lat)
  delete(lon)
  delete(segments)
  delete(geometry)

;---We have to return plot so that attached lines are not lost.
  return(plot)
end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
  dir       = "IND_adm/"
  filenames = "IND_adm" + ispan(0,3,1) + ".shp"
  nfiles    = dimsizes(filenames)

;--- Open workstation.
  wks = gsn_open_wks("ps","shapefiles")

  res               = True
  res@gsnMaximize   = True
  res@gsnDraw       = False
  res@gsnFrame      = False

  res@mpOutlineOn   = False   ; Use outlines from shapefile
  res@mpFillOn      = False  

;---Turn on fancier tickmark labels.
  res@pmTickMarkDisplayMode = "Always"

;---Zoom in on area of interest
  res@mpLimitMode           = "LatLon"
  res@mpMinLatF             =  5
  res@mpMaxLatF             = 37
  res@mpMinLonF             = 65
  res@mpMaxLonF             = 99

;
; Loop through the four given shapefiles, and draw the data from
; the first two on one map, and the other two on their own maps.
;

;---Open first two shapefiles.
  f1 = addfile(dir + filenames(0),"r")
  f2 = addfile(dir + filenames(1),"r")

  res@tiMainString = filenames(0) + "/" + filenames(1)

;---Create a map with a title.
  map = gsn_csm_map(wks,res)

;---Attach polylines to map.
  map = attach_shapefile_polyline(f1,wks,map,"black")
  map = attach_shapefile_polyline(f2,wks,map,"black")

  draw(map)
  frame(wks)

;---Open third shapefile
  f3 = addfile(dir + filenames(2),"r")

;---Create a new map with a different title.
  res@tiMainString = filenames(2)
  map = gsn_csm_map(wks,res)

;---Attach polylines to map.
  map = attach_shapefile_polyline(f3,wks,map,"blue")

  draw(map)
  frame(wks)

;---Open fourth shapefile
  f4 = addfile(dir + filenames(3),"r")

;---Create a new map with a different title.
  res@tiMainString = filenames(3)
  map = gsn_csm_map(wks,res)

;---Attach polylines to map.
  map = attach_shapefile_polyline(f4,wks,map,"brown")

  draw(map)
  frame(wks)

end