<jupyter_start><jupyter_text>Sky with the Moon Up* https://desi.lbl.gov/trac/wiki/CommissioningCommissioningPlanning/commishdataNightSky* https://portal.nersc.gov/project/desi/collab/nightwatch/kpno/20191112/exposures.html* http://desi-www.kpno.noao.edu:8090/nightsum/nightsum-2019-11-12/nightsum.html**Data**  * /global/projecta/projectdirs/desi/spectro/redux/daily/exposures/20191112  * /global/projecta/projectdirs/desi/spectro/redux/daily/exposures/20191112* Exposures 27337-27396**Data Model**  * https://desidatamodel.readthedocs.io/en/latest/DESI_SPECTRO_REDUX/SPECPROD/exposures/NIGHT/EXPID/sky-CAMERA-EXPID.html<NAME><jupyter_code>import os, sys
import numpy as np
import fitsio
import astropy.units as u
from astropy.table import Table, vstack
from astropy.coordinates import SkyCoord
from astropy.io import fits
import desispec.io
import desimodel.io
import speclite
import specsim.atmosphere
import specsim.simulator
from desietcimg.db import DB, Exposures
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context='talk', style='ticks', font_scale=1.1)
%matplotlib inline
topdir = os.getenv('DESI_ROOT')+'/spectro'
gfadir = '/project/projectdirs/desi/users/ameisner/GFA/reduced/v0001'<jupyter_output><empty_output><jupyter_text>Load the GFA filter curve.<jupyter_code>def load_gfa_filter():
    if False:
        rfilt = speclite.filters.load_filter('decam2014-r')
    else:
        if False:
            filtfile = '/global/homes/a/ameisner/ci_throughput/etc/gfa_filter_transmission_DESI-1297.dat'
            filtwave, filtresp, _, _ = np.loadtxt(filtfile, unpack=True)
        else:
            filt = Table.read('/global/homes/a/ameisner/ci_throughput/etc/gfa_throughput-airmass_1.00.fits')
            filtwave, filtresp = filt['LAMBDA_NM'], filt['THROUGHPUT']
            
        filtresp[filtresp < 0] = 0
        filtresp[0] = 0
        filtresp[-1] = 0
        srt = np.argsort(filtwave)
        filtwave, filtresp = filtwave[srt] * 10, filtresp[srt]
        rfilt = speclite.filters.FilterResponse(wavelength=filtwave * u.Angstrom,
                                                response=filtresp,
                                                meta=dict(group_name='gfa', band_name='r'))
        
        return rfilt
rfilt = load_gfa_filter()
plt.plot(rfilt.wavelength, rfilt.response)<jupyter_output><empty_output><jupyter_text>Specify the night and range of GFA and spectrograph exposure IDs.<jupyter_code>night = 20191112
expid_start, expid_end = 27337, 27396
expids = np.arange(expid_end - expid_start + 1) + expid_start<jupyter_output><empty_output><jupyter_text>Initialize access to the online database.https://github.com/desihub/desicmx/blob/master/analysis/gfa/DESI-Online-Database-Tutorial.ipynb<jupyter_code>if not os.path.exists('db.yaml'):
    import getpass
    pw = getpass.getpass(prompt='Enter database password: ')
    with open('db.yaml', 'w') as f:
        print('host: db.replicator.dev-cattle.stable.spin.nersc.org', file=f)
        print('dbname: desi_dev', file=f)
        print('port: 60042', file=f)
        print('user: desi_reader', file=f)
        print(f'password: {pw}', file=f)
    print('Created db.yaml')
db = DB()
ExpInfo = Exposures(db)
for col in ('airmass', 'skyra', 'skydec', 'moonangl', 'moonra', 'moondec'):
    print(col, ExpInfo(27339, col))<jupyter_output>airmass None
skyra None
skydec None
moonangl None
moonra None
moondec None<jupyter_text>Read A. Meisner's reductions to get the sky background in the GFAs vs moon separation.<jupyter_code>def read_gfa(night, expid, raw=False):
    sexpid = '{:08d}'.format(expid)
    rawfname = topdir+'/data/{0}/{1}/gfa-{1}.fits.fz'.format(night, sexpid)
    fname = gfadir+'/{0}/{1}/gfa-{1}_ccds.fits'.format(night, sexpid)    
    if os.path.isfile(fname):
        # get the target-moon separation
        hdr = fitsio.read_header(rawfname, ext=1)
        #print(hdr)
        
        ra, dec, moonra, moondec = hdr['SKYRA'], hdr['SKYDEC'], hdr['MOONRA'], hdr['MOONDEC']
        csky = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
        cmoon = SkyCoord(ra=moonra*u.degree, dec=moondec*u.degree, frame='icrs')
        moonsep = np.array(csky.separation(cmoon).value).astype('f4') # [degree]

        #print('Reading {}'.format(fname))
        data = Table.read(fname)
        #print(data.colnames)
        
        # pack into a table
        out = Table()
        out['sky_mag_ab'] = [np.mean(data['sky_mag_ab'])]
        out['sky_mag_ab_err'] = [np.std(data['sky_mag_ab']) / np.sqrt(len(data))]
        out['moonsep'] = [moonsep]
        #out['airmass'] = [hdr['AIRMASS']]
        
        out = out[np.argsort(out['moonsep'])[::-1]]
        
        return out
    else:
        return None
