;+
; NAME:
;	MONTHS_TO_SEASONS
;
; PURPOSE:
;	This function extracts seasonal data from monthly data.
;
; CATEGORY:
;	Calendar
;
; CALLING SEQUENCE:
;	Result = MONTHS_TO_SEASONS( Data, Season [, Seasonlen] )
;
; INPUTS:
;	Data:  A vector containing the input time series of monthly data.  
;		It can also be an array where the last dimension is time.  The 
;		time dimension must start with a January and end with a 
;		December.  Also see OUTPUT below.
;	Season:  The index (January=0,...) value of the middle month of the 
;		season.
;	Seasonlen:  Optional.  The number of months in a season.  The default 
;		value is "3".  If Seasonlen is even the extra month is added 
;		at the end.
;
; KEYWORD PARAMETERS:
;	ANOMALY:  If set, then the anomaly to the mean over BASEPERIOD is 
;		calculated.  If DESEASONALISE is set, then the anomaly from 
;		the mean annual cycle over BASEPERIOD is calculated.  The 
;		default is to calculate total values.
;	BASEPERIOD:  A vector of [YEARSTART,YEAREND] defining the period over 
;		which to estimate the seasonal cycle when the ANOMALY or 
;		DESEASONALISE keywords are set.  YEARSTART is the first year 
;		index value (month/12) while YEAREND is the last year index 
;		value.  e.g. using [10,19] would use the month indices 120 
;		through 239.
;	CYCLE:  If DESEASONALISE is set, then cycle contains a vector or array 
;		(with calendar month as the last dimension) containing the 
;		calculated seasonal cycle.  Note this is the total seasonal 
;		cycle if ANOMALY is set, or the seasonal cycle anomaly from 
;		the annual mean if ANOMALY is not set.
;	DESEASONALISE:  If set then the seasonal cycle is removed from the 
;		data before the seasonal values are calculated.  This has an 
;		effect if NGOOD is set, because it then removes the bias 
;		arising from missing data in certain months.  It can also 
;		effect the end values if the season overlaps years.
;	NBASEGOOD:  The number of years in the base period required to have 
;		good (non-NaN) values in order for the base value to be 
;		calculated.  The default is 1.
;	NGOOD:  The number of months in a season required to have good 
;		(non-NaN) values in order to calculate a seasonal value.  The 
;		default is NGOOD=Seasonlen.
;
; OUTPUT:
;	Result:  Returns the seasonal time series as a vector or an array, 
;		depending on the dimensions of Data.  The time dimension will 
;		be shorter by a factor of 12.
;	CYCLE:  See above.
;	Data:  If the DEASONALISE and/or ANOMALY options are set then Data is 
;		returned in a deseasonalised and/or anomaly form.
;
; USES:
;	var_type.pro
;
; PROCEDURE:
;	This function calculates seasonal values based on monthly
;	values.
;
; EXAMPLE:
;	Create a random monthly time series ten years long.
;	  x = randomn( seed, 10*12 )
;	Calculate values for the summer (June-August).
;	  result = months_to_seasons( x, 6, 3 )
;
; MODIFICATION HISTORY:
; 	Written by:	Daithi A. Stone (stoned@atm.ox.ac.uk), 2000-08-22.
;	Modified:	DAS, 2000-10-02 (debugged).
;	Modified:	DAS, 2004-12-30 (added BASEPERIOD, CYCLE, 
;			DESEASONALISE, NBASEPERIOD, NGOOD keywords; this 
;			required some altering of main algorithm)
;	Modified:	DAS, 2005-04-25 (can now deal with large number of 
;			spatial points)
;	Modified:	DAS, 2005-08-25 (removed use of constants.pro)
;	Modified:	DAS, 2006-01-06 (altered DEASONALISE keyword to not 
;			include taking anomaly;  added ANOMALY keyword)
;	Modified:	DAS, 2006-02-14 (streamlined some of the code)
;	Modified:	DAS, 2008-03-18 (allowed larger time dimensions through 
;			long integer indices)
;-

;***********************************************************************

FUNCTION MONTHS_TO_SEASONS, $
	Data, $
	Season, $
	Seasonlen, $
	ANOMALY=anomalyopt, DESEASONALISE=deseasonaliseopt, $
	BASEPERIOD=baseperiod, NBASEGOOD=nbasegood, $
	CYCLE=cycle, $
	NGOOD=ngood

;***********************************************************************
; Constants and Variables

; Number of months in a year
mina = 12l
nan = !values.f_nan

; Season length
if not( keyword_set( seasonlen ) ) then seasonlen = 3

; Dimensions of Data
dim = size( data )
; Number of dimensions of Data
ndim = dim[0]
; Number of samples
if ndim eq 1 then begin
  npoint = 1
endif else begin
  npoint = round( product( dim[1:ndim-1] ) )
