E: Ejercicio resuelto: Generacion de una funcion para reducir datos¶
Ejercicio E.1:
Basándonos en el código de los anteriores notebooks crear una función a la que pasándole el directorio donde están nuestras imágenes, bias, darks y flats reduzca la imagen automáticamente.
[1]:
import matplotlib.pyplot as plt
import numpy as np
import glob
from astropy.io import fits
import os
import astroscrappy
[2]:
import astroalign
from astropy.io import fits as ft
def imprimir_info(p, ii):
# Esta función imprime por pantalla la info de la transformación que se aplicará.
# Aunque es innecesaria, sirve para que el usuario sepa que la máquina está haciendo algo.
print("\nAlineando imagen {:}".format(ii))
print("Rotación: {:.2f} grados".format(p.rotation * 180.0 / np.pi))
print("Factor de escala: {:.2f}".format(p.scale))
print("Traslación: (x, y) = ({:.2f}, {:.2f})".format(*p.translation))
def registra_lista(lista, infoAlineado = False):
'''
Función para la alineacion de imagenes a traves de astroalign
lista: tenemos que pasar una lista de la ruta donde se encuentran nuestras imagenes. List of string.
Return: Lista de numpy array con los datos de las imagenes alineadas. List of numpy array.
'''
cantidad=len(lista)
# La primera imagen de la lista será la toma de referencia.
print("\nComenzando la alineación.")
print("\nLa toma de referencia es {:}".format(lista[0]))
blanco=ft.open(lista[0])
img_blanco=blanco[0].data
hdr_blanco=blanco[0].header
blanco.close()
del(lista[0]) # Quito la imagen de referencia del listado
imagenes_alineadas = [img_blanco, ]
for ii in lista:
ff=ft.open(ii)
img_torcida=ff[0].data
hdr_torcida=ff[0].header
ff.close()
p, (pos_img, pos_img_rot) = astroalign.find_transform(img_torcida, img_blanco)
img_aligned = astroalign.register(img_torcida, img_blanco)
hdr_torcida.add_comment("Registrado con Astroalign y PyReduc")
imagenes_alineadas.append(img_aligned[0])
if infoAlineado:
imprimir_info(p, ii)
print("\nAlineado realizado con éxito")
return imagenes_alineadas
Creamos nuestra funcion reducir en donde las entradas de la función serán los directorios donde tenemos guardados los archivos FITS.
[3]:
def reducir(dirBias, dirDarks, dirFlats, dirLights, dirSalida, imaInter = False, infoCosmicray = False):
'''
Función para reducir datos en formato fit. En los directorios solo deben
estar los ficheros particulares de cada caso.
dirBias: directorio donde se encuentran los bias. String.
dirDarks: directorio donde se encuentran los darks. String.
dirFlats: directorio donde se encuentran los bias. String.
dirLights: directorio donde se encuentran los bias. String.
dirSalida: directorio donde irán las imagenes calibradas en formato fit. String.
imaInter: Si queremos guardar las imagenes intermedias pondremos True. Booleano.
infoCosmicray: Si queremos informacion de los rayos cosmicos que ha detectado el algoritmo astroscrappy. Booleano.
Return: Imagen calibrada en forma de matriz. numpy array.
'''
try:
os.mkdir(dirSalida)
except:
print(f"El directorio {dirSalida} ya existe")
bias_list = sorted(glob.glob(f"{dirBias}*"))
darks_list = sorted(glob.glob(f"{dirDarks}*"))
flats_list = sorted(glob.glob(f"{dirFlats}*"))
lights_list = sorted(glob.glob(f"{dirLights}*"))
# Bias
print("Procesando BIAS")
bias_data = []
bias_header = []
for bias in bias_list:
hdul_bias = fits.open(bias)
bias_data.append(hdul_bias[0].data)
bias_header.append(hdul_bias[0].header)
hdul_bias.close()
for r in range(len(bias_list)):
if np.nanstd(bias_data[r]) > np.median(np.nanstd(bias_data))*1.2:
print("Eliminando BIAS: " + bias_list[r])
bias_list.remove(bias_list[r])
superBias = np.nanmedian(bias_data, axis=0)
bias_header[0]['history'] = f"Esta imagen es un stacking de {len(bias_list)} imagenes"
# Darks
print("Procesando DARKs")
darks_data = []
darks_header = []
darks_current = []
for dark in darks_list:
hdul_darks = fits.open(dark)
darks_data.append(hdul_darks[0].data)
darks_header.append(hdul_darks[0].header)
hdul_darks.close()
for dark in range(len(darks_list)):
darks_current.append(darks_data[dark] - superBias)
for i in range(len(darks_list)):
darks_current[i][darks_current[i] < 0] = 0
superDark = np.nanmedian(darks_data, axis=0)
darkCurrent = np.median(darks_current, axis=0)
#Flats
print(f"Procesando FLATs")
flats_data = []
flats_header = []
flat_data_b = []
for flat in flats_list:
hdul_flats = fits.open(flat)
flats_data.append(hdul_flats[0].data)
flats_header.append(hdul_flats[0].header)
hdul_flats.close()
for r in range(len(flats_list)):
if np.nanstd(flats_data[r]) > np.median(np.nanstd(flats_data))*1.1:
print(bias_list[r])
bias_list.remove(bias_list[r])
y_len, x_len = np.shape(flats_data[0])
x_center, y_center = x_len//2, y_len//2
for flat in range(len(flats_list)):
flat_superBias = flats_data[flat] - superBias
flat_norm = flat_superBias / np.median(flat_superBias[y_center-250:y_center+250,x_center-250:x_center+250])
flat_data_b.append(flat_norm)
superFlat = np.nanmedian(flat_data_b, axis=0)
#lights
print(f"Procesando LIGHTs de {flats_header[0]['FILTER']}")
lights_header = []
lights_data = []
FWHM = []
for f in lights_list:
hdul = fits.open(f)
lights_header.append(hdul[0].header)
lights_data.append(hdul[0].data)
FWHM.append(hdul[0].header['FWHM'])
hdul.close()
# print(f"Lights de {lights_header[0]['FILTER']}")
for r in range(len(lights_list)):
if lights_header[r]['FWHM'] > np.median(FWHM)*1.5:
lights_data.remove(lights_data[r])
lights_header.reomve(lights_header[r])
lights_list.remove(lights_list[r])
print(f"Se ha eliminado la imagen {lights_list[r]}")
lights_correction = []
for light in range(len(lights_list)):
light_b =lights_data[light].astype(np.float32) - superBias.astype(np.float32)
light_bd = light_b - darkCurrent.astype(np.float32)
light_bdf = light_bd / superFlat.astype(np.float32)
lights_correction.append(light_bdf)
for hot in range(len(lights_list)):
hot_pixeles = np.where(lights_correction[hot] < 0) # hotpixel será un array con la posicion del los hot pixeles
for i in range(len(hot_pixeles[0])): #
lights_correction[hot][hot_pixeles[0][i]][hot_pixeles[1][i]] = np.median([lights_correction[hot][hot_pixeles[0][i]-1, hot_pixeles[1][0]], lights_correction[hot][hot_pixeles[1][i]+1, hot_pixeles[1][0]]])
lights_correction_cosmicray = []
mask_cosmicray = []
for image in range(len(lights_list)):
mask, cosmicray = astroscrappy.detect_cosmics(lights_correction[image], sigclip=6, sigfrac=15, objlim=5.0, gain=1.0, readnoise=6.5, satlevel=65536.0, niter=4, verbose=infoCosmicray)
lights_correction_cosmicray.append(cosmicray)
mask_cosmicray.append(mask)
imagenes_alineadas = registra_lista(lights_list)
imagenFinal = np.nanmedian(imagenes_alineadas, axis=0)
# Guardamos
if imaInter:
hdu_Bias = fits.PrimaryHDU(data=superBias.astype(np.float32), header=bias_header[0])
hdu_Bias.writeto(dirSalida + "superBias.fit", overwrite=True)
super_dark = fits.PrimaryHDU(data=superDark.astype(np.float32), header=darks_header[0])
super_dark.writeto(dirSalida + "superDark.fit", overwrite=True)
hdu_Dark_current = fits.PrimaryHDU(data=darkCurrent.astype(np.float32))
hdu_Dark_current.writeto(dirSalida + "darkCurrent.fit", overwrite=True)
hdu_FlatH = fits.PrimaryHDU(data=superFlat.astype(np.float32), header=flats_header[0])
hdu_FlatH.writeto(dirSalida + "superFlat.fit",overwrite=True)
hdu_imagenFinal = fits.PrimaryHDU(data=imagenFinal.astype(np.float32), header=lights_header[0])
hdu_imagenFinal.writeto(dirSalida + f"imagenFinal_{lights_header[0]['FILTER']}.fit",overwrite=True)
return imagenFinal
Generamos los directorios y llamamos a la funcion reducir.
[4]:
import os
baseDir = "imagenes/calibracionImagenes/"
dirBias = baseDir + 'bias/'
dirDarks = baseDir + 'darks/'
dirFlatsHa = baseDir + 'flats/ha/'
dirFlatsO3 = baseDir + 'flats/O3/'
dirFlatsS2 = baseDir + 'flats/S2/'
dirLightsHa = baseDir + 'ha/'
dirLightsO3 = baseDir + 'O3/'
dirLightsS2 = baseDir + 'S2/'
dirSalida = 'salidas/salidaCalibracion/sho/'
try:
os.mkdir(dirSalida)
except:
print("Probablemente el directorio ya existe")
ha = reducir(dirBias, dirDarks, dirFlatsHa, dirLightsHa, dirSalida, imaInter = False)
O3 = reducir(dirBias, dirDarks, dirFlatsO3, dirLightsO3, dirSalida, imaInter = False)
S2 = reducir(dirBias, dirDarks, dirFlatsS2, dirLightsS2, dirSalida, imaInter = False)
Probablemente el directorio ya existe
El directorio salidas/salidaCalibracion/sho/ ya existe
Procesando BIAS
Eliminando BIAS: imagenes/calibracionImagenes/bias/Telescope3_BIAS_1x1_0s_-30degC_20200120_0547_000001089.fit
Procesando DARKs
Procesando FLATs
Procesando LIGHTs de H-alpha
Comenzando la alineación.
La toma de referencia es imagenes/calibracionImagenes/ha/NGC2070-20200208@030339-300S-H-alpha.fit
Alineado realizado con éxito
El directorio salidas/salidaCalibracion/sho/ ya existe
Procesando BIAS
Eliminando BIAS: imagenes/calibracionImagenes/bias/Telescope3_BIAS_1x1_0s_-30degC_20200120_0547_000001089.fit
Procesando DARKs
Procesando FLATs
Procesando LIGHTs de Oiii
Comenzando la alineación.
La toma de referencia es imagenes/calibracionImagenes/O3/NGC2070-20200208@035232-300S-Oiii.fit
Alineado realizado con éxito
El directorio salidas/salidaCalibracion/sho/ ya existe
Procesando BIAS
Eliminando BIAS: imagenes/calibracionImagenes/bias/Telescope3_BIAS_1x1_0s_-30degC_20200120_0547_000001089.fit
Procesando DARKs
Procesando FLATs
Procesando LIGHTs de Sii
Comenzando la alineación.
La toma de referencia es imagenes/calibracionImagenes/S2/NGC2070-20200209@020725-600S-Sii.fit
Alineado realizado con éxito