From 5e72b92aa5b2c9529a410963a553d0161a68e6e3 Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Tue, 28 Nov 2017 15:27:01 -0800 Subject: [PATCH 1/8] Start work on sun altitude interpolator --- py/desisurvey/ephemerides.py | 92 +++++++++++++++++++++++++++++++++--- 1 file changed, 86 insertions(+), 6 deletions(-) diff --git a/py/desisurvey/ephemerides.py b/py/desisurvey/ephemerides.py index 916c3404..599d5848 100644 --- a/py/desisurvey/ephemerides.py +++ b/py/desisurvey/ephemerides.py @@ -169,6 +169,16 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25, # tabulating the (ra,dec) of objects to avoid. t_obj = np.linspace(0., 1., num_obj_steps) + # Lookup sun altitude thresholds. + bright_max_alt = ( + config.programs.BRIGHT.max_sun_altitude().to(u.rad).value) + dark_max_alt = ( + config.programs.DARK.max_sun_altitude().to(u.rad).value) + + rng = np.rad2deg(bright_max_alt - dark_max_alt) + print('limits', rng) + max_err = 0. + # Calculate ephmerides for each night. for day_offset in range(num_nights): day = self.start + day_offset * u.day @@ -177,19 +187,34 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25, # Store local noon for this day. row['noon'] = day.mjd # Calculate bright twilight. - mayall.horizon = ( - config.programs.BRIGHT.max_sun_altitude().to(u.rad).value) + mayall.horizon = bright_max_alt row['brightdusk'] = mayall.next_setting( ephem.Sun(), use_center=True) + mjd0 row['brightdawn'] = mayall.next_rising( ephem.Sun(), use_center=True) + mjd0 # Calculate dark / gray twilight. - mayall.horizon = ( - config.programs.DARK.max_sun_altitude().to(u.rad).value) + mayall.horizon = dark_max_alt row['dusk'] = mayall.next_setting( ephem.Sun(), use_center=True) + mjd0 row['dawn'] = mayall.next_rising( ephem.Sun(), use_center=True) + mjd0 + # Calculate the midpoint between BRIGHT/DARK twilight. + mayall.horizon = 0.5 * ( + config.programs.BRIGHT.max_sun_altitude().to(u.rad).value + + config.programs.DARK.max_sun_altitude().to(u.rad).value) + middusk = mayall.next_setting( + ephem.Sun(), use_center=True) + mjd0 + middawn = mayall.next_rising( + ephem.Sun(), use_center=True) + mjd0 + # Compare with linear extrapolations. + middusk_pred = 0.5 * (row['brightdusk'] + row['dusk']) + middawn_pred = 0.5 * (row['brightdawn'] + row['dawn']) + dawn = (row['brightdawn'] - row['dawn']) + dusk = (row['dusk'] - row['brightdusk']) + dusk_err = (middusk_pred - middusk) / dusk * rng + dawn_err = (middusk_pred - middusk) / dawn * rng + #print('dawn/dusk (arcsec)', dusk_err * 3600, dawn_err * 360) + max_err = max(max_err, abs(dusk_err), abs(dawn_err)) # Calculate the moonrise/set for any moon visible tonight. m0 = ephem.Moon() mayall.horizon = '-0:34' # the value that the USNO uses. @@ -212,6 +237,8 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25, row[name + '_ra'][i] = math.degrees(float(model.ra)) row[name + '_dec'][i] = math.degrees(float(model.dec)) + print('max_err', max_err * 3600, 'arcsec') + # Initialize a grid covering each night with 15sec resolution # for tabulating the night program. t_grid = get_grid(step_size=15 * u.s) @@ -291,8 +318,7 @@ def get_moon_illuminated_fraction(self, mjd): Parameters ---------- mjd : float or array - MJD values during a single night where the program should be - tabulated. + MJD values where the illuminated fraction should be tabulated. Returns ------- @@ -423,6 +449,60 @@ def is_full_moon(self, night, num_nights=None): return index - lo in sort_order[-num_nights:] +def get_sun_altitude(row, mjd): + """Return the sun altitude at the specified time for twilight modeling. + + This function only returns values in the range [-10, -20] deg + since values outside this range are either never used for observing + or else considered completely dark. + + The calculation uses linear interpolation of the tabulated times when the + sun passes through the bright and dark maximum altitudes, nominally -13 + and -15 deg, which provides an altitude accurate to better than 30" + between these limits. + + Parameters + ---------- + row : astropy.table.Row + A single row from the ephemerides astropy Table corresponding to the + night in question. + mjd : float + MJD value(s) during a single night where the sun altitude should be + calculated. + + Returns + ------- + astropy.units.Quantity + Sun altitude(s) in degrees at each input mjd value, clipped to + the range [-10, -20]. + """ + # Lookup the bright, dark threshold altitudes. + config = desisurvey.config.Configuration() + bright_max_alt = config.programs.BRIGHT.max_sun_altitude() + dark_max_alt = config.programs.DARK.max_sun_altitude() + + # Validate input timestamp(s) + mjd = np.asarray(mjd) + scalar = mjd.ndim == 0 + mjd = np.atleast1d(mjd) + midnight = 0.5 * (row['dusk'] + row['dawn']) + if np.any(np.abs(mjd - midnight) > 0.5): + raise ValueError('mjd is not associated with the specified night') + + # Use linear interpolation on the dusk timestamps for times before midnight. + sun_alt = np.empty_like(mjd) + at_dusk = mjd < midnight + sun_alt[at_dusk] = bright_max_alt + (bright_max_alt - dark_max_alt) * ( + mjd[at_dusk] - row['brightdusk']) / (row['dusk'] - row['brightdusk']) + + # Use linear interpolation on the dawn timestamps after midnight. + at_dawn = ~at_dusk + sun_alt[at_dawn] = bright_max_alt + (bright_max_alt - dark_max_alt) * ( + mjd[at_dawn] - row['brightdawn']) / (row['dusk'] - row['brightdusk']) + + return sun_alt + + def get_object_interpolator(row, object_name, altaz=False): """Build an interpolator for object location during one night. From 82a851fed432671269241e6f88f36c6c53bab65e Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Tue, 28 Nov 2017 17:43:13 -0800 Subject: [PATCH 2/8] Add Ephemerides method to calculate sun altitude during twilight --- doc/changes.rst | 1 + py/desisurvey/ephemerides.py | 37 ++++----------------- py/desisurvey/test/test_ephemerides.py | 46 ++++++++++++++++++++++---- 3 files changed, 48 insertions(+), 36 deletions(-) diff --git a/doc/changes.rst b/doc/changes.rst index fef015bd..06716030 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -5,6 +5,7 @@ desisurvey change log 0.10.1 (unreleased) ------------------- +* Add Ephemerides method to calculate sun altitude during twilight. * Set the ``EXTNAME`` keyword on the Table returned by ``Progress.get_exposures()``. 0.10.0 (2017-11-09) diff --git a/py/desisurvey/ephemerides.py b/py/desisurvey/ephemerides.py index 599d5848..996e5e74 100644 --- a/py/desisurvey/ephemerides.py +++ b/py/desisurvey/ephemerides.py @@ -175,10 +175,6 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25, dark_max_alt = ( config.programs.DARK.max_sun_altitude().to(u.rad).value) - rng = np.rad2deg(bright_max_alt - dark_max_alt) - print('limits', rng) - max_err = 0. - # Calculate ephmerides for each night. for day_offset in range(num_nights): day = self.start + day_offset * u.day @@ -198,23 +194,6 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25, ephem.Sun(), use_center=True) + mjd0 row['dawn'] = mayall.next_rising( ephem.Sun(), use_center=True) + mjd0 - # Calculate the midpoint between BRIGHT/DARK twilight. - mayall.horizon = 0.5 * ( - config.programs.BRIGHT.max_sun_altitude().to(u.rad).value + - config.programs.DARK.max_sun_altitude().to(u.rad).value) - middusk = mayall.next_setting( - ephem.Sun(), use_center=True) + mjd0 - middawn = mayall.next_rising( - ephem.Sun(), use_center=True) + mjd0 - # Compare with linear extrapolations. - middusk_pred = 0.5 * (row['brightdusk'] + row['dusk']) - middawn_pred = 0.5 * (row['brightdawn'] + row['dawn']) - dawn = (row['brightdawn'] - row['dawn']) - dusk = (row['dusk'] - row['brightdusk']) - dusk_err = (middusk_pred - middusk) / dusk * rng - dawn_err = (middusk_pred - middusk) / dawn * rng - #print('dawn/dusk (arcsec)', dusk_err * 3600, dawn_err * 360) - max_err = max(max_err, abs(dusk_err), abs(dawn_err)) # Calculate the moonrise/set for any moon visible tonight. m0 = ephem.Moon() mayall.horizon = '-0:34' # the value that the USNO uses. @@ -237,8 +216,6 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25, row[name + '_ra'][i] = math.degrees(float(model.ra)) row[name + '_dec'][i] = math.degrees(float(model.dec)) - print('max_err', max_err * 3600, 'arcsec') - # Initialize a grid covering each night with 15sec resolution # for tabulating the night program. t_grid = get_grid(step_size=15 * u.s) @@ -484,23 +461,23 @@ def get_sun_altitude(row, mjd): # Validate input timestamp(s) mjd = np.asarray(mjd) scalar = mjd.ndim == 0 - mjd = np.atleast1d(mjd) + mjd = np.atleast_1d(mjd) midnight = 0.5 * (row['dusk'] + row['dawn']) if np.any(np.abs(mjd - midnight) > 0.5): raise ValueError('mjd is not associated with the specified night') - # Use linear interpolation on the dusk timestamps for times before midnight. - sun_alt = np.empty_like(mjd) + # Use linear interpolation on the dusk data for times before midnight. + sun_alt = np.empty_like(mjd) * u.degree at_dusk = mjd < midnight - sun_alt[at_dusk] = bright_max_alt + (bright_max_alt - dark_max_alt) * ( + sun_alt[at_dusk] = bright_max_alt + (dark_max_alt - bright_max_alt) * ( mjd[at_dusk] - row['brightdusk']) / (row['dusk'] - row['brightdusk']) - # Use linear interpolation on the dawn timestamps after midnight. + # Use linear interpolation on the dawn data for times after midnight. at_dawn = ~at_dusk sun_alt[at_dawn] = bright_max_alt + (bright_max_alt - dark_max_alt) * ( - mjd[at_dawn] - row['brightdawn']) / (row['dusk'] - row['brightdusk']) + mjd[at_dawn] - row['brightdawn']) / (row['brightdawn'] - row['dawn']) - return sun_alt + return np.clip(sun_alt, -20 * u.deg, -10 * u.deg) def get_object_interpolator(row, object_name, altaz=False): diff --git a/py/desisurvey/test/test_ephemerides.py b/py/desisurvey/test/test_ephemerides.py index e4395340..4daddbd8 100644 --- a/py/desisurvey/test/test_ephemerides.py +++ b/py/desisurvey/test/test_ephemerides.py @@ -12,8 +12,8 @@ import astropy.io import desisurvey.config -from desisurvey.ephemerides import Ephemerides, get_grid, get_object_interpolator -from desisurvey.utils import get_date, get_location +from desisurvey.ephemerides import * +from desisurvey import utils class TestEphemerides(unittest.TestCase): @@ -50,6 +50,12 @@ def tearDownClass(cls): # Reset our configuration. desisurvey.config.Configuration.reset() + def setUp(self): + utils.freeze_iers() + + def tearDown(self): + utils._iers_is_frozen = False + def test_getephem(self): """Tabulate one month of ephemerides""" start = datetime.date(2019, 9, 1) @@ -117,7 +123,7 @@ def test_get_grid(self): def test_moon_phase(self): """Verfify moon illuminated fraction for first week of 2020""" ephem = Ephemerides( - get_date('2019-12-31'), get_date('2020-02-02'), + utils.get_date('2019-12-31'), utils.get_date('2020-02-02'), use_cache=False, write_cache=False) for i, jd in enumerate(self.table['jd']): t = Time(jd, format='jd') @@ -128,7 +134,7 @@ def test_moon_phase(self): def test_moon_radec(self): """Verify moon (ra,dec) for first week of 2020""" ephem = Ephemerides( - get_date('2019-12-31'), get_date('2020-02-02'), + utils.get_date('2019-12-31'), utils.get_date('2020-02-02'), use_cache=False, write_cache=False) for i, jd in enumerate(self.table['jd']): t = Time(jd, format='jd') @@ -144,9 +150,9 @@ def test_moon_radec(self): def test_moon_altaz(self): """Verify moon (alt,az) for first week of 2020""" ephem = Ephemerides( - get_date('2019-12-31'), get_date('2020-02-02'), + utils.get_date('2019-12-31'), utils.get_date('2020-02-02'), use_cache=False, write_cache=False) - location = get_location() + location = utils.get_location() for i, jd in enumerate(self.table['jd']): t = Time(jd, format='jd') night = ephem.get_night(t) @@ -160,6 +166,34 @@ def test_moon_altaz(self): sep = truth.separation(calc) self.assertTrue(abs(sep.to(u.deg).value) < 0.3) + def test_sun_alt(self): + """Test interpolated sun altitude""" + start = datetime.date(2019, 9, 1) + stop = datetime.date(2019, 10, 1) + ephem = Ephemerides(start, stop, use_cache=False, write_cache=False) + config = desisurvey.config.Configuration() + balt = config.programs.BRIGHT.max_sun_altitude().to(u.deg).value + dalt = config.programs.DARK.max_sun_altitude().to(u.deg).value + for i in range(ephem.num_nights): + row = ephem.get_row(i) + dusk1, dusk2 = row['brightdusk'], row['dusk'] + dawn1, dawn2 = row['dawn'], row['brightdawn'] + assert (dusk1 < dusk2) and (dawn1 < dawn2) + sun_alt = get_sun_altitude( + row, [dusk1, 0.5 * (dusk1 + dusk2), dusk2, + dawn1, 0.5 * (dawn1 + dawn2), dawn2, + dusk1 - 0.2, dusk2 + 0.2, dawn1 - 0.2, dawn2 + 0.2] + ).to(u.deg).value + assert np.allclose( + sun_alt, [ + balt, 0.5 * (balt + dalt), dalt, + dalt, 0.5 * (dalt + balt), balt, + -10, -20, -20, -10]) + with self.assertRaises(ValueError): + get_sun_altitude(row, dusk1 - 1.) + with self.assertRaises(ValueError): + get_sun_altitude(row, [dusk1, dusk1 - 1.]) + def test_suite(): """Allows testing of only this module with the command:: From d3e62060ca0f0aa5d4ba712c68d53d417bc3752b Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Wed, 29 Nov 2017 11:59:51 -0800 Subject: [PATCH 3/8] Add exposure metadata: moonv,sunalt,sundaz,sunr --- py/desisurvey/progress.py | 43 +++++++++++++++++++---- py/desisurvey/test/test_progress.py | 54 +++++++++++++++++++---------- 2 files changed, 72 insertions(+), 25 deletions(-) diff --git a/py/desisurvey/progress.py b/py/desisurvey/progress.py index a706e5c7..a3782c55 100644 --- a/py/desisurvey/progress.py +++ b/py/desisurvey/progress.py @@ -115,6 +115,20 @@ def __init__(self, restore=None, max_exposures=32): table['moonsep'] = astropy.table.Column( length=num_tiles, shape=(max_exposures,), format='%.1f', description='Moon-tile separation angle in degrees', unit='deg') + table['moonv'] = astropy.table.Column( + length=num_tiles, shape=(max_exposures,), format='%.1f', + description='V-band magnitude of scattered moonlight') + table['sunalt'] = astropy.table.Column( + length=num_tiles, shape=(max_exposures,), format='%.1f', + description='Sun altitude angle in degrees, clipped below -20', + unit='deg') + table['sundaz'] = astropy.table.Column( + length=num_tiles, shape=(max_exposures,), format='%.1f', + description='Sun-tile relative azimuth angle in degrees', + unit='deg') + table['sunr'] = astropy.table.Column( + length=num_tiles, shape=(max_exposures,), format='%.1f', + description='r-band magnitude of twilight scattered sun') # Copy tile data. table['tileid'] = tiles['TILEID'] table['pass'] = tiles['PASS'] @@ -415,7 +429,8 @@ def copy_range(self, mjd_min=None, mjd_max=None): return Progress(restore=table) def add_exposure(self, tile_id, start, exptime, snr2frac, airmass, seeing, - transparency, moonfrac, moonalt, moonsep): + transparency, moonfrac, moonalt, moonsep, moonv, + sunalt, sundaz, sunr): """Add a single exposure to the progress. Parameters @@ -440,6 +455,16 @@ def add_exposure(self, tile_id, start, exptime, snr2frac, airmass, seeing, Moon altitude angle in degrees. moonsep : float Moon-tile separation angle in degrees. + moonv : float + The V-band magnitude of scattered moonlight, or -np.inf if + the moon is below the horizon. + sunalt : float + Sun altitude angle in degrees, clipped below -20. + sundaz : float + Sun-tile relative azimuth angle in degrees. + sunr : float + The r-band magnitude of twilight scattered sun, or -np.inf + if the sun is below -18 deg altitude. """ mjd = start.mjd self.log.info( @@ -478,6 +503,10 @@ def add_exposure(self, tile_id, start, exptime, snr2frac, airmass, seeing, row['moonfrac'][num_exp] = moonfrac row['moonalt'][num_exp] = moonalt row['moonsep'][num_exp] = moonsep + row['moonv'][num_exp] = moonv + row['sunalt'][num_exp] = sunalt + row['sundaz'][num_exp] = sundaz + row['sunr'][num_exp] = sunr # Update this tile's status. row['status'] = 1 if row['snr2frac'].sum() < self.min_snr2 else 2 @@ -486,7 +515,7 @@ def get_exposures(self, start=None, stop=None, tile_fields='tileid,pass,ra,dec,ebmv', exp_fields=('night,mjd,exptime,seeing,transparency,' + 'airmass,moonfrac,moonalt,moonsep,' + - 'program,flavor')): + 'moonv,sunalt,sundaz,sunr,program,flavor')): """Create a table listing exposures in time order. Parameters @@ -515,7 +544,8 @@ def get_exposures(self, start=None, stop=None, Returns ------- astropy.table.Table - Table with the specified columns as uppercase and one row per exposure. + Table with the specified columns as uppercase and one row + per exposure. """ # Get MJD range to show. if start is None: @@ -580,7 +610,8 @@ def get_exposures(self, start=None, stop=None, mjd = table['mjd'].flatten()[order] night = np.empty(len(mjd), dtype=(str, 8)) for i in range(len(mjd)): - night[i] = str(desisurvey.utils.get_date(mjd[i])).replace('-', '') + night[i] = str(desisurvey.utils.get_date(mjd[i]) + ).replace('-', '') output[name.upper()] = astropy.table.Column( night, description='Date at start of night when exposure taken') @@ -605,8 +636,8 @@ def get_exposures(self, start=None, stop=None, program[exppass < 4] = 'DARK' program[exppass == 4] = 'GRAY' - output[name.upper()] = astropy.table.Column(program, - description='Program name') + output[name.upper()] = astropy.table.Column( + program, description='Program name') elif name == 'flavor': flavor = np.empty(len(exppass), dtype=(str, 7)) flavor[:] = 'science' diff --git a/py/desisurvey/test/test_progress.py b/py/desisurvey/test/test_progress.py index b2c19467..d022a683 100644 --- a/py/desisurvey/test/test_progress.py +++ b/py/desisurvey/test/test_progress.py @@ -63,7 +63,8 @@ def test_add_exposures(self): t0 = astropy.time.Time('2020-01-01 07:00') for i, tile_id in enumerate(tiles): p.add_exposure( - tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) self.assertTrue(p.get_tile(tile_id)['snr2frac'][0] == 0.5) last_tile = p.get_tile(tile_id) self.assertTrue(np.array_equal( @@ -82,7 +83,8 @@ def test_restore_status(self): t0 = astropy.time.Time('2020-01-01 07:00') for i, tile_id in enumerate(tiles): p.add_exposure( - tile_id, t0 + i * u.hour, 1e3 * u.s, 2.0, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + i * u.hour, 1e3 * u.s, 2.0, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) good_status = p._table['status'].copy() p._table['status'] = 0 p2 = Progress(p._table) @@ -96,7 +98,8 @@ def test_get_exposures(self): t0 = astropy.time.Time('2020-01-01 07:00') for i, tile_id in enumerate(tiles): p.add_exposure( - tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) explist = p.get_exposures() self.assertEqual(explist.meta['EXTNAME'], 'EXPOSURES') self.assertTrue(np.all(np.diff(explist['MJD']) > 0)) @@ -135,13 +138,16 @@ def test_exposures_incrementing(self): tile_id = t['tileid'][0] t0 = astropy.time.Time('2020-01-01 07:00') t1 = t0 + 1 * u.hour - p.add_exposure(tile_id, t0, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) - p.add_exposure(tile_id, t1, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + p.add_exposure(tile_id, t0, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) + p.add_exposure(tile_id, t1, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) self.assertEqual(p.first_mjd, t0.mjd) self.assertEqual(p.last_mjd, t1.mjd) self.assertEqual(p.num_exp, 2) with self.assertRaises(ValueError): - p.add_exposure(tile_id, t0, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + p.add_exposure(tile_id, t0, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) def test_save_read(self): """Create, save and read a progress table""" @@ -150,7 +156,8 @@ def test_save_read(self): t0 = astropy.time.Time('2020-01-01 07:00') for i, tile_id in enumerate(tiles): p1.add_exposure( - tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) p1.save('p1.fits') p2 = Progress('p1.fits') self.assertEqual(p2.completed(include_partial=True), 5.) @@ -166,7 +173,8 @@ def test_table_ctor(self): t0 = astropy.time.Time('2020-01-01 07:00') for i, tile_id in enumerate(tiles): p1.add_exposure( - tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) p2 = Progress(p1._table) self.assertEqual(p2.completed(include_partial=True), 5.) self.assertEqual(p2.completed(include_partial=False), 0.) @@ -185,11 +193,14 @@ def test_completed_truncates(self): tile_id = p._table['tileid'][0] t0 = astropy.time.Time('2020-01-01 07:00') p.add_exposure( - tile_id, t0 + 1 * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + 1 * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) p.add_exposure( - tile_id, t0 + 2 * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + 2 * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) p.add_exposure( - tile_id, t0 + 3 * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + 3 * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) self.assertEqual(p.completed(include_partial=True), 1.) self.assertEqual(p.completed(include_partial=False), 1.) @@ -210,7 +221,8 @@ def test_completed_only_passes(self): tiles = p._table['tileid'][list(pass1[:n]) + list(pass7[:n])] t0 = astropy.time.Time('2020-01-01 07:00') for tile_id in tiles: - p.add_exposure(tile_id, t0, 1e3 * u.s, 1.5, 1.5, 1.1, 1, 0, 0, 0) + p.add_exposure(tile_id, t0, 1e3 * u.s, 1.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) t0 += 0.1 * u.day self.assertEqual(p.completed(only_passes=(7, 1)), 2 * n) self.assertEqual(p.completed(only_passes=7), n) @@ -230,10 +242,11 @@ def test_max_exposures(self): mjds = 58849. + np.arange(n) tt = astropy.time.Time('2020-01-01 07:00') + np.arange(n) * u.hour for t in tt[:-1]: - p.add_exposure(tile_id, t, 1e3 * u.s, 0.2, 1.5, 1.1, 1, 0, 0, 0) + p.add_exposure(tile_id, t, 1e3 * u.s, 0.2, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) with self.assertRaises(RuntimeError): p.add_exposure(tile_id, tt[-1], 1e3 * u.s, - 0.2, 1.5, 1.1, 1, 0, 0, 0) + 0.2, 1.5, 1.1, 1, 0, 0, 0, 0, 0, 0, 0) def test_summary(self): """Summary contains one row per tile""" @@ -247,9 +260,9 @@ def test_summary(self): for i, t in enumerate(p._table['tileid'][:n]): p.add_exposure( t, t0 + i * u.hour, 1e3 * u.s, 0.25, airmass, seeing, - transp, 0, 0, 0) + transp, 0, 0, 0, 0, 0, 0, 0) p.add_exposure(t, t0 + (i + 0.5) * u.hour, 1e3 * u.s, 0.25, - airmass, seeing, transp, 0, 0, 0) + airmass, seeing, transp, 0, 0, 0, 0, 0, 0, 0) self.assertEqual(len(p.get_summary('observed')), 100) self.assertEqual(len(p.get_summary('completed')), 0) self.assertTrue(np.all(p.get_summary('observed')['nexp'] == 2)) @@ -278,7 +291,8 @@ def test_copy_all(self): t0 = astropy.time.Time('2020-01-01 07:00') for i, tile_id in enumerate(tiles): p1.add_exposure( - tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t0 + i * u.hour, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) p2 = p1.copy_range() self.assertTrue(np.all(np.array(p1._table) == np.array(p2._table))) @@ -289,10 +303,12 @@ def test_copy_some(self): tiles = p1._table['tileid'][:n].data tt = astropy.time.Time('2020-01-01 07:00') + np.arange(n) * u.hour for t, tile_id in zip(tt, tiles): - p1.add_exposure(tile_id, t, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + p1.add_exposure(tile_id, t, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) for t, tile_id in zip(tt, tiles): p1.add_exposure( - tile_id, t + 100 * u.day, 1e3 * u.s, 0.5, 1.5, 1.1, 1, 0, 0, 0) + tile_id, t + 100 * u.day, 1e3 * u.s, 0.5, 1.5, 1.1, 1, + 0, 0, 0, 0, 0, 0, 0) self.assertEqual(p1.completed(), n) # Selects everything. mjd0 = tt[0].mjd From 971da45d247ac0c0722f341e30dee1df658d11f2 Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Thu, 30 Nov 2017 08:54:55 -0800 Subject: [PATCH 4/8] Allow empty tile_fields or exp_fields --- py/desisurvey/progress.py | 4 ++++ py/desisurvey/test/test_progress.py | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/py/desisurvey/progress.py b/py/desisurvey/progress.py index a3782c55..7e54e042 100644 --- a/py/desisurvey/progress.py +++ b/py/desisurvey/progress.py @@ -582,6 +582,8 @@ def get_exposures(self, start=None, stop=None, output = astropy.table.Table() output.meta['EXTNAME'] = 'EXPOSURES' for name in tile_fields.split(','): + if name == '': + continue # handle tile_fields='' name = name.lower() if name == 'index': output[name.upper()] = tile_index @@ -599,6 +601,8 @@ def get_exposures(self, start=None, stop=None, 'Invalid tile field name: {0}.'.format(name)) output[name.upper()] = table[name][tile_index] for name in exp_fields.split(','): + if name == '': + continue # handle exp_fields='' name = name.lower() if name == 'snr2cum': snr2cum = np.cumsum( diff --git a/py/desisurvey/test/test_progress.py b/py/desisurvey/test/test_progress.py index d022a683..19cefe63 100644 --- a/py/desisurvey/test/test_progress.py +++ b/py/desisurvey/test/test_progress.py @@ -120,6 +120,11 @@ def test_get_exposures(self): night = str(desisurvey.utils.get_date(row['MJD'])) self.assertEqual(night, str(desisurvey.utils.get_date(night))) + # Test empty tile_fields or exp_fields. + p.get_exposures(tile_fields='') + p.get_exposures(exp_fields='') + p.get_exposures(tile_fields='', exp_fields='') + #- Test roundtrip to disk expfile = os.path.join(self.tmpdir, 'test-exposures.fits') explist.write(expfile) From 0c80d690db05f374d28458aaf087fd9a8032a6a8 Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Thu, 30 Nov 2017 17:16:41 -0800 Subject: [PATCH 5/8] Implement twilight exposure model --- py/desisurvey/etc.py | 68 +++++++++++++++++++++++++++++----- py/desisurvey/test/test_etc.py | 17 ++++++++- 2 files changed, 75 insertions(+), 10 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index 79f932cb..6f6cdade 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -121,6 +121,62 @@ def airmass_exposure_factor(airmass): return np.power((X / X0), 1.25) +# Linear regression coefficients for converting scattered twilight r-band +# magnitude into an exposure-time correction factor. +_twilightCoefficients = np.array([ + -8.70618444805, 64323409.928312987, 739.2405272968557, + -18763.675456604633, 158813.46020567254]) + +def twilight_exposure_factor(sun_alt, sun_daz, airmass): + """Calculate exposure time factor due to scattered twilight. + + The returned factor is relative to dark conditions when the moon is + below the local horizon. + + This factor is based on a study of SNR for BGS threshold targets. + See ``surveysim:doc/nb/BGS_ETC_Study.ipynb`` for details. + + Parameters + ---------- + sun_alt : float + Altitude angle of the sun above the horizon in degrees. Must be in the + range [-90, -12] deg. Values below -18 deg are considered dark, with + no scattered sun. + sun_daz : float + Relative azimuth angle in degrees between the pointing and the sun. + airmass : float + Airmass used for observing this tile, must be >= 1. + + Returns + ------- + float + Dimensionless factor that exposure time should be increased to + account for increased sky brightness due to scattered twilight. + Will be 1 when the sun is below -18 degrees. + """ + if (sun_alt < -90) or (sun_alt > -12): + raise ValueError('Got invalid sun_alt outside [-90,-12].') + if airmass < 1: + raise ValueError('Got invalid airmass < 1.') + + # No exposure penalty when the sun is far enough below the horizon. + if sun_alt < -18: + return 1. + + # Estimate the altitude angle corresponding to this observing airmass. + # We invert eqn.3 of KS1991 for this (instead of eqn.14). + obs_zenith = np.arcsin(np.sqrt((1 - airmass ** -2) / 0.96)) * u.rad + obj_altitude = 90 * u.deg - obs_zenith.to(u.deg) + + # Calculate scattered twilight r-band brightness. + r = specsim.atmosphere.twilight_surface_brightness( + obj_altitude, sun_alt * u.deg, sun_daz * u.deg).value + + # Evaluate the linear regression model. + X = np.array((1, np.exp(-r), 1/r, 1/r**2, 1/r**3)) + return _twilightCoefficients.dot(X) + + # Linear regression coefficients for converting scattered moon V-band # magnitude into an exposure-time correction factor. _moonCoefficients = np.array([ @@ -140,14 +196,8 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass): This factor is based on a study of SNR for ELG targets and designed to achieve a median SNR of 7 for a typical ELG [OII] doublet at the lower flux limit of 8e-17 erg/(cm2 s A), averaged over the expected ELG target - redshift distribution 0.6 < z < 1.7. - - TODO: - - Check the assumption that exposure time scales with SNR ** -0.5. - - Check if this ELG-based analysis is also valid for BGS targets. - - For details, see the jupyter notebook doc/nb/ScatteredMoon.ipynb in - this package. + redshift distribution 0.6 < z < 1.7. See + ``surveysim:doc/nb/ScatteredMoon.ipynb`` for details. Parameters ---------- @@ -191,7 +241,7 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass): # We invert eqn.3 of KS1991 for this (instead of eqn.14). obs_zenith = np.arcsin(np.sqrt((1 - airmass ** -2) / 0.96)) * u.rad - # Calculate scattered moon V-band brightness at each pixel. + # Calculate scattered moon V-band brightness. V = specsim.atmosphere.krisciunas_schaefer( obs_zenith, moon_zenith, separation_angle, moon_phase, _vband_extinction).value diff --git a/py/desisurvey/test/test_etc.py b/py/desisurvey/test/test_etc.py index a47f8974..d7946509 100644 --- a/py/desisurvey/test/test_etc.py +++ b/py/desisurvey/test/test_etc.py @@ -4,7 +4,7 @@ import astropy.units as u -from desisurvey.etc import exposure_time, moon_exposure_factor +from desisurvey.etc import * class TestExpCalc(unittest.TestCase): @@ -44,6 +44,21 @@ def test_exptime(self): EBV, moon_frac, moon_sep, moon_alt) self.assertGreater(t2, t1) + def test_twilight(self): + # Sun below the horizon + x = twilight_exposure_factor(sun_alt=-20., sun_daz=0., airmass=1.2) + self.assertAlmostEqual(x, 1.0, 4) + + # Exposure times increase as the dawn sun rises, for many daz values. + for daz in (0, 90, 270, -10): + efac = np.empty(7) + for i, alt in enumerate(np.linspace(-18, -12, 7)): + efac[i] = twilight_exposure_factor( + sun_alt=alt, sun_daz=daz, airmass=1.) + self.assertTrue(np.all(efac >= 1.0)) + self.assertTrue(np.all(efac < 1.2)) + self.assertTrue(np.all(np.diff(efac) > 0)) + def test_moon(self): #- Moon below horizon x = moon_exposure_factor( From 9a4f9b57c3b8dd73b371dbb266df2f2b66ebc5fe Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Thu, 30 Nov 2017 17:29:35 -0800 Subject: [PATCH 6/8] Update scattered moon model coefs --- py/desisurvey/etc.py | 10 ++++++---- py/desisurvey/test/test_etc.py | 5 ++--- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index 6f6cdade..5acae377 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -180,8 +180,8 @@ def twilight_exposure_factor(sun_alt, sun_daz, airmass): # Linear regression coefficients for converting scattered moon V-band # magnitude into an exposure-time correction factor. _moonCoefficients = np.array([ - -8.83964463188, -7372368.5041596508, 775.17763895781638, - -20185.959363990656, 174143.69095766739]) + -9.6827757062, -5454732.2391116023, 805.44541757440129, + -20274.598621708101, 170436.98267100245]) # V-band extinction coefficient to use in the scattered moonlight model. # See specsim.atmosphere.krisciunas_schaefer for details. @@ -247,8 +247,10 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass): moon_phase, _vband_extinction).value # Evaluate the linear regression model. - X = np.array((1, np.exp(-V), 1/V, 1/V**2, 1/V**3)) - return _moonCoefficients.dot(X) + X = np.array((1., np.exp(-V), 1/V, 1/V**2, 1/V**3)) + + # Clip result to 1 since fits goes just below 1 at the faint end (V>26). + return np.maximum(1., _moonCoefficients.dot(X)) def exposure_time(program, seeing, transparency, airmass, EBV, diff --git a/py/desisurvey/test/test_etc.py b/py/desisurvey/test/test_etc.py index d7946509..7890afae 100644 --- a/py/desisurvey/test/test_etc.py +++ b/py/desisurvey/test/test_etc.py @@ -65,11 +65,10 @@ def test_moon(self): moon_frac=0.5, moon_sep=60, moon_alt=-30, airmass=1.2) self.assertAlmostEqual(x, 1.0, 4) - #- New moon above horizon still has some impact + #- New moon above horizon still no impact x = moon_exposure_factor( moon_frac=0.0, moon_sep=60, moon_alt=30, airmass=1.2) - self.assertGreater(x, 1.0) - self.assertLess(x, 1.1) + self.assertAlmostEqual(x, 1.0, 4) #- Partial moon takes more time x1 = moon_exposure_factor( From c3aa9aea09299227e210c06cbd007569ff64202a Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Sat, 2 Dec 2017 06:58:49 -0800 Subject: [PATCH 7/8] Use uppercase colnames in surveymovie --- py/desisurvey/scripts/surveymovie.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/py/desisurvey/scripts/surveymovie.py b/py/desisurvey/scripts/surveymovie.py index c66f9d9b..c5e0e9e2 100644 --- a/py/desisurvey/scripts/surveymovie.py +++ b/py/desisurvey/scripts/surveymovie.py @@ -152,15 +152,15 @@ def __init__(self, ephem, progress, start, stop, label, show_scores): start, stop, tile_fields='tileid,index,ra,dec,pass', exp_fields='expid,mjd,night,exptime,snr2cum,seeing,transparency') self.num_exp = len(self.exposures) - self.num_nights = len(np.unique(self.exposures['night'])) + self.num_nights = len(np.unique(self.exposures['NIGHT'])) # Calculate each exposure's LST window. exp_midpt = astropy.time.Time( - self.exposures['mjd'] + self.exposures['exptime'] / 86400., + self.exposures['MJD'] + self.exposures['EXPTIME'] / 86400., format='mjd', location=desisurvey.utils.get_location()) lst_midpt = exp_midpt.sidereal_time('apparent').to(u.deg).value # convert from seconds to degrees. - lst_len = self.exposures['exptime'] / 240. + lst_len = self.exposures['EXPTIME'] / 240. self.lst = np.empty((self.num_exp, 2)) self.lst[:, 0] = wrap(lst_midpt - 0.5 * lst_len) self.lst[:, 1] = wrap(lst_midpt + 0.5 * lst_len) @@ -326,7 +326,7 @@ def init_date(self, date, ephem): hdus.close() # Save index of first exposure on this date. noon = desisurvey.utils.local_noon_on_date(date) - self.iexp0 = np.argmax(self.exposures['mjd'] > noon.mjd) + self.iexp0 = np.argmax(self.exposures['MJD'] > noon.mjd) else: self.warn('Missing scores file: {0}.'.format(scores_name)) # Get interpolator for moon, planet positions during this night. From 3d5d7a3fc7f2555ff83d9b069845916ce567b938 Mon Sep 17 00:00:00 2001 From: David Kirkby Date: Sat, 2 Dec 2017 07:06:41 -0800 Subject: [PATCH 8/8] More colname ucasing --- py/desisurvey/scripts/surveymovie.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/py/desisurvey/scripts/surveymovie.py b/py/desisurvey/scripts/surveymovie.py index c5e0e9e2..a885fcc6 100644 --- a/py/desisurvey/scripts/surveymovie.py +++ b/py/desisurvey/scripts/surveymovie.py @@ -367,9 +367,9 @@ def draw_exposure(self, iexp, nightly): True if a new frame was drawn. """ info = self.exposures[iexp] - mjd = info['mjd'] + mjd = info['MJD'] date = desisurvey.utils.get_date(mjd) - assert date == desisurvey.utils.get_date(info['night']) + assert date == desisurvey.utils.get_date(info['NIGHT']) night = self.ephem.get_night(date) # Initialize status if necessary. if (self.status is None) or nightly: @@ -381,18 +381,18 @@ def draw_exposure(self, iexp, nightly): elif nightly: return False # Update the status for the current exposure. - complete = info['snr2cum'] >= self.config.min_snr2_fraction() - self.status[info['index']] = 2 if complete else 1 + complete = info['SNR2CUM'] >= self.config.min_snr2_fraction() + self.status[info['INDEX']] = 2 if complete else 1 # Update the top-right label. - label = '{} {} #{:06d}'.format(self.label, date, info['expid']) + label = '{} {} #{:06d}'.format(self.label, date, info['EXPID']) if not nightly: label += ' ({:.1f}",{:.2f})'.format( - info['seeing'], info['transparency']) + info['SEEING'], info['TRANSPARENCY']) self.text.set_text(label) if not nightly: # Update current time in program. dt1 = mjd - night['noon'] - dt2 = dt1 + info['exptime'] / 86400. + dt2 = dt1 + info['EXPTIME'] / 86400. self.pline1.set_xdata([dt1, dt1]) self.pline2.set_xdata([dt2, dt2]) if self.show_scores: @@ -415,9 +415,9 @@ def draw_exposure(self, iexp, nightly): sizes[done] = 30. fc[done] = self.completecolor fc[~avail] = self.unavailcolor - if not nightly and (info['pass'] == passnum): + if not nightly and (info['PASS'] == passnum): # Highlight the tile being observed now. - jdx = np.where(self.tileid[sel] == info['tileid'])[0][0] + jdx = np.where(self.tileid[sel] == info['TILEID'])[0][0] fc[jdx] = self.nowcolor scatter.get_sizes()[jdx] = 600. scatter.set_facecolors(fc) @@ -430,7 +430,7 @@ def draw_exposure(self, iexp, nightly): # Update LST lines. x1, x2 = self.lst[iexp] for passnum, (line1, line2) in enumerate(self.lstlines): - ls = '-' if info['pass'] == passnum else '--' + ls = '-' if info['PASS'] == passnum else '--' line1.set_linestyle(ls) line2.set_linestyle(ls) line1.set_xdata([x1, x1]) @@ -485,7 +485,7 @@ def main(args): animator.init_figure(args.nightly) if args.expid is not None: - expid = animator.exposures['expid'] + expid = animator.exposures['EXPID'] assert np.all(expid == expid[0] + np.arange(len(expid))) if (args.expid < expid[0]) or (args.expid > expid[-1]): raise RuntimeError('Requested exposure ID {0} not available.'