import_comcam.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import lsst.daf.butler as dafButler\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "warnings.simplefilter('ignore')\n",
    "plt.set_loglevel('WARNING')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install git+https://github.com/astronomy-commons/lsdb.git@rubin_index\n",
    "%pip install git+https://github.com/astronomy-commons/hats.git@main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hats\n",
    "print(hats.__version__)\n",
    "import lsdb\n",
    "print(lsdb.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load collection\n",
    "repo = \"embargo_new\"\n",
    "collection = \"LSSTComCam/runs/DRP/20241101_20241211/w_2024_50/DM-48128\"\n",
    "instrument = \"LSSTComCam\"\n",
    "\n",
    "butler = dafButler.Butler(repo, collections=collection, instrument=instrument)\n",
    "# older collections\n",
    "butler_49 = dafButler.Butler(repo, collections=\"LSSTComCam/runs/DRP/20241101_20241204/w_2024_49/DM-47988\", instrument=instrument)\n",
    "butler_48 = dafButler.Butler(repo, collections=\"LSSTComCam/runs/DRP/20241101_20241127/w_2024_48/DM-47841\", instrument=instrument)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "bands = [\"u\", \"g\", \"r\", \"i\", \"z\", \"y\"]\n",
    "verbose = False\n",
    "\n",
    "for band in bands:\n",
    "    datasetRefs_dia = butler.query_datasets(\"goodSeeingDiff_differenceExp\", where=f\"band='{band}'\")\n",
    "\n",
    "    if verbose:\n",
    "        print(f\"\\nDataset references for band '{band}':\")\n",
    "        for dr in datasetRefs_dia:\n",
    "            print(dr)\n",
    "\n",
    "    print(f\"Found {len(datasetRefs_dia)} differenceExps for band '{band}'\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forcedsrc_refs = butler.query_datasets(\"forcedSourceTable\")\n",
    "obj_refs = butler.query_datasets(\"objectTable\")\n",
    "\n",
    "print(f\"Number of tracts in forcedsrc_refs (current): {len(forcedsrc_refs)}\")\n",
    "print(f\"Number of tracts in obj_refs (current): {len(obj_refs)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "diaSrc_refs = butler.query_datasets(\"diaSourceTable_tract\")\n",
    "diaObj_refs = butler.query_datasets(\"diaObjectTable_tract\")\n",
    "diaObj_refs_49 = butler_49.query_datasets(\"diaObjectTable_tract\")\n",
    "diaObj_refs_48 = butler_48.query_datasets(\"diaObjectTable_tract\")\n",
    "\n",
    "print(f\"Number of tracts in diaSrc_refs (current): {len(diaSrc_refs)}\")\n",
    "print(f\"Number of tracts in diaObj_refs (current): {len(diaObj_refs)}\")\n",
    "print(f\"Number of tracts in diaObj_refs_49 (weekly 49): {len(diaObj_refs_49)}\")\n",
    "print(f\"Number of tracts in diaObj_refs_48 (weekly 48): {len(diaObj_refs_48)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# List of tract indices to process\n",
    "tract_indices = np.arange(0, len(diaObj_refs))\n",
    "\n",
    "# Initialize lists to hold diaObj and diaSource DataFrames\n",
    "dia_Obj_list = []\n",
    "diaSource_list = []\n",
    "\n",
    "# Loop through the specified tracts and fetch the data\n",
    "for idx in tqdm(tract_indices):\n",
    "    dia_Obj = butler.get('diaObjectTable_tract', dataId=diaObj_refs[idx].dataId)  # Fetch diaObj for this tract\n",
    "    diaSource = butler.get('diaSourceTable_tract', dataId=diaSrc_refs[idx].dataId)  # Fetch diaSource for this tract\n",
    "    \n",
    "    # Append the results to the lists\n",
    "    dia_Obj_list.append(dia_Obj)\n",
    "    diaSource_list.append(diaSource)\n",
    "\n",
    "# Concatenate all diaObj and diaSource DataFrames\n",
    "combined_dia_Obj = pd.concat(dia_Obj_list, ignore_index=False)\n",
    "combined_diaSource = pd.concat(diaSource_list, ignore_index=False)\n",
    "\n",
    "dia_Obj = combined_dia_Obj.copy()\n",
    "diaSource = combined_diaSource.copy()\n",
    "\n",
    "# assuming the lastest one is weekly 50\n",
    "dia_Obj_50 = combined_dia_Obj.copy()\n",
    "diaSource_50 = combined_diaSource.copy()\n",
    "\n",
    "# Compute the histogram of nDiaSources\n",
    "counts, bin_edges = np.histogram(combined_dia_Obj['nDiaSources'], bins=40, range=(0, 40))\n",
    "\n",
    "# Take the logarithm (base 10) of the counts\n",
    "log_counts = np.log10(counts + 1)  # Add 1 to avoid log(0)\n",
    "\n",
    "# Plot the histogram\n",
    "plt.bar(bin_edges[:-1], log_counts, width=np.diff(bin_edges), align='edge', edgecolor='black', alpha=0.7)\n",
    "plt.xlabel('Number of diaSources per diaObject')\n",
    "plt.ylabel('Log10(Number of diaObjects)')\n",
    "plt.title('Log Number of diaObjects per Number of Sources per Object')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_diaSource_50 = lsdb.from_dataframe(\n",
    "    diaSource,\n",
    "    catalog_name=\"diaSource_50\",\n",
    "    catalog_type=\"source\")\n",
    "\n",
    "catalog_diaObj_50 = lsdb.from_dataframe(\n",
    "    dia_Obj,\n",
    "    catalog_name=\"diaObject_50\",\n",
    "    catalog_type=\"object\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save it to disk in HATS format\n",
    "catalog_diaSource_50.to_hats(\"diaSource_50\", overwrite=True)\n",
    "catalog_diaObj_50.to_hats(\"diaObject_50\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forced Sources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Process forcedSourceTable and objectTable (812 tracts)\n",
    "forcedSource_list = []\n",
    "\n",
    "for idx in tqdm(range(len(forcedsrc_refs))):\n",
    "    forcedSource = butler.get('forcedSourceTable', dataId=forcedsrc_refs[idx].dataId)  # Fetch forcedSource\n",
    "    # Append to the lists\n",
    "    forcedSource_list.append(forcedSource)\n",
    "\n",
    "# Concatenate forcedSource and object DataFrames\n",
    "forcedSource = pd.concat(forcedSource_list, ignore_index=False)\n",
    "del forcedSource_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_forcedSource_50 = lsdb.from_dataframe(\n",
    "    forcedSource,\n",
    "    catalog_name=\"forcedSource_50\",\n",
    "    catalog_type=\"source\",\n",
    "    ra_column=\"coord_ra\",\n",
    "    dec_column=\"coord_dec\",\n",
    "    highest_order=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_forcedSource_50.to_hats(\"forcedSource_50\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Process forcedSourceTable and objectTable (812 tracts)\n",
    "object_list = []\n",
    "\n",
    "for idx in tqdm(range(len(forcedsrc_refs))):\n",
    "    obj = butler.get('objectTable', dataId=obj_refs[idx].dataId)  # Fetch object table\n",
    "    # Append to the lists\n",
    "    object_list.append(obj)\n",
    "\n",
    "# Concatenate forcedSource and object DataFrames\n",
    "Objects = pd.concat(object_list, ignore_index=False)\n",
    "del object_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_Obj_50 = lsdb.from_dataframe(\n",
    "    Objects,\n",
    "    catalog_name=\"Obj_50\",\n",
    "    catalog_type=\"object\",\n",
    "    ra_column=\"coord_ra\",\n",
    "    dec_column=\"coord_dec\",\n",
    "    highest_order=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_Obj_50.to_hats(\"Object_50\", overwrite=True)"
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}