; Concepts illustrated:
;   - Overlaying vectors and filled contours on a map
;   - Masking out particular areas in a map
;   - Subsetting a color map
;----------------------------------------------------------------------
; This script shows how to overlay contours and vectors on a map,
; but with the contours limited to specific areas, and the 
; vectors not limited.
;
; The point of this script is to show how to mask contours against a
;geographical boundary, but in a way that allows them to be drawn up to
; the boundary location. This is unlike the shapefile masking examples,
;where grid points are set to missing if they fall outside a boundary,
; and hence you can get blocky features close to the boundary.
;
; The uvt.nc data file can be downloaded from:
;
; http://www.ncl.ucar.edu/Applications/Data/
;
;; This script was written by Yang Zhao (CAMS) 
;; (Chinese Academy of Meteorological Sciences) 
;;  email: 409360946@qq.com    Thanks you!
;;----------------------------------------------------------------------

begin
;;read u,v,t from the data at 500hPa 
        f = addfile("uvt.nc","r")
        level = f->lev
        u = f->U(:,{500},:,:)
        v = f->V(:,{500},:,:)
        t = f->T(:,{500},:,:)

printVarSummary(u)
printVarSummary(v)
printMinMax(u,0)
printMinMax(v,0)

;;create plots;;
        wks = gsn_open_wks("png","Bangladesh") ; send graphics to PNG file
        cmap = read_colormap_file("gui_default")

        res                = True
        res@gsnDraw        = False
        res@gsnFrame       = False 
        res@gsnMaximize    = True 
        res@tmXTOn         = False
        res@tmYROn         = False
        res@sfXArray            = grid_lont        ; Required to tell NCL where to                                          
        res@sfYArray            = grid_latt        ; overlay data on globe.  
        res@gsnLeftString  = ""
        res@gsnRightString = ""

printMinMax(grid_lont,0)
printMinMax(grid_latt,0)

plot = gsn_csm_contour_map(wks,aqual,res)

mkres                = True
mkres@gsnCoordsLat   = grid_latt
mkres@gsnCoordsLon   = grid_lont
;mkres@gsMarkerSizeF = 5      ; may need to adjust this if markers are too large or too small
gsn_coordinates(wks,plot,aqual,mkres)
        
;;set map;;
        mpres                             = res
        mpres@mpDataSetName               = "Earth..4"       ;http://www.ncl.ucar.edu/Document/HLUs/Classes/MapPlotData4_1_earth_4.shtml
        mpres@mpDataBaseVersion           = "MediumRes"
        mpres@mpOutlineOn                 = True
        mpres@mpOutlineSpecifiers         = (/"Eurasia","Bangladesh"/)
        mpres@mpGeophysicalLineThicknessF = 2
        mpres@mpNationalLineThicknessF    = 2
        mpres@mpFillDrawOrder             = "PostDraw"
        mpres@mpFillOn                    = True
        mpres@mpFillAreaSpecifiers        = (/"water",       "land" /)
        mpres@mpSpecifiedFillColors       = (/"deepskyblue2","white"/)
        mpres@mpMaskAreaSpecifiers        = (/"Eurasia","Bangladesh"/)

;;set area;;
        mpres@mpMinLatF                   = 20
        mpres@mpMaxLatF                   = 27 
        mpres@mpMinLonF                   = 88
        mpres@mpMaxLonF                   = 93

;;set contour;;
        cnres                             = res
        cnres@cnFillDrawOrder             = "PreDraw"
        cnres@cnFillOn                    = True
        cnres@cnLinesOn                   = False
        cnres@cnFillPalette               = cmap
        cnres@lbOrientation        = "Vertical"
        
;;set vector;;
        res_vc                            = res
        res_vc@vcGlyphStyle               = "LineArrow"
        res_vc@vcLineArrowThicknessF      = 5
        res_vc@vcMinDistanceF             = 0.01
        res_vc@vcRefLengthF               = 0.03
        res_vc@vcRefAnnoOn               = True
        res_vc@vcRefMagnitudeF           = 40 
        res_vc@vcRefAnnoString1          = "40"
        res_vc@vcRefAnnoSide             = "Top"
        res_vc@vcRefAnnoString2On        = False
        res_vc@vcRefAnnoPerimOn          = False
        res_vc@vcRefAnnoOrthogonalPosF   = -0.12
        res_vc@vcRefAnnoParallelPosF     = 0.999
        res_vc@vcRefAnnoBackgroundColor  = "Purple"
        res_vc@vcVectorDrawOrder         = "PostDraw"
        res_vc@gsnRightString            = "Wind"
  
       
;;plot;;
        map     = gsn_csm_map(wks,mpres)
        contour = gsn_csm_contour(wks,t,cnres)
        vector  = gsn_csm_vector(wks,u,v,res_vc)

;;overlay filled contours and vectors on the map;;
        overlay(map,contour)
        overlay(map,vector)

;;drawing "map" will draw everything: map, contours, and vectors
        draw(map)
        frame(wks)
end