From 951e941e0e2a800e634341e4c3cf9f14157c5aca Mon Sep 17 00:00:00 2001 From: Emir Date: Fri, 26 Jun 2020 14:10:50 +0200 Subject: [PATCH 01/10] add basic reference checking using DAOStarFinder and modify clean_references minimum. --- flows/api/datafiles.py | 1 - flows/load_image.py | 9 ++++ flows/photometry.py | 110 +++++++++++++++++++++++++++++++++++++++-- 3 files changed, 115 insertions(+), 5 deletions(-) diff --git a/flows/api/datafiles.py b/flows/api/datafiles.py index 6cb36ea..a7861f5 100644 --- a/flows/api/datafiles.py +++ b/flows/api/datafiles.py @@ -18,7 +18,6 @@ def get_datafile(fileid): token = config.get('api', 'token', fallback=None) if token is None: raise Exception("No API token has been defined") - r = requests.get('https://flows.phys.au.dk/api/datafiles.php', params={'fileid': fileid}, headers={'Authorization': 'Bearer ' + token}) diff --git a/flows/load_image.py b/flows/load_image.py index 60769b9..39a4f54 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -55,6 +55,7 @@ def load_image(FILENAME): # Specific headers: image.exptime = float(hdr['EXPTIME']) # * u.second + image.nonlin = 1e5 #default value # Timestamp: if origin == 'LCOGT': @@ -66,6 +67,10 @@ def load_image(FILENAME): image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=observatory) image.photfilter = hdr['FILTER'] + #Get non-linear limit + image.nonlin = hdr.get('MAXLIN')#Presumed non-linearity limit from header + image.nonlin = 60000 #From experience, this one is better. + #@TODO: Use actual or some fraction of the non-linearity limit 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 @@ -79,6 +84,10 @@ def load_image(FILENAME): 'i SDSS': 'ip', 'u SDSS': 'up' }.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER']) + #get real non-linear limit + image.nonlin = 80000 #For ALFOSC D, 1x1, 200; the standard for SNe. + #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 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 diff --git a/flows/photometry.py b/flows/photometry.py index 61b11b3..9172694 100644 --- a/flows/photometry.py +++ b/flows/photometry.py @@ -125,6 +125,7 @@ def photometry(fileid, output_folder=None): # Calculate the targets position in the image: target_pixel_pos = image.wcs.all_world2pix([[target['ra'], target['decl']]], 0)[0] + image.target_pixel_pos = target_pixel_pos #Feed this to target so we can use it to clean around it. # Clean out the references: hsize = 10 @@ -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 @@ -231,20 +233,59 @@ def photometry(fileid, output_folder=None): logger.info("FWHM: %f", fwhm) if np.isnan(fwhm): raise Exception("Could not estimate FWHM") + else: + image.fwhm=fwhm #write fwhm to image object + # 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) - + cleanout_references = (len(references) > 6) + logger.info("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) + #daofind_tbl = DAOStarFinder(100, fwhm=fwhm, roundlo=-0.5, roundhi=0.5).find_stars(image.subclean, mask=image.mask) + #Now We use a more struct DAOStarFinder check. + daofind_tbl = DAOStarFinder(fwhm=fwhm, threshold=7.*image.std,exclude_border=True, + sharphi=0.8,sigma_radius=1.1,peakmax=image.nonlin).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 - + + + print(references) + + #Arbitary but we don't want to cut too many references. In that case, go back to less strict version. + #TODO: Change the checking here to a function + min_references=6 + logger.info("Number of references 1: %d", len(references[indx_good])) + if len(references[indx_good])<= min_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 + + logger.info("Number of references 2: %d", len(references[indx_good])) references = references[indx_good] + + + ## Can be used for making cuts based on sharpness or roundness parameters from daofind + calculate_daostar_properties=True + if calculate_daostar_properties: + idx_sources_good=np.zeros(len(daofind_tbl), dtype='bool') + for k, daoref in enumerate(daofind_tbl): + dist = np.sqrt( (daoref['xcentroid'] - references['pixel_column'])**2 + (daoref['ycentroid'] - references['pixel_row'])**2 ) + if np.any(dist <= fwhm/4): + idx_sources_good[k] = True + daoclean = daofind_tbl[idx_sources_good] + print(daoclean) + + #Remove sources within r arcsecond of the target pos from the reference list + #@TODO: Change target with host galaxy and r with a multiple of half-light radius. + references=rm_sources_within(references,target_coord,r=10*u.arcsec) + fig, ax = plt.subplots(1, 1, figsize=(20, 18)) plot_image(image.subclean, ax=ax, scale='log', cbar='right', title=target_name) @@ -574,3 +615,64 @@ def photometry(fileid, output_folder=None): logger.info("Photometry took: %f seconds", toc-tic) return photometry_output + + +def rm_sources_within(references,target_coord,r=10.0*u.arcsec): + '''remove sources within r units of target. + Return references sans the close one.''' + #Create SkyCoord Array + refs=coords.SkyCoord(references['ra'],references['decl'],frame='icrs') + seps=target_coord.separation(refs) #Separation in astropy angles + return references[seps>r] + +#### +# Reference code for removing based on DAOStarFinder parameters. +### + +#def rm_galaxies_jank(image,magnitude_limit=-1.5,strictness=1): +# '''remove galaxies and space junk but using a very janky +# method. Use as basis for writing a proper removal script. ''' +# +# import pandas as pd +# from astropy.stats import sigma_clipped_stats +# #Fail nicely, we don't care if this step doesn't work. +# #Only problem is we won't be able to keyboard interrupt here, which is not ideal. +# #@TODO: Fix naked try except. Probably want to except when DAOStarFinder is fed null arguments +# #Also when image or rows are flat and ic is empty, pandas df cannot be built etc. +# +# try: +# sources = DAOStarFinder(fwhm=FWHM, threshold=7.*image.std,exclude_border=True, +# sharphi=0.8,sigma_radius=1.1,peakmax=image.nonlin).find_stars(image.subclean, mask=image.mask) +# +# +# #index of finite daofind sources +# ic = [i for i, row in enumerate(sources) if np.isfinite(row['mag']) and row['mag'] <= magnitude_limit] +# if strictiness==1: +# return ic +# #I don't even know why this is in pandas, I just copied from my old code +# #Update to just use astropy table and remove the unnecessary for loops. +# elif strictness==2 +# df=pd.DataFrame(index=ic,columns=['s','r1','r2']) +# sp,rp1,rp2=[],[],[] +# for i in ic: +# sp.append(sources[i]['sharpness']) +# rp1.append(sources[i]['roundness1']) +# rp2.append(sources[i]['roundness2']) +# df['s']=sp +# df['r1']=rp1 +# df['r2']=rp2 +# +# df['delta_s']=df.s-df.s.mean() +# df['delta_r1']=df.r1-df.r1.mean() +# df['delta_r2']=df.r2-df.r2.mean() +# +# df['del_r2']=(df.delta_r1**2+df.delta_r2**2)**(1/2) +# +# masked_s=sigma_clip(df['s'],sigma=2.5) +# +# return df[((df.delta_s<-0.1) & (df.del_r2>df.del_r2.median())) | (masked_s.mask==True)].index.values +# except: +# return 'None' + + + \ No newline at end of file From 91630314710b5c4b7e382114d19b65e50788cf14 Mon Sep 17 00:00:00 2001 From: Emir Date: Fri, 26 Jun 2020 14:12:34 +0200 Subject: [PATCH 02/10] add jupyter notebooks to gitignore --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index 1b693bf..177db5a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,10 @@ flows/config.ini flows/casjobs/CasJobs.config +#Jupyter lab & notebooks +*.ipynb +.idea/ + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] From bb61270c164a263251d90ac3d1f6dfc301e3dcb2 Mon Sep 17 00:00:00 2001 From: Emir Date: Fri, 26 Jun 2020 14:27:21 +0200 Subject: [PATCH 03/10] fix dangling print statements --- flows/photometry.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/flows/photometry.py b/flows/photometry.py index 9172694..c65058b 100644 --- a/flows/photometry.py +++ b/flows/photometry.py @@ -240,7 +240,7 @@ 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) > 6) - logger.info("Number of references before cleaning: %d", len(references)) + 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) #Now We use a more struct DAOStarFinder check. @@ -251,14 +251,12 @@ def photometry(fileid, output_folder=None): 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 - - - print(references) + #Arbitary but we don't want to cut too many references. In that case, go back to less strict version. #TODO: Change the checking here to a function min_references=6 - logger.info("Number of references 1: %d", len(references[indx_good])) + logger.debug("Number of references after strict cleaning: %d", len(references[indx_good])) if len(references[indx_good])<= min_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') @@ -267,12 +265,12 @@ def photometry(fileid, output_folder=None): if np.any(dist <= fwhm/4): # Cutoff set somewhat arbitrary indx_good[k] = True - logger.info("Number of references 2: %d", len(references[indx_good])) + logger.debug("Number of references after less strict cleaning: %d", len(references[indx_good])) references = references[indx_good] ## Can be used for making cuts based on sharpness or roundness parameters from daofind - calculate_daostar_properties=True + calculate_daostar_properties=False if calculate_daostar_properties: idx_sources_good=np.zeros(len(daofind_tbl), dtype='bool') for k, daoref in enumerate(daofind_tbl): @@ -280,7 +278,7 @@ def photometry(fileid, output_folder=None): if np.any(dist <= fwhm/4): idx_sources_good[k] = True daoclean = daofind_tbl[idx_sources_good] - print(daoclean) + #Remove sources within r arcsecond of the target pos from the reference list #@TODO: Change target with host galaxy and r with a multiple of half-light radius. From 59c68d2a7aa0d83d0363ddfd48f8c4f12ed0b431 Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 16:30:31 +0200 Subject: [PATCH 04/10] Some of the added functionality is already there Just changing the distance check from pixels to degrees --- flows/photometry.py | 88 ++++++--------------------------------------- 1 file changed, 10 insertions(+), 78 deletions(-) diff --git a/flows/photometry.py b/flows/photometry.py index c65058b..34bc61a 100644 --- a/flows/photometry.py +++ b/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() @@ -131,7 +131,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) @@ -233,10 +234,7 @@ def photometry(fileid, output_folder=None): logger.info("FWHM: %f", fwhm) if np.isnan(fwhm): raise Exception("Could not estimate FWHM") - else: - image.fwhm=fwhm #write fwhm to image object - # 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) > 6) @@ -251,9 +249,9 @@ def photometry(fileid, output_folder=None): 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 - - - #Arbitary but we don't want to cut too many references. In that case, go back to less strict version. + + + #Arbitary but we don't want to cut too many references. In that case, go back to less strict version. #TODO: Change the checking here to a function min_references=6 logger.debug("Number of references after strict cleaning: %d", len(references[indx_good])) @@ -264,11 +262,11 @@ def photometry(fileid, output_folder=None): 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 less strict cleaning: %d", len(references[indx_good])) references = references[indx_good] - - + + ## Can be used for making cuts based on sharpness or roundness parameters from daofind calculate_daostar_properties=False if calculate_daostar_properties: @@ -278,12 +276,7 @@ def photometry(fileid, output_folder=None): if np.any(dist <= fwhm/4): idx_sources_good[k] = True daoclean = daofind_tbl[idx_sources_good] - - - #Remove sources within r arcsecond of the target pos from the reference list - #@TODO: Change target with host galaxy and r with a multiple of half-light radius. - references=rm_sources_within(references,target_coord,r=10*u.arcsec) - + fig, ax = plt.subplots(1, 1, figsize=(20, 18)) plot_image(image.subclean, ax=ax, scale='log', cbar='right', title=target_name) @@ -613,64 +606,3 @@ def photometry(fileid, output_folder=None): logger.info("Photometry took: %f seconds", toc-tic) return photometry_output - - -def rm_sources_within(references,target_coord,r=10.0*u.arcsec): - '''remove sources within r units of target. - Return references sans the close one.''' - #Create SkyCoord Array - refs=coords.SkyCoord(references['ra'],references['decl'],frame='icrs') - seps=target_coord.separation(refs) #Separation in astropy angles - return references[seps>r] - -#### -# Reference code for removing based on DAOStarFinder parameters. -### - -#def rm_galaxies_jank(image,magnitude_limit=-1.5,strictness=1): -# '''remove galaxies and space junk but using a very janky -# method. Use as basis for writing a proper removal script. ''' -# -# import pandas as pd -# from astropy.stats import sigma_clipped_stats -# #Fail nicely, we don't care if this step doesn't work. -# #Only problem is we won't be able to keyboard interrupt here, which is not ideal. -# #@TODO: Fix naked try except. Probably want to except when DAOStarFinder is fed null arguments -# #Also when image or rows are flat and ic is empty, pandas df cannot be built etc. -# -# try: -# sources = DAOStarFinder(fwhm=FWHM, threshold=7.*image.std,exclude_border=True, -# sharphi=0.8,sigma_radius=1.1,peakmax=image.nonlin).find_stars(image.subclean, mask=image.mask) -# -# -# #index of finite daofind sources -# ic = [i for i, row in enumerate(sources) if np.isfinite(row['mag']) and row['mag'] <= magnitude_limit] -# if strictiness==1: -# return ic -# #I don't even know why this is in pandas, I just copied from my old code -# #Update to just use astropy table and remove the unnecessary for loops. -# elif strictness==2 -# df=pd.DataFrame(index=ic,columns=['s','r1','r2']) -# sp,rp1,rp2=[],[],[] -# for i in ic: -# sp.append(sources[i]['sharpness']) -# rp1.append(sources[i]['roundness1']) -# rp2.append(sources[i]['roundness2']) -# df['s']=sp -# df['r1']=rp1 -# df['r2']=rp2 -# -# df['delta_s']=df.s-df.s.mean() -# df['delta_r1']=df.r1-df.r1.mean() -# df['delta_r2']=df.r2-df.r2.mean() -# -# df['del_r2']=(df.delta_r1**2+df.delta_r2**2)**(1/2) -# -# masked_s=sigma_clip(df['s'],sigma=2.5) -# -# return df[((df.delta_s<-0.1) & (df.del_r2>df.del_r2.median())) | (masked_s.mask==True)].index.values -# except: -# return 'None' - - - \ No newline at end of file From f3c82b0e2a1caa32121ac568cee4c3d2cc53050d Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 16:34:40 +0200 Subject: [PATCH 05/10] No need to store target position in image object --- flows/photometry.py | 1 - 1 file changed, 1 deletion(-) diff --git a/flows/photometry.py b/flows/photometry.py index 34bc61a..c3a8b7b 100644 --- a/flows/photometry.py +++ b/flows/photometry.py @@ -125,7 +125,6 @@ def photometry(fileid, output_folder=None): # Calculate the targets position in the image: target_pixel_pos = image.wcs.all_world2pix([[target['ra'], target['decl']]], 0)[0] - image.target_pixel_pos = target_pixel_pos #Feed this to target so we can use it to clean around it. # Clean out the references: hsize = 10 From 3d21a5cac28c803fb3a7597cf78b4bcba3330d4e Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 16:50:04 +0200 Subject: [PATCH 06/10] Cleanup --- flows/api/datafiles.py | 1 + 1 file changed, 1 insertion(+) diff --git a/flows/api/datafiles.py b/flows/api/datafiles.py index a7861f5..6cb36ea 100644 --- a/flows/api/datafiles.py +++ b/flows/api/datafiles.py @@ -18,6 +18,7 @@ def get_datafile(fileid): token = config.get('api', 'token', fallback=None) if token is None: raise Exception("No API token has been defined") + r = requests.get('https://flows.phys.au.dk/api/datafiles.php', params={'fileid': fileid}, headers={'Authorization': 'Bearer ' + token}) From 6e67aa7b55434e707f1a19fb5f94957e06a34b4a Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 16:57:14 +0200 Subject: [PATCH 07/10] Removed unused code-block --- flows/photometry.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/flows/photometry.py b/flows/photometry.py index c3a8b7b..cd0c971 100644 --- a/flows/photometry.py +++ b/flows/photometry.py @@ -265,18 +265,6 @@ def photometry(fileid, output_folder=None): logger.debug("Number of references after less strict cleaning: %d", len(references[indx_good])) references = references[indx_good] - - ## Can be used for making cuts based on sharpness or roundness parameters from daofind - calculate_daostar_properties=False - if calculate_daostar_properties: - idx_sources_good=np.zeros(len(daofind_tbl), dtype='bool') - for k, daoref in enumerate(daofind_tbl): - dist = np.sqrt( (daoref['xcentroid'] - references['pixel_column'])**2 + (daoref['ycentroid'] - references['pixel_row'])**2 ) - if np.any(dist <= fwhm/4): - idx_sources_good[k] = True - daoclean = daofind_tbl[idx_sources_good] - - 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) From fbe24bd8e86c2b41501aeb523d6d49ee4623eeda Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 17:26:22 +0200 Subject: [PATCH 08/10] We should not ignore all Jupyter notebooks --- .gitignore | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/.gitignore b/.gitignore index 177db5a..a0311e0 100644 --- a/.gitignore +++ b/.gitignore @@ -1,10 +1,6 @@ flows/config.ini flows/casjobs/CasJobs.config -#Jupyter lab & notebooks -*.ipynb -.idea/ - # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] @@ -134,3 +130,6 @@ dmypy.json # Pyre type checker .pyre/ + +# JetBrains IDEs +.idea/ From c65d513f3fb9c7319c758439cc7c07a03dce2849 Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 17:47:31 +0200 Subject: [PATCH 09/10] Rewrote DAOStarFind section so we dont repeat ourselves too much --- flows/photometry.py | 51 ++++++++++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/flows/photometry.py b/flows/photometry.py index cd0c971..6d467c7 100644 --- a/flows/photometry.py +++ b/flows/photometry.py @@ -236,35 +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: + 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) - #Now We use a more struct DAOStarFinder check. - daofind_tbl = DAOStarFinder(fwhm=fwhm, threshold=7.*image.std,exclude_border=True, - sharphi=0.8,sigma_radius=1.1,peakmax=image.nonlin).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 - - - #Arbitary but we don't want to cut too many references. In that case, go back to less strict version. - #TODO: Change the checking here to a function - min_references=6 - logger.debug("Number of references after strict cleaning: %d", len(references[indx_good])) - if len(references[indx_good])<= min_references: - daofind_tbl = DAOStarFinder(100, fwhm=fwhm, roundlo=-0.5, roundhi=0.5).find_stars(image.subclean, mask=image.mask) + # 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.nonlin + }, { + 'threshold': 100, # TODO: This should be changed!!! + '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 less strict cleaning: %d", len(references[indx_good])) - references = references[indx_good] + 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) From 75498406dd9559158fc45bac5a9d977b16494198 Mon Sep 17 00:00:00 2001 From: Rasmus Handberg Date: Tue, 30 Jun 2020 18:38:37 +0200 Subject: [PATCH 10/10] Renamed nonlin to peakmax - Changed default peakmax to be undefined (None). - Changed weakest DAOFind threshold to be relative instead of fixed absolute value. --- flows/load_image.py | 20 +++++++++++--------- flows/photometry.py | 4 ++-- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/flows/load_image.py b/flows/load_image.py index 39a4f54..5575f49 100644 --- a/flows/load_image.py +++ b/flows/load_image.py @@ -55,7 +55,7 @@ def load_image(FILENAME): # Specific headers: image.exptime = float(hdr['EXPTIME']) # * u.second - image.nonlin = 1e5 #default value + image.peakmax = None # Maximum value above which data is not to be trusted # Timestamp: if origin == 'LCOGT': @@ -67,10 +67,11 @@ def load_image(FILENAME): image.obstime = Time(hdr['MJD-OBS'], format='mjd', scale='utc', location=observatory) image.photfilter = hdr['FILTER'] - #Get non-linear limit - image.nonlin = hdr.get('MAXLIN')#Presumed non-linearity limit from header - image.nonlin = 60000 #From experience, this one is better. - #@TODO: Use actual or some fraction of the non-linearity limit + + # 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 @@ -84,10 +85,11 @@ def load_image(FILENAME): 'i SDSS': 'ip', 'u SDSS': 'up' }.get(hdr['FILTER'].replace('_', ' '), hdr['FILTER']) - #get real non-linear limit - image.nonlin = 80000 #For ALFOSC D, 1x1, 200; the standard for SNe. - #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 + + # 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 diff --git a/flows/photometry.py b/flows/photometry.py index 6d467c7..7ed7e6a 100644 --- a/flows/photometry.py +++ b/flows/photometry.py @@ -250,9 +250,9 @@ def photometry(fileid, output_folder=None): 'exclude_border': True, 'sharphi': 0.8, 'sigma_radius': 1.1, - 'peakmax': image.nonlin + 'peakmax': image.peakmax }, { - 'threshold': 100, # TODO: This should be changed!!! + 'threshold': 3 * image.std, 'fwhm': fwhm, 'roundlo': -0.5, 'roundhi': 0.5