; Concepts illustrated:
;   - Specifying geolocation_circle parameters
;   - Use geolocation_circle to  calculate circular polygons
;     at different location centers
;   - Using gc_inout to mask out data outside a great circle
;   - Extract a region that spans area of interest
;   - using area_poly_sphere to calculates the area enclosed by an arbitrary polygon on the sphere
;   - use gc_clkwise to verify that the lat and lon are in clockwise order
;   - Use NCL's distributed 1x1 land-sea mask top create a mask file 
;     CAVEAT: landsea_mask uses a 1x1 basemap. TRMM 3B42 is 0.25


   ;---Read in TRMM precipitation data
  
   f           = addfile("all_BoB.nc", "r")
   prc         = f->pre(:,{75:100},{0:30})
   prc_reorder = prc(time|:,latitude|:,longitude|:)
   ;printVarSummary(prc_reorder)
  
  
   value = 1d20
   prc_reorder@missing_value = value                  ; Force change to CF recognized value
   prc_reorder@_FillValue    = value
  
   
   ; TCs rainfall
   ntim = (/16,17,18,19,43,44,45,73,74,75,76,79,80,81,82,83/) ; 1998 TCs period
   tc98 = prc_reorder(ntim,:,:)
   tcavg98 = dim_avg_n_Wrap(tc98,0)
   ;printVarSummary(tcavg98)
 

;---Specify center locations 

 filename = "year_track_98.txt"
 ncols = numAsciiCol(filename)                     ; Calculate the number of columns
 data_ = readAsciiTable(filename,ncols,"float",0)  ; Given the # of columns, we can use readAsciiTable to read this file
 
 slat=data_(:,0)  ; Latitude:  storm centers
 slon=data_(:,1)  ; Longitude: storm centers 
 
  srad      = (/ 500  /)                   ; one or mor radii (great circle distance)
  srad_unit = 1                            ; 0=degrees; 1=km
  N         = 90                           ; # of points; more points nicer 'circle'
  opt_gc    = False

;---Derive circle polygons

  circle    = geolocation_circle(slat, slon, srad, srad_unit, N, opt_gc) ; circle is type list

  ;printVarSummary(circle)                  ; list containing two variables
  
  clat      = circle[0]                    ; For clarity: explicitly extract list elements
  clon      = circle[1]
                                           ; For *each* location and radius
                                           ; rightmost dimension contains polygon lat and lon
  ;printVarSummary(clat)                   ; [location | 15] x [radii | 1] x [circle | 90]
  ;print("------------------")
  ;printVarSummary(clon)                   ; [location | 15] x [radii | 1] x [circle | 90]
  ;print("------------------")


;---For each center and radius: extract extent of each circle [ reduce gc_inout time ]
;---Use NCL's 'reassignment syntax [ := ] to accomodate possiblr changing array sizes

  nCenter   = dimsizes(slat)               ; # of locations
  nRadii    = dimsizes(srad)               ; # of circles at each location
  cdim      = dimsizes(clat)
  nCenter   = cdim(0)                      ; # of centers: same as: nCenter=dimsizes(slat)
  nRadius   = cdim(1)                      ; # of circles: same as: nRadius=dimsizes(srad) 
  nPoints   = cdim(2)                      ; # of points : same as: nPoints=N
 
  do nr=0,nRadius-1
    do nc=0,nCenter-2
       poly_lat := clat(nc,nr,:)            ; clarity; place into separate variable  
       poly_lon := clon(nc,nr,:) 
       
       min_lat  = min(poly_lat)            ; min of current latitude polygon
       max_lat  = max(poly_lat)            ; max       "
       min_lon  = min(poly_lon)            ; min       "    longitude polygon
       max_lon  = max(poly_lon)

;---Extract the desired rectangle of data 

       PRC  := tcavg98({min_lat:max_lat},{min_lon:max_lon})

;---Set points that are outside of the circle of data to missing

       lat2d := conform(PRC, PRC&latitude, 0)          
       lon2d := conform(PRC, PRC&longitude, 1)         
	   ;printVarSummary(lat2d)
	 
       latlon_circle := gc_inout(lat2d,lon2d, poly_lat,poly_lon)
       PRC    = where(latlon_circle,  PRC, PRC@_FillValue)	   
       ;printVarSummary(PRC)
	   
    if (nc.eq.0) then
	 
	                 f_lsm    = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")  

                     lsmtc   = tobyte(landsea_mask(f_lsm->LSMASK,PRC&latitude,PRC&longitude))              
                     lsmtc@long_name  = "Land-Sea Mask: TRMM 3B42"
                     copy_VarCoords(PRC, lsmtc)
	
	                 tc_land_MASK = where(lsmtc.eq.1 .or. lsmtc.eq.2, PRC, PRC@_FillValue) ; Land only
                     copy_VarMeta(PRC, lsmtc)
                     print("tc_landgrid = "+tc_land_MASK)
	
    else					 
                    
		             f_lsm    = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")  

                     lsmtc1   = tobyte(landsea_mask(f_lsm->LSMASK,PRC&latitude,PRC&longitude))              
                     lsmtc1@long_name  = "Land-Sea Mask: TRMM 3B42"
                     copy_VarCoords(PRC, lsmtc1)
	
	                 tc1_land_MASK = where(lsmtc1.eq.1 .or. lsmtc1.eq.2, PRC, PRC@_FillValue) ; Land only
                     copy_VarMeta(PRC, lsmtc1)
                     print("tc1_landgrid = "+tc1_land_MASK)
	 end if

    end do   ; nc
end do      ; nr