Commit a763ef7d authored by PAduverne's avatar PAduverne
Browse files

add the option to add the TELESCOPE keyword in the header - useful for KNC...

add the option to add the TELESCOPE keyword in the header - useful for KNC images. fixed soe typo/minor bug in photometry methods
parent 493681b9
......@@ -8,8 +8,7 @@ Created on Sun Nov 15 00:09:43 2020.
import warnings
import numpy as np
from astropy.stats import (SigmaClip, gaussian_fwhm_to_sigma,
sigma_clipped_stats)
from astropy.stats import (gaussian_fwhm_to_sigma)
from astropy.convolution import Gaussian2DKernel
from astropy import units as u
# from astropy.table import vstack, hstack, join
......@@ -17,11 +16,13 @@ from astropy import units as u
from astropy.wcs import WCS, FITSFixedWarning
from photutils import (detect_sources, deblend_sources, source_properties,
SkyEllipticalAperture, aperture_photometry,
SkyCircularAperture, CircularAperture)
SkyCircularAperture, CircularAperture,
EllipticalAperture)
from photutils.utils import calc_total_error
warnings.filterwarnings("ignore", category=FITSFixedWarning)
def masking(image, coordinate, header, mask_size):
"""
Create a square mask around a position.
......@@ -49,16 +50,17 @@ def masking(image, coordinate, header, mask_size):
xs = coord_image[0][0]
ys = coord_image[0][1]
mask = np.ones_like(image, dtype=bool)
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])
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
def source_detection(data, background_array, background_rms, n_sigma=3.0,
deblend=True, mask=None, threshold_level = 2.,
deblend=True, mask=None, threshold_level=2.,
npixels=5, contrast=0.001):
"""
Detect the sources in an image.
......@@ -93,7 +95,7 @@ def source_detection(data, background_array, background_rms, n_sigma=3.0,
A segmentation image containing the detected sources in the image.
"""
#WARNING data MUST NOT BE background subtracted
# WARNING data MUST NOT BE background subtracted
threshold = background_array + (threshold_level*background_rms)
sigma = n_sigma * gaussian_fwhm_to_sigma
kernel = Gaussian2DKernel(sigma, x_size=3, y_size=3)
......@@ -107,12 +109,13 @@ def source_detection(data, background_array, background_rms, n_sigma=3.0,
if deblend:
if segm is not None:
segm = deblend_sources(data, segm, npixels=npixels,
filter_kernel=kernel,
contrast=contrast)
filter_kernel=kernel,
contrast=contrast)
return segm
def Fixed_Photometry(raw_data, sub_data, segmentation, header, radius,
background_array, background_rms, wcs=None):
background_array, background_rms, wcs=None):
"""
"""
......@@ -125,8 +128,8 @@ def Fixed_Photometry(raw_data, sub_data, segmentation, header, radius,
time = 1.0
error = calc_total_error(raw_data,
background_rms,
effective_gain=time)
background_rms,
effective_gain=time)
cat = source_properties(sub_data, segmentation,
background=background_array,
error=error, wcs=wcs)
......@@ -146,6 +149,7 @@ def Fixed_Photometry(raw_data, sub_data, segmentation, header, radius,
return table
def Kron_Photometry(raw_data, sub_data, segmentation, header,
background_array, background_rms, wcs=None):
"""
......@@ -192,7 +196,7 @@ def Kron_Photometry(raw_data, sub_data, segmentation, header,
cat = source_properties(sub_data, segmentation,
background=background_array,
error=error, wcs=wcs,
kron_params=('correct',2.5,0.001,'exact', 5))
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']
......@@ -213,9 +217,10 @@ def Kron_Photometry(raw_data, sub_data, segmentation, header,
return table
def Isophotal_Photometry(raw_data, sub_data, segmentation, header,
background_array, background_rms,
radius=4, wcs=None):
background_array, background_rms,
radius=4, wcs=None):
""" Compute the instrumental magnitude and several other properties of
a fits file. This is done using aperture photometry
......@@ -244,8 +249,8 @@ def Isophotal_Photometry(raw_data, sub_data, segmentation, header,
time = 1.0
error = calc_total_error(raw_data,
background_rms,
effective_gain=time)
background_rms,
effective_gain=time)
cat = source_properties(sub_data, segmentation,
background=background_array,
error=error, wcs=wcs)
......@@ -263,7 +268,7 @@ def Isophotal_Photometry(raw_data, sub_data, segmentation, header,
b=table['semiminor_axis_sigma'][i]*radius,
theta=table['orientation'][i])
apertures.append(elli)
photometry.append(aperture_photometry(sub_data,
photometry.append(aperture_photometry(sub_data,
elli.to_pixel(wcs))['aperture_sum'][0])
table['aperture'] = apertures
table['aperture_sum'] = photometry
......@@ -273,36 +278,42 @@ def Isophotal_Photometry(raw_data, sub_data, segmentation, header,
return table
def get_source(table, coord_source, wcs, col_name='dist'):
coord_source = wcs.all_world2pix(coord_source, 1)
table[col_name] = np.sqrt((table['xcentroid']/u.pix -
coord_source[0][0])**2
+ (table['ycentroid']/u.pix -
coord_source[0][1])**2)
source = table[table[col_name]==min(table[col_name])]
coord_source[0][0])**2
+ (table['ycentroid']/u.pix -
coord_source[0][1])**2)
source = table[table[col_name] == min(table[col_name])]
return source
def Dophot(data, position, wcs, radius=5.0, aperture=None):
if aperture is None:
position = wcs.all_world2pix(position, 1)
aperture = CircularAperture(position, r=radius)
aperture = CircularAperture(position, r=radius)
aper_photo = aperture_photometry(data, aperture)
else:
aperture = aperture.to_sky(wcs).to_pixel(wcs)
if isinstance(aperture,EllipticalAperture):
aperture.to_sky(wcs).to_pixel(wcs)
else:
aperture = aperture.to_pixel(wcs)
aper_photo = aperture_photometry(data, aperture)
return aper_photo['aperture_sum'][0]
def compute_error(table, a, da_2, db_2):
table['poisson_source_error'] = 2.5/(np.log(10)*\
table['poisson_source_error'] = 2.5/(np.log(10) *
np.sqrt(table['aperture_sum']))
table['poisson_bkg_error'] = 2.5/(np.log(10)*\
table['poisson_bkg_error'] = 2.5/(np.log(10) *
np.sqrt(table['background_sum']))
table['calib_error'] = np.sqrt(table['Magnitude']**2 * da_2 +
db_2 +
table['calib_error'] = np.sqrt(table['Magnitude']**2 * da_2 +
db_2 +
a * table['poisson_source_error']**2)
table['total_error'] = np.sqrt(table['poisson_source_error']**2 +
......@@ -311,9 +322,10 @@ def compute_error(table, a, da_2, db_2):
return table
def compute_snr(table):
table['SNR'] = (table['aperture_sum'])/np.sqrt(table['aperture_sum']+
table['SNR'] = (table['aperture_sum'])/np.sqrt(table['aperture_sum'] +
table['background_sum'])
return table
......@@ -68,7 +68,8 @@ def sanitise_headers(filename):
'DATE-OBS',
'BJD-OBS',
'HJD-OBS',
'JD-OBS']
'JD-OBS',
'TELESCOPE']
hdulist = fits.open(filename)
# Verify and try to fix issue with fits standard
# hdulist.verify("fix")
......@@ -175,8 +176,26 @@ def add_filter(fits_file, band):
with fits.open(fits_file, mode = 'update') as hdul:
hdul[0].header.append(('FILTER', band))
def add_observer(filename, tel):
"""
Add the observer in the fits header
Parameters
----------
filename : fits file
File to be modified..
tel : str
Name opf the observer
Returns
-------
None.
"""
with fits.open(filename, mode = 'update') as hdul:
hdul[0].header.append(('TELESCOPË', tel))
def sanitise_fits(filename):
"""Call function to sanitise fits headers and data"""
sanitise_time(filename)
sanitise_headers(filename)
sanitise_data(filename)
\ No newline at end of file
sanitise_data(filename)
......@@ -8,12 +8,13 @@ Created on Wed Feb 17 22:31:08 2021
import os, os.path
import errno
import glob
from astropy.table import Table
from astropy.io import fits
import importlib
import shutil
import time
import numpy as np
from astropy.table import Table
from astropy.io import fits
from astropy.wcs import WCS
def is_subdir(path, basepath):
""" Checks whether the path is inside basepath """
......@@ -103,7 +104,8 @@ def list_files(paths, pattern=["*.fit", "*.fits", "*.fts"],
if is_psf(f, patterns=['_psf.fits', 'SUB',
'sub', 'conv',
'noisemap', 'clean', 'image',
'badpix', 'dupe', 'background', 'results']):
'badpix', 'dupe', 'background',
'results', '_ref', 'old']):
flag_2keep = False
else:
for text in folder2skip:
......@@ -327,11 +329,36 @@ def get_im_res(rep, image_name, prefix='sub_', hdr=False):
image = fits.open(image)[0].data
return image
def get_observer(header):
if 'TELESCOPE' in header.items():
observer = header['TELESCOPE']
else:
observer = 'GRANDMA'
return observer
def set_results_table(col_names, N_rows):
res_table = Table()
res_table['image_number'] = np.arange(N_rows)
res_table['image_index'] = np.arange(N_rows)
for name in col_names:
res_table[name] = None
return res_table
\ No newline at end of file
return res_table
def get_corner_coords(filename):
"""Get the image coordinates of an image"""
header = fits.getheader(filename)
# Get physical coordinates
Naxis1 = header["NAXIS1"]
Naxis2 = header["NAXIS2"]
pix_coords = [[0, 0, Naxis1, Naxis1], [0, Naxis2, Naxis2, 0]]
# Get physical coordinates
w = WCS(header)
ra, dec = w.all_pix2world(pix_coords[0], pix_coords[1], 1)
return [ra, dec]
\ No newline at end of file
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