Merge pull request #8 from kadrlica/obztak

update from obztak

bin/plot_json

@@ -3,8 +3,8 @@
 Plot json files.
 """
 import argparse
-import pylab as plt
 import os
+import pylab as plt
 import numpy as np
 
 from obztak import get_survey
@@ -25,12 +25,18 @@ if __name__ == "__main__":
                         help='plot chunks of exposures file-by-file.')
     parser.add_argument('-a','--all',action='store_true',
                         help='plot all exposures at once.')
+    parser.add_argument('--airmass',default=None,type=float,
+                        help='Airmass to plot.')
     parser.add_argument('-c','--complete',const=None,nargs='?',action='append',
                         help="fields that have been completed.")
     parser.add_argument('-o','--outfile',
                         help='output file to save')
     args = parser.parse_args()
 
+    # Set the backend for output (doesn't really work here...)
+    #if args.outfile is not None: matplotlib.use('Agg')
+    survey = get_survey()
+
     ext = os.path.splitext(args.outfile)[-1] if args.outfile else None
 
     idx = []
@@ -39,7 +45,7 @@ if __name__ == "__main__":
         fields = fields + FieldArray.read(filename)
         if args.chunk: idx.append(len(fields)-1)
 
-    completed_fields = field_factory(get_survey())
+    completed_fields = field_factory(survey)
     if args.complete is not None:
         for f in args.complete:
             if f is None:
@@ -49,7 +55,14 @@ if __name__ == "__main__":
 
     movie = True if ext in MOVIES else False
 
-    options = dict(airmass=1.4,smash=False)
+    options = dict(airmass=args.airmass,smash=False)
+
+    if args.airmass is not None:
+        options['airmass'] = args.airmass
+    elif survey in ['maglites','pals']:
+        options['airmass'] = 2.0
+    else:
+        options['airmass'] = 1.4
 
     if args.all:
         plotFields(fields[-1],target_fields=fields,completed_fields=completed_fields+fields,options_basemap=options)

bin/plot_nightsum

@@ -31,7 +31,7 @@ if __name__ == "__main__":
     if args.nite:
         date = nite2utc(args.nite)
     else:
-        # Yesterday...
+        # Yesterday...? Would 12*ephem.hour be better?
         date = ephem.Date(ephem.now() - 8*ephem.hour)
     nitestr = utc2nite(date)
 
@@ -49,16 +49,28 @@ if __name__ == "__main__":
     elif survey == 'bliss':
 
         plot_bliss_coverage(fields)
+        plt.suptitle('Coverage (%s)'%nitestr,fontsize=18)
         plt.savefig('nightsum_coverage_%s.png'%nitestr)
 
-        plt.figure()
-        bmap = DECamOrtho(date='2017/02/08 06:00:00')
+        fig,axes = plt.subplots(1,2,figsize=(12,5))
+        plt.sca(axes[0])
+        bmap = DECamOrtho(date='2017/02/08 07:00:00')
+        for b in np.unique(fields['FILTER']):
+            f = fields[fields['FILTER']==b]
+            bmap.draw_focal_planes(f['RA'],f['DEC'],color=COLORS[b],alpha=0.3)
+        bmap.draw_bliss()
+        bmap.draw_galaxy()
+        bmap.draw_des()
+
+        plt.sca(axes[1])
+        bmap = DECamOrtho(date='2017/02/08 19:00:00')
         for b in np.unique(fields['FILTER']):
             f = fields[fields['FILTER']==b]
             bmap.draw_focal_planes(f['RA'],f['DEC'],color=COLORS[b],alpha=0.3)
         bmap.draw_bliss()
         bmap.draw_galaxy()
         bmap.draw_des()
+        plt.suptitle('Coverage (%s)'%nitestr,fontsize=16)
         plt.savefig('nightsum_summary_%s.png'%nitestr)
 
         new = (np.array(map(utc2nite,fields['DATE'])) == nitestr)

obztak/bliss.py

@@ -25,12 +25,12 @@
 PROPID  = '2017A-0260'
 PROPOSER = 'Soares-Santos'
 BANDS = ['g','r','i','z']
-TILINGS = 2
+TILINGS = 4
 
 class BlissSurvey(Survey):
     """ Survey sublcass for BLISS. """
 
-    # 2017A ACTUAL
+    # 2017A SCHEDULED (bliss-windows.csv is actually used)
     nights = [
         ['2017/02/07', 'second'], # phase=90%, set=07:40 (1.5h dark)
         ['2017/02/08', 'second'], # phase=95%, set=08:40 (0.5h dark)
@@ -54,20 +54,34 @@ class BlissSurvey(Survey):
         ['2017/07/15', 'second'], # phase=57%, set=nope
         ]
 
+    extra_night = [
+        ['2017/08/05', 'second'], # phase=99%, set=06:40 (4h dark)
+        ['2017/08/06', 'full'],   # phase=100%, set=nope
+        ['2017/08/07', 'full'],   # phase=99%, set=nope
+    ]
+
     alfredo_nights = [
         ['2017/03/06', 'full'],  # level=9
         ['2017/03/07','second'], # level=10
         ['2017/06/21','first'],
         ]
 
     def prepare_fields(self, infile=None, outfile=None, mode='bliss_rotate', plot=True, smcnod=False):
-        """ Create the list of fields to be targeted by this survey.
+        """Create the list of fields to be targeted by the BLISS survey.
+
+        Selection starts with 3 regions:
+        - P9 - Planet 9 region above DES footprint (priority 1)
+        - LIGO - Region targeted based on LIGO sensitivity maps (and good for MW)
+        - Alfredo - Overlap with Alfredo's eRosita companion survey.
+
+        Fields that have been previously covered by DECam are removed
+        from the LIGO and Alfredo fooptrint regions.
 
         Parameters:
         -----------
         infile : File containing all possible field locations.
         outfile: Output file of selected fields
-        mode   : Mode for dithering: 'smash_dither', 'smash_rotate', 'decam_dither', 'none'
+        mode   : Mode for dithering. default: 'bliss_rotate'
         plot   : Create an output plot of selected fields.
 
         Returns:
@@ -160,22 +174,26 @@ def dither(ra,dec,dx,dy):
         p9 = self.planet9(fields['RA'],fields['DEC'])
         ligo = self.ligo_mw(fields['RA'],fields['DEC'])
         alfredo = self.alfredo(fields['RA'],fields['DEC'])
+        p9v2 = self.planet9v2(fields['RA'],fields['DEC'])
 
         fields['PRIORITY'][p9] = 1
         fields['PRIORITY'][ligo] = 2
         fields['PRIORITY'][alfredo] = 3
 
-        sel = (p9 | ligo | alfredo)
+        #sel = (p9 | ligo | alfredo)
+        sel = (p9 | ligo | alfredo | p9v2)
 
         # Apply telescope constraints
         sel &= (fields['DEC'] > constants.SOUTHERN_REACH)
 
         # Apply covered fields
-        sel &= self.uncovered(fields['RA'],fields['DEC'],fields['FILTER'])[0]
+        #sel &= self.uncovered(fields['RA'],fields['DEC'],fields['FILTER'])[0]
+        # Apply covered fields (but take everything in P9 region)
+        uncovered = self.uncovered(fields['RA'],fields['DEC'],fields['FILTER'])[0]
+        sel &= ( (p9v2) | (p9 & (fields['RA'] > 180)) | uncovered )
 
         fields = fields[sel]
 
-
         logging.info("Number of target fields: %d"%len(fields))
         logging.debug("Unique priorities: ",np.unique(fields['PRIORITY']))
 
@@ -219,6 +237,9 @@ def dither(ra,dec,dx,dy):
             if not sys.flags.interactive:
                 plt.show(block=True)
 
+        # Just 3rd tiling
+        #fields = fields[fields['TILING'] == 3]
+
         if outfile: fields.write(outfile)
 
         return fields
@@ -264,12 +285,28 @@ def ligo_mw(ra,dec):
     @staticmethod
     def planet9(ra,dec):
         """The high-probability region for Planet 9"""
-        sel = ((ra > 305)|(ra < 15)) & (dec > -40) & (dec < -30) # v0
+        sel = ((ra > 305)|(ra < 15)) & (dec > -40) & (dec < -30) # v0, v4
         #sel = ((ra > 305)|(ra < 15)) & (dec > -40) & (dec < -25) # v1
         #sel = ((ra > 305)|(ra < 15)) & (dec > -35) & (dec < -25) # v2
         #sel = ((ra > 305)|(ra < 15)) & (dec > -40) & (dec < -28) # v3
         return sel
 
+    @staticmethod
+    def planet9v2(ra,dec):
+        """The other high-probability region for Planet 9"""
+        from matplotlib.path import Path
+        data = np.genfromtxt(fileio.get_datafile('blissII-poly.txt'),
+                             names=['RA','DEC','POLY'])
+        poly = data[data['POLY'] == 1]
+        path = Path(zip(poly['RA'],poly['DEC']))
+        sel = path.contains_points(np.vstack([ra,dec]).T)
+
+        poly = data[data['POLY'] == 4]
+        path = Path(zip(poly['RA'],poly['DEC']))
+        sel |= path.contains_points(np.vstack([ra,dec]).T)
+
+        return sel
+
     @staticmethod
     def alfredo(ra,dec):
         """Alfredo's eRosita survey extension"""
