Commit 62a32692 authored by PAduverne's avatar PAduverne
Browse files

add the main script for photometry. add the calibration scripts and method....

add the main script for photometry. add the calibration scripts and method. Add the saving of the results
parent e29a8483
......@@ -5,7 +5,6 @@ Created on Sun Nov 15 19:48:07 2020.
@author: duverne
"""
"""Muphoten main module for background treatment."""
import os
from photutils import (MeanBackground, MedianBackground,
ModeEstimatorBackground, MMMBackground,
......@@ -15,11 +14,14 @@ from astropy.stats import SigmaClip
from muphoten.utils import createFITS
from muphoten.plotting import plot_background
"""Muphoten main module for background treatment."""
def background_estimation(data, n_sigma=3.0, box_size=(30, 30),
bkg_estimator='sex',
filter_size=(3, 3), mask=None):
"""
Estimation of the background of a fits image with photutils background2D method.
Estimation of the background of a fits image with photutils\
background2D method.
Parameters
----------
......@@ -105,41 +107,84 @@ def subtract(data, bkg):
def sub_background(data, header, image_name, n_sigma=3.0, box_size=(30, 30),
bkg_estimator='sex', filter_size=(3, 3), mask=None,
path2save=None, png=False):
"""
Perform the background subtraction in a fits image.
Possible to save the background and the backgrounds subtracted image
as a fits images.
Parameters
----------
data : numpy array
Image as a array.
header : astropy header object
Header of the image.
image_name : str
name of the image to save the cleaned image and the background.
n_sigma : float, optional
sigma level from the median. The default is 3.0.
box_size : tuple size 2, optional
Size of the boxe in which the background is estimated.
The default is (30, 30).
bkg_estimator : str, optional
Estimator name. The default is 'sex'.
filter_size : tuple size 2, optional
Size of the filter. The default is (3, 3).
mask : boolean array, optional
Array used to mask pixels if some regions
of the image have to be excludes from
the estimation of the background.
The default is None.
path2save : str, optional
Path where the are saved. The default is None.
png : Bool, optional
Used to save the background as a png file. The default is False.
Returns
-------
data_cleaned : numpy array
Background subtracted image.
background_array : numpy array
Background array.
background_rms : numpy array
2D Standard devitation of the image.
"""
# Computing the background with photutils
bkg_obj = background_estimation(data,
n_sigma = n_sigma,
box_size = box_size,
bkg_estimator=bkg_estimator,
filter_size=(3, 3), mask=mask)
# Convert the Background object into a numpy array
background_array, background_rms = bkg_to_array(bkg_obj, rms=True)
# Subtracting the background to the image
data_cleaned = data - background_array
# Save clean image if required
if path2save is not None:
clean_filename = os.path.join(path2save,
'clean_' + image_name + '.fits')
createFITS(data_cleaned, header, filename=clean_filename)
# Save background image
background_filename = os.path.join(path2save,
'background/background_' +
image_name)
createFITS(background_array, header,
filename=background_filename+'.fits')
# Save RMS image
RMS_filename = os.path.join(path2save,
'background/rms_' + image_name)
createFITS(background_rms, header,
filename=RMS_filename+'.fits')
if png:
plot_background(data, background_array, dif=False,
save=True, path=background_filename+'.png')
return data_cleaned, background_array, background_rms
# Computing the background with photutils
bkg_obj = background_estimation(data,
n_sigma = n_sigma,
box_size = box_size,
bkg_estimator=bkg_estimator,
filter_size=filter_size,
mask=mask)
# Convert the Background object into a numpy array
background_array, background_rms = bkg_to_array(bkg_obj, rms=True)
# Subtracting the background to the image
data_cleaned = data - background_array
# Save clean image if required
if path2save is not None:
clean_filename = os.path.join(path2save,
'clean_' + image_name + '.fits')
createFITS(data_cleaned, header, filename=clean_filename)
# Save background image
background_filename = os.path.join(path2save,
'background/background_' +
image_name)
createFITS(background_array, header,
filename=background_filename+'.fits')
# Save RMS image
RMS_filename = os.path.join(path2save,
'background/rms_' + image_name)
createFITS(background_rms, header,
filename=RMS_filename+'.fits')
if png:
plot_background(data, background_array, dif=False,
save=True, path=background_filename+'.png')
return data_cleaned, background_array, background_rms
This diff is collapsed.
......@@ -6,31 +6,38 @@ Created on Mon Apr 19 16:58:12 2021
@author: duverne, IJClab, Orsay, duverne@protonmail.com
"""
"""Muphoten main script for photometry."""
import argparse
import sys
import os
import time
import warnings
from astropy.io import fits
from astropy.wcs import WCS, FITSFixedWarning
import time
from astropy.table import Table
import hjson
warnings.filterwarnings("ignore", category=FITSFixedWarning)
warnings.simplefilter(action="ignore", category=FutureWarning)
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,
getCata, CataList, list_files, load_config,
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, catalog_choice, xmatch, joining,
source_coordinates, from_PS2Johnson)
from muphoten.photometry import (source_detection, Kron_Photometry,
Fixed_Photometry, Isophotal_Photometry)
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)
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."""
......@@ -59,24 +66,24 @@ def main():
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=False,
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",
......@@ -102,7 +109,7 @@ def main():
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
......@@ -112,25 +119,57 @@ def main():
print('Loading configuration for Muphoten,')
config = load_config(args.telescope)
mu =config['muphoten']['conf']
mu = config['muphoten']['conf']
with open(mu) as json_file:
mu_config = 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()
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_magnitude',
'transient_calibrated_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)
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')
......@@ -143,13 +182,16 @@ def main():
# Getting different useful informations : band, MJD,
# catalog for calibration.
mjd = header.get('MJD-OBS')
results['time'][i] = header.get('MJD-OBS')
if args.catalog:
catalog = getCata(args.catalog)
catalog, system = getCata(args.catalog)
cata_name = args.catalog
band = get_filter(header)
else:
catalog = catalog_choice(header)
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()
......@@ -164,10 +206,10 @@ def main():
data_clean, background,\
rms = 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'],
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)
......@@ -184,21 +226,22 @@ def main():
start_photo_time = time.time()
# Detection of the sources & segmentation
segmentation = source_detection(data, background, rms,
n_sigma=mu_config['det_n_sigma'],
deblend=mu_config['det_deblend'],
threshold_level=mu_config['det_treshold'],
npixels=mu_config['det_npixels'],
contrast=mu_config['contrast'])
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, background)
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]:
......@@ -211,19 +254,106 @@ def main():
wcs=wcs_data, path2save=calib_rep)
if args.photo == phot_list[2]:
print('Performing Fixed size Photometry')
radius = mu_config['fixed_factor'] * FWHM_psf[0]
radius = mu_conf['fixed_factor'] * FWHM_psf[0]
photometry_table = Fixed_Photometry(data, data_clean, segmentation,
header, radius, background,
header, radius, background,
rms, wcs=wcs_data)
plot_image(data_clean, sky_ap=photometry_table['aperture'],
wcs=wcs_data, path2save=calib_rep)
coord_sources = source_coordinates(photometry_table)
crossmatch = xmatch(coord_sources, catalog, 5.0)
crossmatch = from_PS2Johnson(band, crossmatch)
photometry_table = joining(photometry_table, crossmatch)
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['starresults.show_in_browser(jsviewer=True)_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_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
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:
results['transient_detection'][i] = 1
# Computing the number of counts in a aperture of % pixels to have
# an upper limit on the magnitude of the transient.
flux = Dophot(data_sub,
transients,
wcs_sub,
radius=5.0)
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 photmetry'
# raise ValueError('{} is not an avalaible photometry'
# 'in Muphoten'.format(args.photo))
return 0
......
......@@ -6,7 +6,6 @@ Created on Sun Apr 18 18:11:45 2021
@author: duverne, IJClab, Orsay, duverne@protonmail.com
"""
"""Muphoten script for background substraction."""
import argparse
import sys
import os
......@@ -14,11 +13,13 @@ import warnings
from astropy.io import fits
import hjson
warnings.simplefilter(action="ignore", category=FutureWarning)
from muphoten.background import sub_background
from muphoten.utils import (list_files, get_filename, make_results_dir,
getTel, load_config, mkdir_p)
getTel, load_config)
warnings.simplefilter(action="ignore", category=FutureWarning)
"""Muphoten script for background substraction."""
def main():
"""Muphoten background subtraction."""
......@@ -35,7 +36,7 @@ def main():
required=True,
type=str,
help="Path where the images are stored. ")
parser.add_argument("--telescope",
dest="telescope",
choices=telescope_list,
......@@ -56,7 +57,7 @@ def main():
required=False,
action="store_true",
help="Skip already processed files")
parser.add_argument("--save-png",
dest="save_png",
action="store_true",
......@@ -69,7 +70,7 @@ def main():
config = load_config(args.telescope)
mu =config['muphoten']['conf']
with open(mu) as json_file:
mu_config = hjson.load(json_file)
mu_config = hjson.load(json_file)
# Loading the image to analyse
images = list_files(args.images)[0]
......
......@@ -6,11 +6,12 @@ Created on Sun Nov 15 00:09:43 2020.
@author: duverne
"""
import numpy as np
import warnings
import numpy as np
from astropy.stats import (SigmaClip, gaussian_fwhm_to_sigma,
sigma_clipped_stats)
from astropy.convolution import Gaussian2DKernel
from astropy import units as u
# from astropy.table import vstack, hstack, join
# from astropy import units as u, coordinates as coord
from astropy.wcs import WCS, FITSFixedWarning
......@@ -18,11 +19,10 @@ from photutils import (detect_sources, deblend_sources, source_properties,
SkyEllipticalAperture, aperture_photometry,
SkyCircularAperture, CircularAperture)
from photutils.utils import calc_total_error
from astropy import units as u
warnings.filterwarnings("ignore", category=FITSFixedWarning)
def masking_the_source(image, coordinate, header, mask_size):
def masking(image, coordinate, header, mask_size):
"""
Create a square mask around a position.
......@@ -44,12 +44,16 @@ def masking_the_source(image, coordinate, header, mask_size):
"""
w = WCS(header)
coord_image = w.all_world2pix([coordinate], 1)
coord_image = w.all_world2pix(coordinate, 1)
shape = np.shape(image)
xs = coord_image[0][0]
ys = coord_image[0][1]
mask = np.ones_like(image, dtype=bool)
mask[(int(ys) - int(mask_size)):(int(ys) + int(mask_size)),
(int(xs) - int(mask_size)):(int(xs) + int(mask_size))] = False
yo=max(int(ys) - int(mask_size), 0)
y1=min(int(ys) + int(mask_size), shape[1])
xo=max(int(xs) - int(mask_size), 0)
x1=min(int(xs) + int(mask_size), shape[0])
mask[yo:y1, xo:x1] = False
return mask
......@@ -101,12 +105,10 @@ def source_detection(data, background_array, background_rms, n_sigma=3.0,
segm = detect_sources(data, threshold, npixels=npixels,
filter_kernel=kernel)
if deblend:
segm_deblend = deblend_sources(data, segm, npixels=npixels,
if segm is not None:
segm = deblend_sources(data, segm, npixels=npixels,
filter_kernel=kernel,
contrast=contrast)
return segm_deblend
return segm
def Fixed_Photometry(raw_data, sub_data, segmentation, header, radius,
......@@ -135,7 +137,7 @@ def Fixed_Photometry(raw_data, sub_data, segmentation, header, radius,
table['ycentroid'].info.format = '.4f'
apertures = SkyCircularAperture(table['sky_centroid'], r=radius*u.pix)
photometry = aperture_photometry(sub_data, apertures.to_pixel(wcs))
table['aperture'] = apertures
table['aperture_sum'] = photometry['aperture_sum']
table['Magnitude'] = -2.5 * np.log10(table['aperture_sum'])
......@@ -170,7 +172,8 @@ def Kron_Photometry(raw_data, sub_data, segmentation, header,
Returns
-------
table : astropy table
Astropy table with all the informations about the sources detected in data.
Astropy table with all the informations about the
sources detected in data.
"""
# WARNING data MUST NOT BE BG substracted here
......@@ -192,8 +195,7 @@ def Kron_Photometry(raw_data, sub_data, segmentation, header,
kron_params=('correct',2.5,0.001,'exact', 5))
columns = ['id', 'xcentroid', 'ycentroid', 'source_sum', 'area',
'background_at_centroid', 'background_mean', 'background_sum',
'sky_centroid', 'kron_aperture', 'kron_flux', 'kron_radius'
]
'sky_centroid', 'kron_aperture', 'kron_flux', 'kron_radius']
table = cat.to_table(columns=columns)
table['xcentroid'].info.format = '.4f' # optional format
table['ycentroid'].info.format = '.4f'
......@@ -202,8 +204,9 @@ def Kron_Photometry(raw_data, sub_data, segmentation, header,