;*************************************************
; polyg_9.ncl
;
; Concepts illustrated:
;   - Drawing the meteorological subdivisions of India
;   - Paneling two plots on a page
;   - Attaching polylines to a map plot
;   - Attaching filled polygons to a map plot
;   - Changing the color and thickness of polylines
;   - Zooming in on India on a cylindrical equidistant map
;
;*************************************************

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

begin
     wks  = gsn_open_wks("ps","polyg")            ; open a ps file
     gsn_define_colormap(wks,"amwg")

     res                             = True
     res@mpDataBaseVersion           = "Ncarg4_1"	  ; higher resolution basemap
     res@mpFillOn                    = False        ; turn off gray fill
     res@mpGeophysicalLineColor      = "Black"      ; color of cont. outlines
     res@mpGeophysicalLineThicknessF = 1.5          ; thickness of outlines
     res@mpOutlineOn = True

     res@mpMaxLatF                   = 38           ; choose subregion to plot          
     res@mpMinLatF                   = 6
     res@mpMaxLonF                   = 100
     res@mpMinLonF                   = 65
     res@gsnFrame = False
     res@gsnDraw = False
     plot = new(2,graphic)
     plot(0) = gsn_csm_map(wks,res)                 ; draw blank map for met division outlines
     plot(1) = gsn_csm_map(wks,res)  			  ; draw blank map for met division filling
	
     resp                  = True                      ; polyline mods desired
     resp@gsLineColor      = "red"                     ; color of lines
     resp@gsLineThicknessF = 2.0                       ; thickness of lines
;===========================================================================================
;     Read in ascii files. Each ascii file contains the lat/lon coordinates for 1 
;     Meteorological Subdivision. The ascii files are of different sizes, and are formatted:
;     like this:
;
;     36 draw polyf 88.8499 27.2099 89 27.36 88.9499 27.4099 88.9499 27.51 88.9 27.6599 
;
;     where the first number is the station number (36 in above snippet), the second
;     number is the longitude of the first data point (88.8499), the third number is the
;     latitude of the first data point (27.2099), the fourth number is the longitude of 
;     the second data point (89.), and so on...
;===========================================================================================
     dum = new(2909,"graphic")	; need to create a dummy array to house each polyline,
				; in this case, there will be 2909 lines drawn. The
				; number 2909 was calculated previously. 
     dum2 = dum				; for 2nd panel of panel plot											
     pgres = True
     colorarr = random_uniform(2,14,34)		; set up some random colors to fill the met divisions with.
									; for second panel
									
     cntr = 0					; polyline counter	
     ccntr= 0					; color counter
     nfils = systemfunc("ls ../../DrawIndMetSub/india*.asc")
     do gg = 0,dimsizes(nfils)-1						; loop thru all ascii files
        arrT = asciiread(nfils(gg),(/-1/),"float")		; read in an ascii file
        numpts = dimsizes(arrT(1::2))					; number of data points
        print(numpts+" points for station "+arrT(0))
        finarr = new((/2,numpts/),"float")				; create new array to house coordinates
        finarr(0,:) = (/ arrT(2::2) /)					; put latitudes in finarr(0,:)
        finarr(1,:) = (/ arrT(1::2) /)					; put longitudes in finarr(1,:)
        delete(arrT)
        delete(numpts)
      
        do hh = 0,dimsizes(finarr(0,:))-2		; draw polylines for subdivision borders for 1st panel
           dum(cntr) = gsn_add_polyline(wks,plot(0),finarr(1,hh:hh+1),finarr(0,hh:hh+1),resp)
     	   cntr = cntr+1
        end do	
	   
	pgres@gsFillColor = colorarr(ccntr)		; draw filled polygons for 2nd panel
	dum2(ccntr) = gsn_add_polygon(wks,plot(1),finarr(1,:),finarr(0,:),pgres)
	ccntr = ccntr+1
	delete(finarr)
     end do							
     gsn_panel(wks,plot,(/2,1/),False)
end