diff --git a/potpyri/instruments/MOSFIRE.py b/potpyri/instruments/MOSFIRE.py index 6d0f8b1..88f7170 100755 --- a/potpyri/instruments/MOSFIRE.py +++ b/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 diff --git a/potpyri/primitives/absphot.py b/potpyri/primitives/absphot.py index 6b32084..5430ae9 100755 --- a/potpyri/primitives/absphot.py +++ b/potpyri/primitives/absphot.py @@ -107,13 +107,20 @@ 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: @@ -121,6 +128,23 @@ def get_catalog(self, coords, catalog, filt, log=None): 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) diff --git a/potpyri/primitives/image_procs.py b/potpyri/primitives/image_procs.py index 18481f7..27c49fc 100755 --- a/potpyri/primitives/image_procs.py +++ b/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 diff --git a/potpyri/primitives/photometry.py b/potpyri/primitives/photometry.py index 23f91bd..8695391 100755 --- a/potpyri/primitives/photometry.py +++ b/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} diff --git a/potpyri/primitives/solve_wcs.py b/potpyri/primitives/solve_wcs.py index 89f7502..ab387e9 100755 --- a/potpyri/primitives/solve_wcs.py +++ b/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)