; 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


;---Read in TRMM precipitation data
;1998
 f   = addfile("all_BoB.nc", "r")
 prc = f->pre(:,{75:100},{0:30})
 prc_reorder = prc(time|:,latitude|:,longitude|:)

  
  value = 1d20
  prc_reorder@missing_value = value                  ; Force change to CF recognized value
  prc_reorder@_FillValue    = value
  
 ; 1998 
 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-1
       poly_lat = clat(nc,nr,:)            ; clarity; place into separate variable  
       poly_lon = clon(nc,nr,:) 
       ;print(poly_lat)
	 	   

       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)          ; [36] x [38]
       lon2d := conform(PRC, PRC&longitude, 1)         ; [36] x [38]
	   ;printVarSummary(lat2d)
	 
       latlon_circle := gc_inout(lat2d,lon2d, poly_lat,poly_lon)
       PRC    = where(latlon_circle,  PRC, PRC@_FillValue)
       ;printVarSummary(PRC)
	   
    end do   ; nc
 end do      ; nr
    
	
	; Land Sea mask for TCs lat and lon
	
	f_lsm    = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")  

    lsmTCs   = tobyte(landsea_mask(f_lsm->LSMASK,lat2d,lon2d))              ;   - 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
    lsmTCs@long_name  = "Land-Sea Mask: TRMM 3B42"
    copy_VarCoords(lat2d, lsmTCs)
	;printVarSummary(lsmTCs)           ; [36] x [40]
	;print(lsmTCs)
	
	    dims   = dimsizes(lsmTCs)
	    PRClat = dims(0)
	    PRClon = dims(1)
	  
	TCs_lat_MASK = where(PRClat.eq.1 .or. PRClat.eq.2, PRClat, PRC@_FillValue) ; Land only
    copy_VarCoords(PRClat, TCs_lat_MASK)
	;print(TCs_lat_MASK)


    TCs_lon_MASK = where(PRClon.eq.1 .or. PRClon.eq.2, PRClon, PRC@_FillValue) ; Land only
    copy_VarCoords(PRClon, TCs_lon_MASK)
	;print(TCs_lon_MASK)
	 
	 
	; Calculate area
	rsph = 6371    ; Earth radius in Km   ; return area will be Km
    
	
	; TCs latitude and longitude
	;TCslat = ndtooned(lat2d)
    ;TCslon = ndtooned(lon2d)
	
	;TCslat1=TCslat(::-1)
	;TCslon1=TCslon(::-1)

    qorder     = gc_clkwise(TCs_lat_MASK, TCs_lon_MASK)            ; verify if lat and lon are in clockwise
    TCs_area = area_poly_sphere(TCs_lat_MASK, TCs_lon_MASK, rsph)
    print ("TCs_area="+TCs_area+"    qorder="+qorder)