Commit 629d372a authored by PAduverne's avatar PAduverne
Browse files

Final version of the photometry script

parent 27ac4f10
......@@ -910,9 +910,9 @@ def calibration(source, catalog, cata_phot_syst, stars, wcs_data,
coefs, V = np.ma.polyfit(source['Magnitude'], source[band],
1, cov=True)
source['Calibrated Mag'] = pl.polyval(source['Magnitude'],
source['Calibrated_Mag'] = pl.polyval(source['Magnitude'],
np.flip(coefs))
source['error magnitude'] = source['Calibrated Mag'] - source[band]
source['error magnitude'] = source['Calibrated_Mag'] - source[band]
da_2, db_2 = V[0,0], V[1,1]
return source, coefs, da_2, db_2, star
......
......@@ -9,11 +9,11 @@ Created on Mon Apr 19 16:58:12 2021
import argparse
import sys
import os
import time
# import time
import warnings
from astropy.io import fits
from astropy.wcs import WCS, FITSFixedWarning
from astropy.table import Table
from astropy import units as u
import hjson
import numpy as np
import numpy.polynomial.polynomial as pl
......@@ -39,9 +39,10 @@ warnings.simplefilter(action="ignore", category=FutureWarning)
"""Muphoten main script for photometry."""
def main():
"""Console script for muphoten."""
path_muphoten = getpath()
# path_muphoten = getpath()
telescope_list = getTel()
catalog_list = CataList()
phot_list = PhotList()
......@@ -60,6 +61,11 @@ def main():
type=str,
help="Path to images.")
parser.add_argument("--output-file",
required=True,
type=str,
help="Name of the output file.")
parser.add_argument("--telescope",
dest="telescope",
choices=telescope_list,
......@@ -121,12 +127,12 @@ def main():
config = load_config(args.telescope)
mu = config['muphoten']['conf']
with open(mu) as json_file:
mu_conf = hjson.load(json_file)
mu_conf = hjson.load(json_file)
# Preparing the files for analysis
images = list_files(args.images, get_subdirs=False)
print('{} images to calibrate'.format(len(images)))
glob_start = time.time()
# glob_start = time.time()
results_file = os.path.join(args.images, 'results')
col_name = ['filename',
'telescope',
......@@ -138,8 +144,8 @@ def main():
'transient_sky_centroid',
'transient_flux',
'transient_background_sum',
'transient_ins_magnitude',
'transient_magnitude',
'transient_calibrated_magnitude',
'transient_SNR',
'transient_error_bkg',
'transient_error',
......@@ -168,9 +174,9 @@ def main():
image_rep = os.path.dirname(image)
calib_rep = os.path.join(image_rep, 'calibration')
print(calib_rep)
sub_im = os.path.join(image_rep, 'subtraction/sub_'+
image.split('/')[-1])
start_im_time = time.time() # Computing time for one image processing
sub_im = os.path.join(image_rep, 'subtraction/sub_' +
image.split('/')[-1])
# start_im_time = time.time() # Computing time for one image processing
print('Processing image n° {} : {}'.format(i, image_name))
print('Sanitising header')
sanitise_fits(image)
......@@ -194,7 +200,7 @@ def main():
print('Subtracting background in image : {}'.format(image_name))
# Getting the backgroud, rms and bkg subtracted images
start_bkg_time = time.time()
# start_bkg_time = time.time()
if args.do_clean:
if args.save_bkg:
background_rep = make_sub_result_rep(image_rep,
......@@ -204,14 +210,14 @@ def main():
path2save = None
data_clean, background,\
rms = sub_background(data, header,
image_name,
n_sigma=mu_conf['bkg_n_sigma'],
box_size=mu_conf['bkg_box'],
bkg_estimator=mu_conf['bkg_estimator'],
filter_size=mu_conf['bkg_filter_size'],
path2save=path2save,
png=args.save_bkg)
rms = sub_background(data, header,
image_name,
n_sigma=mu_conf['bkg_n_sigma'],
box_size=mu_conf['bkg_box'],
bkg_estimator=mu_conf['bkg_estimator'],
filter_size=mu_conf['bkg_filter_size'],
path2save=path2save,
png=args.save_bkg)
else:
background_rep = os.path.join(image_rep, 'background')
......@@ -221,9 +227,9 @@ def main():
background = get_im_res(background_rep, image_name,
prefix='background_')
rms = get_im_res(background_rep, image_name, prefix='rms_')
bkg_time = time.time() - start_bkg_time
# bkg_time = time.time() - start_bkg_time
start_photo_time = time.time()
# start_photo_time = time.time()
# Detection of the sources & segmentation
segmentation = source_detection(data, background, rms,
n_sigma=mu_conf['det_n_sigma'],
......@@ -237,7 +243,7 @@ def main():
FWHM_psf = psfex(image, config, verbose="QUIET")
results['psf'][i] = FWHM_psf
print(FWHM_psf[0])
if args.photo == phot_list[0]: # Kron photometry
if args.photo == phot_list[0]: # Kron photometry
print('Performing Kron Photometry')
photometry_table = Kron_Photometry(data, data_clean, segmentation,
header, background, rms,
......@@ -263,16 +269,16 @@ def main():
print('Performing calibration using {}'.format(cata_name))
photometry_table, coefs,\
da_2, db_2, star_ref = calibration(source=photometry_table,
catalog=catalog,
cata_phot_syst=system,
stars=stars,
wcs_data=wcs_data,
n_sigma=mu_conf['clipping_calib'],
cut_mag=mu_conf['cut_mag'],
band=band)
da_2, db_2, star_ref = calibration(source=photometry_table,
catalog=catalog,
cata_phot_syst=system,
stars=stars,
wcs_data=wcs_data,
n_sigma=mu_conf['clipping_calib'],
cut_mag=mu_conf['cut_mag'],
band=band)
plot_calib(band, photometry_table, coefs, path2save=calib_rep)
photometry_table.write(calib_rep +'/calibration.dat',
photometry_table.write(calib_rep + '/calibration.dat',
format='ascii.commented_header',
overwrite=True)
......@@ -287,33 +293,33 @@ def main():
star_ref = compute_snr(star_ref)
results['star_ref_tot_error'][i] = star_ref['total_error'][0]
results['star_ref_magnitude_cata'] = star_ref[band][0]
results['starresults.show_in_browser(jsviewer=True)_ref_error_cata'][i] = star_ref['error'][0]
results['star_ref_error_cata'][i] = star_ref['error'][0]
results['star_ref_flux'][i] = star_ref['aperture_sum'][0]
results['star_ref_Magnitude'][i] = star_ref['Magnitude'][0]
results['star_ref_magnitude'][i] = star_ref['Calibrated_Mag'][0]
results['star_ref_bkg_sum'][i] = star_ref['background_sum'][0]
results['star_ref_SNR'][i] = star_ref['SNR'][0]
results['star_ref_distance'][i] = star_ref['dist'][0]
if star_ref['dist'][0]:
results['star_detection'][i] = 0 # Reference star detected
if star_ref['dist'][0] < 5:
results['star_detection'][i] = 0 # Reference star detected
else:
results['star_detection'][i] = 1 # Reference star not detected
results['star_detection'][i] = 1 # Reference star not detected
# Detecting the transient in the subtracted image.
print('Treating subtracted image')
hdu_sub = fits.open(sub_im)[0]
data_sub, header_sub = hdu_sub.data, hdu_sub.header
wcs_sub = WCS(header_sub)
# start_sub_time = time.time()
mask=(data_sub == 1e-30) | (data_sub == 0)
mask = (data_sub == 1e-30) | (data_sub == 0)
data_clean_sub, background_sub,\
rms_sub = sub_background(data_sub, header_sub,
image_name,
n_sigma=mu_conf['bkg_n_sigma_sub'],
box_size=mu_conf['bkg_box_sub'],
bkg_estimator=mu_conf['bkg_estimator'],
filter_size=mu_conf['bkg_filter_size_sub'],
mask=mask)
rms_sub = sub_background(data_sub, header_sub,
image_name,
n_sigma=mu_conf['bkg_n_sigma_sub'],
box_size=mu_conf['bkg_box_sub'],
bkg_estimator=mu_conf['bkg_estimator'],
filter_size=mu_conf['bkg_filter_size_sub'],
mask=mask)
# Creating a mask centered on the transient to avoid the source
# detection in all the sub image. Otherwise, there could issues
# because of artefacts in the sub image.
......@@ -329,32 +335,78 @@ def main():
mask=mask)
# If the transient is not detected in the image
if seg_sub is None:
dist = 30.
else:
if args.photo == phot_list[0]: # Kron photometry
print('Performing Kron Photometry in sub image')
photometry_sub = Kron_Photometry(data_sub, data_clean_sub,
seg_sub,
header_sub, background_sub,
rms_sub, wcs_sub)
if args.photo == phot_list[1]:
print('Performing Isophotal Photometry in sub image')
photometry_sub = Isophotal_Photometry(data_sub,
data_clean_sub,
seg_sub, header_sub,
background_sub,
rms_sub, wcs_sub)
if args.photo == phot_list[2]:
print('Performing Fixed size Photometry in sub image')
radius = mu_conf['fixed_factor'] * FWHM_psf[0]
photometry_sub = Fixed_Photometry(data_sub,
data_clean_sub,
seg_sub,
header_sub, radius,
background_sub,
rms_sub, wcs_sub)
transient = get_source(photometry_sub, transients, wcs_sub)
dist = transient['dist'][0] / u.pix
if dist < 5.: # Considering that the transient has been detected
results['transient_detection'][i] = 0
bkg_sum = Dophot(data_sub, transients,
wcs_sub, aperture=transient['aperture'])
flux = transient['aperture_sum'][0]
transient['background_sum'] = bkg_sum
transient = compute_error(transient, coefs[0], da_2, db_2)
transient = compute_snr(transient)
results['transient_sky_centroid'] = transient['sky_centroid']
results['transient_flux'][i] = flux
results['transient_ins_magnitude'][i] = transient['Magnitude'][0]
results['transient_magnitude'][i] = pl.polyval(transient['Magnitude'][0],
np.flip(coefs))
results['transient_error_tot'][i] = transient['total_error'][0]
results['transient_SNR'][i] = transient['SNR'][0]
results['transient_background_sum'][i] = bkg_sum
results['transient_error_bkg'][i] = 2.5 / \
(np.log(10)*np.sqrt(bkg_sum))
results['transient_error'][i] = 2.5/(np.log(10)*np.sqrt(flux))
results['transient_error_calib'][i] = transient['calib_error'][0]
results['transient_error_tot'][i] = transient['total_error'][0]
else:
results['transient_detection'][i] = 1
# Computing the number of counts in a aperture of % pixels to have
# Computing the number of counts in a aperture of 5 pixels to have
# an upper limit on the magnitude of the transient.
flux = Dophot(data_sub,
transients,
wcs_sub,
radius=5.0)
flux = Dophot(data_sub, transients, wcs_sub, radius=5.0)
bkg_sum = Dophot(background, transients, wcs_data, radius=5.0)
mag = -2.5*np.log10(flux)
results['transient_flux'][i] = abs(flux)
results['transient_magnitude'][i] = -2.5*np.log10(flux)
results['transient_detection'][i] = 1
results['transient_background_sum'][i] = Dophot(background,
transients,
wcs_data,
radius=5.0)
results['transient_detection'][i] =
results.show_in_browser(jsviewer=True)
# background_object_sub = back.background_estimation(data_sub, n_sigma=3.0,
# box_size=(30,30),
# bkg_estimator='sex',
# filter_size=(3, 3),
# mask=mask)
# background_array_sub, background_rms_sub = back.bkg_to_array(background_object_sub)
# data_cleaned_sub = data_sub - background_array_sub
# if args.photo not in phot_list:
# raise ValueError('{} is not an avalaible photometry'
# 'in Muphoten'.format(args.photo))
results['transient_ins_magnitude'][i] = mag
results['transient_magnitude'][i] = pl.polyval(mag,
np.flip(coefs))
results['transient_background_sum'][i] = bkg_sum
results['transient_SNR'][i] = flux / np.sqrt(flux + bkg_sum)
results['transient_error_bkg'][i] = 2.5 / \
(np.log(10)*np.sqrt(bkg_sum))
results['transient_error'][i] = 2.5/(np.log(10)*np.sqrt(flux))
results['transient_error_calib'][i] = np.sqrt(mag**2 * da_2 +
db_2 +
coefs[0]*results['transient_error'][i]**2)
results['transient_error_tot'][i] = np.sqrt(results['transient_error'][i]**2 +
results['transient_error_bkg'][i]**2 +
results['transient_error_calib'][i]**2)
results.write(results_file + '/' + args.output_name,
format='ascii.commented_header')
return 0
......
......@@ -283,11 +283,15 @@ def get_source(table, coord_source, wcs, col_name='dist'):
source = table[table[col_name]==min(table[col_name])]
return source
def Dophot(data, position, wcs, radius=5.0):
def Dophot(data, position, wcs, radius=5.0, aperture=None):
position = wcs.all_world2pix(position, 1)
aperture = CircularAperture(position, r=radius)
aper_photo = aperture_photometry(data, aperture)
if aperture is None:
position = wcs.all_world2pix(position, 1)
aperture = CircularAperture(position, r=radius)
aper_photo = aperture_photometry(data, aperture)
else:
aperture = aperture.to_sky(wcs).to_pixel(wcs)
aper_photo = aperture_photometry(data, aperture)
return aper_photo['aperture_sum'][0]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment