From c8ae670ea7a6bdd4d77972960a915f4180dcd273 Mon Sep 17 00:00:00 2001 From: Yael-II <154762804+Yael-II@users.noreply.github.com> Date: Tue, 3 Dec 2024 23:05:59 +0100 Subject: [PATCH] Begining of data reduction --- Source/FIDR.py | 144 +++++++++++++++++++++++++++++ Source/astrobs_v1_toolbox.py | 9 +- Source/data_reduction_tools.py | 159 +++++++++++++++++++++++++++++++++ requirements.txt | 8 ++ 4 files changed, 317 insertions(+), 3 deletions(-) create mode 100644 Source/FIDR.py create mode 100644 Source/data_reduction_tools.py diff --git a/Source/FIDR.py b/Source/FIDR.py new file mode 100644 index 0000000..e036ed9 --- /dev/null +++ b/Source/FIDR.py @@ -0,0 +1,144 @@ +import os +import numpy as np +import data_reduction_tools as tools +import astroalign as aa + +TYPE = ["image", "spectroscopy"][1] +target = "NGC_40" +IN_DIR = "./Input/" +OUT_DIR = "./Output/" +LIGHT_DIR = IN_DIR + "Light/" +DARK_DIR = IN_DIR + "Dark/" +FLAT_DIR = IN_DIR + "Flat/" +BIAS_DIR = IN_DIR + "Bias/" +THAR_DIR = IN_DIR + "ThAr/" + +def spectro_reduction(): + # BIAS + filelist = tools.list_fits(BIAS_DIR) + tools.stack(filelist, BIAS_DIR, "master_bias.fits") + + # DARK + #filelist = tools.list_fits(DARK_DIR) + #tools.debias(filelist, DARK_DIR) + #filelist = ["debiased_" + fname for fname in filelist \ + # if not "debiased_" in fname] + #tools.stack(filelist, DARK_DIR, "master_debiased_dark.fits") + + # FLAT + filelist = tools.list_fits(FLAT_DIR) + tools.debias(filelist, FLAT_DIR) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.normalize_spectrum_flat(filelist) + filelist = ["normalized_" + fname for fname in filelist \ + if not "normalized_" in fname] + tools.stack(filelist, FLAT_DIR, "master_normalized_debiased_flat.fits") + + # LIGHT + filelist = tools.list_fits(LIGHT_DIR) + tools.debias(filelist, LIGHT_DIR) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.stack(filelist, LIGHT_DIR, "master_debiased_light.fits") + tools.reduction_operation(target + ".fits", + master_light = "master_debiased_light.fits", + master_flat = "master_normalized_debiased_flat.fits") + # MAIN + # tools.plot_spectrum("master_debiased_flat.fits", FLAT_DIR) + # tools.plot_spectrum("reduced_" + target + ".fits") + return 0 + +def spectro_calibration(): + # BIAS + filelist = tools.list_fits(BIAS_DIR) + tools.stack(filelist, BIAS_DIR, "master_bias.fits") + + # DARK + #filelist = tools.list_fits(DARK_DIR) + #tools.debias(filelist, DARK_DIR) + #filelist = ["debiased_" + fname for fname in filelist \ + # if not "debiased_" in fname] + #tools.stack(filelist, DARK_DIR, "master_debiased_dark.fits") + + # FLAT + filelist = tools.list_fits(FLAT_DIR) + tools.debias(filelist, FLAT_DIR) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.normalize_spectrum_flat(filelist) + filelist = ["normalized_" + fname for fname in filelist \ + if not "normalized_" in fname] + tools.stack(filelist, FLAT_DIR, "master_normalized_debiased_flat.fits") + + # ThAr + filelist = tools.list_fits(THAR_DIR) + tools.debias(filelist, THAR_DIR) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.stack(filelist, THAR_DIR, "master_debiased_ThAr.fits") + tools.reduction_operation("ThAr.fits", + light_directory = THAR_DIR, + master_light = "master_debiased_ThAr.fits", + master_flat = "master_normalized_debiased_flat.fits") + # MAIN + tools.plot_spectrum("reduced_ThAr.fits") + return 0 + +def image_reduction(): + # BIAS + filelist = tools.list_fits(BIAS_DIR) + tools.stack(filelist, BIAS_DIR, "master_bias.fits") + + # DARK + filelist = tools.list_fits(DARK_DIR) + tools.debias(filelist, DARK_DIR) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.stack(filelist, DARK_DIR, "master_debiased_dark.fits") + + # FLAT + filter_list = [d for d in os.listdir(LIGHT_DIR) if not d[0] == "."] + for filter_ in filter_list: + filter_dir = filter_ + "/" + filelist = tools.list_fits(FLAT_DIR + filter_dir) + tools.debias(filelist, FLAT_DIR + filter_dir) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.normalize_image_flat(filelist, FLAT_DIR + filter_dir) + filelist = ["normalized_" + fname for fname in filelist \ + if not "normalized_" in fname] + tools.stack(filelist, FLAT_DIR + filter_dir, + "master_normalized_debiased_flat_{}.fits".format(filter_)) + + # LIGHT + for filter_ in filter_list: + filter_dir = filter_ + "/" + filelist = tools.list_fits(LIGHT_DIR + filter_dir) + tools.debias(filelist, LIGHT_DIR + filter_dir) + filelist = ["debiased_" + fname for fname in filelist \ + if not "debiased_" in fname] + tools.stack(filelist, + LIGHT_DIR + filter_dir, + "master_debiased_light_{}.fits".format(filter_)) + tools.reduction_operation(target + "_" + filter_ + ".fits", + master_light = "master_debiased_light_{}.fits".format(filter_), + light_directory = LIGHT_DIR + filter_dir, + master_dark = "master_debiased_dark.fits", + master_flat = "master_normalized_" \ + + "debiased_flat_{}.fits".format(filter_), + flat_directory = FLAT_DIR + filter_dir) + + + return 0 + +def main(): + if TYPE == "image": + image_reduction() + if TYPE == "spectroscopy": + spectro_reduction() + spectro_calibration() + return 0 + +if __name__ == "__main__": + main() diff --git a/Source/astrobs_v1_toolbox.py b/Source/astrobs_v1_toolbox.py index c1808e9..20d66bb 100644 --- a/Source/astrobs_v1_toolbox.py +++ b/Source/astrobs_v1_toolbox.py @@ -347,7 +347,7 @@ def add_simbad(table: Table, ra = coord.Angle(obj["ra"][i], unit=u.degree) dec = coord.Angle(obj["dec"][i], unit=u.degree) except KeyError: - ra = coord.Angle(obj["RA"][i], unit=u.degree) + ra = coord.Angle(obj["RA"][i], unit=u.hourangle) dec = coord.Angle(obj["DEC"][i], unit=u.degree) ra_str = ra.to_string(unit=u.hourangle, sep=":", @@ -727,6 +727,8 @@ def make_plot(table: Table, ra = coord.Angle(table["ra"][i], unit=u.hourangle).degree if ra > window_west: ra -= 360 + if ra < window_east: + ra += 360 ra_before = ra + a_before ra_after = ra - a_after @@ -748,12 +750,13 @@ def make_plot(table: Table, horizontalalignment="right", verticalalignment="center") ra_text = table["ra"][i] dec_text = table["dec"][i] - ax.text(time_before, i, " {}{}".format(ra_text, dec_text)) + ax.text(time_before, i, " {}{}".format(ra_text, dec_text), + horizontalalignment="left", verticalalignment="center") ax.set_xlabel("Observation date and time (UTC)") ax.set_yticks(range(N), table["seq"]) - plt.show(block=False) + plt.show(block=True) # obs_time = np.arange(sun_set, sun_rise, dtype="datetime64[m]") return 0 diff --git a/Source/data_reduction_tools.py b/Source/data_reduction_tools.py new file mode 100644 index 0000000..8305edf --- /dev/null +++ b/Source/data_reduction_tools.py @@ -0,0 +1,159 @@ +import os +import numpy as np +import matplotlib.pyplot as plt +from astropy.io import fits + + +IN_DIR = "./Input/" +OUT_DIR = "./Output/" +LIGHT_DIR = IN_DIR + "Light/" +DARK_DIR = IN_DIR + "Dark/" +FLAT_DIR = IN_DIR + "Flat/" +BIAS_DIR = IN_DIR + "Bias/" +THAR_DIR = IN_DIR + "ThAr/" + +PLT_STYLE = "YII_light_1" + +def list_fits(directory: str = LIGHT_DIR): + """ + Returns the list of files in the directory + @params: + - directory: the directory containing the fits + @output: + - list of fits files in the directory + """ + filelist = [fname for fname in os.listdir(directory) if ".fit" in fname] + # works for fit and fits + return filelist + +def stack(filelist: list, + in_directory: str, + out_name: str, + out_directory: str = None, + method:str = "median"): + in_data = [] + if out_directory == None: + out_directory = in_directory + N = len(filelist) + for file in filelist: + if not "master_" in file: + in_data.append(fits.getdata(in_directory + file)) + in_data = np.array(in_data) + if method == "mean": + out_data = np.mean(in_data, axis=0) + if method == "median": + out_data = np.median(in_data, axis=0) + header = fits.getheader(in_directory+file) + header["history"] = "stacking with {} files".format(N) + fits.writeto(out_directory \ + + out_name, out_data, header, overwrite=True) + return 0 + +def debias(filelist: list, + in_directory: str, + master_bias: str = "master_bias.fits", + bias_directory: str = BIAS_DIR, + out_directory: str = None): + if out_directory == None: + out_directory = in_directory + for file in filelist: + if not "debiased_" in file: + in_data = fits.getdata(in_directory + file) + header = fits.getheader(in_directory + file) + bias = fits.getdata(bias_directory+master_bias) + out_data = in_data - bias + header["history"] = "debiased with master bias" + fits.writeto(out_directory + "debiased_" + file, \ + out_data, header, \ + overwrite=True) + return 0 + +def normalize_image_flat(filelist: str, + in_directory: str = FLAT_DIR, + out_directory: str = None): + if out_directory == None: + out_directory = in_directory + + for file in filelist: + if not "normalized_" in file: + data = fits.getdata(in_directory + file) + head = fits.getheader(in_directory + file) + normalized_data = data/np.median(data) + fits.writeto(out_directory + "normalized_" + file, + normalized_data, overwrite=True) + return 0 + + +def normalize_spectrum_flat(filelist: str, + in_directory: str = FLAT_DIR, + out_directory: str = None, + fit_degree: int = 3, + d_lim: int = 10000): + if out_directory == None: + out_directory = in_directory + + poly = np.polynomial.chebyshev + for file in filelist: + if not "normalized_" in file: + data = fits.getdata(in_directory + file).flatten() + head = fits.getheader(in_directory + file) + pixel = list(range(len(data))) + mask = np.ones_like(data) + mask[np.argwhere(data < d_lim)] = 0 + fit_coefs = poly.chebfit(pixel, data, fit_degree, w=mask) + fit_data = poly.chebval(pixel, fit_coefs) + normalized_data = data/fit_data + normalized_data[np.argwhere(normalized_data<0.5)] = np.nan + fits.writeto(out_directory + "normalized_" + file, + normalized_data, overwrite=True) + return 0 + +def reduction_operation(target: str, + master_light: list, + light_directory: str = LIGHT_DIR, + master_dark: str = None, + dark_directory: str = DARK_DIR, + master_flat: str = None, + flat_directory: str = FLAT_DIR, + out_directory: str = OUT_DIR): + + data = fits.getdata(light_directory + master_light) + head = fits.getheader(light_directory + master_light) + if master_dark != None: + dark = fits.getdata(dark_directory + master_dark) + data = data - dark + head["history"] = "removed dark with {}".format(master_dark) + if master_flat != None: + flat = fits.getdata(flat_directory + master_flat) + data = data / flat + head["history"] = "flatten with {}".format(master_flat) + fits.writeto(out_directory + "reduced_" + target, + data, head, overwrite=True) + + + return 0 + +def plot_spectrum(reduced: str, + directory: str = OUT_DIR): + if PLT_STYLE in plt.style.available: plt.style.use(PLT_STYLE) + data = fits.getdata(directory + reduced).flatten() + head = fits.getheader(directory + reduced) + pixel = list(range(len(data))) + + plt.plot(pixel, data) + plt.show(block=True) + return 0 + +def find_peaks_highest(data: list, + N: int = 10): + """Find the N highest peaks in the data""" + peaks = [] + mask = np. + for i in range(N): + None + return 0 + + + + + diff --git a/requirements.txt b/requirements.txt index 2bb6e43..707dcf8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ +astroalign==2.6.1 astropy==7.0.0 astropy-iers-data==0.2024.11.25.0.34.48 astroquery==0.4.7 @@ -9,13 +10,16 @@ cycler==0.12.1 fonttools==4.55.0 html5lib==1.1 idna==3.10 +imageio==2.36.1 jaraco.classes==3.4.0 jaraco.context==6.0.1 jaraco.functools==4.1.0 keyring==25.5.0 kiwisolver==1.4.7 +lazy_loader==0.4 matplotlib==3.9.3 more-itertools==10.5.0 +networkx==3.4.2 numpy==2.1.3 packaging==24.2 pillow==11.0.0 @@ -25,7 +29,11 @@ python-dateutil==2.9.0.post0 pyvo==1.6 PyYAML==6.0.2 requests==2.32.3 +scikit-image==0.24.0 +scipy==1.14.1 +sep-pjw==1.3.7 six==1.16.0 soupsieve==2.6 +tifffile==2024.9.20 urllib3==2.2.3 webencodings==0.5.1