%pylab inline
# Make default figure size bigger
matplotlib.rcParams['savefig.dpi'] = 150
matplotlib.rcParams['figure.figsize'] = (4,3)

import FPIprocess
from matplotlib import dates
import pytz
import fpiinfo
from datetime import datetime, timedelta
from scipy import interpolate

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['datetime']
`%matplotlib` prevents importing * from pylab and numpy

Characterizing FPI Etalon Gap Fluctuations as a Wind Error

Update: 29 Jun 2015

Brian Harding

Navigate using the arrow keys

A few months ago, we implemented a technique to estimate the "calibration error" of our FPI wind measurements. This is in contrast to "fit error", also known as statistical error, which is caused by noise in the CCD. The "calibration error" is caused by the fact that the etalon gap drifts over time. Ideally, we measure this drift by taking laser calibration exposures often enough, but we have not always captured the history of this drift perfectly, especially in the early years. Environmental factors like etalon heaters, air conditioning, and ambient temperature fluctuations can cause rapid etalon gap fluctuations.

The first technique, implemented a few months ago, used a simple technique to account for this extra error by increasing the error bars proportional to how fast the etalon gap was changing:

instr_name = 'minime05'
year = 2013
doy = 7

wlim = ((-120,100))

lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']

direction = r['direction']

tsky = array([w for (w,d) in zip(r['sky_times'], direction) if d=='Zenith'])
vertwind = array([w for (w,d) in zip(r['LOSwind'], direction) if d=='Zenith'])
sigma_vertwind = array([w for (w,d) in zip(r['sigma_LOSwind'], direction) if d=='Zenith'])
sigma_fit_vertwind = array([w for (w,d) in zip(r['sigma_fit_LOSwind'], direction) if d=='Zenith'])

tlaser = r['laser_times']
gap = r['laser_value']['t']
sigma_gap = r['laser_stderr']['t']

figure(figsize=(7,5))

