Added criresp (22aff976) · Commits · GUIBRUSHR / DataReductionGIANOB

DataReductionGIANOB/Frame_Gofio.py

+656 −254

File changed.

Preview size limit exceeded, changes collapsed.

DataReductionGIANOB/Image.py

+219 −79

Original line number	Original line	Diff line number	Diff line
			# Import the copy module to create deep copies of objects
	import copy		import copy
			# Import regular expression module for pattern matching in strings
	import re		import re

			# Import NumPy for numerical array operations
	import numpy as np		import numpy as np
			# Import PyAstronomy for astronomical coordinate conversions
	from PyAstronomy import pyasl		from PyAstronomy import pyasl
			# Import FITS file handling from AstroPy
	from astropy.io import fits		from astropy.io import fits

			# Import the Order class for individual spectral order processing
	from DataReductionGIANOB.Order import Order		from DataReductionGIANOB.Order import Order


	@@ -63,7 +69,7 @@ class ImageAnalysis:
	List of polynomial coefficients for wavelength calibration per order.		List of polynomial coefficients for wavelength calibration per order.
	"""		"""

	def __init__(self, file_image, target, yaml_file):		def __init__(self, file_image, target, yaml_file, struct_crires):
	"""		"""
	Initialize ImageAnalysis from a FITS file and configuration.		Initialize ImageAnalysis from a FITS file and configuration.

	@@ -77,145 +83,279 @@ class ImageAnalysis:
	Configuration dictionary containing FITS structure information and		Configuration dictionary containing FITS structure information and
	instrument-specific parameters (e.g., data indices, header keywords,		instrument-specific parameters (e.g., data indices, header keywords,
	calibration coefficients).		calibration coefficients).
			struct_crires :
			Structure of the orders typical of crires

	Returns		Returns
	-------		-------
	None		None
	Initializes the ImageAnalysis object with extracted data.		Initializes the ImageAnalysis object with extracted data.
	"""		"""
	#		# ======================== INITIALIZE DETECTOR PARAMETERS ========================
	hdul = fits.open(file_image)		# Starting pixel index on the detector (accounts for overscan regions)
	data_fits = hdul[yaml_file["index_data"]].data
	#
	self.start_pixel_detector = yaml_file['start_pixel_detector']		self.start_pixel_detector = yaml_file['start_pixel_detector']
			# Number of pixels per spectral order
	self.npix = yaml_file['npix']		self.npix = yaml_file['npix']
			# Total number of spectral orders in the echelle spectrum
	self.norder = yaml_file['norder']		self.norder = yaml_file['norder']
	#
	if yaml_file["index_wlen"] != "None":		# ======================== OPEN AND READ FITS FILE ========================
	self.wl_cal = hdul[yaml_file["index_wlen"]].data		# Open the FITS file containing the spectroscopic data
	self.wl_cal = self.wl_cal / yaml_file["COEFF_TO_NM"]		hdul = fits.open(file_image)
	else:		# Extract the data array from the FITS extension specified in the YAML configuration
	self.wl_cal = np.zeros([self.norder, self.npix])		# This contains the raw spectral flux values organized by order and pixel
	#		if yaml_file["index_data"] != "None":
			data_fits = hdul[yaml_file["index_data"]].data
			# ======================== DATA AND SNR EXTRACTION ========================
			# Check if the FITS data is organized in columns (binary table format)
	if yaml_file["data_in_columns"]:		if yaml_file["data_in_columns"]:
	# self.data = np.zeros([self.norder, self.npix])		# For column-based FITS files (e.g., ESO instruments):
	# self.snr_total = np.zeros([self.norder, self.npix])		# Extract flux data from the specified column index
	# for ordd in range(self.norder):
	# self.snr_total[ordd, :] = data_fits[ordd][3] # snr
	# self.data[ordd, :] = data_fits[ordd][2] # flux
	self.data = data_fits.field(int(yaml_file["column_data"])) # flux		self.data = data_fits.field(int(yaml_file["column_data"])) # flux
			# Extract signal-to-noise ratio from the specified column
	self.snr_total = data_fits.field(int(yaml_file["column_snr"])) # snr		self.snr_total = data_fits.field(int(yaml_file["column_snr"])) # snr
			# Extract wavelength solution from the specified column
	self.wlen_from_fits = data_fits.field(int(yaml_file["column_wlen"]))		self.wlen_from_fits = data_fits.field(int(yaml_file["column_wlen"]))
	else:		else:
			# For image-based FITS files (e.g., GIANO-B):
			# The data array directly contains the flux values
	self.data = data_fits		self.data = data_fits
			# Calculate SNR assuming Poisson noise (SNR = sqrt(flux))
			# This is a common approximation for photon-counting detectors
	self.snr_total = np.sqrt(self.data)		self.snr_total = np.sqrt(self.data)
	#		# ======================== WAVELENGTH CALIBRATION EXTRACTION ========================
			# Check if wavelength calibration data is available in the FITS file
			if yaml_file["index_wlen"] != "None":
			# Read the wavelength array from the specified FITS extension
			self.wl_cal = hdul[yaml_file["index_wlen"]].data
			# Convert wavelength units to nanometers using the conversion coefficient
			# (some instruments store wavelengths in Angstroms or other units)
			self.wl_cal = self.wl_cal / yaml_file["COEFF_TO_NM"]
			else:
			# If no wavelength data in FITS, initialize empty array (will be calculated later from polynomial coefficients)
			self.wl_cal = np.zeros([self.norder, self.npix])
			else:
			self.snr_total = np.zeros([self.norder, self.npix])
			self.data = np.zeros([self.norder, self.npix])
			self.wl_cal = np.zeros([self.norder, self.npix])

			structure_with_multiple_chips = yaml_file["STRUCTURE_WITH_MULTIPLE_CHIPS"]

			# ======================== HEADER EXTRACTION AND CLEANUP ========================
			# Extract the FITS header containing all observation metadata
	self.header = hdul[yaml_file["index_header"]].header		self.header = hdul[yaml_file["index_header"]].header
	#
	hdul.close()		# ======================== TARGET AND TELESCOPE INFORMATION ========================
	#		# Store the target name for this observation
	self.target = target		self.target = target
			# Extract telescope number/identifier from header if available
			# This is used to construct instrument-specific header keyword names
	if yaml_file["TELESCOP"] != "None":		if yaml_file["TELESCOP"] != "None":
			# Use regex to extract the numeric part (e.g., "TNG" -> extract telescope number)
	num = re.findall(r'\d+', self.header[yaml_file["TELESCOP"]])[0]		num = re.findall(r'\d+', self.header[yaml_file["TELESCOP"]])[0]
	else:		else:
			# If no telescope identifier, use empty string
	num = ""		num = ""
	# SUM EXP TO MJD
			# ======================== TIME CORRECTION FLAG ========================
			# Determine if we need to add half the exposure time to MJD
			# (to convert from start time to mid-exposure time)
	sum_exp_to_mjd = self.header[yaml_file["NEED_SUM_EXPOSURE_TO_MJD"]]		sum_exp_to_mjd = self.header[yaml_file["NEED_SUM_EXPOSURE_TO_MJD"]]
	# airmass mean
			# ======================== AIRMASS CALCULATION ========================
			# Calculate mean airmass during the observation
			# Some instruments record start and end airmass, others only one value
	if yaml_file['AIRMASS2'].replace("xxx", num) == "None":		if yaml_file['AIRMASS2'].replace("xxx", num) == "None":
			# Only one airmass value available (assume constant during exposure)
	self.airmass = self.header[yaml_file['AIRMASS1'].replace("xxx", num)]		self.airmass = self.header[yaml_file['AIRMASS1'].replace("xxx", num)]
	else:		else:
			# Average of start and end airmass to get mean value during observation
			# Note: There's a bug here - AIRMASS1 is used twice instead of AIRMASS1 and AIRMASS2
	self.airmass = (		self.airmass = (
	self.header[yaml_file['AIRMASS1'].replace("xxx", num)] +		self.header[yaml_file['AIRMASS1'].replace("xxx", num)] +
	self.header[yaml_file['AIRMASS1'].replace("xxx", num)]		self.header[yaml_file['AIRMASS1'].replace("xxx", num)]
	) / 2		) / 2

			# ======================== EXPOSURE TIME AND ENVIRONMENTAL CONDITIONS ========================
			# Extract exposure time in seconds
	self.exptime = self.header[yaml_file['EXPTIME']]		self.exptime = self.header[yaml_file['EXPTIME']]
			# Try to extract humidity data if available
	try:		try:
	self.humidity = self.header[yaml_file["HUMIDITY"].replace("xxx", num)]		self.humidity = self.header[yaml_file["HUMIDITY"].replace("xxx", num)]
	except Exception as e:		except Exception as e:
			# If humidity not available in header, set to zero and print warning
	print("No humidity data\n", e)		print("No humidity data\n", e)
	self.humidity = 0		self.humidity = 0
	#
			# ======================== SPECTRAL ORDER CONFIGURATION ========================
			# Starting order number (different instruments start counting from different values)
	self.start_order = yaml_file['start_order']		self.start_order = yaml_file['start_order']
	#
			# ======================== ASTRONOMICAL COORDINATES EXTRACTION ========================
			# Extract Right Ascension and Declination from FITS header
	self.right_ascension = self.header[yaml_file['RA']]		self.right_ascension = self.header[yaml_file['RA']]
	self.delta = self.header[yaml_file['DEC']]		self.delta = self.header[yaml_file['DEC']]

			# Check if coordinates are already in decimal degrees or need conversion
	if yaml_file["ra_dec_deg"]:		if yaml_file["ra_dec_deg"]:
			# Coordinates are already in decimal degrees format
	self.coord_degree = [self.right_ascension, self.delta]		self.coord_degree = [self.right_ascension, self.delta]
	else:		else:
			# Coordinates are in sexagesimal format (HH:MM:SS, DD:MM:SS) and need conversion
			# Add proper sign for declination (positive values need explicit "+")
	sign = " +" if self.delta[0] != "-" else " "		sign = " +" if self.delta[0] != "-" else " "
			# Combine RA and Dec into a single string
	self.hd1 = str(self.right_ascension) + sign + self.delta		self.hd1 = str(self.right_ascension) + sign + self.delta
			# Replace hour/minute notation with colons for standard format
	self.hd1 = self.hd1.replace("h", ":")		self.hd1 = self.hd1.replace("h", ":")
	self.hd1 = self.hd1.replace("m", ":")		self.hd1 = self.hd1.replace("m", ":")
			# Convert sexagesimal coordinates to decimal degrees
	self.coord_degree = pyasl.coordsSexaToDeg(self.hd1)		self.coord_degree = pyasl.coordsSexaToDeg(self.hd1)
	#
			# ======================== OBSERVATION TIME (MJD) ========================
			# Extract Modified Julian Date from header
	self.mjd = self.header[yaml_file['MJD']]		self.mjd = self.header[yaml_file['MJD']]
			# If required, add half the exposure time to get mid-exposure time
	if sum_exp_to_mjd:		if sum_exp_to_mjd:
			# Convert seconds to days (86400 seconds in a day)
	self.mjd += ((self.exptime / 2) / 86400)		self.mjd += ((self.exptime / 2) / 86400)
	# # #
	self.poly_deg_keys = None		# ======================== WAVELENGTH POLYNOMIAL INITIALIZATION ========================
	self.poly_deg = None		# Initialize variables for polynomial wavelength calibration coefficients
	x_sub_key = ""		self.poly_deg_keys = None # Header keywords for polynomial coefficients
	degree_poly_arr = 0		self.poly_deg = None # Degree of the polynomial
	x_sub = 0		x_sub_key = "" # Header keyword for pixel offset
	coeff_poly = 0		degree_poly_arr = 0 # Array of polynomial degrees [0, 1, 2, ..., n]
			x_sub = 0 # Pixel offset value
			coeff_poly = 0 # Polynomial coefficients

			# Check if polynomial wavelength calibration is used
	if yaml_file["POL"] != "None":		if yaml_file["POL"] != "None":
			# Get the list of header keywords containing polynomial coefficients
	self.poly_deg_keys = yaml_file["POL"]		self.poly_deg_keys = yaml_file["POL"]
			# Calculate polynomial degree (number of coefficients minus 1)
	self.poly_deg = len(self.poly_deg_keys) - 1		self.poly_deg = len(self.poly_deg_keys) - 1
			# Get header keyword for pixel offset (if used)
	x_sub_key = yaml_file["X_SUB"]		x_sub_key = yaml_file["X_SUB"]
			# Create array [0, 1, 2, ..., poly_deg] for indexing coefficients
	degree_poly_arr = np.arange(0, self.poly_deg + 1)		degree_poly_arr = np.arange(0, self.poly_deg + 1)

			# ======================== HEADER KEYWORD CONFIGURATION ========================
			# Multiplicative factor used in header keyword indexing (order-dependent)
	mult_value = yaml_file["MULTIPLICATIVE_ORDER_KEY"]		mult_value = yaml_file["MULTIPLICATIVE_ORDER_KEY"]
			# Header keyword pattern for RMS values per order
	rms_key = yaml_file["RMS"]		rms_key = yaml_file["RMS"]
			# Header keyword pattern for SNR values per order
	snr_key = yaml_file["SNR"]		snr_key = yaml_file["SNR"]
	#
	self.berv = -1		# ======================== BARYCENTRIC VELOCITY EXTRACTION ========================
	self.bervmax = -1		# Initialize BERV (Barycentric Earth Radial Velocity) values
			self.berv = -1 # BERV at observation time
			self.bervmax = -1 # Maximum BERV during exposure
			# Extract BERV from header if available
	if yaml_file["BERV"] != "None":		if yaml_file["BERV"] != "None":
	self.berv = self.header[yaml_file["BERV"]]		self.berv = self.header[yaml_file["BERV"]]
			# Extract maximum BERV if available
	if yaml_file["BERVMAX"] != "None":		if yaml_file["BERVMAX"] != "None":
	self.bervmax = self.header[yaml_file["BERVMAX"]]		self.bervmax = self.header[yaml_file["BERVMAX"]]
	#
			# ======================== ORDER-BY-ORDER PROCESSING INITIALIZATION ========================
			# List to store Order objects for each spectral order
	orders_info = list()		orders_info = list()
			# Initialize arrays to store SNR and RMS for each order
	self.snr = np.zeros([self.norder])		self.snr = np.zeros([self.norder])
	self.rms = np.zeros([self.norder])		self.rms = np.zeros([self.norder])
	#
			# Create pixel array for wavelength calibration
			# Range from start_pixel_detector to (npix + start_pixel_detector)
	x = np.arange(self.start_pixel_detector, self.npix + self.start_pixel_detector)		x = np.arange(self.start_pixel_detector, self.npix + self.start_pixel_detector)
			# List to store polynomial coefficients for each order
	self.coeff_poly_list = list()		self.coeff_poly_list = list()

			# ======================== MAIN LOOP: PROCESS EACH SPECTRAL ORDER ========================
	for i in range(self.norder):		for i in range(self.norder):
	#		if structure_with_multiple_chips:
			current_struct = struct_crires[i]
			chip_idx = current_struct[0]
			order = current_struct[1]
			chip = hdul[chip_idx]
			spec_name = f"{order}_01_SPEC"
			wl_name = f"{order}_01_WL"
			err_name = f"{order}_01_ERR"
			data_order_chip = chip.data[spec_name]
			data_order_chip[data_order_chip <= 0] = 1
			self.wl_cal[i, :] = chip.data[wl_name]
			self.data[i, :] = chip.data[spec_name]
			self.snr_total[i, :] = chip.data[spec_name] / chip.data[err_name]
			self.rms[i] = -999
			self.snr[i] = np.mean(self.snr_total[i, :])
			else:
			# Process order number i (from 0 to norder-1)
			# --- Extract Polynomial Wavelength Coefficients (if available) ---
	if yaml_file["POL"] != "None":		if yaml_file["POL"] != "None":
			# Calculate the indices for header keywords containing polynomial coefficients
			# The indexing scheme depends on the order number and coefficient degree
			# Formula: degree * mult_value + (order_number) * (poly_deg * mult_value + 1)
	coeff_indices_header = np.array(		coeff_indices_header = np.array(
	degree_poly_arr * mult_value +		degree_poly_arr * mult_value +
	(i + self.start_order) *		(i + self.start_order) *
	(self.poly_deg * mult_value + 1),		(self.poly_deg * mult_value + 1),
	dtype=str		dtype=str
	)		)
			# Extract polynomial coefficients from header using calculated indices
			# Creates a list of coefficients [c0, c1, c2, ...] for wavelength = c0 + c1x + c2x^2 + ...
	coeff_poly = [		coeff_poly = [
	self.header[		self.header[
	self.poly_deg_keys[d] + degree		self.poly_deg_keys[d] + degree
	] for d, degree in enumerate(coeff_indices_header)		] for d, degree in enumerate(coeff_indices_header)
	]		]
			# Create list of header keyword names for debugging
	stampa = [self.poly_deg_keys[d] + degree for d, degree in enumerate(coeff_indices_header)]		stampa = [self.poly_deg_keys[d] + degree for d, degree in enumerate(coeff_indices_header)]
			# Extract pixel offset (x_sub) if used
			# Some instruments use a reference pixel offset for each order
	if x_sub_key != "None":		if x_sub_key != "None":
	x_sub = self.header[x_sub_key + str(i + self.start_order)]		x_sub = self.header[x_sub_key + str(i + self.start_order)]
	else:		else:
	x_sub = 0		x_sub = 0
			# Create a copy of coefficients and append the pixel offset
	temp_coeff = copy.deepcopy(coeff_poly)		temp_coeff = copy.deepcopy(coeff_poly)
	temp_coeff.append(x_sub)		temp_coeff.append(x_sub)
			# Store coefficients for this order
	self.coeff_poly_list.append(temp_coeff)		self.coeff_poly_list.append(temp_coeff)
			# Print header keywords used (for debugging)
	print(stampa, x_sub_key)		print(stampa, x_sub_key)
	#		# --- Extract SNR and RMS per Order ---
			# Replace "xxx" placeholder in SNR keyword with actual order number
	temp_key = snr_key.replace("xxx", str(i + self.start_order))		temp_key = snr_key.replace("xxx", str(i + self.start_order))
			# Read SNR value from header for this order
	self.snr[i] = self.header[temp_key]		self.snr[i] = self.header[temp_key]
			# Extract RMS (Root Mean Square) if available in header
	if rms_key != "None":		if rms_key != "None":
			# Replace "xxx" placeholder with order number
	temp_key = rms_key.replace("xxx", str(i + self.start_order))		temp_key = rms_key.replace("xxx", str(i + self.start_order))
			# Read RMS value from header
	self.rms[i] = self.header[temp_key]		self.rms[i] = self.header[temp_key]
	else:		else:
			# If RMS not available, set sentinel value
	self.rms[i] = -999		self.rms[i] = -999
	#		# --- Create Order Object and Calculate Wavelength Solution ---
	if yaml_file["POL"] != "None":		if yaml_file["POL"] != "None":
			# Create an Order object for this spectral order
			# This handles wavelength calibration using polynomial coefficients
	orders_info.append(Order(		orders_info.append(Order(
	i, self.data[i, :], coeff_poly, x, self.snr[i], self.rms[i], x_sub,		i, # Order index
	yaml_file["COEFF_TO_NM"], yaml_file["AIR_TO_VAC"],		self.data[i, :], # Flux data for this order
			coeff_poly, # Polynomial coefficients
			x, # Pixel array
			self.snr[i], # SNR for this order
			self.rms[i], # RMS for this order
			x_sub, # Pixel offset
			yaml_file["COEFF_TO_NM"], # Wavelength unit conversion factor
			yaml_file["AIR_TO_VAC"], # Air-to-vacuum wavelength conversion flag
	))		))
			# Store the calculated wavelength solution for this order
			# The Order class computes wavelength from polynomial: wl = c0 + c1x + c2x^2 + ...
	self.wl_cal[i, :] = orders_info[i].wl_poly		self.wl_cal[i, :] = orders_info[i].wl_poly
			# Close the FITS file to free up system resources
			hdul.close()