@@ -283,6 +320,19 @@ def alfredo(ra,dec):
 
     @staticmethod
     def uncovered(ra,dec,band):
+        """
+        Determine which fields haven't been previously covered by DECam
+
+        Parameters:
+        -----------
+        ra  : right ascension of field
+        dec : declination of field
+        band: observing band
+
+        Returns:
+        --------
+        sel, frac : selection of fields and coverage fraction
+        """
         import healpy as hp
         #dirname = '/home/s1/kadrlica/projects/bliss/v0/data'
         #basename = 'decam_coverage_90s_%s_n1024.fits.gz'
@@ -298,11 +348,11 @@ def uncovered(ra,dec,band):
             idx = (band==b)
             filename1 = os.path.join(dirname,basename1%b)
             logging.info("Reading %s..."%os.path.basename(filename1))
-            skymap1 = hp.read_map(filename1,verbose=True)
+            skymap1 = hp.read_map(filename1,verbose=False)
 
             filename2 = os.path.join(dirname,basename2%b)
             logging.info("Reading %s..."%os.path.basename(filename2))
-            skymap2 = hp.read_map(filename2,verbose=True)
+            skymap2 = hp.read_map(filename2,verbose=False)
 
             # Apply minimum exposure time selection
             #skymap = (skymap > 30)
