DM

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\n@author: Pierre Guiot\\nAdaptation en notebook: Raphaël Wicker, Cateline Lantz\\n'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "@author: Pierre Guiot\n",
    "Adaptation en notebook: Raphaël Wicker, Cateline Lantz\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "import glob\n",
    "from astropy.io import fits # Pour ouvrir et manipuler les données\n",
    "from scipy.ndimage import rotate # Pour la rotation des images\n",
    "from scipy.optimize import curve_fit # Pour la régression linéaire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture des données\n",
    "\n",
    "Les prochaines lignes permettent d'accéder aux images de science, de Dark, de calib et potentiellement de Flats, pour les lire et récupérer les données sous forme d'array numpy.\n",
    "\n",
    "Les chemins doivent prendre la forme `chemin/vers/les/fichiers/*`, l'astérique servant à signifier que vous récupérez tous les fichiers à l'intérieur du dossier qui vous intéresse.\n",
    "\n",
    "Dans la cellule ci-dessous, vous pourrez trouver trois boucles `for` (une pour chaque type de donnée, quatre si vous utilisez le flat), ainsi que des tableaux `XXX_cube`:\n",
    "\n",
    "* Les boucles `for` parcourent tous les fichiers présents dans le chemin fourni, en les triant dans l'ordre alphanumérique, et en fournissant une variable \"compteur\" et une variable \"nom de fichier\" à la boucle `for`.\n",
    "* Les tableaux `XXX_cube` contiennent 3 dimensions selon lesquelles sont stockées les données : (n,l,c) où n est le numéro de l'image dans la série d'acquisitions, l la dimension des lignes de l'image et c la dimension des colonnes de l'image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "science_path = \" \" # Chemin vers les données de science\n",
    "dark_path = \" \" # Chemin vers les dark\n",
    "calib_path = \" \" # Chemin vers les calibs\n",
    "# Si utilisation des flats, rajouter le chemin\n",
    "\n",
    "science_cube = []\n",
    "dark_cube = []\n",
    "calib_cube = []\n",
    "\n",
    "for i, file in enumerate(sorted(glob.glob(science_path))): \n",
    "    hdu = fits.open(file)\n",
    "    image = hdu[0].data\n",
    "    science_cube.append(image)\n",
    "\n",
    "science_cube = np.asarray(science_cube, dtype=np.float64)\n",
    "\n",
    "for i, file in enumerate(sorted(glob.glob(dark_path))):\n",
    "    hdu = fits.open(file)\n",
    "    image = hdu[0].data\n",
    "    dark_cube.append(image)\n",
    "\n",
    "dark_cube = np.asarray(dark_cube, dtype=np.float64)\n",
    "\n",
    "for i, file in enumerate(sorted(glob.glob(calib_path))):\n",
    "    hdu = fits.open(file)\n",
    "    image = hdu[0].data\n",
    "    calib_cube.append(image)\n",
    "\n",
    "calib_cube = np.asarray(calib_cube, dtype=np.float64)\n",
    "# Si utilisation des flats, créer un flat_cube sur le modèle des blocs de code précédents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation\n",
    "\n",
    "Afin d'améliorer la visualisation des contrastes des images, vous pourrez spécifier des mots-clef `vmin` et `vmax` à l'intérieur des `plt.imshow()`.\n",
    "Sans spécifier ces seuils vous risquez de vous retrouver avec des images entièrement blanches, noires, ou saturées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "plt.figure(1)\n",
    "plt.clf()\n",
    "plt.imshow(science_cube[0], interpolation=None, origin='lower', cmap='Greys_r') #vmin = , vmax = \n",
    "plt.title('First science image')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.figure(2)\n",
    "plt.clf()\n",
    "plt.imshow(dark_cube[0], interpolation=None, origin='lower', cmap='Greys_r')\n",
    "plt.title('First Dark image')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.figure(3)\n",
    "plt.clf()\n",
    "plt.imshow(np.mean(dark_cube, axis=0), interpolation=None, origin='lower', cmap='Greys_r')\n",
    "plt.title('Stacked master dark image')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.figure(4)\n",
    "plt.clf()\n",
    "plt.imshow(calib_cube[0], interpolation=None, origin='lower', cmap='Greys_r')\n",
    "plt.title('First Calib image')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.figure(5)\n",
    "plt.clf()\n",
    "plt.imshow(np.mean(calib_cube, axis=0), interpolation=None, origin='lower', cmap='Greys_r')\n",
    "plt.title('Stacked Calib image')\n",
    "plt.colorbar()\n",
    "\n",
    "# Si utilisation des flats, les visualiser sur le modèle des blocs de code précédents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Réduction des données\n",
    "\n",
    "Dans cette section, le but est de réduire les données brutes de science, pour parvenir à des données exploitables scientifiquement. Une image réduite $I_r$ est la soustraction de l'image brute $I_b$ et de la moyenne des darks, ou Master Dark, MD : $I_t = I_b - MD$. En cas d'utilisation des Flats, il faut également diviser le tout par la moyenne des Flats, le Master Flat MF : $I_r = \\frac{I_b - MD}{MF}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(science_cube.shape[0]):\n",
    "    science_cube[i,:,:] = science_cube[i,:,:] - np.mean(dark_cube, axis=0) #Ajouter le master flat si nécessaire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "plt.figure(6)\n",
    "plt.clf()\n",
    "plt.imshow(science_cube[0,:,:], interpolation=None, origin='lower', cmap='Greys_r') #vmin = , vmax = \n",
    "# Essayer de modifier le premier indice du cube pour tester.\n",
    "plt.title('First science image with master dark subtracted')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rotation des spectres\n",
    "\n",
    "Le plus souvent, les spectres que vous obtiendrez en sortie de la réduction seront inclinés. Etant donné que votre but sera d'analyser des spectres, cela sera plus simple s'il suffit de récupérer une ligne dans un tableau numpy pour faire ce travail. Il faudra donc faire tourner l'image selon un angle que vous aurez déterminé pour que le spectre soit horizontal.\n",
    "\n",
    "Pour pouvoir étalonner facilement les spectres, il faut faire cette démarche avec l'image de calibration également.\n",
    "\n",
    "L'angle optimal de rotation peut être déterminé à l'aide de DS9. Attention l'angle optimal pour la rotation par python est l'opposé de l'angle donné par DS9 $\\theta_{opt, python} = - \\theta_{opt, DS9}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Données de science ###\n",
    "\n",
    "rotated_cube = [] \n",
    "rotation_angle = # Angle en degrés à déterminer à l'aide de DS9 pour mettre le spectre horizontal.\n",
    "for i in range(science_cube.shape[0]):\n",
    "    rotated_cube.append(rotate(science_cube[i,:,:], rotation_angle))\n",
    "rotated_cube = np.asarray(rotated_cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectre de référence pour l'étalonnage spectral\n",
    "Extraction du spectre de référence de la lampe à l'Argon-Néon : trouver l'indice de la ligne centrale du spectre selon laquelle on va effectuer une coupe pour extraire le spectre 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Images de calibration ###\n",
    "image = calib_cube[0,:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "plt.figure(8)\n",
    "plt.clf()\n",
    "plt.imshow(rotate(image,rotation_angle), interpolation=None, origin='lower', cmap='Greys_r') #vmin = ,  vmax =\n",
    "plt.title('First calib image rotated')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_cut = # rentrer l'indice de la ligne centrale du spectre\n",
    "ar_ne_spectrum = rotate(image,rotation_angle)[pixel_cut,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "### Données de science ###\n",
    "\n",
    "plt.figure(7)\n",
    "plt.clf()\n",
    "plt.imshow(rotated_cube[0,:,:], interpolation=None, origin='lower', cmap='Greys_r') #vmin = ,  vmax =\n",
    "plt.title('First science image rotated')\n",
    "plt.colorbar()\n",
    "\n",
    "### Images de calibration ###\n",
    "\n",
    "plt.figure(9)\n",
    "plt.clf()\n",
    "plt.plot(ar_ne_spectrum)\n",
    "plt.xlabel('Array index')\n",
    "plt.ylabel('Arbitrary intensity')\n",
    "plt.title('Ar Ne calibration lamp spectrum')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Etalonnage spectral à l'aide la lampe ArNe\n",
    "\n",
    "Dans cette section, le but sera d'étalonner votre spectre, en reliant les positions des pics en pixels aux longueurs d'onde correspondantes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Définition de la zone du spectre\n",
    "Ici le but est de trouver dans quelle zone de l'image se trouve le signal significatif du spectre : en effectuant une coupe verticale (dimension spatiale de l'image, perpendiculaire au spectre) dans l'image, on repère la zone d'intérêt et les zones qui ne sont que du bruit. Un critère quantitatif permettant de déterminer ces limites pourrait être l'écart-type.\n",
    "\n",
    "Notez bien que selon le sujet que vous traitez, vous pouvez avoir besoin de l'information contenue le long de la direction spatiale. Dans ce cas, il ne faut pas utiliser de spectres moyennés verticalement sur plusieurs pixels, mais à des spectres extraits sur une seule ligne."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Coupe verticale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "plt.figure(10)\n",
    "plt.clf()\n",
    "plt.plot(rotated_cube[0,:,2210]) # Indice de la colonne où la coupe est effectuée. À modifier\n",
    "plt.ylim(0,1000) # À modifier\n",
    "plt.ylabel('Arbitrary intensity')\n",
    "plt.xlabel('Array index in px')\n",
    "plt.grid()\n",
    "plt.title('Spatial cut in Atik One images')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extraction du spectre 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "average_pixels= # Pixel central autour duquel les 20 pixels voisins sont moyennés pour \"lisser\" le spectre extrait.\n",
    "\n",
    "#======================================================================================================\n",
    "# À ne faire que si la dimension spatiale de vos données n'a pas d'importance ! Sinon, ne pas moyenner.\n",
    "extract_avg_spectrum = np.mean(rotated_cube[0,average_pixels:average_pixels+20,:], axis=0)\n",
    "#======================================================================================================\n",
    "\n",
    "central_pix= # Pixel central où le spectre 1D est extrait (sans moyenne). À modifier\n",
    "plt.figure(11)\n",
    "plt.clf()\n",
    "plt.plot(extract_avg_spectrum, 'k-', lw=1, label='spatially averaged')\n",
    "plt.plot(rotated_cube[0,central_pix,:], 'r-', label='central cut')\n",
    "plt.ylabel('Arbitrary intensity')\n",
    "plt.xlabel('Arrays index, not wavelengths yet')\n",
    "plt.title(' ') # À modifier\n",
    "plt.grid()\n",
    "plt.ylim()\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correspondance pixels/longueurs d'onde\n",
    "Utiliser le spectre de référence de l'ArNe présenté dans le polycopié pour repérer les raies et faire la correspondance entre leur longueur d'onde lue sur ce spectre et les indices en pixels lus sur vos graphiques précédents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "plt.figure(12)\n",
    "plt.clf()\n",
    "plt.plot(ar_ne_spectrum)\n",
    "plt.xlabel('Array index')\n",
    "plt.ylabel('Arbitrary intensity')\n",
    "plt.title('Ar Ne calibration lamp spectrum')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixels_peaks = np.array([ ]) # A remplir\n",
    "peaks_lambda = np.array([ ]) # A remplir\n",
    "def affine(x, a=1, b=1):\n",
    "    return a*x + b\n",
    "\n",
    "calib_wvl = curve_fit(affine, pixels_peaks, peaks_lambda, p0=[2, 50])\n",
    "\n",
    "calib_wavelengths = affine(range(len(extract_avg_spectrum)), a=calib_wvl[0][0], b=calib_wvl[0][1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation des spectres étalonnés"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "\n",
    "plt.figure(13)\n",
    "plt.clf()\n",
    "plt.plot(calib_wavelengths, extract_avg_spectrum, 'k-', lw=1, label='spatially averaged')\n",
    "plt.plot(calib_wavelengths, rotated_cube[0,central_pix,:], 'r-', label='central cut')\n",
    "plt.ylabel('Arbitrary intensity')\n",
    "plt.xlabel('Wavelengths in nm')\n",
    "plt.title('Reduced and calibrated spectra')\n",
    "plt.grid()\n",
    "plt.legend(loc='best')\n",
    "plt.ylim(0, 2000) # A modifier\n",
    "\n",
    "\n",
    "print(\"Traitement terminé, place à l'analyse !\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}