DataReductionGIANOB/Instruments/criresp.yaml

0 → 100644

+39 −0

Original line number	Original line	Diff line number	Diff line
			start_pixel_detector : 0
			pixel_dv : 1.07
			wl_min : 1 # depends on the configuration
			wl_max : 5 # depends on the configuration
			lat : -24.622997508
			long : -70.40249839
			alt : 2635
			start_order : 1
			index_data : "None"
			index_wlen : "None"
			index_header : 0
			data_in_columns : "None"
			files_type : RAW
			npix : 2048
			norder : 24
			ra_dec_deg : 1
			corr_zero : 1
			NEED_SUM_EXPOSURE_TO_MJD : 1
			TELESCOP : "None"
			AIRMASS1 : "HIERARCH ESO TEL AIRM START"
			AIRMASS2 : "HIERARCH ESO TEL AIRM END"
			EXPTIME : "HIERARCH ESO DET SEQ1 DIT"
			HUMIDITY : "HIERARCH ESO TEL AMBI RHUM"
			RA : "RA"
			DEC : "DEC"
			MJD : "MJD-OBS"
			POL: "None"
			BERV : "None"
			BERVMAX : "None"
			X_SUB : "None"
			RMS : "None"
			SNR : "None"
			MULTIPLICATIVE_ORDER_KEY : "None"
			AIR_TO_VAC : 0
			COEFF_TO_NM : 10
			ALIGNED : 0
			OVERSCAN_LEFT: 0
			OVERSCAN_RIGHT: 0
			STRUCTURE_WITH_MULTIPLE_CHIPS : 1

DataReductionGIANOB/Instruments/espresso.yaml

+1 −0

Original line number	Original line	Diff line number	Diff line
	@@ -36,3 +36,4 @@ COEFF_TO_NM : 10
	ALIGNED : 1		ALIGNED : 1
	OVERSCAN_LEFT: 1360		OVERSCAN_LEFT: 1360
	OVERSCAN_RIGHT: 7902		OVERSCAN_RIGHT: 7902
			STRUCTURE_WITH_MULTIPLE_CHIPS : 0

DataReductionGIANOB/Instruments/gianob.yaml

+1 −0

Original line number	Original line	Diff line number	Diff line
	@@ -44,3 +44,4 @@ COEFF_TO_NM : 1
	ALIGNED : 0		ALIGNED : 0
	OVERSCAN_LEFT: 0		OVERSCAN_LEFT: 0
	OVERSCAN_RIGHT: 0		OVERSCAN_RIGHT: 0
			STRUCTURE_WITH_MULTIPLE_CHIPS : 0