@@ -367,6 +417,12 @@ def query(cls, **kwargs):
         FROM exposure where propid = '%(propid)s' and exptime > 89
         and discard = False and delivered = True and flavor = 'object'
         and object like '%(object_fmt)s%%'
+        -- i-band AOS failures 'sqrt(pow(qc_fwhm,2)-pow(dimm2see,2)) > 0.9'
+        and id not in (652771,652794,652795,652796,652797,652799,652800,652803,652804,652806,652807,652808,652809,652810,652812,652813,652814,652815,652816,652817,652818,652819,652820,652821,652822,652823,652824,652825,652826,652827,652828,652829,652830,652831,652832,652833,652834,652835,652836,652837,652838,652839)
+        -- z-band AOS failers 'sqrt(pow(qc_fwhm,2)-pow(dimm2see,2)) > 0.7'
+        and id not in (652692,652693,652694,652695,652702,652703,652704,652705,652706,652707,652709,652752,652753,652760,652762,652763,652764,652765,652773,652774,652775,652776,652777,652781,652782,652784,652785,652786,652787,652788,652789,652790,652791,652792,652801,652802,652811)
+        -- t_eff values (careful about nulls)
+        and not (qc_teff < 0.1 and date < '2017/08/07 12:00:00')
         ORDER BY utc_beg %(limit)s
         """%kwargs
         return query