Clean references

.gitignore

@@ -130,3 +130,6 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# JetBrains IDEs
+.idea/

flows/load_image.py

@@ -55,6 +55,7 @@ def load_image(FILENAME):
 
 		# Specific headers:
 		image.exptime = float(hdr['EXPTIME']) # * u.second
+		image.peakmax = None # Maximum value above which data is not to be trusted
 
 		# Timestamp:
 		if origin == 'LCOGT':
@@ -67,6 +68,11 @@ def load_image(FILENAME):
 
 			image.photfilter = hdr['FILTER']
 
+			# Get non-linear limit
+			# TODO: Use actual or some fraction of the non-linearity limit
+			#image.peakmax = hdr.get('MAXLIN') # Presumed non-linearity limit from header
+			image.peakmax = 60000 #From experience, this one is better.
+
 		elif hdr.get('TELESCOP') == 'NOT' and hdr.get('INSTRUME') in ('ALFOSC FASU', 'ALFOSC_FASU') and hdr.get('OBS_MODE') == 'IMAGING':
 			image.site = api.get_site(5) # Hard-coded the siteid for NOT
 			image.obstime = Time(hdr['DATE-AVG'], format='isot', scale='utc', location=image.site['EarthLocation'])
@@ -80,6 +86,11 @@ def load_image(FILENAME):
 				'u SDSS': 'up'
 			}.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER'])
 
+			# Get non-linear limit
+			# Obtained from http://www.not.iac.es/instruments/detectors/CCD14/LED-linearity/20181026-200-1x1.pdf
+			# TODO: grab these from a table for all detector setups of ALFOSC
+			image.peakmax = 80000 # For ALFOSC D, 1x1, 200; the standard for SNe.
+
 		elif hdr.get('FPA.TELESCOPE') == 'PS1' and hdr.get('FPA.INSTRUMENT') == 'GPC1':
 			image.site = api.get_site(6) # Hard-coded the siteid for Pan-STARRS1
 			image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=image.site['EarthLocation'])

flows/photometry.py

@@ -54,7 +54,7 @@ def photometry(fileid, output_folder=None):
 
 	# Settings:
 	ref_mag_limit = 22 # Lower limit on reference target brightness
-	ref_target_dist_limit = 30 # Reference star must be further than this away to be included
+	ref_target_dist_limit = 10 * u.arcsec # Reference star must be further than this away to be included
 
 	logger = logging.getLogger(__name__)
 	tic = default_timer()
@@ -130,7 +130,8 @@ def photometry(fileid, output_folder=None):
 	hsize = 10
 	x = references['pixel_column']
 	y = references['pixel_row']
-	references = references[(np.sqrt((x - target_pixel_pos[0])**2 + (y - target_pixel_pos[1])**2) > ref_target_dist_limit)
+	refs_coord = coords.SkyCoord(ra=references['ra'], dec=references['decl'], unit='deg', frame='icrs')
+	references = references[(target_coord.separation(refs_coord) > ref_target_dist_limit)
 		& (x > hsize) & (x < (image.shape[1] - 1 - hsize))
 		& (y > hsize) & (y < (image.shape[0] - 1 - hsize))]
 	# 		& (references[ref_filter] < ref_mag_limit)
@@ -161,6 +162,7 @@ def photometry(fileid, output_folder=None):
 		exclude_percentile=50.0
 		)
 	image.background = bkg.background
+	image.std = bkg.background_rms_median
 
 	# Create background-subtracted image:
 	image.subclean = image.clean - image.background
@@ -234,18 +236,48 @@ def photometry(fileid, output_folder=None):
 
 	# Use DAOStarFinder to search the image for stars, and only use reference-stars where a
 	# star was actually detected close to the references-star coordinate:
-	cleanout_references = (len(references) > 20)
-
+	min_references = 6
+	cleanout_references = (len(references) > 6)
+	logger.debug("Number of references before cleaning: %d", len(references))
 	if cleanout_references:
-		daofind_tbl = DAOStarFinder(100, fwhm=fwhm, roundlo=-0.5, roundhi=0.5).find_stars(image.subclean, mask=image.mask)
-		indx_good = np.zeros(len(references), dtype='bool')
-		for k, ref in enumerate(references):
-			dist = np.sqrt( (daofind_tbl['xcentroid'] - ref['pixel_column'])**2 + (daofind_tbl['ycentroid'] - ref['pixel_row'])**2 )
-			if np.any(dist <= fwhm/4): # Cutoff set somewhat arbitrary
-				indx_good[k] = True
-
-		references = references[indx_good]
-
+		# Arguments for DAOStarFind, starting with the strictest and ending with the
+		# least strict settings to try:
+		# We will stop with the first set that yield more than the minimum number
+		# of reference stars.
+		daofind_args = [{
+			'threshold': 7 * image.std,
+			'fwhm': fwhm,
+			'exclude_border': True,
+			'sharphi': 0.8,
+			'sigma_radius': 1.1,
+			'peakmax': image.peakmax
+		}, {
+			'threshold': 3 * image.std,
+			'fwhm': fwhm,
+			'roundlo': -0.5,
+			'roundhi': 0.5
+		}]
+
+		# Loop through argument sets for DAOStarFind:
+		for kwargs in daofind_args:
+			# Run DAOStarFind with the given arguments:
+			daofind_tbl = DAOStarFinder(**kwargs).find_stars(image.subclean, mask=image.mask)
+
+			# Match the found stars with the catalog references:
+			indx_good = np.zeros(len(references), dtype='bool')
+			for k, ref in enumerate(references):
+				dist = np.sqrt( (daofind_tbl['xcentroid'] - ref['pixel_column'])**2 + (daofind_tbl['ycentroid'] - ref['pixel_row'])**2 )
+				if np.any(dist <= fwhm/4): # Cutoff set somewhat arbitrary
+					indx_good[k] = True
+
+			logger.debug("Number of references after cleaning: %d", np.sum(indx_good))
+			if np.sum(indx_good) >= min_references:
+				references = references[indx_good]
+				break
+
+	logger.debug("Number of references after cleaning: %d", len(references))
+
+	# Create plot of target and reference star positions:
 	fig, ax = plt.subplots(1, 1, figsize=(20, 18))
 	plot_image(image.subclean, ax=ax, scale='log', cbar='right', title=target_name)
 	ax.scatter(references['pixel_column'], references['pixel_row'], c='r', marker='o', alpha=0.6)