endelse
; Integer type for point counters
ipoint0 = 0
if var_type( npoint ) eq 3 then ipoint0 = long( ipoint0 )
; Number of months in Data
nmonth = dim[ndim]
; Number of seasons (years)
nyear = nmonth / mina

; Vector of the season's month index values
index = indgen( seasonlen ) - ( seasonlen - 1 ) / 2 + season

; Number of non-missing months required to calculate a seasonal value
if not( keyword_set( ngood ) ) then ngood = seasonlen

; The default base period to use as the reference when deseasonalising or 
; taking the anomaly
if not( keyword_set( baseperiod ) ) then baseperiod = [ 0, nyear - 1 ]
; The number of years in this period
nbase = baseperiod[1] - baseperiod[0] + 1
; The default number of good years required for calculating a reference
if not( keyword_set( nbasegood ) ) then nbasegood = 1

; Reform Data to space-time format.
; Take the transpose too because this speeds things up later
indata = transpose( reform( data, npoint, nmonth ) )

;***********************************************************************
; Options to Deseasonalise and Take Anomalies

; If deseasonalising has been requested
if keyword_set( deseasonaliseopt ) then begin
  ; Initialise the array containing the reference seasonal cycle.
  ; Note we are doing all months in the season because this ends up being 
  ; easier and more flexible (for instance for when seasonlen>12).
  cycle = fltarr( mina, npoint )
  ; Initialise an index of a calendar month in all of the years
  indexall = indgen( nyear ) * mina
  ; Initialise an index of a calendar month in the base years
  indexbase = ( indgen( nbase ) + baseperiod[0] ) * mina
  ; Iterate through points
  for i = ipoint0, npoint - 1 do begin
    ; Iterate through months in a season
    for j = 0, mina - 1 do begin
      ; Copy the calendar month we want
      temp = indata[j+indexbase,i]
      ; Find the good values
      id = where( finite( temp ) eq 1, nid )
      ; If we do not have enough good values then record this season as missing
      if nid lt nbasegood then begin
        cycle[j,i] = nan
      ; If we do have enough good values then record their mean
      endif else begin
        cycle[j,i] = mean( temp[id] )
      endelse
    endfor
  endfor
  ; Iterate through points
  for i = ipoint0, npoint - 1 do begin
    ; Take the anomaly from the annual mean unless we want anomaly output
    if not( keyword_set( anomalyopt ) ) then begin
      cycle[*,i] = cycle[*,i] - mean( cycle[*,i], nan=1 )
    endif
    ; Iterate through months in a season
    for j = 0, mina - 1 do begin
      ; Subtract this seasonal mean from this calendar month in all years
      indata[j+indexall,i] = indata[j+indexall,i] - cycle[j,i]
    endfor
  endfor
endif

; If we want to take the anomaly from the annual mean rather than from the 
; seasonal mean
if keyword_set( anomalyopt ) $
    and not( keyword_set( deseasonaliseopt ) ) then begin
  ; Initialise an index months in the base years
  indexbase = indgen( nbase * mina ) + baseperiod[0] * mina
  ; Iterate through points
  for i = ipoint0, npoint - 1 do begin
    ; Copy the climatology data
    temp = indata[indexbase,i]
    ; Find the good values
    id = where( finite( temp ) eq 1, nid )
    ; If we do not have enough good values then record this point as missing
    if nid lt nbasegood * mina then begin
      indata[*,i] = nan
    ; If we do have enough good values then record their mean
    endif else begin
      indata[*,i] = indata[*,i] - mean( temp, nan=1 )
    endelse
  endfor
endif

;***********************************************************************
; Convert to Seasonal Data

; Initialise the output array (later reformed to the dimensions of Data)
result = nan * fltarr( npoint, nyear )

; Iterate through years
for j = 0, nyear - 1 do begin
  ; Determine which months in our season are within our time series
  id = where( ( index + mina * j ge 0 ) $
      and ( index + mina * j le nmonth - 1 ) )
  ; Iterate through points
  for i = ipoint0, npoint - 1 do begin
    ; Copy the values for this point and year
    temp = indata[j*mina+index[id],i]
    ; Determine where we have good values
    idgood = where( finite( temp ) eq 1, nidgood )
    ; If we have enough good values then calculate the seasonal mean
    if nidgood ge ngood then result[i,j] = mean( temp[idgood] )
  endfor
endfor

;***********************************************************************
; Post-processing

; Reform various output to original dimensions
if ndim eq 1 then begin
  result = reform( result )
  if keyword_set( cycle ) then cycle = reform( cycle )
endif else begin
  result = reform( result, [dim[1:ndim-1],nyear] )
  if keyword_set( cycle ) then begin
    cycle = reform( transpose( cycle ), [dim[1:ndim-1],mina] )
  endif
endelse

;***********************************************************************
; The End

return, result
END