Commit b578958e authored by PAduverne's avatar PAduverne
Browse files

change the script names

parent ff14160e
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 19 16:58:12 2021
@author: duverne, IJClab, Orsay, duverne@protonmail.com
"""
import argparse
import sys
import os
# import time
import warnings
from astropy.io import fits
from astropy.wcs import WCS, FITSFixedWarning
from astropy import units as u
import hjson
import numpy as np
import numpy.polynomial.polynomial as pl
from muphoten.background import sub_background
from muphoten.utils import (load_coord_file, getpath,
getTel, make_results_dir, get_filename,
list_files, load_config, set_results_table,
make_sub_result_rep, get_im_res, PhotList)
from muphoten.sanitise import sanitise_fits
from muphoten.catalog import (get_filter, getCata, CataList, calibration,
catalog_choice, )
from muphoten.photometry import (source_detection, Kron_Photometry, masking,
Fixed_Photometry, Isophotal_Photometry,
Dophot, get_source, compute_error,
compute_snr)
from muphoten.psfex import psfex
from muphoten.plot_image import (plot_sementation_image,
plot_image, plot_calib)
warnings.filterwarnings("ignore", category=FITSFixedWarning)
warnings.simplefilter(action="ignore", category=FutureWarning)
"""Muphoten main script for photometry."""
def main():
"""Console script for muphoten."""
# path_muphoten = getpath()
telescope_list = getTel()
catalog_list = CataList()
phot_list = PhotList()
parser = argparse.ArgumentParser(description="Performing photometry of\
transient objects.")
parser.add_argument("--coord",
dest="coord",
required=True,
type=str,
help="Coordinates file's path.")
parser.add_argument("--images",
dest="images",
required=True,
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,
required=True,
type=str,
help="Telescope that acquired the images.")
parser.add_argument("--catalog",
dest="catalog",
choices=catalog_list,
required=False,
type=str,
help="Catalog to use for the calibration.")
parser.add_argument("--photo",
dest="photo",
choices=phot_list,
required=True,
type=str,
help="Type of photometry. Use muphoten configuration \
file to modify the radius of the isophote or the \
factor used for the fixed apertures.\
Default is iso.")
parser.add_argument("--keep-old",
"--keep",
dest="keep",
required=False,
action="store_true",
help="Keep previous results")
parser.add_argument("--skip-processed",
"--skip",
dest="skip",
required=False,
action="store_true",
help="Skip already processed files")
parser.add_argument("--Do-clean",
dest="do_clean",
action="store_true",
help="Perform the background subbtraction.\
Default : False")
parser.add_argument("--save-bkg",
dest="save_bkg",
action="store_true",
help="Save the background images as fits and png.\
Default : False")
args = parser.parse_args()
# Loading the coordinates of the analysed objects
transients, stars = load_coord_file(args.coord)
print('Computing lightcurve for {} transients'.format(len(transients)))
print('Using {} stars to veto poor quality images'.format(len(stars)))
print('Loading configuration for Muphoten,')
config = load_config(args.telescope)
mu = config['muphoten']['conf']
with open(mu) as 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()
results_file = os.path.join(args.images, 'results')
col_name = ['filename',
'telescope',
'filter',
'time',
'psf',
'a', 'da',
'b', 'db',
'transient_sky_centroid',
'transient_flux',
'transient_background_sum',
'transient_ins_magnitude',
'transient_magnitude',
'transient_SNR',
'transient_error_bkg',
'transient_error',
'transient_error_calib',
'transient_error_tot',
'transient_detection',
'star_ref_magnitude',
'star_ref_magnitude_cata',
'star_ref_tot_error',
'star_ref_error_cata',
'star_ref_distance',
'star_ref_flux',
'star_ref_bkg_sum',
'star_ref_SNR',
'star_detection']
results = set_results_table(col_name, len(images))
results['telescope'] = args.telescope
for i, image in enumerate(images):
image = make_results_dir(image,
directory='calibration',
keep=args.keep,
skip=args.skip)
image_name = get_filename(image)
results['filename'][i] = image_name
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
print('Processing image n° {} : {}'.format(i, image_name))
print('Sanitising header')
sanitise_fits(image)
print(image)
# Load the data
hdu = fits.open(image)[0]
data, header = hdu.data, hdu.header
wcs_data = WCS(header)
# Getting different useful informations : band, MJD,
# catalog for calibration.
results['time'][i] = header.get('MJD-OBS')
if args.catalog:
catalog, system = getCata(args.catalog)
cata_name = args.catalog
band = get_filter(header)
else:
catalog, system, cata_name = catalog_choice(header)
band = get_filter(header)
results['filter'][i] = band
print('Subtracting background in image : {}'.format(image_name))
# Getting the backgroud, rms and bkg subtracted images
# start_bkg_time = time.time()
if args.do_clean:
if args.save_bkg:
background_rep = make_sub_result_rep(image_rep,
directory='background')
path2save = image_rep
else:
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)
else:
background_rep = os.path.join(image_rep, 'background')
data_clean = get_im_res(image_rep, image_name,
prefix='clean_')
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
# start_photo_time = time.time()
# Detection of the sources & segmentation
segmentation = source_detection(data, background, rms,
n_sigma=mu_conf['det_n_sigma'],
deblend=mu_conf['det_deblend'],
threshold_level=mu_conf['det_treshold'],
npixels=mu_conf['det_npixels'],
contrast=mu_conf['contrast'])
plot_sementation_image(segmentation, path2save=calib_rep)
# Computing PSF
print('Computing PSF for image n° {}'.format(i))
FWHM_psf = psfex(image, config, verbose="QUIET")
results['psf'][i] = FWHM_psf
print(FWHM_psf[0])
if args.photo == phot_list[0]: # Kron photometry
print('Performing Kron Photometry')
photometry_table = Kron_Photometry(data, data_clean, segmentation,
header, background, rms,
wcs=wcs_data)
plot_image(data_clean, ap=photometry_table['kron_aperture'],
wcs=wcs_data, path2save=calib_rep)
if args.photo == phot_list[1]:
print('Performing Isophotal Photometry')
photometry_table = Isophotal_Photometry(data, data_clean,
segmentation, header,
background, rms,
wcs=wcs_data)
plot_image(data_clean, sky_ap=photometry_table['aperture'],
wcs=wcs_data, path2save=calib_rep)
if args.photo == phot_list[2]:
print('Performing Fixed size Photometry')
radius = mu_conf['fixed_factor'] * FWHM_psf[0]
photometry_table = Fixed_Photometry(data, data_clean, segmentation,
header, radius, background,
rms, wcs=wcs_data)
plot_image(data_clean, sky_ap=photometry_table['aperture'],
wcs=wcs_data, path2save=calib_rep)
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)
plot_calib(band, photometry_table, coefs, path2save=calib_rep)
photometry_table.write(calib_rep + '/calibration.dat',
format='ascii.commented_header',
overwrite=True)
# Saving the calibration results.
results['a'][i] = coefs[0]
results['da'][i] = np.sqrt(da_2)
results['b'][i] = coefs[1]
results['db'][i] = np.sqrt(db_2)
# Saving the results for the reference star.
star_ref = compute_error(star_ref, coefs[0], da_2, db_2)
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['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['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] < 5:
results['star_detection'][i] = 0 # Reference star detected
else:
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)
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)
# 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.
mask = masking(data_clean_sub, transients, header_sub, 150)
seg_sub = source_detection(data_sub,
background_sub,
rms_sub,
n_sigma=mu_conf['det_n_sigma_sub'],
deblend=mu_conf['det_deblend_sub'],
threshold_level=mu_conf['det_treshold_sub'],
npixels=mu_conf['det_npixels_sub'],
contrast=mu_conf['contrast_sub'],
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 5 pixels to have
# an upper limit on the magnitude of the transient.
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_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
if __name__ == "__main__":
sys.exit(main()) # pragma: no cover
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 18 18:11:45 2021
@author: duverne, IJClab, Orsay, duverne@protonmail.com
"""
import argparse
import sys
import os
import warnings
from astropy.io import fits
import hjson
from muphoten.background import sub_background
from muphoten.utils import (list_files, get_filename, make_results_dir,
getTel, load_config)
warnings.simplefilter(action="ignore", category=FutureWarning)
"""Muphoten script for background substraction."""
def main():
"""Muphoten background subtraction."""
telescope_list = getTel()
parser = argparse.ArgumentParser(description="Subtracts background in fits\
images using photutils. Can ave the\
background subtracted images or the\
background image.")
parser.add_argument("--images",
dest="images",
required=True,
type=str,
help="Path where the images are stored. ")
parser.add_argument("--telescope",
dest="telescope",
choices=telescope_list,
required=True,
type=str,
help="Alias for the available telescopes.")
parser.add_argument("--keep-old",
"--keep",
dest="keep",
required=False,
action="store_true",
help="Keep previous results")
parser.add_argument("--skip-processed",
"--skip",
dest="skip",
required=False,
action="store_true",
help="Skip already processed files")
parser.add_argument("--save-png",
dest="save_png",
action="store_true",
help="Save the background images as fits and png.\
Default : False")
args = parser.parse_args()
print('Loading configuration for Muphoten,')
config = load_config(args.telescope)
mu =config['muphoten']['conf']
with open(mu) as json_file:
mu_config = hjson.load(json_file)
# Loading the image to analyse
images = list_files(args.images)[0]
print('Processing {} images'.format(len(images)))
print(images)
for i, image in enumerate(images):
image = make_results_dir(image,
directory='background',
keep=args.keep,
skip=args.skip)
image_name = get_filename(image)
print('Subtracting background in image n° {} : {}'.format(i,
image_name))
hdu = fits.open(image)[0]
data, header = hdu.data, hdu.header
path2save=os.path.dirname(image)
sub_background(data, header, image_name,
n_sigma=mu_config['bkg_n_sigma'],
box_size=mu_config['bkg_box'],
bkg_estimator=mu_config['bkg_estimator'],
filter_size=mu_config['bkg_filter_size'],
path2save=path2save, png=args.save_png)
return 0
if __name__ == "__main__":
sys.exit(main()) # pragma: no cover
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