#data = read_gfa(night, 27390)
#data
def gfa_sky(night):
    out = []
    for expid in expids:
        dd = read_gfa(night, expid)
        if dd is not None:
            out.append(dd)
    return vstack(out)
#Table(fitsio.read(os.getenv('DESI_ROOT')+'/users/ameisner/GFA/files/skymags-prelim.fits')).colnames
Table(fitsio.read(os.getenv('DESI_ROOT')+'/users/ameisner/GFA/files/skymags-prelim.with_airmass.fits')).colnames<jupyter_output><empty_output><jupyter_text>Read Aaron's updated GFA analysis outputs.<jupyter_code>def gfa_sky_updated():
    #gfafile = os.getenv('DESI_ROOT')+'/users/ameisner/GFA/files/skymags-prelim.fits'
    gfafile = os.getenv('DESI_ROOT')+'/users/ameisner/GFA/files/skymags-prelim.with_airmass.fits'

    cat = Table(fitsio.read(gfafile, lower=True))
    cat = cat[(cat['expid'] >= expid_start) * (cat['expid'] <= expid_end)]

    ra, dec, moonra, moondec = cat['skyra'], cat['skydec'], cat['ra_moon_deg'], cat['dec_moon_deg']
    #print(moonra, moondec)
    
    csky = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
    cmoon = SkyCoord(ra=moonra*u.degree, dec=moondec*u.degree, frame='icrs')
    moonsep = np.array(csky.separation(cmoon).value) # [degree]
    
    out = Table()
    out['camera'] = cat['extname']
    out['expid'] = cat['expid']
    out['airmass'] = cat['airmass']
    out['ra'] = cat['skyra']
    out['dec'] = cat['skydec']
    out['moon_phase'] = cat['moon_phase']
    out['sky_mag_ab'] = cat['skymag_median_top_camera'].astype('f4')
    #out['sky_mag_ab_err'] = [np.std(data['sky_mag_ab']) / np.sqrt(len(data))]
    out['moonsep'] = moonsep.astype('f4')
    
    out = out[np.argsort(out['moonsep'])[::-1]]
    
    return out    
if False:
    gfa = gfa_sky(night)
else:
    gfa = gfa_sky_updated()
gfa
fig, ax = plt.subplots(figsize=(8, 6))
for cam, mark in zip(set(gfa['camera']), ('s', 'o')):
    ww = gfa['camera'] == cam
    ax.scatter(gfa['moonsep'][ww], gfa['sky_mag_ab'][ww], 
               label=cam, marker=mark, s=80, alpha=0.7)
    #ax.invert_yaxis()
    #ax.set_ylim(17.5, 14.5)
    ax.set_xlabel('Target-Moon Separation (degree)')
    ax.set_ylabel('Sky Brightness (AB mag)')
ax.legend()<jupyter_output><empty_output><jupyter_text>Read the spectroscopic reductions to get the sky spectra vs moon separation.<jupyter_code>skymodel = specsim.simulator.Simulator('desi').atmosphere
skymodel.airmass = 1.1
skymodel.moon.moon_phase = 0.05
skymodel.moon.moon_zenith = 40 * u.deg
for moonsep in (45, 10):
    skymodel.moon.separation_angle = moonsep * u.deg
    skymodelflux = (skymodel.surface_brightness * np.pi * (0.75 * u.arcsec) ** 2).to(
        u.erg / (u.Angstrom * u.cm ** 2 * u.s))
    plt.plot(skymodel._wavelength, skymodelflux.value)
def read_spec(night, expid):
    sexpid = '{:08d}'.format(expid)
    datadir = topdir+'/redux/daily/exposures/{}/{}/'.format(night, sexpid)
    if os.path.isdir(datadir):
        #fr = desispec.io.read_frame('{}/frame-r3-{}.fits'.format(datadir, sexpid))
        sp = desispec.io.read_frame('{}/sframe-r3-{}.fits'.format(datadir, sexpid))
        sky = desispec.io.read_sky('{}/sky-r3-{}.fits'.format(datadir, sexpid))
        return sp, sky
        #return fr, sp, sky
    else:
        return None, None    
