Source code for lsst.sims.maf.utils.obs2sqlite

import numpy as np
import pandas as pd
from lsst.sims.skybrightness import stupidFast_RaDec2AltAz, SkyModel
import lsst.sims.skybrightness_pre as sb
from lsst.sims.utils import raDec2Hpid, m5_flat_sed, Site
import healpy as hp
import sqlite3
import ephem
import sys

__all__ = ['obs2sqlite']


class mjd2night(object):
    def __init__(self, mjd_start=59580.035):
        self.site = Site(name='LSST')
        self.obs = ephem.Observer()
        self.obs.lat = self.site.latitude_rad
        self.obs.lon = self.site.longitude_rad
        self.obs.elevation = self.site.height
        self.mjd = mjd_start
        self.sun = ephem.Sun()
        self.generate_sunsets()

    def generate_sunsets(self, nyears=13, day_pad=50):
        """
        Generate the sunrise times for LSST so we can label nights by MJD
        """

        # Set observatory horizon to zero
        doff = ephem.Date(0)-ephem.Date('1858/11/17')

        self.obs.horizon = 0.

        # Swipe dates to match sims_skybrightness_pre365
        mjd_start = self.mjd
        mjd_end = np.arange(59560, 59560+365.25*nyears+day_pad+366, 366).max()
        step = 0.25
        mjds = np.arange(mjd_start, mjd_end+step, step)
        setting = mjds*0.

        # Stupid Dublin Julian Date
        djds = mjds - doff
        sun = ephem.Sun()

        for i, (mjd, djd) in enumerate(zip(mjds, djds)):
            sun.compute(djd)
            setting[i] = self.obs.previous_setting(sun, start=djd, use_center=True)
        setting = setting + doff

        # zomg, round off crazy floating point precision issues
        setting_rough = np.round(setting*100.)
        u, indx = np.unique(setting_rough, return_index=True)
        self.setting_sun_mjds = setting[indx]
        left = np.searchsorted(self.setting_sun_mjds, mjd_start)
        self.setting_sun_mjds = self.setting_sun_mjds[left:]

    def __call__(self, mjd):
        """
        Convert an mjd to a night integer.
        """
        return np.searchsorted(self.setting_sun_mjds, mjd)


[docs]def obs2sqlite(observations_in, location='LSST', outfile='observations.sqlite', slewtime_limit=5., full_sky=False, radians=True): """ Utility to take an array of observations and dump it to a sqlite file, filling in useful columns along the way. observations_in: numpy array with at least columns of ra : RA in degrees dec : dec in degrees mjd : MJD in day filter : string with the filter name exptime : the exposure time in seconds slewtime_limit : float Consider all slewtimes larger than this to be closed-dome time not part of a slew. """ # Set the location to be LSST if location == 'LSST': telescope = Site('LSST') # Check that we have the columns we need needed_cols = ['ra', 'dec', 'mjd', 'filter'] in_cols = observations_in.dtype.names for col in needed_cols: if needed_cols not in in_cols: ValueError('%s column not found in observtion array' % col) n_obs = observations_in.size sm = None # make sure they are in order by MJD observations_in.sort(order='mjd') # Take all the columns that are in the input and add any missing names = ['filter', 'ra', 'dec', 'mjd', 'exptime', 'alt', 'az', 'skybrightness', 'seeing', 'night', 'slewtime', 'fivesigmadepth', 'airmass', 'sunAlt', 'moonAlt'] types = ['|S1'] types.extend([float]*(len(names)-1)) observations = np.zeros(n_obs, dtype=list(zip(names, types))) # copy over the ones we have for col in in_cols: observations[col] = observations_in[col] # convert output to be in degrees like expected if radians: observations['ra'] = np.degrees(observations['ra']) observations['dec'] = np.degrees(observations['dec']) if 'exptime' not in in_cols: observations['exptime'] = 30. # Fill in the slewtime. Note that filterchange time gets included in slewtimes if 'slewtime' not in in_cols: # Assume MJD is midpoint of exposures mjd_sec = observations_in['mjd']*24.*3600. observations['slewtime'][1:] = mjd_sec[1:]-mjd_sec[0:-1] - observations['exptime'][0:-1]*0.5 - observations['exptime'][1:]*0.5 closed = np.where(observations['slewtime'] > slewtime_limit*60.) observations['slewtime'][closed] = 0. # Let's just use the stupid-fast to get alt-az if 'alt' not in in_cols: alt, az = stupidFast_RaDec2AltAz(np.radians(observations['ra']), np.radians(observations['dec']), telescope.latitude_rad, telescope.longitude_rad, observations['mjd']) observations['alt'] = np.degrees(alt) observations['az'] = np.degrees(az) # Fill in the airmass if 'airmass' not in in_cols: observations['airmass'] = 1./np.cos(np.pi/2. - np.radians(observations['alt'])) # Fill in the seeing if 'seeing' not in in_cols: # XXX just fill in a dummy val observations['seeing'] = 0.8 if 'night' not in in_cols: m2n = mjd2night() observations['night'] = m2n(observations['mjd']) # Sky Brightness if 'skybrightness' not in in_cols: if full_sky: sm = SkyModel(mags=True) for i, obs in enumerate(observations): sm.setRaDecMjd(obs['ra'], obs['dec'], obs['mjd'], degrees=True) observations['skybrightness'][i] = sm.returnMags()[obs['filter']] else: # Let's try using the pre-computed sky brighntesses sm = sb.SkyModelPre(preload=False) full = sm.returnMags(observations['mjd'][0]) nside = hp.npix2nside(full['r'].size) imax = float(np.size(observations)) for i, obs in enumerate(observations): indx = raDec2Hpid(nside, obs['ra'], obs['dec']) observations['skybrightness'][i] = sm.returnMags(obs['mjd'], indx=[indx])[obs['filter']] sunMoon = sm.returnSunMoon(obs['mjd']) observations['sunAlt'][i] = sunMoon['sunAlt'] observations['moonAlt'][i] = sunMoon['moonAlt'] progress = i/imax*100 text = "\rprogress = %.2f%%"%progress sys.stdout.write(text) sys.stdout.flush() observations['sunAlt'] = np.degrees(observations['sunAlt']) observations['moonAlt'] = np.degrees(observations['moonAlt']) # 5-sigma depth for fn in np.unique(observations['filter']): good = np.where(observations['filter'] == fn) observations['fivesigmadepth'][good] = m5_flat_sed(fn, observations['skybrightness'][good], observations['seeing'][good], observations['exptime'][good], observations['airmass'][good]) conn = sqlite3.connect(outfile) df = pd.DataFrame(observations) df.to_sql('observations', conn)