subplot(3,1,1)
errorbar(tlaser,gap, yerr=sigma_gap, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

subplot(3,1,2)
errorbar(tsky,vertwind, yerr=sigma_fit_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('Nominal Sky Measurement')
ylabel('Vertical Wind [m/s]')

subplot(3,1,3)
errorbar(tsky,vertwind, yerr=sigma_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With Calibration Error')
ylabel('Vertical Wind [m/s]')

tight_layout()

But sometimes, this technique seemed to be a bit too aggressive with adding extra error, especially during times when the etalon gap slope is large.

instr_name = 'minime01'
year = 2012
doy = 85

wlim = ((-120,100))

lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']

direction = r['direction']

tsky = array([w for (w,d) in zip(r['sky_times'], direction) if d=='Zenith'])
vertwind = array([w for (w,d) in zip(r['LOSwind'], direction) if d=='Zenith'])
sigma_vertwind = array([w for (w,d) in zip(r['sigma_LOSwind'], direction) if d=='Zenith'])
sigma_fit_vertwind = array([w for (w,d) in zip(r['sigma_fit_LOSwind'], direction) if d=='Zenith'])

tlaser = r['laser_times']
gap = r['laser_value']['t']
sigma_gap = r['laser_stderr']['t']

figure(figsize=(7,5))

subplot(3,1,1)
errorbar(tlaser,gap, yerr=sigma_gap, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

subplot(3,1,2)
errorbar(tsky,vertwind, yerr=sigma_fit_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('Nominal Sky Measurement')
ylabel('Vertical Wind [m/s]')

subplot(3,1,3)
errorbar(tsky,vertwind, yerr=sigma_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With Calibration Error')
ylabel('Vertical Wind [m/s]')

tight_layout()

A large slope is no problem, as long as the variation is captured. What's really a problem is the curvature. So instead of a first-order calibration error, let's try a second-order:

def calc_wind_calib_error(dt, gaps, my_dt):
    '''
    Calculate the error induced by the inaccurate (not imprecise!) knowledge
    of the etalon gap, at the time indicated.
    INPUTS:
        dt: array of floats, number of seconds since first laser exposure, 
            for each laser exposure. The first element should be 0.
        gaps: array of floats, the estimated gap at each laser time
        my_dt: number of seconds since first laser exposure: the time you want
               to calculate the calibration error at.
    OUTPUTS:
        wind_cal_error: m/s. This is nan if it's not straddled by laser 
                        exposures.
    HISTORY:
        2015 Feb 02: Written by Brian J. Harding
        2015 Jun 30: Rewritten by Brian J. Harding to use gap temporal curvature 
                     instead of gradient as proxy for calibration error.
    '''

    c = 3e8
    #dt = np.array([(laser_time - laser_times[0]).total_seconds() for laser_time in laser_times])

    if my_dt < dt[0] or my_dt > dt[-1]: # not bounded by laser images
        return np.nan
    elif my_dt < dt[1]: # extrapolate from second laser image
        my_dt = dt[1]
    elif my_dt > dt[-2]: # extrapolate from second-to-last laser image
        my_dt = dt[-2]
    
    # Calculate the cal_error at each laser image and linearly interpolate to my_dt
    cal_error_vec = zeros(len(dt))
    for k in range(1,len(dt)-1):
        # Extract gap before, during, and after this image
        gt0 = dt[k-1]
        gt1 = dt[k]
        gt2 = dt[k+1]
        g0 = gaps[k-1]
        g1 = gaps[k]
        g2 = gaps[k+1]
        # Estimate uncertainty due to curvature
        g1interp = (gt1-gt0) * (g2-g0)/(gt2-gt0) + g0
        g1diff = g1 - g1interp
        # Convert to an uncertainty in wind
        cal_error = abs(c/g1*g1diff)
        cal_error_vec[k] = cal_error
    
    # linearly interpolate to my_dt
    sfit = interpolate.interp1d(dt[1:-1], cal_error_vec[1:-1])
    wind_cal_error = sfit(my_dt)
    return wind_cal_error

instr_name = 'minime01'
year = 2012
doy = 85

wlim = ((-120,100))

lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']

direction = r['direction']

tsky = array([w for (w,d) in zip(r['sky_times'], direction) if d=='Zenith'])
vertwind = array([w for (w,d) in zip(r['LOSwind'], direction) if d=='Zenith'])
sigma_vertwind = array([w for (w,d) in zip(r['sigma_LOSwind'], direction) if d=='Zenith'])
sigma_fit_vertwind = array([w for (w,d) in zip(r['sigma_fit_LOSwind'], direction) if d=='Zenith'])

tlaser = r['laser_times']
gap = r['laser_value']['t']
sigma_gap = r['laser_stderr']['t']

dt = array([(t - tlaser[0]).total_seconds() for t in tlaser])
dtsky = array([(t - tlaser[0]).total_seconds() for t in tsky])
sigma_cal2_vertwind = zeros(len(vertwind))
for j in range(len(dtsky)):
    cal_err = calc_wind_calib_error(dt, gap, dtsky[j])
    sigma_cal2_vertwind[j] = cal_err
    
    
    
figure(figsize=(7,5*4./3))

subplot(4,1,1)
errorbar(tlaser,gap, yerr=sigma_gap, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

subplot(4,1,2)
errorbar(tsky,vertwind, yerr=sigma_fit_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('Nominal Sky Measurement')
ylabel('Vertical Wind [m/s]')

subplot(4,1,3)
errorbar(tsky,vertwind, yerr=sigma_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With Calibration Error')
ylabel('Vertical Wind [m/s]')

subplot(4,1,4)
sigma2_vertwind = sqrt(sigma_fit_vertwind**2 + sigma_cal2_vertwind**2)
errorbar(tsky,vertwind, yerr=sigma2_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With New Calibration Error')
ylabel('Vertical Wind [m/s]')

tight_layout()

instr_name = 'minime05'
year = 2013
doy = 7

wlim = ((-120,100))

lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']

direction = r['direction']

tsky = array([w for (w,d) in zip(r['sky_times'], direction) if d=='Zenith'])
vertwind = array([w for (w,d) in zip(r['LOSwind'], direction) if d=='Zenith'])
sigma_vertwind = array([w for (w,d) in zip(r['sigma_LOSwind'], direction) if d=='Zenith'])
sigma_fit_vertwind = array([w for (w,d) in zip(r['sigma_fit_LOSwind'], direction) if d=='Zenith'])

tlaser = r['laser_times']
gap = r['laser_value']['t']
sigma_gap = r['laser_stderr']['t']

dt = array([(t - tlaser[0]).total_seconds() for t in tlaser])
dtsky = array([(t - tlaser[0]).total_seconds() for t in tsky])
sigma_cal2_vertwind = zeros(len(vertwind))
for j in range(len(dtsky)):
    cal_err = calc_wind_calib_error(dt, gap, dtsky[j])
    sigma_cal2_vertwind[j] = cal_err
    
    
    
figure(figsize=(7,5*4./3))

subplot(4,1,1)
errorbar(tlaser,gap, yerr=sigma_gap, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

subplot(4,1,2)
errorbar(tsky,vertwind, yerr=sigma_fit_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('Nominal Sky Measurement')
ylabel('Vertical Wind [m/s]')

subplot(4,1,3)
errorbar(tsky,vertwind, yerr=sigma_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With Calibration Error')
ylabel('Vertical Wind [m/s]')

subplot(4,1,4)
sigma2_vertwind = sqrt(sigma_fit_vertwind**2 + sigma_cal2_vertwind**2)
errorbar(tsky,vertwind, yerr=sigma2_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With New Calibration Error')
ylabel('Vertical Wind [m/s]')

tight_layout()

instr_name = 'minime05'
year = 2014
doy = 333

wlim = ((-120,100))

lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']

direction = r['direction']

tsky = array([w for (w,d) in zip(r['sky_times'], direction) if d=='Zenith'])
vertwind = array([w for (w,d) in zip(r['LOSwind'], direction) if d=='Zenith'])
sigma_vertwind = array([w for (w,d) in zip(r['sigma_LOSwind'], direction) if d=='Zenith'])
sigma_fit_vertwind = array([w for (w,d) in zip(r['sigma_fit_LOSwind'], direction) if d=='Zenith'])

tlaser = r['laser_times']
gap = r['laser_value']['t']
sigma_gap = r['laser_stderr']['t']

dt = array([(t - tlaser[0]).total_seconds() for t in tlaser])
dtsky = array([(t - tlaser[0]).total_seconds() for t in tsky])
sigma_cal2_vertwind = zeros(len(vertwind))
for j in range(len(dtsky)):
    cal_err = calc_wind_calib_error(dt, gap, dtsky[j])
    sigma_cal2_vertwind[j] = cal_err
    
    
    
figure(figsize=(7,5*4./3))

subplot(4,1,1)
errorbar(tlaser,gap, yerr=sigma_gap, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

subplot(4,1,2)
errorbar(tsky,vertwind, yerr=sigma_fit_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('Nominal Sky Measurement')
ylabel('Vertical Wind [m/s]')

subplot(4,1,3)
errorbar(tsky,vertwind, yerr=sigma_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With Calibration Error')
ylabel('Vertical Wind [m/s]')

subplot(4,1,4)
sigma2_vertwind = sqrt(sigma_fit_vertwind**2 + sigma_cal2_vertwind**2)
errorbar(tsky,vertwind, yerr=sigma2_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With New Calibration Error')
ylabel('Vertical Wind [m/s]')

tight_layout()

New idea: Use difference between linear interpolation and spline fit

instr_name = 'minime05'
year = 2013
doy = 7

instr_name = 'minime01'
year = 2012
doy = 85

instr_name = 'minime05'
year = 2014
doy = 333




lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']




sigma_gap_multiplier = 10.0

laser_times = r['laser_times']
gaps = r['laser_value']['t']
sigma_gaps = r['laser_stderr']['t']
weight_gaps = 1./(sigma_gap_multiplier * sigma_gaps)
s = 1.0*len(gaps) # Expected value of chi^2

dt = array([(t - laser_times[0]).total_seconds() for t in laser_times])
sfit = interpolate.UnivariateSpline(dt, gaps, w=weight_gaps, s=s)

dtvec = linspace(dt[0],dt[-1],5000)
tvec = [laser_times[0] + timedelta(seconds=dti) for dti in dtvec]
gap_interp = sfit(dtvec)


figure(figsize=(6,2))
errorbar(laser_times, gaps, yerr=sigma_gaps, fmt='k.-')
plot(tvec, gap_interp, 'r-')
plot()
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((laser_times[0], laser_times[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

<matplotlib.text.Text at 0x7efe55f74350>

def calc_wind_calib_error(dt, gaps, sigma_gaps, my_dt):
    '''
    Calculate the error induced by the inaccurate (not imprecise!) knowledge
    of the etalon gap, at the time indicated.
    INPUTS:
        dt: array of floats, number of seconds since first laser exposure, 
            for each laser exposure. The first element should be 0.
        gaps: array of floats, the estimated gap at each laser time
        sigma_gaps: array of floats, statistical error of gaps
        my_dt: number of seconds since first laser exposure: the time you want
               to calculate the calibration error at.
    OUTPUTS:
        wind_cal_error: m/s. This is nan if it's not straddled by laser 
                        exposures.
    HISTORY:
        2015 Feb 02: Written by Brian J. Harding
        2015 Jun 30: Rewritten by Brian J. Harding to use gap temporal curvature 
                     instead of gradient as proxy for calibration error.
    '''
    #### THE IDEA ####
    # Use the difference between linear interpolation and spline interpolation
    # as the metric for calibration error.
    ##################
    
    c = 3e8
    #dt = np.array([(laser_time - laser_times[0]).total_seconds() for laser_time in laser_times])

    if my_dt < dt[0] or my_dt > dt[-1]: # not bounded by laser images
        return np.nan
    elif my_dt < dt[1]: # extrapolate from second laser image
        my_dt = dt[1]
    elif my_dt > dt[-2]: # extrapolate from second-to-last laser image
        my_dt = dt[-2]
        
    # Spline fit of gaps
    sigma_gap_multiplier = 20.0
    weight_gaps = 1./(sigma_gap_multiplier * sigma_gaps)
    s = 1.0*len(gaps) # Expected value of chi^2
    sfit = interpolate.UnivariateSpline(dt, gaps, w=weight_gaps, s=s)
    
    # Calculate the cal_error at each laser image and linearly interpolate to my_dt
    cal_error_vec = zeros(len(dt))
    for k in range(len(dt)):
        # Calculate diff b/t measured and spline interp
        g1diff = sfit(dt[k]) - gaps[k]
        cal_error = abs(c/g1*g1diff)
        cal_error_vec[k] = cal_error
    
    # linearly interpolate to my_dt
    sfit = interpolate.interp1d(dt, cal_error_vec)
    wind_cal_error = sfit(my_dt)
    return wind_cal_error

instr_name = 'minime05'
year = 2013
doy = 7

instr_name = 'minime01'
year = 2012
doy = 85

instr_name = 'minime05'
year = 2014
doy = 333

wlim = ((-120,100))

lev0 = FPIprocess.load_level0(instr_name, year, doy)
r = lev0['FPI_Results']

direction = r['direction']

tsky = array([w for (w,d) in zip(r['sky_times'], direction) if d=='Zenith'])
vertwind = array([w for (w,d) in zip(r['LOSwind'], direction) if d=='Zenith'])
sigma_vertwind = array([w for (w,d) in zip(r['sigma_LOSwind'], direction) if d=='Zenith'])
sigma_fit_vertwind = array([w for (w,d) in zip(r['sigma_fit_LOSwind'], direction) if d=='Zenith'])

tlaser = r['laser_times']
gap = r['laser_value']['t']
sigma_gap = r['laser_stderr']['t']

dt = array([(t - tlaser[0]).total_seconds() for t in tlaser])
dtsky = array([(t - tlaser[0]).total_seconds() for t in tsky])
sigma_cal2_vertwind = zeros(len(vertwind))
for j in range(len(dtsky)):
    cal_err = calc_wind_calib_error(dt, gap, sigma_gap, dtsky[j])
    sigma_cal2_vertwind[j] = cal_err
    
    
    
figure(figsize=(7,5*4./3))

subplot(4,1,1)
errorbar(tlaser,gap, yerr=sigma_gap, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
title('Laser Calibration')
ylabel('Etalon Gap [m]')

subplot(4,1,2)
errorbar(tsky,vertwind, yerr=sigma_fit_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('Nominal Sky Measurement')
ylabel('Vertical Wind [m/s]')

subplot(4,1,3)
errorbar(tsky,vertwind, yerr=sigma_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With Calibration Error')
ylabel('Vertical Wind [m/s]')

subplot(4,1,4)
sigma2_vertwind = sqrt(sigma_fit_vertwind**2 + sigma_cal2_vertwind**2)
errorbar(tsky,vertwind, yerr=sigma2_vertwind, fmt='k.-')
gca().xaxis.set_major_formatter(dates.DateFormatter('%H'))
grid()
xlim((tlaser[0], tlaser[-1]))
ylim(wlim)
title('With New Calibration Error')
ylabel('Vertical Wind [m/s]')

tight_layout()

import subprocess
import shutil
nb_name = 'FPI_Gap_Misknowledge_Error_Update.ipynb' # There's probably a way to get this programmatically
cmd = ['ipython', 'nbconvert', '--to', 'slides', nb_name, \
       '--reveal-prefix', '"http://cdn.jsdelivr.net/reveal.js/2.6.2"',\
      '--template','output_toggle']
# The --template command points to a custom template which opens and closes the input code cells with a click
subprocess.call(cmd)
html_name = '%s.slides.html' % (nb_name.split('.')[0])
shutil.copy(html_name, '/home/bhardin2/public_html/slideshows/')