diff --git a/Projects/Bulge_exploration/Explore_Bulge_Part1A_IngestDataSet.ipynb b/Projects/Bulge_exploration/Explore_Bulge_Part1A_IngestDataSet.ipynb new file mode 100644 index 0000000..0e5d827 --- /dev/null +++ b/Projects/Bulge_exploration/Explore_Bulge_Part1A_IngestDataSet.ipynb @@ -0,0 +1,213 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bulge data ingestion\n", + "\n", + "In this notebook we show how to ingest the raw DECam images. We are using one of the five fields of the bulge data set from [Saha et al. (2019)](https://arxiv.org/abs/1902.05637). The field has been observed from Apr. 1 - 3, 2015, with 43 visits and 62 detectors. The raw data is available on `/project/stack-club/course_data/decam-bulge-rawdata`. \n", + "\n", + "The raw images occupy around 50 GB of disk space.\n", + "\n", + "**Ingesting and processing this one field takes several days and around 300 GB of disk space (42 visits x 62 CCDs).**\n", + "\n", + "We have created another Jupyter notebook `Explore_Bulge_singleVisit.ipynb`, which only ingest a single visit and a single CCD speeding up the process of ingestion to a few minutes. The whole ingested data set can be found `/project/stack-club/course_data/DECAM_BULGE`.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import sqlite3\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setting up the directories." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Raw data directory\n", + "DATA_DIR = '/project/stack-club/course_data/decam-bulge-rawdata/'\n", + "\n", + "# Repo directory\n", + "REPO_DIR = '/project/stack-club/course_data/DECAM_BULGE/' \n", + "CALIB_DIR = REPO_DIR + \"CALIB\"\n", + "RERUN_DIR = REPO_DIR + \"rerun\"\n", + "\n", + "# Create the REPO_DIR and add the mapper.\n", + "! mkdir -p {REPO_DIR}\n", + "! echo \"lsst.obs.decam.DecamMapper\" > {REPO_DIR+\"_mapper\"}\n", + "! mkdir -p {CALIB_DIR}\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Following the standard procedure, we start by linking the reference catalogs to the respository directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "args = REPO_DIR+'/ref_cats'\n", + "!mkdir {args}\n", + "\n", + "args = '/datasets/refcats/htm/v1/sdss-dr9-fink-v5b/ '+REPO_DIR+'/ref_cats/sdss'\n", + "! ln -s {args}\n", + "\n", + "args = '/datasets/refcats/htm/v1/ps1_pv3_3pi_20170110/ '+REPO_DIR+'/ref_cats/ps1_pv3_3pi_20170110'\n", + "! ln -s {args}\n", + "\n", + "args = '/datasets/refcats/htm/v1/gaia_DR1_v1/ '+REPO_DIR+'/ref_cats/gaia'\n", + "! ln -s {args}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## We have to ingest several calibration images like crosstalk, defects, flats, biases and the images itself:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting crosstalk\n", + "args = REPO_DIR + ' --calib ' + CALIB_DIR + ' ' + DATA_DIR + 'DECam/crosstalk/'\n", + "! ingestCuratedCalibs.py {args}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting defects\n", + "args = REPO_DIR + ' --calib ' + CALIB_DIR + ' ' + DATA_DIR + 'DECam/defects/'\n", + "! ingestCuratedCalibs.py {args}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting the raw images\n", + "args = REPO_DIR + ' ' + DATA_DIR + \"raw/c4d_*.fits.fz --mode=link\"\n", + "! ingestImages.py {args} " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting flat and bias images\n", + "args = REPO_DIR + ' --calib ' + CALIB_DIR + ' ' + DATA_DIR + \\\n", + " \"cal/c4d_150403*.fits.fz --mode=link --validity 999 \"\n", + "! ingestCalibs.py {args}\n", + "\n", + "#Works only if I ingest one night of cal data????" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run `processCcd.py` in command line!!\n", + "\n", + "`processCcd.py` will several image calibration and analsis steps, including astrometric solution, source detection and photometry.\n", + "\n", + "\n", + "This will take a lot of time and space. It will also output the log and fill-up the notebook with data, therefore we run this command in a terminal.\n" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "'''\n", + "#args = REPO_DIR + \" --calib \" + CALIB_DIR + \" --rerun \" + RERUN_DIR + \\\n", + "# \" --id \" \n", + "#! processCcd.py {args}\n", + "#print(args)\n", + "'''" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "processCcd.py /project/stack-club/course_data/DECAM_BULGE/ --calib /project/stack-club/course_data/DECAM_BULGE/CALIB --rerun /project/stack-club/course_data/DECAM_BULGE/rerun --id > /home/mrabus/process_log.txt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## After ingesting and processing the images we can continue with Part 2." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.7.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Projects/Bulge_exploration/Explore_Bulge_Part1B_IngestSingleVisit.ipynb b/Projects/Bulge_exploration/Explore_Bulge_Part1B_IngestSingleVisit.ipynb new file mode 100644 index 0000000..c7b129c --- /dev/null +++ b/Projects/Bulge_exploration/Explore_Bulge_Part1B_IngestSingleVisit.ipynb @@ -0,0 +1,380 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bulge data small ingestion\n", + "\n", + "In this notebook we show how to ingest the raw DECam images. We are using one of the five fields of the bulge data set from [Saha et al. (2019)](https://arxiv.org/abs/1902.05637). The field has been observed from Apr. 1 - 3, 2015, with 43 visits and 62 detectors. The raw data is available on `/project/stack-club/course_data/decam-bulge-rawdata`. Here we will ingest only one visit and one CCD per night. . \n", + "\n", + "The raw images occupy around 50 GB of disk space.\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import sqlite3\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "import lsst.daf.persistence as dafPersist\n", + "import lsst.afw.table as afwTable\n", + "import lsst.afw.display as afwDisplay\n", + "\n", + "afwDisplay.setDefaultBackend('matplotlib') \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setting up the directories." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Raw data directory\n", + "DATA_DIR = '/project/stack-club/course_data/decam-bulge-rawdata/'\n", + "\n", + "# Repo directory, we will create a Butler repository in the home directory\n", + "REPO_DIR = '/home/mrabus/DATA/DECAM_BULGE_SingleVisit/' \n", + "CALIB_DIR = REPO_DIR + \"CALIB\"\n", + "RERUN_DIR = REPO_DIR + \"rerun\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#if repo directory does not exist create it\n", + "\n", + "if not os.path.isdir(REPO_DIR):\n", + " # Create the REPO_DIR and add the mapper.\n", + " ! mkdir -p {REPO_DIR}\n", + " ! echo \"lsst.obs.decam.DecamMapper\" > {REPO_DIR+\"_mapper\"}\n", + " ! mkdir -p {CALIB_DIR}\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Following the standard procedure, we start by linking the reference catalogs to the respository directory.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#args = '/project/shared/decamBulge/saha2/ref_cats ' + REPO_DIR\n", + "#! ln -s {args}\n", + "\n", + "args = REPO_DIR+'/ref_cats'\n", + "!mkdir {args}\n", + "\n", + "args = '/datasets/refcats/htm/v1/sdss-dr9-fink-v5b/ '+REPO_DIR+'/ref_cats/sdss'\n", + "! ln -s {args}\n", + "\n", + "args = '/datasets/refcats/htm/v1/ps1_pv3_3pi_20170110/ '+REPO_DIR+'/ref_cats/ps1_pv3_3pi_20170110'\n", + "! ln -s {args}\n", + "\n", + "args = '/datasets/refcats/htm/v1/gaia_DR1_v1/ '+REPO_DIR+'/ref_cats/gaia'\n", + "! ln -s {args}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### We have to ingest several calibration images like crosstalk, defects, flats, biases and the images itself:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting crosstalk\n", + "args = REPO_DIR + ' --calib ' + CALIB_DIR + ' ' + DATA_DIR + 'DECam/crosstalk/'\n", + "! ingestCuratedCalibs.py {args}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting defects\n", + "args = REPO_DIR + ' --calib ' + CALIB_DIR + ' ' + DATA_DIR + 'DECam/defects/'\n", + "! ingestCuratedCalibs.py {args}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# ingesting raw DECam images\n", + "args = REPO_DIR + ' ' + DATA_DIR + \"raw/c4d_*.fits.fz --mode=link\"\n", + "! ingestImages.py {args} " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#ingesting flat and bias images\n", + "args = REPO_DIR + ' --calib ' + CALIB_DIR + ' ' + DATA_DIR + \\\n", + " \"cal/c4d_150403*.fits.fz --mode=link --validity 999 \"\n", + "! ingestCalibs.py {args}\n", + "\n", + "#Works only if I ingest one night of cal data????" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Processing single visit and one CCD for each night!\n", + "\n", + "Now we are processing a single visit on each night:\n", + "\n", + "night 150402, visit = 427320, CCD = 1\n", + "\n", + "night 150403, visit = 427940, CCD = 1\n", + "\n", + "night 150404, visit = 428626, CCD = 1\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "args = REPO_DIR + \" --calib \" + CALIB_DIR + \" --rerun \" + RERUN_DIR + \\\n", + " \" --id visit=427320 ccdnum=1 \" \n", + "! processCcd.py {args}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "args = REPO_DIR + \" --calib \" + CALIB_DIR + \" --rerun \" + RERUN_DIR + \\\n", + " \" --id visit=427940 ccdnum=1 \" \n", + "! processCcd.py {args}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "args = REPO_DIR + \" --calib \" + CALIB_DIR + \" --rerun \" + RERUN_DIR + \\\n", + " \" --id visit=428626 ccdnum=1 \" \n", + "! processCcd.py {args}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## After processing the images we can start with the butler and displaying images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# point Butler to the current rerun directory\n", + "butler = dafPersist.Butler(RERUN_DIR)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "good_visits = []\n", + "\n", + "# show available data\n", + "metadata = butler.queryMetadata('calexp',['visit','ccd','filter'])\n", + "for dataset in metadata:\n", + " dataId = {'visit': int(dataset[0]), 'ccd': int(dataset[1]), 'filter':dataset[2]}\n", + " if butler.datasetExists('src', dataId=dataId):\n", + " print(dataId)\n", + " \n", + " good_visits.append(dataId['visit'])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have three visits on three different nights.\n", + "Next, plot for each visit calexp." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "CCD = 1\n", + "\n", + "for ii,visit in enumerate(good_visits):\n", + "\n", + " calexp = butler.get('calexp', visit=int(visit), ccd=CCD)\n", + "\n", + " plt.figure(ii)\n", + " display = afwDisplay.Display(frame=ii, backend='matplotlib')\n", + " display.scale(\"linear\", \"zscale\")\n", + " display.mtv(calexp[500:1500,2500:3000])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "All three visits on the different nights correspond to the same pointing. We see from the images, that for the last night the seeing was not good." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Create a source catalog dictionary for the three visits and print for each visit number of detected sources\n", + "\n", + "srcCatalog = {}\n", + "for visit in good_visits:\n", + "\n", + " dataId={'visit': visit, 'ccd': 1} #This dataId exists\n", + " srcCatalog[visit] = butler.get('src', dataId=dataId)\n", + " print( 'visit {}, number of detections {}'.format(visit, len(srcCatalog[visit])) )\n", + " print( 'visit {}, median effective PSF area {}'.format(visit, np.median(srcCatalog[visit]['base_PsfFlux_area'])) )\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Because for the last night the seeing was not good, the number of detected objects is lower and the effective PSF area is larger.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "srcCatalog[427320].asAstropy().to_pandas().head()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot for each visit the detected source \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "markers = ['+', 'x', '2']\n", + "\n", + "plt.figure(figsize=(10,10))\n", + "\n", + "for ii, visit in enumerate(good_visits):\n", + "\n", + " plt.scatter(srcCatalog[visit]['coord_ra'], srcCatalog[visit]['coord_dec'], marker=markers[ii], label='{}'.format(visit))\n", + "\n", + "plt.legend()\n", + "plt.xlabel('RA')\n", + "plt.ylabel('Dec')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It looks very crowdy. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.7.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/Projects/Bulge_exploration/Explore_Bulge_Part2_ExploreData.ipynb b/Projects/Bulge_exploration/Explore_Bulge_Part2_ExploreData.ipynb new file mode 100644 index 0000000..37e21f3 --- /dev/null +++ b/Projects/Bulge_exploration/Explore_Bulge_Part2_ExploreData.ipynb @@ -0,0 +1,392 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Bulge data exploration\n", + "\n", + "In Part A, we showed how to ingest raw DECam images and to process them. A processed data set is available at `/project/stack-club/course_data/DECAM_BULGE`. In this notebook we will use the Butler to explore the processed data. \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Make plots available to the notebook\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy.table import Table\n", + "from astropy.visualization import hist\n", + "import matplotlib.pyplot as plt\n", + "import lsst.daf.persistence as dafPersist\n", + "import lsst.afw.table as afwTable\n", + "import lsst.afw.display as afwDisplay\n", + "\n", + "afwDisplay.setDefaultBackend('matplotlib') \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Repo directory and rerun directory\n", + "REPO_DIR = '/project/stack-club/course_data/DECAM_BULGE/' \n", + "RERUN_DIR = REPO_DIR + \"rerun\"\n", + "\n", + "# Directory where we save our data, like pandas data frames, tables, light curves ...\n", + "parquet_save_path = '/home/mrabus/DATA/'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## After processing the images we can start with the butler and displaying images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Create Butler with the rerun directory of the processed DECam Bulge data \n", + "butler = dafPersist.Butler(RERUN_DIR)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## We create a two pandas data frames and a list of cal. exposures for CCD1.\n", + "\n", + "We have go through the metadata and validate for each dataId that it exits. If the dataset exits we write the coordinates min., max., and center based on the source catalog into a pandas data frame. All CCDs for a certain pointing should have the same time, i.e. you point all CCDs at once to the field, therefore, we only need one data frame which associates the visit to the time of observation. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "coordinate_list = []\n", + "calexp_list = []\n", + "visit_date_list = []\n", + "metadata = butler.queryMetadata('calexp',['visit','ccd','filter'])\n", + "# Iterate to metadata\n", + "for dataset in metadata:\n", + " dataId = {'visit': int(dataset[0]), 'ccd': int(dataset[1]), 'filter':dataset[2]}\n", + " #Check if the data set has a source catalog\n", + " if butler.datasetExists('src', dataId=dataId):\n", + " srcCatalog = butler.get('src', dataId=dataId).asAstropy() # get the source catalog\n", + " # get the minimum and maximum RA/DEC\n", + " raMax = srcCatalog['coord_ra'].max()\n", + " raMin = srcCatalog['coord_ra'].min()\n", + " decMax = srcCatalog['coord_dec'].max()\n", + " decMin = srcCatalog['coord_dec'].min()\n", + " #Calculate the center \n", + " raCenter = 0.5*(raMax + raMin)\n", + " decCenter = 0.5*(decMax + decMin)\n", + " # Number of detected sources\n", + " nr_detected_sources = len(srcCatalog)\n", + " # Get the median effective PSF area\n", + " medianPSFarea = np.median(srcCatalog['base_PsfFlux_area'])\n", + " #Append to the list.\n", + " coordinate_list.append([int(dataset[0]), int(dataset[1]), dataset[2], raCenter, decCenter, \n", + " raMin, raMax, decMin, decMax, nr_detected_sources, medianPSFarea])\n", + " # for ccd1, create a list with visit and time of observation. (Should be the same for all other CCDs)\n", + " if int(dataset[1]) == 1:\n", + " # get the calexp for CCD1\n", + " calexp = butler.get('calexp', visit=int(dataset[0]), ccd=1)\n", + " # Append the calexp in the list\n", + " calexp_list.append( calexp )\n", + " #Get visit info, to extract time of observation\n", + " exp_visit_info = calexp.getInfo().getVisitInfo()\n", + " visit_date = exp_visit_info.getDate()\n", + " visit_date_list.append( [int(dataset[0]), visit_date.toPython()] )\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Write the lists to a panda data frame \n", + "df_valid_visists = pd.DataFrame(coordinate_list, columns=['visit', 'ccd', 'DECAM_filter', 'ra_center', 'dec_center', \n", + " 'ra_min', 'ra_max', 'dec_min', 'dec_max', 'nr_detected_sources', 'median_effPSF_area'])\n", + "df_visit_date = pd.DataFrame(visit_date_list, columns=['visit', 'timestamp'])\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Save the pandas data frame as parque on disk\n", + "df_valid_visists.to_parquet( os.path.join(parquet_save_path,'df_valid_visits.parquet.gzip'), compression='gzip')\n", + "df_visit_date.to_parquet( os.path.join(parquet_save_path,'df_visits_date.parquet.gzip'), compression='gzip')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Read the paque files\n", + "df_valid_visists = pd.read_parquet( os.path.join(parquet_save_path,'df_valid_visits.parquet.gzip'), engine='fastparquet' )\n", + "df_visit_data = pd.read_parquet( os.path.join(parquet_save_path,'df_visits_date.parquet.gzip'), engine='fastparquet' )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Show the visit time of observations\n", + "df_visit_data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Show the beginning of the data frame which has all valid visits and coordinates in it.\n", + "df_valid_visists.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Get unique visit IDs\n", + "df_valid_visists.visit.unique()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Query all ccd 1 in the pandas data frame\n", + "valid_visit_ccd1 = df_valid_visists.query('ccd == 1')\n", + "valid_visit_ccd1.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that all pointing centers for CCD 1 are in the same field." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Print the standard deviation of the pointing centers for each visit and for CCD1 in arcsec\n", + "print(valid_visit_ccd1['ra_center'].std()*u.deg.to(u.arcsec),valid_visit_ccd1['dec_center'].std()*u.deg.to(u.arcsec))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#plot the first ten images of CCD1 to inspect visually the pointing:\n", + "\n", + "for ii,calexp in enumerate(calexp_list[:10]):\n", + "\n", + " plt.figure(ii)\n", + " display = afwDisplay.Display(frame=ii, backend='matplotlib')\n", + " display.scale(\"linear\", \"zscale\")\n", + " #display only a small region of the calexp.\n", + " display.mtv(calexp[500:1500,2500:3000])\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Sort valid visit from ccd1 \n", + "valid_visit_ccd1 = valid_visit_ccd1.sort_values(by=['nr_detected_sources'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "valid_visit_ccd1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "unique_visits = valid_visit_ccd1.visit.unique()\n", + "ccd = 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataId={'visit': int(unique_visits[0]), 'ccd': ccd}\n", + "srcCatalog1 = butler.get('src', dataId=dataId).asAstropy().to_pandas()\n", + "dataId={'visit': int(unique_visits[1]), 'ccd': ccd}\n", + "srcCatalog2 = butler.get('src', dataId=dataId).asAstropy().to_pandas()\n", + "\n", + "srcCatalog1 = srcCatalog1.sort_values(by=['id'])\n", + "srcCatalog2 = srcCatalog2.sort_values(by=['id'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "srcCatalog1.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "srcCatalog2.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "all_dist = np.array([])\n", + "#create master star list\n", + "\n", + "\n", + "dataId={'visit': int(unique_visits[0]), 'ccd': ccd}\n", + "\n", + "srcCatalog = butler.get('src', dataId=dataId) # get the source catalog for the first visit and make this the master star list\n", + "master_coordinates = SkyCoord(srcCatalog['coord_ra']*u.deg, srcCatalog['coord_dec']*u.deg)\n", + "master_starID = srcCatalog['id']\n", + "\n", + "coord_table = Table([master_starID, master_coordinates], names=('id', int(unique_visits[0])))\n", + "\n", + "for visit in unique_visits[1:]:\n", + " dataId={'visit': int(visit), 'ccd': ccd}\n", + " srcCatalog = butler.get('src', dataId=dataId)\n", + " coordinates = SkyCoord(srcCatalog['coord_ra']*u.deg, srcCatalog['coord_dec']*u.deg)\n", + " idx, d2d, d3d = master_coordinates.match_to_catalog_sky(coordinates)\n", + " coord_table.add_column(coordinates[idx], name=int(visit))\n", + " coord_table[f'distance {int(visit)}'] = d2d.arcsec*u.arcsec\n", + " all_dist = np.append(all_dist,d2d.arcsec)\n", + " print('visit: {} max. dist. {:.3f} arcsec std. dist. {:.3f} arcsec'.format(visit, np.max(d2d.arcsec), np.std(d2d.arcsec)))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot histogram of distances\n", + "fig, ax = plt.subplots(1, 1, figsize=(10, 4))\n", + "hist(all_dist, bins='freedman', ax=ax, histtype='stepfilled',\n", + " alpha=0.75, density=True)\n", + "ax.set_xlabel('distance [arcsec]')\n", + "ax.set_ylabel('Density(distance)')\n", + "ax.set_xlim(-0.001,0.01)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "TODO:\n", + "\n", + "- Check astrometry\n", + "- Make light curves\n", + "- create co-add image\n", + "- run ap-pipe\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.7.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}