;----------------------------------------------------------------------
; shapefiles_11.ncl
;
; Concepts illustrated:
;   - Reading shapefiles
;   - Creating a mask variable using outlines from a shapefile
;   - Attaching markers to a map
;   - Attaching polylines to a map plot
;   - Masking a data array based on a geographical area obtained from a shapefile
;----------------------------------------------------------------------
; This script creates a new data mask based on the outline of the
; United States.  It then draws two plots: the original data,
; and the data with the USA mask.
;
; The "gz_2010_us_020_00_5m.shp" shapefile is from
; http://www.census.gov/geo/www/cob/
;
; The "nationalp010g.shp" shapefile is from:
; http://www.nationalatlas.gov/atlasftp-1m.html#nationp
;
; Try both of these files to see if they work for your purposes.
; The national file is larger, and hence takes longer to process
; (206 wall clock seconds on a Mac versus 17 wall clock seconds.)
;
; Once you create the mask, you can set CREATE_MASK to False so
; it doesn't get created again when you run this script.
;
;----------------------------------------------------------------------
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 "shapefile_mask_data.ncl"

;----------------------------------------------------------------------
; This is the main code
;----------------------------------------------------------------------
begin
;---Name of shapefile containing USA outlines
   shp_fname = "GSAs_simplified.shp"
  
;---Rough area we are interested in. Everything outside this will be masked.
   minlon =  -08.0
   maxlon =  45.0
   minlat =  29.00
   maxlat =  48.0

;---Read in zonal winds
   a = addfile("dox_med_1.4-101.nc","r")
   dox = a->dox
   ;printVarSummary(dox)
   dox_dims = dimsizes(dox)
   ntim = dox_dims(0)
   nlev = dox_dims(1)
   nlat = dox_dims(2)
   nlon = dox_dims(3)

;---Create a new mask using a shapefile of USA
   opt             = True
   opt@return_mask = True
;  opt@debug       = True
   opt@minlon      = minlon     ; Makes the shapefile masking
   opt@maxlon      = maxlon     ; go faster.
   opt@minlat      = minlat
   opt@maxlat      = maxlat
   opt@shape_var   = "SECT_COD"
   opt@shape_names = ("GSA01")
   pmask2d  = shapefile_mask_data(dox(0,0,:,:),shp_fname,opt)

;---Mask original 4D variable against a 4D pmask array
   pmask4d    = conform_dims(dox_dims,pmask2d,(/2,3/))
   dox_masked = mask(dox,pmask4d,1)   

   copy_VarMeta(dox,dox_masked)       ; copy attributes and coordinate arrays
   dox_masked!1 = "lev"               ; rename three of the dimensions
   dox_masked!2 = "lat"
   dox_masked!3 = "lon"

   printVarSummary(dox)
   printVarSummary(dox_masked)
   printMinMax(dox,0)
   printMinMax(dox_masked,0)

   wks = gsn_open_wks("eps","plot_dox01_0-209")

   cmap = read_colormap_file("cmp_b2r")
  
   res                       = True     ; plot mods desired
   res@gsnMaximize           = True    
   res@gsnDraw               = False
   res@gsnFrame              = False


   res@cnFillOn              = True     ; turn on color fill
   res@cnFillPalette         = cmap(::-1,:)   ; reverse color map
   res@cnLinesOn             = False    ; turn off contour lines
   res@cnLineLabelsOn        = False    ; turn off contour line labels
   res@lbLabelBarOn          = False    ; will turn on in panel
   res@cnLevelSelectionMode  = "ManualLevels"	; manually set the contour levels with the following 3 resources
   res@cnMinLevelValF        = 0.			; set the minimum contour level
   res@cnMaxLevelValF        = 300.			; set the maximum contour level
   res@cnLevelSpacingF       = 10	

   res@mpOutlineOn           = False
   res@mpDataBaseVersion     = "MediumRes"
   res@mpMinLonF             = minlon
   res@mpMaxLonF             = maxlon
   res@mpMinLatF             = minlat
   res@mpMaxLatF             = maxlat
   res@mpCenterLonF          = (minlon+maxlon)/2.
   res@pmTickMarkDisplayMode = "Always"           ; nicer tickmarks
   res@pmTitleZone           = 4                  ; moves title down

   res@gsnAddCyclic          =  False             ; Regional data

