Merging recent updates to dev branch

potpyri/instruments/MOSFIRE.py

@@ -33,7 +33,7 @@ def __init__(self):
         self.min_exptime = 1.0
 
         # Run dark/bias/flat calibration?
-        self.dark = False
+        self.dark = True
         self.bias = False
         self.flat = True
 

potpyri/primitives/absphot.py

@@ -107,20 +107,44 @@ def get_catalog(self, coords, catalog, filt, log=None):
 
         seps = med_coord.separation(coords)
         max_sep = np.max(seps.to(u.deg).value)
+
+        cols = [cat_ra, cat_dec, cat_mag, cat_err]
+        # Add Kron mag if the catalog is PS1
+        if cat_ID=='II/349':
+            cols.append(f'{filt}Kmag')
 
-        vizier = Vizier(columns=[cat_ra, cat_dec, cat_mag, cat_err])
+        vizier = Vizier(columns=cols)
         vizier.ROW_LIMIT = -1
         if log: 
             log.info(f'Getting {catalog} catalog with ID {cat_ID} in filt {filt}')
             log.info(f'Querying around {coord_ra}, {coord_dec} deg')
-        cat = vizier.query_region(med_coord, width=1.2*max_sep*u.degree, 
+        Vizier.clear_cache()
+        width = np.max([2.0 * max_sep, 0.5])
+        cat = vizier.query_region(med_coord, width=width*u.degree, 
             catalog=cat_ID)
 
         if len(cat)>0:
             cat = cat[0]
             cat = cat[~np.isnan(cat[cat_mag])]
             cat = cat[cat[cat_err]>0.]
 
+            if cat_ID=='II/349':
+                if log:
+                    log.info('Cutting on Kron magnitudes')
+                else:
+                    print('Cutting on Kron magnitudes')
+
+                nsources = len(cat)
+                cat_kron = f'{filt}Kmag'
+                mask = cat[cat_mag]-cat[cat_kron] < 0.1
+                cat = cat[mask]
+                nkron = len(cat)
+
+                if log:
+                    log.info(f'Cut catalog from {nsources} to {nkron}')
+                else:
+                    print(f'Cut catalog from {nsources} to {nkron}')
+
             cat.rename_column(cat_ra, 'ra')
             cat.rename_column(cat_dec, 'dec')
             cat.rename_column(cat_mag, 'mag')
@@ -182,10 +206,10 @@ def find_zeropoint(self, cmpfile, tel, match_radius=2.5*u.arcsec,
 
             cat, catalog, cat_ID = self.get_catalog(coords, catalog, filt, log=log)
 
-        min_mag = self.get_minmag(filt)
-        cat = cat[cat['mag']>min_mag]
-
         if cat:
+            min_mag = self.get_minmag(filt)
+            cat = cat[cat['mag']>min_mag]
+
             coords_cat = SkyCoord(cat['ra'], cat['dec'], unit='deg')
 
             idx, d2, d3 = coords_cat.match_to_catalog_sky(coords)

potpyri/primitives/image_procs.py

@@ -27,6 +27,7 @@
 from ccdproc import cosmicray_lacosmic
 
 # Internal dependencies
+from potpyri.primitives import photometry
 from potpyri.primitives import solve_wcs
 from potpyri.utils import utilities
 
@@ -138,8 +139,11 @@ def align_images(reduced_files, paths, tel, binn, use_wcs=None, fieldcenter=None
 
             if log:
                 log.info('Rejecting the following images for high astrometric dispersion:')
-                for i,m in enumerate(mask):
-                    if not m: log.info(solved_images[i])
+                if np.all(mask):
+                    log.info('No images rejected')
+                else:
+                    for i,m in enumerate(mask):
+                        if not m: log.info(solved_images[i])
 
             solved_images = np.array(solved_images)[mask]
     else:
@@ -451,72 +455,20 @@ def stack_data(stacking_data, tel, masks, errors, mem_limit=8.0e9, log=None):
         else:
             raise Exception('Could not get stacking method for these images')
 
-    weights = []
-    for i,error in enumerate(errors):
-        ivar = 1./error.data**2
-        mask = masks[i]
-        ivar[mask.data.astype(bool)]=0.
-        weights.append(ivar)
-    weights=np.array(weights)
-
-    new_data = []
     exptimes = []
     for i,stk in enumerate(stacking_data):
         mask = masks[i].data.astype(bool)
         stacking_data[i].data[mask] = np.nan
         exptimes.append(float(tel.get_exptime(stk.header)))
-
+    
     exptimes = np.array(exptimes)
     scale = 1./exptimes
 
-    all_data = np.array([s.data for s in stacking_data])
-
     if len(scale)==1:
         stack_method='average'
 
-    # Determine if intermediate stacks are needed
-    # Get size of individual frame, account for data, noise, and mask
-    exdata = stacking_data[0].data
-    size = 64.0/8.0 * exdata.shape[0]*exdata.shape[1] * 2 + exdata.shape[0]*exdata.shape[1]
-
-    if log:
-        log.info(f'Image size is {size}')
-    else:
-        print(f'Image size is {size}')
-
-    # Get maximum size of chunk
-    max_chunk = int(np.floor(mem_limit / 4.0 / (size)))
-    if max_chunk==0: max_chunk=1
-
-    if len(stacking_data)<=max_chunk:
-        sci_med = combine(stacking_data, weights=weights, scale=scale,
+    sci_med = combine(stacking_data, scale=scale,
             method=stack_method, mem_limit=mem_limit)
-    else:
-        nimgs = len(stacking_data) * 1.0
-        nchunks = int(np.ceil(nimgs/(1.0*max_chunk)))
-
-        if log:
-            log.info(f'Splitting stacking into {nchunks} chunks')
-        else:
-            print(f'Splitting stacking into {nchunks} chunks')
-
-        inter_med = []
-        for i in np.arange(nchunks):
-            if log:
-                log.info(f'Stacking chunk {i+1}/{nchunks}')
-            else:
-                print(f'Stacking chunk {i+1}/{nchunks}')
-            start = i*max_chunk
-            end = (i+1)*max_chunk
-            if end>len(stacking_data): end=None
-
-            sci_med = combine(stacking_data[start:end], 
-                weights=weights[start:end], scale=scale[start:end],
-                method=stack_method, mem_limit=mem_limit)
-
-            inter_med.append(sci_med)
-
-        sci_med = combine(inter_med, method=stack_method)
 
     sci_med = sci_med.to_hdu()
 
@@ -602,9 +554,9 @@ def mask_satellites(images, filenames, log=None):
 
         os.remove(tmpfile)
 
-def create_mask(science_data, saturation, rdnoise, sigclip=3.5, 
-    sigfrac=0.2, objlim=4.5, niter=6, outpath='', grow=0, satellites=True,
-    cosmic_ray=True, log=None):
+def create_mask(science_data, saturation, rdnoise, sigclip=3.0, 
+    sigfrac=0.10, objlim=4.0, niter=6, outpath='', grow=0, cosmic_ray=True, 
+    fsmode='convolve', cleantype='medmask', log=None):
 
     t_start = time.time()
 
@@ -616,6 +568,18 @@ def create_mask(science_data, saturation, rdnoise, sigclip=3.5,
     hdu = fits.open(science_data)
     data = np.ascontiguousarray(hdu[0].data.astype(np.float32))
 
+    # Estimate FWHM size for LACosmic
+    table = photometry.run_sextractor(science_data, log=log)
+    if table is not None:
+        # Clip fwhm_stars by fwhm
+        fwhm, meanfwhm, stdfwhm = sigma_clipped_stats(table['FWHM_IMAGE'])
+        if log: 
+            log.info(f'Using FWHM for cosmic rays: {fwhm}')
+        else:
+            print(f'Using FWHM for cosmic rays: {fwhm}')
+    else:
+        fwhm = 3.5
+
     # Astroscrappy requires added sky background, so add this value back
     # Set the sky background to some nominal value if it is too low for CR rej
     skybkg = hdu[0].header['SKYBKG']
@@ -656,10 +620,15 @@ def create_mask(science_data, saturation, rdnoise, sigclip=3.5,
     mask_cr = np.zeros(data.shape)
     mask_cr = mask_cr.astype(np.uint8)
 
+    psfsize = int(np.round(2.5*fwhm))
+    if psfsize%2==0: psfsize+=1
+
     if cosmic_ray:
         newdata, mask_cr = cosmicray_lacosmic(data,
             readnoise=rdnoise, satlevel=saturation, verbose=True,
-            sigclip=sigclip, sigfrac=sigfrac, objlim=objlim, niter=niter)
+            sigclip=sigclip, sigfrac=sigfrac, objlim=objlim, niter=niter,
+            psffwhm=fwhm, psfsize=psfsize, fsmode=fsmode,
+            cleantype=cleantype)
 
         mask_cr = mask_cr.astype(np.uint8) #set cosmic ray mask type
     else:
@@ -692,7 +661,7 @@ def create_mask(science_data, saturation, rdnoise, sigclip=3.5,
     nsat = np.sum(mask & 4 == 4)
     ngrow = np.sum(mask & 8 == 8)
 
-    mask_hdu = fits.PrimaryHDU(mask) #create mask Primary HDU
+    mask_hdu = fits.ImageHDU(mask) #create mask Primary HDU
     mask_hdu.header['VER'] = (__version__, 
         'Version of image procedures used.') #name of log file
     mask_hdu.header['USE'] = 'Complex mask using additive flags.'#header comment

potpyri/primitives/photometry.py

@@ -520,7 +520,7 @@ def photloop(stack, phot_sn_min=3.0, phot_sn_max=40.0, fwhm_init=5.0, log=None):
     epsf = None ; fwhm = None
     while signal_to_noise > phot_sn_min:
 
-        if log: log.info(f'Trying photometry with final S/N={signal_to_noise}')
+        if log: log.info(f'Trying PSF generation with S/N={signal_to_noise}')
         star_param = {'snthresh_psf': signal_to_noise*2.0,
                       'fwhm_init': fwhm_init,
                       'snthresh_final': signal_to_noise}

potpyri/primitives/solve_wcs.py

@@ -335,7 +335,7 @@ def solve_astrometry(file, tel, binn, paths, radius=0.5, replace=True,
         return(False)
 
 def align_to_gaia(file, tel, radius=0.5, max_search_radius=5.0*u.arcsec,
-    save_centroids=False, log=None):
+    save_centroids=False, min_gaia_match=5, log=None):
 
     cat = get_gaia_catalog(file, log=log)
 
@@ -379,7 +379,7 @@ def align_to_gaia(file, tel, radius=0.5, max_search_radius=5.0*u.arcsec,
     else:
         print(f'Found {len(cat)} stars in the image')
 
-    if len(cat)<7:
+    if len(cat)<min_gaia_match:
         if log:
             log.info(f'Too few stars in {file}.  Skipping Gaia alignment...')
         else:
@@ -472,6 +472,20 @@ def align_to_gaia(file, tel, radius=0.5, max_search_radius=5.0*u.arcsec,
         else:
             print(f'Matched to {len(cat_coords_match)} coordinates')
 
+    if len(cat_coords_match)<min_gaia_match:
+        if log:
+            log.info(f'Too few matched stars in {file}.  Skipping Gaia alignment...')
+        else:
+            print(f'Too few matched stars in {file}.  Skipping Gaia alignment...')
+
+        # Set these values to some default so pipeline won't break
+        hdu[0].header['RADISP']=(1.0, 'Dispersion in R.A. of WCS [Arcsec]')
+        hdu[0].header['DEDISP']=(1.0, 'Dispersion in Decl. of WCS [Arcsec]')
+
+        hdu.writeto(file, overwrite=True, output_verify='silentfix')
+
+        return(True)
+
     # Bootstrap WCS solution and get center pixel
     ra_cent = [] ; de_cent = []
     coords = SkyCoord(cat_coords_match['RA_ICRS'], cat_coords_match['DE_ICRS'], 

potpyri/primitives/sort_files.py

@@ -282,14 +282,14 @@ def sort_files(files, file_list, tel, paths, incl_bad=False, log=None):
                 file_type = 'BIAS'
                 moved_path = paths['raw']
                 bias_num += 1
-            elif is_science(hdr, tel):
-                file_type = 'SCIENCE'
-                moved_path = paths['raw']
-                sci_num += 1
             elif is_dark(hdr, tel):
                 file_type = 'DARK'
                 moved_path = paths['raw']
                 dark_num += 1
+            elif is_science(hdr, tel):
+                file_type = 'SCIENCE'
+                moved_path = paths['raw']
+                sci_num += 1
             else:
                 file_type = 'BAD'
                 moved_path = paths['bad']

potpyri/utils/options.py

@@ -53,7 +53,7 @@ def init_options():
         help='''Minimum signal-to-noise to try in photometry loop.''')
     params.add_argument('--phot-sn-max',
         type=float,
-        default=40.0,
+        default=20.0,
         help='''Maximum signal-to-noise to try in photometry loop.''')
     params.add_argument('--fwhm-init',
         type=float,