#sp, sky = read_spec(night, 27339)
def spec_sky():
    import astropy.time
    import astropy.coordinates
  
    # hack!
    from astropy.utils.iers import conf
    from astropy.utils import iers
    conf.auto_max_age = None  
    iers.Conf.iers_auto_url.set('ftp://cddis.gsfc.nasa.gov/pub/products/iers/finals2000A.all')    
    
    moonphase = 0.99 # from elog
    loc = astropy.coordinates.EarthLocation.of_site('Kitt Peak')
    
    # skymodel
    skymodel = specsim.simulator.Simulator('desi').atmosphere

    _rfilt = speclite.filters.FilterSequence([rfilt])
    rand = np.random.RandomState(seed=1)
    allmoonsep, meansky, stdsky, skymodelmag = [], [], [], []
    specsky, specwave = [], []
    for expid in expids:
        sp, sky = read_spec(night, expid)
        if sky is not None:
            # blarg! there's no metadata in the headers (or in the database), so
            # choose the RA, Dec, and airmass from the nearest GFA exposure.
            this = np.argmin(np.abs(gfa['expid'] - expid))
            ra, dec = gfa['ra'][this], gfa['dec'][this]
            moonphase, airmass = gfa['moon_phase'][this], gfa['airmass'][this], 
            
            # get the object-moon separation!!
            #moonsep.append(sky.header['PROGRAM'][:3])
            hdr = sky.header
            date, obstime = hdr['DATE-OBS'], hdr['TIME-OBS']
            mjd, exptime = hdr['MJD-OBS'], hdr['EXPTIME']
            
            time = astropy.time.Time(mjd, format='mjd')
            #time = astropy.time.Time(date, format='isot', scale='utc')
            moonpos = astropy.coordinates.get_moon(time, loc)
            moonra, moondec = moonpos.ra, moonpos.dec
            
            csky = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
            cmoon = SkyCoord(ra=moonra, dec=moondec, frame='icrs')
            moonsep = csky.separation(cmoon)
            allmoonsep.append(moonsep.value)
            
            moon_altaz = moonpos.transform_to(astropy.coordinates.AltAz(obstime=time, location=loc))
            moon_az = moon_altaz.az.value * u.degree
            moon_zenith = (90. - moon_altaz.alt.value) * u.degree
            
            # get the model sky brightness
            skymodel.airmass = airmass
            skymodel.moon.moon_phase = 1 - moonphase
            skymodel.moon.moon_zenith = moon_zenith
            skymodel.moon.separation_angle = moonsep
            print(skymodel.airmass, skymodel.moon.moon_phase, skymodel.moon.moon_zenith.value, 
                  skymodel.moon.separation_angle.value)
            
            skymodelflux = (skymodel.surface_brightness * np.pi * (0.75 * u.arcsec) ** 2).to(
                u.erg / (u.Angstrom * u.cm ** 2 * u.s))
            pad_skymodelflux, pad_skymodelwave = _rfilt.pad_spectrum(
                skymodelflux.value, skymodel._wavelength)
            
            #plt.plot(pad_skymodelwave, pad_skymodelflux)
            #import pdb ; pdb.set_trace()
            
            skymodel_abmags = _rfilt.get_ab_magnitudes(pad_skymodelflux, pad_skymodelwave)[rfilt.name]
            skymodelmag.append(skymodel_abmags)

            # now the data: convolve the spectrum with the r-band filter curve
            keep = sp.fibermap['OBJTYPE'] == 'SKY'
            padflux, padwave = _rfilt.pad_spectrum(sky.flux[keep, :], sky.wave, method='edge')
            abmags = _rfilt.get_ab_magnitudes(padflux, padwave)[rfilt.name]
            
            meansky.append(np.mean(abmags))
            stdsky.append(np.std(abmags) / np.sqrt(len(abmags)))
            
            # get the median spectra for this moon separation
            specwave.append(sky.wave)
            specsky.append(np.percentile(sky.flux[keep, :], axis=0, q=50))
            #q25sky = np.percentile(sky.flux[keep, :], axis=0, q=25)
            #q75sky = np.percentile(sky.flux[keep, :], axis=0, q=75)
            #these = rand.choice(keep, size=20, replace=False)
            #[plt.plot(sky.wave, sky.flux[ii, :], alpha=1.0) for ii in these]
            #plt.plot(sky.wave, medsky, alpha=0.7, color='k')            
            
    out = Table()
    out['moonsep'] = np.hstack(allmoonsep).astype('f4')
    out['meansky'] = np.hstack(meansky).astype('f4')
    out['stdsky'] = np.hstack(stdsky).astype('f4')
    out['specwave'] = specwave
    out['specsky'] = specsky
    out['skymodel'] = np.hstack(skymodelmag).astype('f4')
    
    out = out[np.argsort(out['moonsep'])[::-1]]
    
    return out
