Observing Hour Angle Assignments

This notebook generates plots for DESI-doc-3060: "Observing Hour Angle Assignments for the DESI Survey".


In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
import astropy.table
import astropy.units as u

In [3]:
import desisurvey.config
import desisurvey.optimize
import desisurvey.old.schedule
import desisurvey.plots

Initialize a configuration object for finding files in $DESISURVEY:


In [4]:
config = desisurvey.config.Configuration()

Load scheduler with tabulated data needed to create an optimizer below. Assumes that $DESISURVEY is set and contains the necessary file.


In [5]:
scheduler = desisurvey.old.schedule.Scheduler()

Load optimized HA assignments calculated with flat initialization using:

surveyinit --init flat --save flat.fits
surveyinit --init flat --save flat0.fits --max-cycles 0

The first command takes ~25 minutes to run. The second command is much faster and just saves the initial HA assignments for plotting below.


In [6]:
flat = astropy.table.Table.read(config.get_path('flat.fits'))
flat0 = astropy.table.Table.read(config.get_path('flat0.fits'))

Load optimized HA assignments calculated with HA=0 initialization using:

surveyinit --init zero --save zero.fits
surveyinit --init zero --save zero0.fits --max-cycles 0

In [7]:
zero = astropy.table.Table.read(config.get_path('zero.fits'))
zero0 = astropy.table.Table.read(config.get_path('zero0.fits'))

Sky Plots of HA Assignments


In [8]:
def HA_plot(what, save=None):
    desisurvey.plots.plot_sky_passes(what['RA'], what['DEC'], what['PASS'], what['HA'],
                                 clip_lo=-35, clip_hi=+35, cmap='seismic',
                                 label='Hour Angle [deg]', save=save);

In [9]:
HA_plot(flat0) #, save='../tex/note3060/HA_flat0.png');



In [10]:
HA_plot(flat) #, save='../tex/note3060/HA_flat.png');



In [11]:
HA_plot(zero) #, save='../tex/note3060/HA_zero.png');


Histograms of LST Allocation


In [12]:
def LST_plot(what, nbins=192, save=None):
    
    tileid = what['TILEID']
    passnum = what['PASS']
    ha = what['HA']
    texp = what['OBSTIME']
    passes = dict(DARK=(0, 3), GRAY=(4, 4), BRIGHT=(5, 7))
    
    fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)

    for i, program in enumerate(('DARK', 'GRAY', 'BRIGHT')):
        sel = (passnum >= passes[program][0]) & (passnum <= passes[program][1])
        opt = desisurvey.optimize.Optimizer(
            scheduler, program, tileid[sel], init='array',
            initial_ha=ha[sel], nbins=nbins, stretch=1.0,
            smoothing_radius=10*u.deg)
        
        # Compare predicted (with no stretch) and tabulated exposure times.
        texp_pred, _ = opt.get_exptime(opt.ha)
        texp_pred *= 24. * 3600. / 360. # convert from deg to sec
        ratio = texp[sel] / texp_pred
        print('{0} stretch = {1:.3f} +/- {2:.3f}, loss = {3:.3f}%'
              .format(program, np.median(ratio), np.std(ratio),
                      1e2 * opt.loss_history[-1]))
        assert np.std(ratio) < 0.005
        
        # Update stretch and recalculate plan.
        opt.stretch = np.median(ratio)
        opt.plan_tiles = opt.get_plan(opt.ha)
        opt.use_plan()
        
        plan_sum = opt.plan_hist.sum()
        avail_sum = opt.lst_hist_sum
        margin = (avail_sum - plan_sum) / plan_sum
        print('{} plan uses {:.1f}h with {:.1f}h avail ({:.1f}% margin).'
         .format(program, plan_sum, avail_sum, 1e2 * margin))

        # Plot results for this program.
        ax = axes[i]
        # Plot available time.
        ax.hist(opt.lst_centers, bins=opt.lst_edges, weights=opt.lst_hist,
                histtype='stepfilled', fc=(1,0.7,0.7), ec='none', label='Avail')
        # Plot planned LST usage.
        ax.hist(opt.lst_centers, bins=opt.lst_edges, weights=opt.plan_hist,
                histtype='stepfilled', fc=(0,0,1,0.25), ec='none', label='Plan')
        # Plot rescaled plan to better visualize MSE.
        rescale = opt.lst_hist_sum / opt.plan_hist.sum()
        ax.hist(opt.lst_centers, bins=opt.lst_edges, weights=rescale * opt.plan_hist,
                histtype='step', ec='b', lw=1, alpha=0.75)

        ax.legend(loc='upper left')
        ax.set_ylabel('{0} Hours / LST bin'.format(program), fontsize='large')
        ax.set_xlim(opt.lst_edges[0], opt.lst_edges[-1])
        
        
    ax.set_xlabel('Local Sidereal Time [deg]', fontsize='large')
    plt.subplots_adjust(hspace=0.07, left=0.05, right=0.98, top=0.98, bottom=0.05)
    if save:
        plt.savefig(save)

In [13]:
LST_plot(flat0) #, save='../tex/note3060/LST_flat0.png')


INFO:optimize.py:171:__init__: DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
DARK stretch = 1.400 +/- 0.000, loss = 2.961%
DARK plan uses 4801.6h with 5319.0h avail (10.8% margin).
INFO:optimize.py:171:__init__: GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
GRAY stretch = 1.500 +/- 0.000, loss = 3.551%
GRAY plan uses 1289.8h with 1357.8h avail (5.3% margin).
INFO:optimize.py:171:__init__: BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
INFO:DESI:BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
BRIGHT stretch = 2.000 +/- 0.000, loss = 4.140%
BRIGHT plan uses 1560.9h with 1899.8h avail (21.7% margin).

