;https://www.ncl.ucar.edu/Applications/HDF.shtml
;	  https://www.ncl.ucar.edu/Applications/Scripts/smap_l3_1.ncl
;*************************************************************** 
;----------------------------------------------------------------------
  ;load "BlAqGrYeOrReVi200_modified.ncl"
;----------------------------------------------------------------------
; Iran Provinces shapefile
undef("add_outlines_to_map")
procedure add_outlines_to_map(wks,plot_smap)
local fnames, colors, lnres, n
begin
;
; These two shapefiles contain waterways and administrative 
; geographic info. We will add them to the existing map using
; polylines.
;
;  fnames = (/"TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3","../2shapefile/iran_province"/) + ".shp"
  fnames = (/"/home/taghizade/96_8_26/shp_ostan/iran_province"/) + ".shp"

;  colors = (/"LightBlue", "Tan"/)    ; water, administrative
  colors = (/"Tan"/)    ; water, administrative

  lnres                  = True      ; resources for polylines
  lnres@gsLineThicknessF = 2.0       ; 2x as thick

; Ehsan: Following loop is unnecessary
  do n=0,dimsizes(fnames)-1
    dumstr = unique_string("primitive")
    lnres@gsLineColor = colors(n)
    plot_smap@$dumstr$ = gsn_add_shapefile_polylines(wks, plot_smap, fnames(n), lnres)
  end do
end

;----------------------------------------------------------------------
; Add markers for cities to given map plot.
;
; 22 synoptic stations in NW
undef("add_places_to_map")
procedure add_places_to_map(wks,plot)
local f, names, lat, lon, txres, mkres
begin
;---Read shapefile with "places" information.
  f = addfile("/home/taghizade/96_8_26/shp_stations_NW_22/stations_NW_22.shp","r")

  names      = f->icao_code;name;station_id;
  lon        = f->x
  lat        = f->y
  num_places = dimsizes(names)

;---Set up text resources to label random areas of interest.
  txres               = True
  txres@txFontHeightF = 0.015
  txres@txJust        = "topright"

;---Set up marker resources to mark random areas of interest.
  mkres               = True
  mkres@gsMarkerIndex = 16
  mkres@gsMarkerColor = "brown"
  mkres@gsMarkerSizeF = 7.

;---Areas we want to put a label and a marker.
  names_of_interest = (/"ARAL","OITU","AZUH","AZTA","OITT",\
                        "ARAD","AZTS","OITR","OITM","AZTY",\
                        "AZTY","OIGG","OIGG","Zürich","OIGG",\
                        "ZAZN","OIIK","OINR","OINN", \
						"OICS","OIII","OICC","OIHH"/)

;
; Loop through each of the "name" places in the shapefile, and see if
; it's on our list of ones we want to put a label and marker for.
;
; Also, make sure we haven't already encountered this place. The
; shapefile seems to contain duplicate entries (Zurich has two entries,
; for example).
;
; If you find this looping code to be too slow, you can use 
; gsn_polymarker and gsn_text instead.
;
  names_of_interest@_FillValue = "missing"
  do i=0,num_places-1
    if(any(names(i).eq.names_of_interest)) then
      ii = ind(names_of_interest.eq.names(i))
      names_of_interest(ii) = "missing"
      dumstr = unique_string("primitive")

;---Attach the marker to the map.
      plot@$dumstr$ = gsn_add_polymarker(wks, plot, lon(i), lat(i), mkres)
      dumstr = unique_string("primitive")
;

;---Attach the text string to the map. 
      plot@$dumstr$ = gsn_add_text(wks, plot, "  " + names(i), \
                                    lon(i), lat(i), txres)
      delete(ii)
    end if
  end do
end

;----------------------------------------------------------------------
; Add markers for Azarshahr and Ajabshir to given map plot.
undef("add_AzarAjab_to_map")
procedure add_AzarAjab_to_map(wks,plot)
local f, names, lat, lon, txres, mkres
begin
;---Read shapefile with "places" information.
  f = addfile("/home/taghizade/96_8_26/shp_Azarshahr_Ajabshir/Azarshahr_Ajabshir.shp","r")

  names      = f->name
  lon        = f->x
  lat        = f->y
  num_places = dimsizes(names)

