;
;Copyright 1996-2013 United States Government as represented by the
;Administrator of the National Aeronautics and Space Administration.
;All Rights Reserved.
; This software may be used, copied, or redistributed as long as it is not
; sold and this copyright notice is reproduced on each copy made. This
; routine is provided as is without any express or implied warranties
; whatsoever.
;
pro cmb_setemptybin2fill,b,b_nbin,fillval
ii=where( b_nbin eq 0)
if ii[0] ne -1 then b[ii] = fillval
print,'in cmb_setemptybin2fill fillval=',fillval
end
pro cmb_fill_array,b,xb,b_nbin
; fills empty bins by linear interpolating across them.
; input/output
; b[nc,nt] a two dimensional array of binned data, trailing dimesnion is time, if scaler b[1,nt]
; xb[nt] is center time of the binned data
; exptrapoled bin values set to value of closest valid bin value
; b_nbin[nc,nt] is the number of data points used in each bin to binning.
ii=where(b_nbin eq 0)
if ii[0] eq -1 then return
si0 = size(b,/structure)
ndims = si0.n_dimensions
;help,si0,/str
if ndims eq 1 then begin
iz = where(b_nbin eq 0) ;find bins with no data
iv = where(b_nbin ne 0)
if iz[0] eq -1 or n_elements(iv) lt 2 then return
y = b
yz = interpol(b[iv],xb[iv],xb[iz])
y[iz] =yz
b[*]=y
xr = [min(xb[iv],i0),max(xb[iv],i1)]
i=where(xb lt xr[0]) & if i[0] ne -1 then b[i]= b[iv[i0]]
i=where(xb gt xr[1]) & if i[0] ne -1 then b[i]= b[iv[i1]]
return
endif else begin
for ic = 0l,si0.dimensions[0]-1 do begin
iz = where(b_nbin[ic,*] eq 0) ;find bins with no data
iv = where(b_nbin[ic,*] ne 0)
if iz[0] eq -1 or n_elements(iv) lt 2 then goto, next
y = b[ic,*]
yz = interpol(b[ic,iv],xb[iv],xb[iz])
y[iz] =yz
b[ic,*]=y
xr = [min(xb[iv],i0),max(xb[iv],i1)]
i=where(xb lt xr[0]) & if i[0] ne -1 then b[ic,i]= b[ic,iv[i0]]
i=where(xb gt xr[1]) & if i[0] ne -1 then b[ic,i]= b[ic,iv[i1]]
next:
endfor
endelse
end
pro cmb_timebin_array,tin00,serin0,t0,tf,dtbin $
,tout,serout,fillval, seroutuncertainty, missing_variance_flag $
,serout_flag=serout_flag, fill_empty_bins=fill_empty_bins $
,autobad=autobad,sigmul=sigmul $
,VALIDMIN=VALIDMIN, VALIDMAX=VALIDMAX
;input
; t0-start time of binning
; tf-end time of binning
; dtbin time increment
; note: t0, tf, dtbin, tin0 must be in the same time units, i.e. Julian days, NSSDC epoch
; serin0 input time series, series of scalars, vectors or matrices of serin0 = serin0[ .,.,.,nt] nt is the number of timesteps in series
; fillval - fill value
;output
; tout - centered times of binned data ; tout will be in the same time units as tin.
; serout - averaged binned timeseries
; seroutuncertainty - standard deviaton in each bin
; missing_variance_flag - value to indicate a missing sample variance. RY (08/08/2019)
; serout_flag=serout_flag -number of samples per bin
a = serin0
;print,'min/max of a', minmax(a,fill=fillval) ;diagnostic
tin0 = tin00
npts = n_elements(tin0)
;help, fillval, validmin, validmax
cmb_filter_time_series,a,fillval=fillval, validmax=validmax, validmin=validmin ; set values outside of valdimin/max to fillval
;print, 'cmb_timebin_array' & help, a,fillval,validmin,validmax
ii=where( finite(a) eq 0 ) & if ii[0] ne -1 then a[ii]=fillval
if keyword_set(autobad) then begin
a = cmb_collapse(a, ndimsorg = ndimsorg0)
ig = cmb_autobad(a,sigmul,fill=fillval, /filter_by_fill)
ii=where(ig eq 0)
if ii[0] ne -1 then a[ii]=fillval
a = cmb_collapse(a, ndimsorg = ndimsorg0,/revert)
endif
if dtbin le 0 then begin ;don't bin data, return full resolution
tout = tin0
serout = a
serout_flag = ''
return
endif
dt = dtbin
tout = cmb_timesofbin(t0,tf,dt)
ntb = n_elements(tout)
; below if added SAB 6/19/2017
if npts eq 1 then BEGIN ; Ron Yurow pointed out an error is only 1 time instance in timeseries
print,'***** WARNING ONLY 1 TIME STEP ******'
serout = serin0 + intarr(n_elements(tout) )
serout_flag= serout*0
i = where( tin0 ge tout[0:ntb-2] and tin0 lt tout[1:*])
if i[0] ne -1 then serout_flag[i] = 1
Return
endif
si0 = size(a,/structure)
if si0.n_dimensions ne 1 then begin ;matrix scalar timeseries
ndin = si0.dimensions[0:si0.n_dimensions-1]
nt = ndin[si0.n_dimensions-1] ;time steps
na = 1l & for i=0,si0.n_dimensions-2 do na = na*ndin[i] ;number of elements in matrix at a timestep
nd1 = [na,nt]
a = reform(a,nd1,/overwrite) ; collaspe matrix times series into vector series
sia = size(a,/structure)
nc = sia.dimensions[0]
ndout = ndin & ndout[si0.n_dimensions-1] = ntb ; form internal dimensions of output array
b = dblarr(na,ntb)*a[0] ;force variable type of b to be the same as a.
endif else begin ;scalar time series
a = transpose(a)
b = transpose(dblarr(ntb)*a[0])
ndout = ntb
nc=1
endelse
b_nbin = b
b2 = b ;*1d0 ;SAB 9/25/2017 compute uncertainty
iabmap = floor((tin0-t0)/dt)
bin_variance = b
;convert matrix time series into vector time series
for ic = 0l, nc-1 do begin
;ig = reform(where( abs(a[ic,*]) lt 0.99*abs(fillval)))
ig = reform(where( abs(a[ic,*]) ne abs(fillval)))
;help,ig & stop
if ig[0] ne -1 then begin
; This variables are used for the calculation of the variance using
; Welford method
; Ron Yurow (July 24, 2018)
m = DBLARR (ntb)
s = DBLARR (ntb)
prev = DBLARR (ntb)
for jpt=0l,n_elements(ig)-1 do begin
ipt = ig[jpt]
if iabmap[ipt] ge 0 and iabmap[ipt] lt ntb then begin
b[ic,iabmap[ipt]]= b[ic,iabmap[ipt]] + a[ic,ipt]
; This method is no longer used for computing variance
; Ron Yurow (July 24, 2018)
; b2[ic,iabmap[ipt]]= b2[ic,iabmap[ipt]] + a[ic,ipt]^2 ;SAB 9/25/2017 compute uncertainty
b_nbin[ic,iabmap[ipt]] = b_nbin[ic,iabmap[ipt]] + 1
; Use the Welford method to compute the variance. The Welford method is described in the
; following pseudocode
; variance(samples):
; M := 0
; S := 0
; for k from 1 to N:
; x := samples[k]
; oldM := M
; M := M + (x-M)/k
; S := S + (x-M)*(x-oldM)
; return S/(N-1)
; The inner part of the loop is impmlemented below.
; Ron Yurow (July 24, 2018)
prev [iabmap[ipt]] = m [iabmap[ipt]]
m [iabmap[ipt]] = m [iabmap[ipt]] + (a[ic,ipt] - m [iabmap[ipt]]) / b_nbin [ic,iabmap[ipt]]
s [iabmap[ipt]] = s [iabmap[ipt]] + (a[ic,ipt] - m [iabmap[ipt]]) * (a [ic,ipt] - prev [iabmap[ipt]])
endif
endfor
; Do the actual variance calculation based on the terms previously calculated.
; Ron Yurow (July 24, 2018)
nz = WHERE (b_nbin [ic, *] gt 0)
bin_variance [ic, nz] = s [nz] / (reform (b_nbin[ic, *]) - 1) [nz]
endif
; if ic eq 5 then stop
endfor
;average
ii=where(b_nbin ne 0)
if ii[0] ne -1 then b[ii] = b[ii]/b_nbin[ii] ; compute mean of each bin
; compute uncertainty of the mean
ii=where(b_nbin gt 2) ;
;if ii[0] ne -1 then b2[ii] = $
; sqrt((b2[ii]/b_nbin[ii] - b[ii]^2) $
; *(b_nbin[ii]/(b_nbin[ii]-1))/b_nbin[ii]) ; uncertainty of the mean ;SAB 9/25/2017 compute uncertainty
if ii[0] ne -1 then BEGIN
; No longer needeed, as variance is now calculated using the Welford method.
; Ron Yurow (July 24, 2018)
; sigma2 = (b2[ii] - b_nbin[ii]*b[ii]^2)/(b_nbin[ii]-1) ; standard deviation
; b2[ii] = sqrt( sigma2 ) ;/b_nbin[ii]) ; comment out uncertainty of mean computatin SAB 3/28/2018
b2 [ii] = sqrt (bin_variance [ii]) ; Calculate standard deviation based on variance.
endif
ii=where(b_nbin le 2)
; Changed by Ron Yurow (08/08/2019)
; if ii[0] ne -1 then b2[ii] = fillval ;max(b2) changed SAB 3/28/2018
if ii[0] ne -1 then b2[ii] = missing_variance_flag
; end compute uncertainty of the mean
;fill in empty bins
if keyword_set(fill_empty_bins) then cmb_fill_array,b,tout,b_nbin else cmb_setemptybin2fill,b,b_nbin,fillval
serout = reform(b,ndout,/overwrite)
seroutuncertainty = reform(b2,ndout,/overwrite)
serout_flag = reform(b_nbin,ndout,/overwrite)
if (where(serout_flag eq 0))(0) eq -1 then emptybin=0 else emptybin = n_elements(where(serout_flag eq 0))
print,'Fraction of missing binned points = ', float(emptybin)/n_elements(serout_flag)
end