anu_astr4004_2024_week7a_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3c5366df",
   "metadata": {},
   "source": [
    "# ANU ASTR4004 2024 - Week 7 (17+19 September 2024)\n",
    "\n",
    "Author: Dr Sven Buder (sven.buder@anu.edu.au)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a203675",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Time-yourself\" data-toc-modified-id=\"Time-yourself-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Time yourself</a></span></li><li><span><a href=\"#Preprocessing-of-Data-Files\" data-toc-modified-id=\"Preprocessing-of-Data-Files-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Preprocessing of Data Files</a></span><ul class=\"toc-item\"><li><span><a href=\"#Creating-data-directories-and-downloading-files\" data-toc-modified-id=\"Creating-data-directories-and-downloading-files-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Creating data directories and downloading files</a></span></li><li><span><a href=\"#Reading-CSV-(Comma-Separated-Values)-Files\" data-toc-modified-id=\"Reading-CSV-(Comma-Separated-Values)-Files-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Reading CSV (Comma Separated Values) Files</a></span></li><li><span><a href=\"#FITS-Images-(the-I-in-Flexible-Image-Transport-System)\" data-toc-modified-id=\"FITS-Images-(the-I-in-Flexible-Image-Transport-System)-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>FITS Images (the I in Flexible Image Transport System)</a></span><ul class=\"toc-item\"><li><span><a href=\"#Metadata:-store-useful-information-about-your-tables\" data-toc-modified-id=\"Metadata:-store-useful-information-about-your-tables-2.3.1\"><span class=\"toc-item-num\">2.3.1&nbsp;&nbsp;</span>Metadata: store useful information about your tables</a></span></li><li><span><a href=\"#The-world-coordinate-system-at-your-fingertips\" data-toc-modified-id=\"The-world-coordinate-system-at-your-fingertips-2.3.2\"><span class=\"toc-item-num\">2.3.2&nbsp;&nbsp;</span>The world coordinate system at your fingertips</a></span></li></ul></li></ul></li><li><span><a href=\"#Working-with-Data-Catalogues\" data-toc-modified-id=\"Working-with-Data-Catalogues-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Working with Data Catalogues</a></span><ul class=\"toc-item\"><li><span><a href=\"#Downloading-Data-via-ADQL/SQL\" data-toc-modified-id=\"Downloading-Data-via-ADQL/SQL-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Downloading Data via ADQL/SQL</a></span></li><li><span><a href=\"#Joining-2-catalogues\" data-toc-modified-id=\"Joining-2-catalogues-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Joining 2 catalogues</a></span></li><li><span><a href=\"#Cleaning-Catalogues\" data-toc-modified-id=\"Cleaning-Catalogues-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Cleaning Catalogues</a></span></li><li><span><a href=\"#Scaling-entries-(e.g.-for-fitting)\" data-toc-modified-id=\"Scaling-entries-(e.g.-for-fitting)-3.4\"><span class=\"toc-item-num\">3.4&nbsp;&nbsp;</span>Scaling entries (e.g. for fitting)</a></span></li><li><span><a href=\"#Visualising-targets-with-ipyaladin\" data-toc-modified-id=\"Visualising-targets-with-ipyaladin-3.5\"><span class=\"toc-item-num\">3.5&nbsp;&nbsp;</span>Visualising targets with ipyaladin</a></span></li></ul></li><li><span><a href=\"#Interpolation\" data-toc-modified-id=\"Interpolation-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Interpolation</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b966a0f0",
   "metadata": {},
   "source": [
    "If you want to get an automatic Table of Contents as above and codefolding, you can install notebook extensions for jupyter:\n",
    "\n",
    "```bash\n",
    "pip install jupyter_contrib_nbextensions\n",
    "jupyter contrib nbextension install --user\n",
    "```\n",
    "\n",
    "Then navigate to the \"Nbextensions\" bar when you start Jupyter and look for \"Table of Contents (2)\" as well as \"Codefolding\"  \n",
    "Note: You might have to tweak a few setting to avoid numbering Notebook headers like the one above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a44ccb03",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    %matplotlib inline\n",
    "    %config InlineBackend.figure_format='retina'\n",
    "except:\n",
    "    pass\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Make the size and fonts larger for this presentation\n",
    "plt.rcParams['figure.figsize'] = (10, 8)\n",
    "plt.rcParams['font.size'] = 16\n",
    "plt.rcParams['lines.linewidth'] = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15617de8",
   "metadata": {},
   "source": [
    "## Time yourself"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc612f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import wraps\n",
    "import time\n",
    "\n",
    "def timeit(func):\n",
    "    \"\"\"Prints the runtime of decorated functions\"\"\"\n",
    "    @wraps(func)\n",
    "    def wrapper_timer(*args, **kwargs):\n",
    "        start_time = time.perf_counter()\n",
    "        value = func(*args, **kwargs)\n",
    "        end_time = time.perf_counter()\n",
    "        run_time = end_time - start_time\n",
    "        print(f\"Finished {func.__name__!r} in {run_time:.4f} secs\")\n",
    "        return value\n",
    "    return wrapper_timer\n",
    "\n",
    "@timeit\n",
    "def waste_some_time(num_times):\n",
    "    for _ in range(num_times):\n",
    "        sum([i**2 for i in range(10000)])\n",
    "\n",
    "waste_some_time(999)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d310717",
   "metadata": {},
   "source": [
    "## Preprocessing of Data Files\n",
    "\n",
    "Different researchers (and industries) prefer different data formats and even in astronomy you will encounter most of the following files:\n",
    "\n",
    "| Name   | File Ending | Advantages                                                                 | Encounter | How to Read with Python                       |\n",
    "|--------|-------------|---------------------------------------------------------------------------|-|----------------------------------------------------|\n",
    "| Excel  | .xlsx       | User-friendly,<br>supports formatting,<br>widely used in business.      | Collaboration spreadsheets | `pandas.read_excel('file.xlsx')` |\n",
    "| CSV    | .csv        | Lightweight,<br>human-readable,<br>easily imported/exported,<br>universal format.    | CDS VizieR | `pandas.read_csv('file.csv')`                           |\n",
    "| FITS   | .fits       | Optimized for astronomy data,<br>supports metadata,<br>handles large datasets.    | Astronomy Survey Images and Catalogues | `table = astropy.table.Table.read('file.fits)`<br>`table.to_pandas()`     |\n",
    "| HDF5   | .hdf5  | Efficient storage for large, complex datasets,<br>supports metadata,<br>fast I/O. | Astronomy Simulations | `pandas.read_hdf('file.h5', key='data')`                |\n",
    "| Xarray | .nc         | Handles labeled n-dimensional arrays,<br>integrates with Dask for big data,<br>Well tested on TB-scale in climate modelling.    | Future Simulations/Observations? | `xarray.load_dataarray('file.nc')` |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce229270",
   "metadata": {},
   "source": [
    "The above table gives you a very incomplete idea of the vast amount of table formats out there. As our measurement and simulation data increases, we will move downwards. The last two file types allow you to only read in columns/rows of files or only their metadata - great if you only care about a small part in a terra-byte large file!\n",
    "\n",
    "Currently, most astronomers are *stuck* in the lower middle of the chart - including myself. Because I work a lot with `FITS` files, I usually work with`astropy.table`.  \n",
    "\n",
    "In terms of python packages, `pandas` is widely used in astronomy and data science. \n",
    "Switching between them is not super-difficult, so I will try to show you examples of both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ba01a79",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from astropy.table import Table # this works well for FITS data catalogues\n",
    "from astropy.io import fits # this is your more agnostic way to work for FITS images"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7841f6df",
   "metadata": {},
   "source": [
    "### Creating data directories and downloading files \n",
    "\n",
    "You might need to download a file (if you have not already).\n",
    "\n",
    "The `os` package comes in handy, because it allows you to check if a file or directory exists, and even create a directory if necessary!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd3f7741",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# Check if 'data' directory exists, if not, create it\n",
    "if not os.path.exists('data'):\n",
    "    os.makedirs('data')\n",
    "\n",
    "# You can also check if a specific file exists\n",
    "file_path = 'data/reduced_TAN_C14.fits'\n",
    "\n",
    "if not os.path.exists(file_path):\n",
    "    try:\n",
    "        # Download the file from a URL, if it does not exist\n",
    "        import wget\n",
    "        wget.download(\n",
    "            url='http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',\n",
    "            out=file_path\n",
    "        )\n",
    "    except Exception as e:\n",
    "        raise RuntimeError(f\"Could not download FITS file: {e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61b06fea",
   "metadata": {},
   "source": [
    "### Reading CSV (Comma Separated Values) Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dba84d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Option 1: Pandas\n",
    "galah_data = pd.read_csv('data/galah_dr3_allstar_m67_lite.csv')\n",
    "galah_data[:5] # let's just display the first 5 rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5733ec1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Option 2: Astropy\n",
    "galah_data = Table.read('data/galah_dr3_allstar_m67_lite.csv')\n",
    "galah_data[:5] # let's just display the first 5 rows"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd56d35",
   "metadata": {},
   "source": [
    "### FITS Images (the I in Flexible Image Transport System)\n",
    "\n",
    "also check out the astropy tutorials: https://learn.astropy.org/tutorials/FITS-cubes.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "572e8964",
   "metadata": {},
   "source": [
    "#### Metadata: store useful information about your tables\n",
    "\n",
    "Here, we will work with a FITS file from the Radio astronomy survey HI4PI.  \n",
    "The specific region of the sky that was observed in this tile (C14) is the Small Magellanic Cloud (SMC)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3851ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FITS Images\n",
    "fits_file = fits.open('data/reduced_TAN_C14.fits')\n",
    "print(fits_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e43e67c7",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "#fits_file is now a list of objects, each one corresponding to a different extension\n",
    "print('This file has {0} extensions.'.format(len(fits_file)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbca721f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#we can also just dump information for this fits file\n",
    "print(fits_file.info())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce2dd5ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print the header for the Primary Extension\n",
    "print(fits_file[0].header.items)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "890e8503",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Try-except: your helper if you think something might fail\n",
    "try:\n",
    "    # try to read out data and header\n",
    "    fits_data = fits_file[0].data\n",
    "    fits_header = fits_file[0].header\n",
    "\n",
    "except:\n",
    "    # this might have failed, if you have closed the fits_file and execute this cell out of order\n",
    "    fits_file = fits.open('data/reduced_TAN_C14.fits')\n",
    "    fits_data = fits_file[0].data\n",
    "    fits_header = fits_file[0].header\n",
    "    \n",
    "# don't forget to close the fits_file to free up memory\n",
    "fits_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91412aa1",
   "metadata": {},
   "source": [
    "#### The world coordinate system at your fingertips"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdc0b7c1",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# The pixels are reported in Galactic longitude and latitude\n",
    "# Imagine you want to translate to Right Ascension and Declination?\n",
    "\n",
    "# GLON and GLAT come in AXIS1 and AXIS2\n",
    "naxis1 = fits_header['NAXIS1']  # Number of pixels along the 1st axis\n",
    "naxis2 = fits_header['NAXIS2']  # Number of pixels along the 2nd axis\n",
    "\n",
    "# It is good practice to report FITS measurements in an agreed reference frame, \n",
    "# the World Coordinate System (WCS)\n",
    "from astropy.wcs import WCS\n",
    "wcs = WCS(fits_header, naxis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad8fb31f",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Imagine you want to just zoom in onto the SMC: \n",
    "import astropy.units as u\n",
    "\n",
    "# Define the desired latitude and longitude range\n",
    "lat_range = [-46, -40] * u.deg  # GLAT range\n",
    "lon_range = [306, 295] * u.deg  # GLON range\n",
    "\n",
    "# Using the WCS object to convert world coordinates (GLON, GLAT) to pixel coordinates\n",
    "\n",
    "# Convert (GLON_min, GLAT_min) and (GLON_max, GLAT_max) to pixel coordinates\n",
    "glon_min_pix, glat_min_pix = wcs.world_to_pixel_values(lon_range[0], lat_range[0])  # Lower-left corner\n",
    "glon_max_pix, glat_max_pix = wcs.world_to_pixel_values(lon_range[1], lat_range[1])  # Upper-right corner\n",
    "\n",
    "# Convert to integer pixel coordinates for slicing\n",
    "glon_min_pix = int(np.floor(glon_min_pix))\n",
    "glon_max_pix = int(np.ceil(glon_max_pix))\n",
    "glat_min_pix = int(np.floor(glat_min_pix))\n",
    "glat_max_pix = int(np.ceil(glat_max_pix))\n",
    "\n",
    "sub_wcs = wcs.deepcopy()  # Make a copy of the original WCS\n",
    "sub_wcs.wcs.crpix[0] -= glon_min_pix  # Adjust reference pixel for GLON\n",
    "sub_wcs.wcs.crpix[1] -= glat_min_pix  # Adjust reference pixel for GLAT\n",
    "\n",
    "sub_data_cube = fits_data[:, glat_min_pix:glat_max_pix, glon_min_pix:glon_max_pix] * u.K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26b8dea4",
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Let's estimate the 0th and 1st Moments, i.e. Intensity Map and Intensity-weighted Velocity Map \n",
    "\n",
    "cdelt3 = fits_header['CDELT3'] * u. m / u.s  # Velocity increment per channel (e.g., km/s)\n",
    "crval3 = fits_header['CRVAL3'] * u. m / u.s  # Velocity value at the reference pixel\n",
    "crpix3 = fits_header['CRPIX3']  # Reference pixel (1-based index in FITS)\n",
    "\n",
    "# Create the velocity axis (1D array for velocity channels)\n",
    "n_channels = sub_data_cube.shape[0]  # Number of velocity channels\n",
    "velocity_axis = ((np.arange(n_channels) - (crpix3 - 1)) * cdelt3 + crval3)\n",
    "\n",
    "# Select velocity between -300 and 300 km/s\n",
    "velocity_select = (velocity_axis > -300 * u.km/u.s) & (velocity_axis < 302 * u.km/u.s)\n",
    "\n",
    "# Calculate the zeroth moment: integral of I delta_v = sum(I * cdelt3)\n",
    "intensity_sum = np.nansum(sub_data_cube[velocity_select,:,:], axis=0) * cdelt3\n",
    "\n",
    "# Calculate the first moment (intensity-weighted mean velocity)\n",
    "# i.e. weighted sum of velocities / sum of intensities\n",
    "\n",
    "# Weighted sum of velocities = sum(I * velocity-center * cdelt3)\n",
    "# Note: velocity center has to be offset by half a bin from velocity_axis\n",
    "weighted_velocity_sum = np.nansum(sub_data_cube[velocity_select,:,:] * (velocity_axis[velocity_select, np.newaxis, np.newaxis] + 0.5 * cdelt3) * cdelt3, axis=0)\n",
    "\n",
    "# Calculate the first moment (intensity-weighted mean velocity)\n",
    "intensity_weighted_mean_velocity = np.where(intensity_sum > 0, weighted_velocity_sum / intensity_sum, np.nan)\n",
    "\n",
    "# Convert Intensity sum to a Column Density assuming optically thin media\n",
    "# Following Eq. (2) of this paper: https://www.aanda.org/articles/aa/pdf/2016/10/aa29178-16.pdf\n",
    "hi_column_density = (intensity_sum * 1.823 * 10**18 / (u.cm * u.cm) * u.s / u.K / u.km).to(1/u.cm**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eb5b684",
   "metadata": {},
   "source": [
    "If you want to read more about this, check out the HI4PI paper: https://ui.adsabs.harvard.edu/abs/2016A%26A...594A.116H/abstract  \n",
    "(or ask Naomi McClure-Griffiths here at ANU, one of the co-authors, and lead of the Parkes Radio observations of this amazing project)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b056dd39",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Let's do some fun plotting! \n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(10, 5),subplot_kw={'projection': sub_wcs})\n",
    "\n",
    "# image = ax.imshow(data_slice, origin='lower', cmap='RdBu_r', aspect='auto')\n",
    "image = ax.imshow(intensity_weighted_mean_velocity.to(u.km/u.s).value, origin='lower', cmap='RdBu_r', vmin=0, vmax=200)\n",
    "\n",
    "# Add a colorbar\n",
    "cbar = plt.colorbar(image, pad = 0.15)\n",
    "cbar.set_label(r'Line-of-Sight Velocity / $\\mathrm{km\\,s^{-1}}$', size=16)\n",
    "\n",
    "# Add axes labels\n",
    "ax.set_xlabel(\"Galactic Longitude / deg\", fontsize=16)\n",
    "ax.set_ylabel(\"Galactic Latitude / deg\", fontsize=16)\n",
    "ax.set_ylim(ax.get_ylim()[::-1])\n",
    "\n",
    "# Overplot column density contours \n",
    "levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22)  # Define contour levels to use\n",
    "ax.contour(hi_column_density.value, cmap='Greys_r', alpha=0.5, levels=levels)\n",
    "\n",
    "# Overlay set of RA/Dec Axes\n",
    "overlay = ax.get_coords_overlay('fk5')\n",
    "overlay.grid(color='white', ls='dotted', lw=2)\n",
    "overlay[0].set_axislabel('Right Ascension (J2000) / deg', fontsize=16)\n",
    "overlay[1].set_axislabel('Declination (J2000) / deg', fontsize=16)\n",
    "\n",
    "contour_legend = matplotlib.lines.Line2D([0], [0], color='gray', lw=2, label='HI Column Density Contours')\n",
    "ax.legend(handles=[contour_legend], loc='lower left', fontsize=12)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('plot_examples/figures/column_density.png',dpi=200,bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "# This is adjusted from a tutorial of the astropy collaboration: https://learn.astropy.org/tutorials/FITS-cubes.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae537f30",
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Let's do some fun plotting! \n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(10, 6),subplot_kw={'projection': sub_wcs})\n",
    "\n",
    "# image = ax.imshow(data_slice, origin='lower', cmap='RdBu_r', aspect='auto')\n",
    "image = ax.imshow(intensity_weighted_mean_velocity.to(u.km/u.s).value, origin='lower', cmap='RdBu_r', vmin=0, vmax=200)\n",
    "\n",
    "# Add a colorbar\n",
    "cbar = plt.colorbar(image, pad = 0.05)\n",
    "cbar.set_label(r'Line-of-Sight Velocity / $\\mathrm{km\\,s^{-1}}$', size=16)\n",
    "\n",
    "# Add axes labels\n",
    "ax.set_xlabel(\"Galactic Longitude / deg\", fontsize=16)\n",
    "ax.set_ylabel(\"Galactic Latitude / deg\", fontsize=16)\n",
    "ax.set_ylim(ax.get_ylim()[::-1])\n",
    "\n",
    "# Overplot column density contours \n",
    "# levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22)  # Define contour levels to use\n",
    "levels = (0.1, 0.5, 1, 3, 5, 7, 10)  # Define contour levels to use\n",
    "contour = ax.contour((hi_column_density / 10**21).value, cmap='Greys_r', alpha=0.5, levels=levels)\n",
    "# Add a colorbar for the column density contours\n",
    "cbar2 = plt.colorbar(contour, ax=ax, orientation='horizontal')\n",
    "cbar2.set_label(r'HI Column Density / $10^{21}\\,\\mathrm{cm^{-2}}$', size=16)\n",
    "\n",
    "# Overlay set of RA/Dec Axes\n",
    "overlay = ax.get_coords_overlay('fk5')\n",
    "overlay.grid(color='white', ls='dotted', lw=2)\n",
    "overlay[0].set_axislabel('Right Ascension (J2000) / deg', fontsize=16)\n",
    "overlay[1].set_axislabel('Declination (J2000) / deg', fontsize=16)\n",
    "\n",
    "# contour_legend = matplotlib.lines.Line2D([0], [0], color='gray', lw=2, label='HI Column Density Contours')\n",
    "# ax.legend(handles=[contour_legend], loc='lower left', fontsize=12)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('plot_examples/figures/column_density2.png',dpi=200,bbox_inches='tight')\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "# This is adjusted from a tutorial of the astropy collaboration: https://learn.astropy.org/tutorials/FITS-cubes.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fc044ce",
   "metadata": {},
   "source": [
    "## Working with Data Catalogues"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ce955dc",
   "metadata": {},
   "source": [
    "### Downloading Data via ADQL/SQL\n",
    "\n",
    "Remember how we read in a data table from the GALAH survey above with spectroscopic measurements?\n",
    "\n",
    "The table has the identifier of the 1.8 billion source large Gaia DR3 catalogue in it.\n",
    "\n",
    "Downloading the full Gaia DR3 catalogue is not useful, but we can use the `astroquery` package to perform a query with the Astronomical Data Query Language (ADQL), an extension of the Structured Query Language (SQL) to include functions when querying data.\n",
    "\n",
    "Let's try to download the matches in Gaia DR3!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b85c78f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.gaia import Gaia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "835bd336",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract only the Gaia DR3 source IDs for crossmatching\n",
    "galah_data['source_id'] = galah_data['dr3_source_id']\n",
    "gaia_source_ids = galah_data['dr3_source_id'].tolist()\n",
    "\n",
    "# Convert the source IDs to an Astropy table to use in the query (to not upload too much data)\n",
    "# with the actual identifier in Gaia DR3 (source_id)\n",
    "gaia_source_ids_table = Table([gaia_source_ids], names=['source_id'])\n",
    "\n",
    "# Define and execute the ADQL query to crossmatch with Gaia DR3\n",
    "query = f\"\"\"\n",
    "SELECT * \n",
    "FROM gaiadr3.gaia_source AS gaia\n",
    "JOIN TAP_UPLOAD.t1 AS galah\n",
    "ON gaia.source_id = galah.source_id\n",
    "\"\"\"\n",
    "\n",
    "# Upload the source_id table for crossmatching\n",
    "job = Gaia.launch_job_async(query=query, upload_resource=gaia_source_ids_table, upload_table_name=\"t1\")\n",
    "gaiadr3_match = job.get_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4920baf",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gaiadr3_match[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7601917e",
   "metadata": {},
   "source": [
    "### Joining 2 catalogues\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80bd7e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fix missing keyword\n",
    "galah_data['source_id'] = galah_data['dr3_source_id']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dd00cd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import join\n",
    "\n",
    "gaia_dr3_galah = join(gaiadr3_match, galah_data, keys='source_id')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c13db2d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_dr3_galah"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13740dd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cmd_and_kiel(data, colormap_left = 'snr_c2_iraf', colormap_left_label = 'GALAH DR3 SNR CCD2'):\n",
    "\n",
    "    # compare Gaia DR3 and GALAH DR3 measurements\n",
    "    f, gs = plt.subplots(1, 2, figsize=(10,4))\n",
    "\n",
    "    # Left panel (Color-Magnitude Diagram, CMD)\n",
    "    ax = gs[0]\n",
    "    ax.text(0.05, 0.95, 'a)', transform=ax.transAxes, fontsize=14, ha = 'left', va='top')\n",
    "\n",
    "    sc = ax.scatter(\n",
    "        data['bp_rp'],\n",
    "        data['phot_g_mean_mag'] + 5 * np.log10(data['parallax']/100.), \n",
    "        c=data[colormap_left],\n",
    "        cmap='viridis', s=10\n",
    "    )\n",
    "    ax.invert_yaxis()\n",
    "    ax.set_xlabel(r'$G_\\mathrm{BP} - G_\\mathrm{RP}~/~\\mathrm{mag}$')\n",
    "    ax.set_ylabel(r'$G~/~\\mathrm{mag}$')\n",
    "    ax.set_title(r'$Gaia$ DR3')\n",
    "\n",
    "    # Adding a colorbar for SNR in the left panel\n",
    "    cbar = plt.colorbar(sc, ax=ax)\n",
    "    cbar.set_label(colormap_left_label)\n",
    "\n",
    "    # Right panel (Teff-logg diagram)\n",
    "    ax = gs[1]\n",
    "    ax.text(0.05, 0.95, 'b)', transform=ax.transAxes, fontsize=14, ha = 'left', va='top')\n",
    "    sc2 = ax.scatter(\n",
    "        data['teff'],\n",
    "        data['logg'], \n",
    "        c=data['fe_h'],\n",
    "        cmap='plasma', s=10\n",
    "    )\n",
    "    ax.invert_xaxis()\n",
    "    ax.invert_yaxis()\n",
    "    ax.set_xlabel(r'$T_\\mathrm{eff}~/~\\mathrm{K}$')\n",
    "    ax.set_ylabel(r'$\\log (g~/~\\mathrm{cm\\,s^{-2}}$')\n",
    "    ax.set_title('GALAH DR3')\n",
    "\n",
    "    # Adding a colorbar for [Fe/H] in the right panel\n",
    "    cbar2 = plt.colorbar(sc2, ax=ax)\n",
    "    cbar2.set_label('[Fe/H]')\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    plt.close()\n",
    "    \n",
    "plot_cmd_and_kiel(data=gaia_dr3_galah)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4ee40e4",
   "metadata": {},
   "source": [
    "### Cleaning Catalogues"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "715a68fc",
   "metadata": {},
   "source": [
    "Quite often, catalogues do not have all measurements.\n",
    "\n",
    "Just above, you have seen the error message\n",
    "```\n",
    "RuntimeWarning: invalid value encountered in log10\n",
    "gaia_dr3_galah['phot_g_mean_mag'] + 5 * np.log10(gaia_dr3_galah['parallax']/100.),\n",
    "```\n",
    "\n",
    "This message is expected if np.log10() is applied to a value that is not positive (this could be either negative value or value that is \"Not a Value\" aka NaN). Let's get a first idea:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0d73a96",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_dr3_galah['parallax']/100."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e71482c9",
   "metadata": {},
   "source": [
    "We already see negative values, but we can also imagine values that are *Not a Number* aka NaN or they might have a bitmask \"flag\" that indicates their quality.\n",
    "\n",
    "Negative parallax measurements are true measurements! They just tell us that the source is quite far away. Later in this course how we can use our **prior** knowledge that distances from us have to be positive to still extract something useful out of these measures.\n",
    "\n",
    "For now, we will simply identify these measurements and not use them to avoid error messages.\n",
    "For a research paper, this would be a selection cut that has needs to be documented for reproducability!\n",
    "\n",
    "**NaN entries**  \n",
    "\n",
    "You can identify NaN entries with the check, e.g. effective temperature $T_\\mathrm{eff}$.\n",
    "If you have checked if a value is NaN, you can also invert the result with a `~` (switch True <-> False):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87511ff1",
   "metadata": {},
   "outputs": [],
   "source": [
    "check_if_parallax_values_finite = np.isfinite(gaia_dr3_galah['parallax'])\n",
    "check_if_parallax_values_nan = np.isnan(gaia_dr3_galah['parallax'])\n",
    "check_if_parallax_values_not_fintie = ~np.isfinite(gaia_dr3_galah['parallax'])\n",
    "check_if_parallax_positive = gaia_dr3_galah['parallax'] > 0\n",
    "\n",
    "# you can also use the np.where function:\n",
    "where_parallax_not_positive = np.where(~check_if_parallax_positive==True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50f24e6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_dr3_galah['parallax'][where_parallax_not_positive]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6440ac5f",
   "metadata": {},
   "source": [
    "**bitmask flags**\n",
    "\n",
    "A bitmask is a way of storing multiple boolean (True/False) values in a single integer by representing each condition with a different bit in the binary representation of the number. For example, you might use bitmasks in a catalog to encode several flags in a compact form.\n",
    "\n",
    "You can for example imagine a bitmask flag `0111`, which would add up to `8*0 + 4*1 + 2*1 + 1*1 = 7`. So with just 1 number, you can check for 4 different things!\n",
    "\n",
    "You can also check if just a specific bitmask is raised:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47c78a53",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_bit_raised(flag, bit):\n",
    "    \"\"\"\n",
    "    Check if *flag* has the value for a specific *bit* raised\n",
    "    \"\"\"\n",
    "    return (flag & bit) != 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41bc1245",
   "metadata": {},
   "source": [
    "In GALAH DR3, one of these quality flags is `flag_sp`, which includes a lot of details about quality checks and extra information about a star, for example, if we think it is a binary star.\n",
    "\n",
    "According to the documentation (Table 4 from https://ui.adsabs.harvard.edu/abs/2021MNRAS.506..150B)\n",
    "that would mean raised flags 32 (spectroscopic binary) or 64 (photometric binary).\n",
    "\n",
    "We can find the stars that are flagged as binaries with an `|` check, which means `or`.\n",
    "You could also find stars that are spectroscopic *and* photometric binaries by replacing `|` with `&`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dde211a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "check_if_binary = (\n",
    "    is_bit_raised(gaia_dr3_galah['flag_sp'], 32) &\n",
    "    is_bit_raised(gaia_dr3_galah['flag_sp'], 64)\n",
    ")\n",
    "\n",
    "gaia_dr3_galah['binary'] = check_if_binary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44fc4a14",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cmd_and_kiel(\n",
    "    gaia_dr3_galah[check_if_parallax_positive],\n",
    "    colormap_left='binary',\n",
    "    colormap_left_label='Binary True/False'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a4db36d",
   "metadata": {},
   "source": [
    "### Scaling entries (e.g. for fitting)\n",
    "\n",
    "You can imagine that any type of optimisation of a function has to start with an initial value.\n",
    "In some cases, this value might be off by quite a lot and the fitting algorithm might be stuck, or take a lot of time to get close to the final result.\n",
    "\n",
    "There are 2 ways to improve your situation:\n",
    "1. Start with a better initial guess  \n",
    "2. Rescale your values, so that they are between -1 and 1 or scaled to their standard deviation  \n",
    "\n",
    "Take a look at: https://scikit-learn.org/stable/modules/preprocessing.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0fbf899",
   "metadata": {},
   "source": [
    "### Visualising targets with ipyaladin\n",
    "\n",
    "You can find ipyaladin at https://github.com/cds-astro/ipyaladin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13fec866",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from ipyaladin import Aladin\n",
    "aladin = Aladin()\n",
    "aladin\n",
    "\n",
    "# Initialize Aladin and set the view to M67 in 2MASS J color image\n",
    "aladin = Aladin(target='M67', survey='2MASS-J', fov=0.5)  # Adjust the field of view as needed\n",
    "\n",
    "# Display Aladin widget\n",
    "aladin"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29858391",
   "metadata": {},
   "source": [
    "## Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d38d4fe7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import CubicSpline\n",
    "from scipy.interpolate import interp1d\n",
    "\n",
    "# sample sparse data points from a sine function\n",
    "sample_x = np.linspace(-np.pi, np.pi, 7)\n",
    "sample_y = np.sin(sample_x)\n",
    "plt.plot(sample_x, sample_y, 'ko')\n",
    "\n",
    "# the true sine function\n",
    "z = np.linspace(np.min(sample_x), np.max(sample_x), 100)\n",
    "plt.plot(z, np.sin(z), label='real')\n",
    "\n",
    "# linear interpolation\n",
    "f1 = interp1d(sample_x, sample_y)\n",
    "plt.plot(z, f1(z), label='interp1d')\n",
    "\n",
    "# cubic spline interpolation\n",
    "f2 = CubicSpline(sample_x, sample_y)\n",
    "plt.plot(z, f2(z), label='CubicSpline')\n",
    "\n",
    "# add legend\n",
    "plt.legend();\n",
    "plt.grid(linestyle='--');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9136e90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's plot the residual\n",
    "plt.plot(z, f1(z)-np.sin(z), label='interp1d residuals', color='orange')\n",
    "plt.plot(z, f2(z)-np.sin(z), label='CubicSpline residuals', color='green')\n",
    "plt.legend();\n",
    "plt.grid(linestyle='--');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59baded6",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.interpolate import RegularGridInterpolator\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# Example: 3D linear interpolation\n",
    "\n",
    "# Define the grid points in each dimension\n",
    "x = np.linspace(0, 10, 5)  # 5 points in the x direction\n",
    "y = np.linspace(0, 10, 5)  # 5 points in the y direction\n",
    "z = np.linspace(0, 10, 5)  # 5 points in the z direction\n",
    "\n",
    "# Create a 3D grid of values (you can replace this with actual data)\n",
    "# For demonstration, we'll just use the sum of coordinates for simplicity\n",
    "values = np.zeros((len(x), len(y), len(z)))\n",
    "\n",
    "for i in range(len(x)):\n",
    "    for j in range(len(y)):\n",
    "        for k in range(len(z)):\n",
    "            values[i, j, k] = x[i] + y[j] + z[k]\n",
    "\n",
    "# Create the interpolator object\n",
    "interpolator = RegularGridInterpolator((x, y, z), values)\n",
    "\n",
    "# Generate a finer 3D grid for visualization using interpolation\n",
    "x_fine = np.linspace(0, 10, 30)\n",
    "y_fine = np.linspace(0, 10, 30)\n",
    "z_fine = np.linspace(0, 10, 30)\n",
    "\n",
    "x_fine_mesh, y_fine_mesh, z_fine_mesh = np.meshgrid(x_fine, y_fine, z_fine)\n",
    "\n",
    "# Prepare points for interpolation\n",
    "fine_points = np.array([x_fine_mesh.ravel(), y_fine_mesh.ravel(), z_fine_mesh.ravel()]).T\n",
    "fine_values = interpolator(fine_points).reshape(x_fine_mesh.shape)\n",
    "\n",
    "# Plotting\n",
    "fig = plt.figure(figsize=(10, 8))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# Plot the original grid points in red\n",
    "x_grid, y_grid, z_grid = np.meshgrid(x, y, z, indexing='ij')\n",
    "ax.scatter(x_grid, y_grid, z_grid, color='red', s=50, label='Original Grid Points')\n",
    "\n",
    "# Plot interpolated values as a scatter plot on the fine grid\n",
    "scatter = ax.scatter(x_fine_mesh, y_fine_mesh, z_fine_mesh, c=fine_values, cmap='viridis', alpha=0.7)\n",
    "\n",
    "# Add color bar for interpolated values\n",
    "colorbar = fig.colorbar(scatter, ax=ax, pad=0.1, label='Interpolated Value')\n",
    "\n",
    "# Labels and title\n",
    "ax.set_title('3D Linear Interpolation')\n",
    "ax.set_xlabel('X')\n",
    "ax.set_ylabel('Y')\n",
    "ax.set_zlabel('Z')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08d2baa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.interpolate import LinearNDInterpolator\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# Example: 3D interpolation on a non-regular grid\n",
    "\n",
    "# Generate random non-regular grid points in 3D space\n",
    "np.random.seed(42)\n",
    "num_points = 30  # Number of scattered data points\n",
    "x = np.random.uniform(0, 10, num_points)\n",
    "y = np.random.uniform(0, 10, num_points)\n",
    "z = np.random.uniform(0, 10, num_points)\n",
    "\n",
    "# Values at the random grid points (for example, we'll use the sum of the coordinates)\n",
    "values = x + y + z\n",
    "\n",
    "# Create the LinearNDInterpolator object for scattered data\n",
    "interpolator = LinearNDInterpolator(list(zip(x, y, z)), values)\n",
    "\n",
    "# Generate a finer grid for visualization (in 3D space)\n",
    "x_fine = np.linspace(0, 10, 30)\n",
    "y_fine = np.linspace(0, 10, 30)\n",
    "z_fine = np.linspace(0, 10, 30)\n",
    "\n",
    "# Meshgrid for finer points\n",
    "x_fine_mesh, y_fine_mesh, z_fine_mesh = np.meshgrid(x_fine, y_fine, z_fine)\n",
    "\n",
    "# Prepare the fine points for interpolation\n",
    "fine_points = np.array([x_fine_mesh.ravel(), y_fine_mesh.ravel(), z_fine_mesh.ravel()]).T\n",
    "\n",
    "# Perform interpolation at the fine points\n",
    "fine_values = interpolator(fine_points)\n",
    "\n",
    "# Reshape the results to match the fine grid shape\n",
    "fine_values = fine_values.reshape(x_fine_mesh.shape)\n",
    "\n",
    "# Plotting\n",
    "fig = plt.figure(figsize=(10, 8))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "# Plot the original scattered grid points in red\n",
    "ax.scatter(x, y, z, color='red', s=50, label='Original Scattered Points')\n",
    "\n",
    "# Plot the interpolated values as a scatter plot on the fine grid\n",
    "# Mask invalid (NaN) values that could occur where interpolation is not defined\n",
    "valid_mask = ~np.isnan(fine_values)\n",
    "scatter = ax.scatter(\n",
    "    x_fine_mesh[valid_mask],\n",
    "    y_fine_mesh[valid_mask],\n",
    "    z_fine_mesh[valid_mask],\n",
    "    c=fine_values[valid_mask],\n",
    "    cmap='viridis',\n",
    "    alpha=0.7,\n",
    ")\n",
    "\n",
    "# Add a color bar to represent the interpolated values\n",
    "colorbar = fig.colorbar(scatter, ax=ax, pad=0.1, label='Interpolated Value')\n",
    "\n",
    "# Labels and title\n",
    "ax.set_title('3D Interpolation on a Non-Regular Grid')\n",
    "ax.set_xlabel('X')\n",
    "ax.set_ylabel('Y')\n",
    "ax.set_zlabel('Z')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc7f19",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "279.59375px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}