;---Set up text resources to label random areas of interest.
  txres               = False
  ;txres@txFontHeightF = 0.015
  ;txres@txJust        = "topright"

;---Set up marker resources to mark random areas of interest.
  mkres               = True
  mkres@gsMarkerIndex = 16
  mkres@gsMarkerColor = "black"
  mkres@gsMarkerSizeF = 7.

  do i=0,num_places-1
      dumstr = unique_string("primitive")

;---Attach the marker to the map.
      plot@$dumstr$ = gsn_add_polymarker(wks, plot, lon(i), lat(i), mkres)
      dumstr = unique_string("primitive")
;
;---Attach the text string to the map. 
      plot@$dumstr$ = gsn_add_text(wks, plot, " ", \
                                    lon(i), lat(i), txres)
  end do
end

;------------------------------------------------------------------------
;
; This procedure adds a line to a plot, making sure that it is returned
; to a unique variable name, and that this variable is retained even
; outside this procedure call.
;
procedure add_polyline(wks,plot,xval,yval)
local plres,dum,i, str;, y
begin
  ;plres = True
  ;plres@gsLineColor = line_color   ; Set the line color.
  plres                  = True                      ; polyline mods desired
  plres@gsLineColor      = "red"                     ; color of lines
  plres@gsLineThicknessF = 4.0                       ; thickness of lines
  ;plres@gsLineLabelString= "test"                    ; adds a line label string

; create array of dummy graphic variables. This is required, b/c each line
; must be associated with a unique dummy variable.

  dum = new(4,graphic)  

; draw each line separately. Each line must contain two points.
  do i = 0 , 3
    dum(i)=gsn_add_polyline(wks,plot,xval(i:i+1),yval(i:i+1),plres)      
  end do 

  ;str1 = unique_string("polyline")   ; "unique_string" will return a unique
                                    ; string every time it is called from
                                    ;  within a single NCL session.
  str2 = unique_string("dum")
;
; You can then use this unique string as an attribute variable name
; that gets attached to the plot variable. This ensures that this
; value will live for the duration of the script.
;
  ;plot@$str1$ = gsn_add_polyline(wks, plot, xval, yval, plres)
;https://www.ncl.ucar.edu/Support/talk_archives/2012/0899.html
  plot@$str2$ = dum
end
;-----------------------------------------------------------------------------------
; smap_l3_1.ncl 
; 
; Concepts illustrated: 
;   - Reading a SMAP HDF5 level 3 file with groups
;   - Use 'direct' syntax to access variable within groups
;   - Manually adding _FillValue to latitude and longitude
;   - Plot
;***************************************************************
; These library 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" 
;**************************************************************
;
;;===============================================================
; SMAP values are provided on the global cylindrical EASE-Grid 2.0.
; Each grid cell has a nominal area of approximately 36 x 36 km2 
; regardless of longitude and latitude. Using this projection, all 
; global data arrays have dimensions of 406 rows and 964 columns. 
;;===============================================================
begin
;---Read specified h5 file
  diri = "/home/taghizade/96_8_26/SPL2SMP_E/"
  fili = "SMAP_L2_SM_P_E_11986_D_20170430T030839_R15060_001.h5"
  pthi = diri + fili
  f    = addfile(pthi, "r")
    
;--Set group (begin and end with /); and desired variable
  grp_smrd  = "/Soil_Moisture_Retrieval_Data/"
  varName   = "soil_moisture"
  var_path  = grp_smrd + varName
  sm        = f->$var_path$
  printVarSummary(sm)
  printMinMax(sm, 0)

  latName   = "latitude"
  lat_path  = grp_smrd + latName
  lat2d     = f->$lat_path$
 ;printVarSummary(lat2d)
 ;printMinMax(lat2d, 0)

  lonName   = "longitude"
  lon_path  = grp_smrd + lonName
  lon2d     = f->$lon_path$
 ;printVarSummary(lon2d)
 ;printMinMax(lon2d, 0)

;---Manually add a _FillValue to lat/lon; not sure why 
;---_FillValue not associated with the variable
  lat2d@_FillValue = -9999.0
  lon2d@_FillValue = -9999.0
 ;printMinMax(lat2d, 0)
 ;printMinMax(lon2d, 0)