spec = spec_sky()
spec
fig, ax = plt.subplots(figsize=(8, 6))

ref = spec['moonsep'].argmax()
#print(ref, spec['specsky'][ref])
for ss in spec:
    #print(ss['specwave'], ss['specsky'] / spec['specsky'][ref])
    ax.plot(ss['specwave'], ss['specsky'],# / spec['specsky'][0],
            label=r'{:.0f}$^{{\circ}}$'.format(ss['moonsep']))
ax.plot(rfilt.wavelength, rfilt.response / np.max(rfilt.response) * np.median(spec['specsky'][-1]),
        color='k', ls='--', lw=2)

ax.set_xlim(5500, 7800)
#ax.set_yscale('log')
ax.set_xlabel('Wavelength ($\AA$)')
ax.set_ylabel('Sky Spectra (counts / $\AA$)')
ax.legend(ncol=5, fontsize=12)
#ax.set_title('Night {}'.format(night))
fig, ax = plt.subplots(figsize=(8, 6))

for cam, mark in zip(set(gfa['camera']), ('s', 'o')):
    ww = gfa['camera'] == cam
    ax.plot(gfa['moonsep'][ww], gfa['sky_mag_ab'][ww], '{}-'.format(mark),
            alpha=0.8, markersize=10, label='GFA-{}'.format(cam))

ax.plot(spec['moonsep'], spec['skymodel'], 'o--', color='k', label='Model Sky')
ax.invert_xaxis()
#ax.set_ylim(17.5, 14.5)
ax.set_xlabel('Target-Moon Separation (degree)')
ax.set_ylabel('$r$-band Sky Brightness (AB mag)')
ax.legend()
fig, ax = plt.subplots(figsize=(10, 8))

specref = spec['meansky'][spec['moonsep'].argmax()]

for cam, mark in zip(set(gfa['camera']), ('s', 'o')):
    ww = gfa['camera'] == cam
    gfaref = gfa['sky_mag_ab'][gfa['moonsep'].argmax()]
    ax.plot(gfa['moonsep'][ww], gfa['sky_mag_ab'][ww] - gfaref, '{}-'.format(mark),
            alpha=0.8, markersize=10, 
            label='GFA-{}'.format(cam))#, marker=mark, s=80)
#ax.scatter(gfa['moonsep'], gfa['sky_mag_ab'] - gfaref, label='GFAs')

ax.plot(spec['moonsep'], spec['meansky'] - specref, '^--', 
        alpha=0.8, markersize=13, label='Spectra')

modelref = spec['skymodel'][spec['moonsep'].argmax()]
ax.plot(spec['moonsep'], spec['skymodel'] - modelref, 'p-.', 
        alpha=0.8, markersize=13, label='Sky Model')

ax.invert_xaxis()
#ax.set_ylim(17.5, 14.5)
ax.set_xlabel('Target-Moon Separation (degree)')
ax.set_ylabel('Relative $r$-band Sky Brightness (AB mag)')
ax.legend(loc='lower left')
#ax.set_title('Night {}'.format(night))<jupyter_output><empty_output><jupyter_text>Playing around<jupyter_code>import specsim.simulator
desi = specsim.simulator.Simulator('desi', num_fibers=1)
desi.instrument.fiberloss_method = 'table'
desi.simulate()
plt.plot(desi.atmosphere._wavelength.value, desi.atmosphere.surface_brightness.value)
plt.xlim(5500, 8500)
from astropy.coordinates import EarthLocation
loc = EarthLocation.of_site('Kitt Peak')
loc

time = astropy.time.Time(t['DATE'], format='jd')

moon_position = astropy.coordinates.get_moon(time, loc)
moon_ra = moon_position.ra.value
moon_dec = moon_position.dec.value
moon_position_altaz = moon_position.transform_to(astropy.coordinates.AltAz(obstime=time, location=location))
moon_alt = moon_position_altaz.alt.value
moon_az = moon_position_altaz.az.value
EarthLocation.from_geodetic(lat='-30d10m10.78s', lon='-70d48m23.49s', height=2241.4*u.m)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1613486629634414, 'fraction_numerical': 0.01808551866683891, 'mean_word_length': 3.3394058295964126, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 34, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '37035', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6104, 'n_tokens_neox': 5812, 'n_words': 925, 'path': 'analysis/specsky/sky-with-moon.ipynb', 'repo_name': 'tanveerkarim/desicmx', 'size': '15490'}