This notebook generates plots for DESI-doc-3060: "Observing Hour Angle Assignments for the DESI Survey".
In [1]:
%pylab inline
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 0The 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'))
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');
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')
In [14]:
LST_plot(flat) #, save='../tex/note3060/LST_flat.png')
In [15]:
LST_plot(zero0) #, save='../tex/note3060/LST_zero0.png')
In [16]:
LST_plot(zero) #, save='../tex/note3060/LST_zero.png')
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()