;---Be sure to use same levels for both plots
   res@gsnRightString = "mmol/m3"
   res@gsnLeftString  = "Dissolved_Molecular_Oxygen"

   lnres = True
   pres = True
   pres@gsnMaximize       = True    
   pres@gsnPanelLabelBar  = True

;---Comment out the plotting for now because it's time consuming.
;;   do nt=0,dox_dims(0)-1
;;     res@tiMainString = "Original data, nt = " + nt
;;     plot_orig = gsn_csm_contour_map(wks,dox(nt,0,:,:),res)
;;     res@tiMainString = "Masked data, nt = " + nt
;;     plot_mask = gsn_csm_contour_map(wks,dox_masked(nt,0,:,:),res)
;;     id_orig = gsn_add_shapefile_polylines(wks,plot_orig,shp_fname,lnres)
;;     id_mask = gsn_add_shapefile_polylines(wks,plot_mask,shp_fname,lnres)
;;
;;     gsn_panel(wks,(/plot_orig,plot_mask/),(/2,1/),pres)
;;   end do

;---Open file for writing
   diro = "./"                     ; Output directory
   filo = "dox_med_1.4-101.nc"
   system("/bin/rm -f " + diro + filo)    ; remove if exists
   foutput  = addfile (diro + filo, "c")  ; open output file
   setfileoption(foutput,"DefineMode",True)

;---Write global file attributes
   fAtt               = True            ; assign file attributes
   fAtt@title         = "NCL Efficient Approach to netCDF Creation"
   fAtt@source_file   =  "dox_med_1.4-101.nc"
   fAtt@Conventions   = "None"
   fAtt@creation_date = systemfunc ("date")

   fileattdef( foutput, fAtt )            ; copy file attributes

;---Write file dimension names and sizes
   dimNames = (/"time","lev", "lat", "lon"/)
   dimSizes = (/ -1   ,  nlev,  nlat,  nlon /)
   dimUnlim = (/ True , False, False, False/)
   filedimdef(foutput,dimNames,dimSizes,dimUnlim)

   filevardef(foutput, "time" ,typeof(dox_masked&time),getvardims(dox_masked&time))
   filevardef(foutput, "lev"  ,typeof(dox_masked&lev),getvardims(dox_masked&lev))
   filevardef(foutput, "lat"  ,typeof(dox_masked&lat),getvardims(dox_masked&lat))
   filevardef(foutput, "lon"  ,typeof(dox_masked&lon),getvardims(dox_masked&lon))
   filevardef(foutput, "dox",typeof(dox_masked),getvardims(dox_masked))
                                                              ; different from name on script
;---Write variable attributes
   filevarattdef(foutput,"dox",dox_masked)                           ; copy T attributes
   filevarattdef(foutput,"time" ,dox_masked&time)         ; copy time attributes
   filevarattdef(foutput,"lev"  ,dox_masked&lev)          ; copy lev attributes
   filevarattdef(foutput,"lat"  ,dox_masked&lat)          ; copy lat attributes
   filevarattdef(foutput,"lon"  ,dox_masked&lon)          ; copy lon attributes

   setfileoption(foutput,"DefineMode",False)

;---Write variable values
   foutput->time   = (/dox_masked&time/)
   foutput->lev    = (/dox_masked&lev/)
   foutput->lat    = (/dox_masked&lat/)
   foutput->lon    = (/dox_masked&lon/)
   foutput->dox = (/dox_masked/)
delete(dox)
end