;---Sample plot options
  pltDir   = "./"
  pltType  = "png"
  pltName  = "smap_l3_1_2_SMAP_L2_SM_P_E_11986_D_20170430T030839_R15060_001"
  pltTitle =  fili

;---
  pltPath = pltDir+pltName
  wks = gsn_open_wks(pltType,pltPath)

  res                  = True              ; Plot modes desired.
  res@gsnMaximize      = True              ; Maximize plot
  res@gsnAddCyclic     = False

  res@cnFillOn         = True              ; color plot desired
  res@cnLinesOn        = False             ; turn off contour lines
  res@cnLineLabelsOn   = False             ; turn off contour labels
  res@cnFillMode       = "RasterFill"      ; turn raster on      
  res@cnFillPalette    = "BlAqGrYeOrReVi200"

  if (varName.eq."soil_moisture") then
      res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
      res@cnMinLevelValF       = 0.05               ; set min contour level
      res@cnMaxLevelValF       = 0.95               ; set max contour level
      res@cnLevelSpacingF      = 0.05               ; set contour spacing
  end if

;START------------------------------------------------------------------------------
       res@gsnDraw               = False    ; don't draw
       res@gsnFrame              = False    ; don't advance frame
       res@mpDataSetName         = "Earth..4"   ; 1 through 4
	   res@mpOutlineBoundarySets = "National"
       res@mpFillOn              = False
       ;res@cnFillOn              = True     ; turn on color fill
       ;res@cnLinesOn             = False       ; Turn off contour lines
       ;res@cnLineLabelsOn        = False       ; Turn off contour lines
	   ;res@cnRasterSmoothingOn  = True
       res@mpMinLatF             = 33.3;min(prc&lat)-1;min(prc&latitude)
       res@mpMaxLatF             = 41.3;max(prc&lat)+1;max(prc&latitude)
       res@mpMinLonF             = 40.0;min(prc&lon)-1;min(prc&longitude)
       res@mpMaxLonF             = 52.5;max(prc&lon)+1;max(prc&longitude)
       res@mpCenterLonF          = 0.5*(41.0 + 51.5)
;END--------------------------------------------------------------------------------

;---Resources for plotting original (source) data
  res@sfXArray         = lon2d
  res@sfYArray         = lat2d
  res@trGridType       = "TriangularMesh"  

  res@tiMainString     = fili
  res@gsnLeftString    = "SMAP: Soil Moisture"    ; default long_name is too long
  plot_smap = gsn_csm_contour_map(wks,sm,res)

;START------------------------------------------------------------------------------
       add_outlines_to_map(wks,plot_smap)
       add_places_to_map(wks,plot_smap)
       add_AzarAjab_to_map(wks,plot_smap)
       figure_strings = "figure_strings";filesc(24:29)
       resP = True
       resP@gsnMaximize         = True                ; make ps/eps/pdf large [no effect x11]
	   resP@gsnPanelFigureStrings = figure_strings
	   resP@gsnPanelFigureStringsFontHeightF = 0.01
	   resP@wkPSResolution = 3000
       ;resP@gsnPaperOrientation = "Portrait"          ; force portrait
       ;*resP@gsnPanelLabelBar    = True                ; add common colorbar
       ;*resP@lbLabelAutoStride   = True
       ;*resP@lbLabelFontHeightF  = 0.0150              ; change font size
       resP@txString        = "GPM Level 3 IMERG Late Daily 0.1 x 0.1 Degree Precipitation";prc@long_name;"APHRODITE: RU: "+ymdPlot
  ;gsn_panel(wks,plot_smap,(/1,1/),resP)               ; now draw as one plot
  
;************************************************
; create points for box
;************************************************
; add the box
;
       ypts = (/ 36.0, 39.0, 39.0, 36.0, 36.0 /)
       xpts = (/ 44.0, 44.0, 48.0, 48.0, 44.0 /)
;
       add_polyline(wks,plot_smap,xpts,ypts)
;************************************************
; label the box with additional text
;************************************************
  ;tres       =  True
  ;tres@txFontHeightF = 0.02
  ;tres@txBackgroundFillColor = 0
  ;tres@txFontOpacityF = 2.0
  ;tres@txString = "sample"
  ;gsn_text(wks,plot_smap,"sample",45.51,37.3,tres)

  draw(plot_smap)
  frame(wks)
end
;END--------------------------------------------------------------------------------