In [14]:
LST_plot(flat) #, save='../tex/note3060/LST_flat.png')


INFO:optimize.py:171:__init__: DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
DARK stretch = 1.400 +/- 0.000, loss = 2.603%
DARK plan uses 4784.9h with 5319.0h avail (11.2% margin).
INFO:optimize.py:171:__init__: GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
GRAY stretch = 1.500 +/- 0.000, loss = 3.118%
GRAY plan uses 1284.4h with 1357.8h avail (5.7% margin).
INFO:optimize.py:171:__init__: BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
INFO:DESI:BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
BRIGHT stretch = 2.000 +/- 0.000, loss = 3.613%
BRIGHT plan uses 1553.0h with 1899.8h avail (22.3% margin).

In [15]:
LST_plot(zero0) #, save='../tex/note3060/LST_zero0.png')


INFO:optimize.py:171:__init__: DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
DARK stretch = 1.400 +/- 0.000, loss = 0.000%
DARK plan uses 4663.5h with 5319.0h avail (14.1% margin).
INFO:optimize.py:171:__init__: GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
GRAY stretch = 1.500 +/- 0.000, loss = 0.000%
GRAY plan uses 1245.6h with 1357.8h avail (9.0% margin).
INFO:optimize.py:171:__init__: BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
INFO:DESI:BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
BRIGHT stretch = 2.000 +/- 0.000, loss = 0.000%
BRIGHT plan uses 1498.8h with 1899.8h avail (26.8% margin).

In [16]:
LST_plot(zero) #, save='../tex/note3060/LST_zero.png')


INFO:optimize.py:171:__init__: DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:DARK: 5319.0h for 8038 tiles (texp_nom 1000.0 s, stretch 1.000).
DARK stretch = 1.400 +/- 0.000, loss = 3.329%
DARK plan uses 4818.8h with 5319.0h avail (10.4% margin).
INFO:optimize.py:171:__init__: GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
INFO:DESI:GRAY: 1357.8h for 2005 tiles (texp_nom 1000.0 s, stretch 1.000).
WARNING:optimize.py:258:__init__: Clipped 3 HA assignments to airmass limits.
WARNING:DESI:Clipped 3 HA assignments to airmass limits.
WARNING:optimize.py:260:__init__: Max clip is 0.1 deg for tile 23665.
WARNING:DESI:Max clip is 0.1 deg for tile 23665.
GRAY stretch = 1.500 +/- 0.000, loss = 3.898%
GRAY plan uses 1294.1h with 1357.8h avail (4.9% margin).
INFO:optimize.py:171:__init__: BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
INFO:DESI:BRIGHT: 1899.8h for 6028 tiles (texp_nom 300.0 s, stretch 1.000).
BRIGHT stretch = 2.000 +/- 0.000, loss = 3.915%
BRIGHT plan uses 1557.5h with 1899.8h avail (22.0% margin).

Exposure Time Factors


In [17]:
import surveysim.weather

In [18]:
import desisurvey.etc

In [19]:
def etc_factors():
    # Look up simulated weather from a surveysim run.
    w = surveysim.weather.Weather(restore='weather_123.fits')
    transp = w._table['transparency']
    seeing = w._table['seeing']
    f1 = desisurvey.etc.transparency_exposure_factor(transp)
    f2 = desisurvey.etc.seeing_exposure_factor(seeing)
    
    # Calculate empirical exposure time factor for transparency.
    F1data = 1. / np.mean(1. / f1)
    # Calculate predicted factor from mean, variance of transparency.
    mu1, var1 = np.mean(transp), np.var(transp)
    g1 = 0.
    F1pred = 1. / mu1

    # Calculate empirical exposure time factor for seeing.
    F2data = 1. / np.mean(1. / f2)
    # Calculate predicted factor from mean, variance of seeing.
    a,b,c = 4.6, -1.55, 1.15
    x0 = 1.1
    f0 = a + b * x0 + c * x0 ** 2
    mu2, var2 = np.mean(seeing), np.var(seeing)
    fmu = a + b * mu2 + c * mu2 ** 2
    g2 = (b + 2 * c * mu2) ** 2 / fmu ** 2 - c / fmu
    F2pred = (fmu / f0) / (1 + g2 * var2)
    
    # Print summary info for the note.
    print('Seeing & ${:.2f}$" & ${:.2f}$" & ${:.2f}$" & {:.3f} & {:.3f} \\\\'
          .format(np.median(seeing), mu2, np.sqrt(var2), g2, F2pred))
    print('Transparency & ${:.2f}$ & ${:.2f}$ & ${:.2f}$ & {:.3f} & {:.3f} \\\\'
          .format(np.median(transp), mu1, np.var(var1), g1, F1pred))
    
    print('Seeing pred {:.3f} data {:.3f}'.format(F2pred, F2data))
    print('Transp pred {:.3f} data {:.3f}'.format(F1pred, F1data))
    
    F12data = 1. / np.mean(1. / (f1 * f2))
    print('Combined F12 = {:.3f} , F1 * F2 = {:.3f}'
          .format(F12data, F1data * F2data))
    
etc_factors()


Seeing & $1.19$" & $1.41$" & $0.76$" & -0.113 & 1.175 \\
Transparency & $1.00$ & $0.93$ & $0.00$ & 0.000 & 1.078 \\
Seeing pred 1.175 data 1.118
Transp pred 1.078 data 1.078
Combined F12 = 1.204 , F1 * F2 = 1.205