{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make Your Own Reduction Scrypt (photometry)\n",
    "\n",
    "*Work in Progress*\n",
    "\n",
    "ASTROPOP is not a reduction script by itself, but a library containing (almost) everything you need to create your reduction script by yourself.\n",
    "\n",
    "In this guide, we will follow a standard reduction procedure using the ASTROPOP modules and you can follow it to perform your own reduction.\n",
    "\n",
    "If you want to use this notebook directly in your Jupyter, you can download it at [this link](https://raw.githubusercontent.com/sparc4-dev/astropop/main/docs/ipynb/diy_reduction_script.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Organize Your Data\n",
    "\n",
    "The first step is to organize your data. *Do not forget to backup all the raw data before starting your reduction!*\n",
    "\n",
    "Here we will perform a reduction using just bias and flatfield corrections, in a single band. No dark frame subtraction will be used, so we do not have to set it.\n",
    "\n",
    "First, we will have to create the lists containing the image names we need. You can do it manually, using [glob](https://docs.python.org/3/library/glob.html) package, or the built-in file organizer of ASTROPOP. This last option allow filtering using header keys, but may be slow and memory consuming for very large sets (with thousands of images). For didatic reasons, we will to it with ASTROPOP file organizer.\n",
    "\n",
    "It is convenient to set now some directory names for raw, reduced and temporary files, as in the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from astropy import units as u\n",
    "\n",
    "# astropop used modules\n",
    "from astropop.image import imcombine, processing, imarith\n",
    "from astropop.framedata import read_framedata\n",
    "from astropop.image.register import register_framedata_list\n",
    "from astropop.photometry import background, starfind\n",
    "from astropop.logger import logger\n",
    "logger.setLevel('INFO')\n",
    "\n",
    "# Directories where we will storefiles\n",
    "base_dir = os.path.expanduser('~/astropop-tutorial/photometry-example')\n",
    "raw_dir = os.path.join(base_dir, 'raw')  # Directory containing the rawdata\n",
    "reduced_dir = os.path.join(base_dir, 'reduced')  # Directory to store processed data\n",
    "tmp_dir = os.path.join(base_dir, 'tmp')  # Directory to store temporary files\n",
    "\n",
    "# create the dir if not exists\n",
    "os.makedirs(raw_dir, exist_ok=True)  \n",
    "os.makedirs(reduced_dir, exist_ok=True)\n",
    "os.makedirs(tmp_dir, exist_ok=True)\n",
    "\n",
    "from matplotlib import pyplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will a set of images of HD5980 star taken in Observatório do Pico dos Dias (OPD/LNA) and available in [github/sparc4-dev/astropop-data](https://github.com/sparc4-dev/astropop-data) repository.\n",
    "\n",
    "To download these data in the raw folder, follow the script below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downloading the data. If you already downloaded it, you can skip this cell\n",
    "import urllib.request\n",
    "import os\n",
    "\n",
    "with urllib.request.urlopen('https://raw.githubusercontent.com/sparc4-dev/astropop-data/main/raw_images/opd_ixon_bc_hd5980/filelist.txt') as f:\n",
    "    for line in f.readlines():\n",
    "        url = line.decode('UTF-8')\n",
    "        filename = os.path.join(raw_dir, url.split('/')[-1].split('?')[0])\n",
    "        urllib.request.urlretrieve(url, filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `astropop.file_collection.FitsFileGroup` we will create a dynamic library of fits files. This library allows iterating over its files, filtering by header keywords, accessing file names or direct access to a given index. In this case, we will access only the file names via `FitsFileGroup.files` property.\n",
    "\n",
    "For the `filtered` method, the arguments is a dictionary containing all keywords you want to filter in the header. Only equality is allowed. In this example, we will filter our files by the `obstype` key, which can have values `BIAS`, `FLAT` and `SCIENCE`. For science, if you have multiple stars, you can use a second key, `object`, to filter the desired images. This is the standard of our observations, but your observatory may contian a different standard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropop.file_collection import FitsFileGroup\n",
    "\n",
    "main_fg = FitsFileGroup(location=raw_dir, fits_ext=['.fits'], ext=0)\n",
    "\n",
    "print(f'Total files: {len(main_fg)}')\n",
    "\n",
    "# Filter files by header keywords\n",
    "bias_fg = main_fg.filtered({'obstype': 'BIAS'})  # OBSTYPE='BIAS'\n",
    "print(f'Bias files: {len(bias_fg)}')\n",
    "print(bias_fg.files)\n",
    "flat_fg = main_fg.filtered({'obstype': 'FLAT'})  # OBSTYPE='FLAT'\n",
    "print(f'Flat files: {len(flat_fg)}')\n",
    "print(flat_fg.files)\n",
    "star_fg = main_fg.filtered({'obstype': 'SCIENCE', 'object': 'HD5980'})  # OBSTYPE='SCIENCE' and OBJECT='HD5980'\n",
    "print(f'Star files: {len(star_fg)}')\n",
    "print(star_fg.files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Master Calibration Frames\n",
    "\n",
    "Before process the science images, we need to create the `master_bias` and `master_flat` calibration frames.\n",
    "\n",
    "### Master Bias\n",
    "\n",
    "Bias files will be corrected by gain, cosmic removal and after they will be median combined to generate the `master_bias` file.\n",
    "\n",
    "The first step is to read the fits files to `FrameData` instances. `use_memmap_backend` feature will be used to save RAM. `FitsFileGroup`s already have a feature to read and iterate over the `FrameData` created instances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if you already have the master bias and the master flat, you can skip the next two cells\n",
    "bias_frames = list(bias_fg.framedata(unit='adu', use_memmap_backend=True))\n",
    "\n",
    "# lest extract gain from the first image\n",
    "gain = float(bias_frames[0].header['GAIN'])*u.electron/u.adu  # using quantities is better for safety\n",
    "print('gain:', gain)\n",
    "\n",
    "# Perform calibrations\n",
    "for i, frame in enumerate(bias_frames):\n",
    "    print(f'processing frame {i+1} from {len(bias_frames)}')\n",
    "    processing.cosmics_lacosmic(frame, inplace=True)\n",
    "    processing.gain_correct(frame, gain, inplace=True)\n",
    "\n",
    "# combine\n",
    "master_bias = imcombine(bias_frames, method='median')\n",
    "\n",
    "# save to disk\n",
    "mbias_name = os.path.join(reduced_dir, 'master_bias.fits')\n",
    "master_bias.write(mbias_name, overwrite=True, no_fits_standard_units=True)\n",
    "master_bias.meta['astropop master_type'] = 'bias'\n",
    "print('master_bias statistics', master_bias.statistics())\n",
    "\n",
    "del bias_frames  # remove tmp frames from memory and disk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Master Flat\n",
    "\n",
    "In addition to cosmic/gain correction, flat frames must be subtracted by the master bias and normalized after the combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_frames = list(flat_fg.framedata(unit='adu', use_memmap_backend=True))\n",
    "\n",
    "for i, frame in enumerate(flat_frames):\n",
    "    print(f'processing frame {i+1} from {len(flat_frames)}')\n",
    "    processing.cosmics_lacosmic(frame, inplace=True)\n",
    "    processing.gain_correct(frame, gain, inplace=True)\n",
    "    processing.subtract_bias(frame, master_bias, inplace=True)\n",
    "\n",
    "# combine and normalize\n",
    "master_flat = imcombine(flat_frames, method='median')\n",
    "mean_value = master_flat.mean()\n",
    "print('flat mean value:', mean_value)\n",
    "master_flat = imarith(master_flat, mean_value, '/')\n",
    "master_bias.meta['astropop master_type'] = 'flat'\n",
    "\n",
    "# save to disk\n",
    "mflat_name = os.path.join(reduced_dir, 'master_flat.fits')\n",
    "master_flat.write(mflat_name, overwrite=True)\n",
    "print('master_flat statistics', master_flat.statistics())\n",
    "\n",
    "del flat_frames  # remove tmp frames from memory and disk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic CCD Image Processing\n",
    "\n",
    "Now, with our master images created, its time to performe the processing of our science images.\n",
    "\n",
    "The steps are:\n",
    "- cosmic ray removal\n",
    "- gain correct\n",
    "- master_bias subtraction\n",
    "- division by normalized flat\n",
    "\n",
    "We already have the master images loaded into the memory. But, considering you have the master frames already created, lets load them from the disk using `read_framedata` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of framedata of our science frames\n",
    "star_frames = list(star_fg.framedata(unit='adu', use_memmap_backend=True))\n",
    "\n",
    "# needed parameters for calibration\n",
    "gain = float(star_frames[0].header['GAIN'])*u.electron/u.adu  # using quantities is better for safety\n",
    "print('gain:', gain)\n",
    "mbias_name = os.path.join(reduced_dir, 'master_bias.fits')\n",
    "mflat_name = os.path.join(reduced_dir, 'master_flat.fits')\n",
    "master_bias = read_framedata(mbias_name)\n",
    "print('bias: ', master_bias.statistics())\n",
    "master_flat = read_framedata(mflat_name)\n",
    "print('flat:', master_flat.statistics())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# process the science frames and save individual files calibrated files to disk\n",
    "for i, frame in enumerate(star_frames):\n",
    "    print(f'processing frame {i+1} from {len(star_frames)}')\n",
    "    processing.cosmics_lacosmic(frame, inplace=True)\n",
    "    processing.gain_correct(frame, gain, inplace=True)\n",
    "    processing.subtract_bias(frame, master_bias, inplace=True)\n",
    "    processing.flat_correct(frame, master_flat, inplace=True)\n",
    "    \n",
    "for frame, name in zip(star_frames, star_fg.files):\n",
    "    basename = os.path.basename(name)\n",
    "    frame.write(os.path.join(reduced_dir, basename), overwrite=True, no_fits_standard_units=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Align and Stack\n",
    "\n",
    "To get a better static photometry, it is good to have several images ofthe field, align them and stack to reduce the noise. This process of align images is called [image registration](https://en.wikipedia.org/wiki/Image_registration).\n",
    "\n",
    "Astropop has two algorithms of image registration:\n",
    "- `cross-correlation`: that computes the phase difference in fourier space using the cross-correlation of the images. This method allow register any images. But it may have less precision and is vulnerable for structures of bad rows/columns in the image.\n",
    "- `asterism-matching`: it uses groups of 3 detected stars in each image to estimate the transform between the by similarity. You can only apply the method when you have good stars in the field. It supports translation and rotation and it's errors depend on the errors of the centroid detection.\n",
    "\n",
    "As we are working with star filed images, lets use the asterism matching. Astropop directly preform registration in `FrameData` container lists!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "registered_frames = register_framedata_list(star_frames, algorithm='asterism-matching', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To combine the images we will use the `imcombine` function of astropop. In our case we will perform the sum of all images, excluding masked pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined = imcombine(registered_frames, method='sum')\n",
    "pyplot.figure()\n",
    "pyplot.imshow(combined.data, vmin=combined.median().value, vmax=np.percentile(combined.data, 99.5), origin='lower')\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Source Detection and Photometry\n",
    "\n",
    "We will perform a very standard aperture photometry reduction using commot astropop default routines. Another routines are available and read the documentation is recommended for customizations.\n",
    "\n",
    "To detect background, we will use the `background` function of astropop. As our image has a big nebula that covers a big part of the image, we will use a space-varying background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(combined.data)\n",
    "bkg, rms = background(data, global_bkg=False)\n",
    "vmin=np.median(combined.data)\n",
    "vmax=np.percentile(combined.data, 99.5)\n",
    "\n",
    "fig, ax = pyplot.subplots(1, 3, sharex=True, sharey=True, figsize=(8, 4))\n",
    "ax[0].imshow(combined.data, vmin=vmin, vmax=vmax, origin='lower')\n",
    "ax[0].set_title('Data')\n",
    "ax[1].imshow(bkg, vmin=vmin, vmax=vmax, origin='lower')\n",
    "ax[1].set_title('Background')\n",
    "ax[2].imshow(data-bkg, vmin=0, vmax=vmax, origin='lower')\n",
    "ax[2].set_title('Data-Background')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets identify the sources in the image. For this we will use the `starfind` function. It will perform the identification of punctual sources only. It has an advantage against `daofind` beacuse it operates blindly, with no necessity of meassure the FWHM previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sources = starfind(data, threshold=10, background=bkg, noise=rms)\n",
    "print('Total found stars:', len(sources))\n",
    "pyplot.figure()\n",
    "pyplot.imshow(data-bkg, vmin=0, vmax=np.percentile(combined.data, 95), origin='lower')\n",
    "pyplot.plot(sources['x'], sources['y'], 'w+', ms=5)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Photometry Calibration with Online Catalogs"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "162041a95dc008afa7a541b7cce9944c08c53eb30e22b146ef152cba0893111e"
  },
  "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}