diff --git a/DP1/100_How_to_Use_RSP_Tools/105_Image_reprocessing/105_3_Forced_Photometry.ipynb b/DP1/100_How_to_Use_RSP_Tools/105_Image_reprocessing/105_3_Forced_Photometry.ipynb index 0e6b7543..3f74e6ae 100644 --- a/DP1/100_How_to_Use_RSP_Tools/105_Image_reprocessing/105_3_Forced_Photometry.ipynb +++ b/DP1/100_How_to_Use_RSP_Tools/105_Image_reprocessing/105_3_Forced_Photometry.ipynb @@ -10,7 +10,7 @@ "id": "325aa8a5-92fd-4913-a471-ad1617343be6", "metadata": {}, "source": [ - "# 105.3. Forced Photometry\n", + "# 105.3. Forced photometry\n", "\n", "
\n", "\n", @@ -22,8 +22,9 @@ "Data Release: Data Preview 1
\n", "Container Size: large
\n", "LSST Science Pipelines version: r29.2.0
\n", - "Last verified to run: 2025-11-03
\n", - "Repository: github.com/lsst/tutorial-notebooks
" + "Last verified to run: 2025-12-15
\n", + "Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)\\\n", + "DOI: [10.11578/rubin/dc.20250909.20](https://doi.org/10.11578/rubin/dc.20250909.20)" ] }, { @@ -31,7 +32,7 @@ "id": "ae760c59-bf94-45b3-9f8e-b4551b22e193", "metadata": {}, "source": [ - "**Learning objective:** Measure sources from an image via forced photometry using IDs and coordinates.\n", + "**Learning objective:** Make forced measurements on an image for a set of coordinates.\n", "\n", "**LSST data products:** `deep_coadd`, `object` table\n", "\n", @@ -42,9 +43,7 @@ "Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.\n", "\n", "**Get Support:**\n", - "Everyone is encouraged to ask questions or raise issues in the \n", - "Support Category \n", - "of the Rubin Community Forum.\n", + "Everyone is encouraged to ask questions or raise issues in the [Support Category](https://community.lsst.org/c/support) of the Rubin Community Forum.\n", "Rubin staff will respond to all questions posted there." ] }, @@ -56,9 +55,16 @@ "source": [ "## 1. Introduction\n", "\n", - "This notebook performs forced photometry using the Data Preview 1 (DP1) data products on an input image, given a set of input IDs and RA/Dec coordinates. The `ForcedMeasurementDriverTask` serves as a convenience function, allowing users to measure sources without running a full-scale `pipetask` command or dealing with additional pipeline setup. This task can measure fluxes, shapes, and other properties at known source positions (e.g., from a pre-existing catalog) on a given image, without running source detection. (Note that for brevity we disable shape measurement in this notebook.)\n", + "In general, \"forced\" photometry refers to a measurement made at a fixed coordinate in an image, regardless of whether an above-threshold (SNR $>5$) source was detected at that location in that particular image.\n", + "Forced photometry uses the point-spread-function (PSF) of the image, and measures the flux assuming the shape and size of the PSF at the provided coordinates.\n", "\n", - "**Related tutorials:** 100-level Detect and measure sources" + "This notebook performs forced photometry using the Data Preview 1 (DP1) data products on an input image, given a set of input coordinates (RA/Dec).\n", + "\n", + "The `ForcedMeasurementDriverTask` serves as a convenience function, allowing users to measure sources without running a full-scale `pipetask` command or dealing with additional pipeline setup.\n", + "This task can measure fluxes, shapes, and other properties at known source positions (e.g., from a pre-existing catalog) on a given image, without running source detection.\n", + "For brevity, shape measurement is disabled in this tutorial.\n", + "\n", + "**Related tutorials:** The 100-level tutorial on how to detect and measure sources. The 200-level tutorials on the PSF." ] }, { @@ -74,8 +80,9 @@ "id": "ccb4fe61-835a-47ac-a691-e086ce766e79", "metadata": {}, "source": [ - "Import `numpy`, a fundamental package for scientific computing with arrays in Python\n", - "(numpy.org). From the `lsst` package, we import `Butler` to access the LSST data repository, which allows us to retrieve images and catalogs using data IDs. Also import `ForcedMeasurementDriverConfig` and `ForcedMeasurementDriverTask`, which are mid-level config and task classes used to set up and execute forced photometry." + "Import standard python packages `numpy`, `matplotlib`, and `astropy`.\n", + "\n", + "From the `lsst` package, import the TAP service (`get_tap_service`) and the `Butler` for data access, `afw_display` to view images, and `ForcedMeasurementDriverConfig` and `ForcedMeasurementDriverTask`: mid-level configuration and task classes used to set up and execute forced photometry." ] }, { @@ -88,14 +95,14 @@ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import astropy.units as u\n", + "from astropy.io import ascii\n", + "\n", + "from lsst.rsp import get_tap_service\n", "from lsst.daf.butler import Butler\n", "from lsst.pipe.tasks.measurementDriver import (\n", " ForcedMeasurementDriverConfig,\n", " ForcedMeasurementDriverTask,\n", - ")\n", - "from lsst.rsp import get_tap_service\n", - "from lsst.utils.plotting import set_rubin_plotstyle\n", - "set_rubin_plotstyle()" + ")" ] }, { @@ -106,40 +113,6 @@ "### 1.2. Define parameters and functions" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "8e2ccc6d-c56f-4481-8e4b-83370d0cb2f1", - "metadata": {}, - "outputs": [], - "source": [ - "def fix_flux_columns(table):\n", - " \"\"\"Rename the columns in an output table from the ForcedMeasurementDriverTask.\n", - "\n", - " This task is necessary to rename \"instFlux\" columns to \"flux\" and units from\n", - " \"ct\" to \"nJy\".\n", - "\n", - " Parameters\n", - " ----------\n", - " table : Astropy Table\n", - " An Astropy table resulting from executing the ForcedMeasurementDriver\n", - "\n", - " Outputs\n", - " -------\n", - " table : Astropy Table\n", - " Transformed version of the input table, with columns and units updated.\n", - " \"\"\"\n", - " \n", - " cols_with_instflux = [col for col in np.array(table.colnames) if 'instFlux' in col]\n", - " \n", - " for col in cols_with_instflux:\n", - " new_colname = str.replace(col, 'instFlux', 'flux')\n", - " table.rename_column(col, new_colname)\n", - " table[new_colname].unit = u.nJy\n", - "\n", - " return table" - ] - }, { "cell_type": "markdown", "id": "d4edf616-4ed4-4bbb-94c8-5b1fefaafd24", @@ -178,14 +151,63 @@ "assert service is not None" ] }, + { + "cell_type": "markdown", + "id": "a61d47d5-32d0-4005-a6f2-bc93d17cf489", + "metadata": {}, + "source": [ + "Create a function called `fix_flux_columns` to rename output columns from the forced measurement task that are incorrectly labeled as instrument fluxes, when they are in fact not instrumental (i.e., in units of counts) but are calibrated fluxes in units of nJy.\n", + "This mistake will be fixed in future versions of the LSST Science Pipelines." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e2ccc6d-c56f-4481-8e4b-83370d0cb2f1", + "metadata": {}, + "outputs": [], + "source": [ + "def fix_flux_columns(table):\n", + " \"\"\"Rename the columns in an output table from the ForcedMeasurementDriverTask.\n", + "\n", + " This task is necessary to rename \"instFlux\" columns to \"flux\" and units from\n", + " \"ct\" to \"nJy\".\n", + "\n", + " Parameters\n", + " ----------\n", + " table : Astropy Table\n", + " An Astropy table resulting from executing the ForcedMeasurementDriver\n", + "\n", + " Outputs\n", + " -------\n", + " table : Astropy Table\n", + " Transformed version of the input table, with columns and units updated.\n", + " \"\"\"\n", + "\n", + " cols_with_instflux = [col for col in np.array(table.colnames) if 'instFlux' in col]\n", + "\n", + " for col in cols_with_instflux:\n", + " new_colname = str.replace(col, 'instFlux', 'flux')\n", + " table.rename_column(col, new_colname)\n", + " table[new_colname].unit = u.nJy\n", + "\n", + " return table" + ] + }, { "cell_type": "markdown", "id": "16eeb034-6728-406e-bfc8-7b875c3a9416", "metadata": {}, "source": [ - "## 2. Load an image and table\n", + "## 2. Measure objects in a deep coadd\n", + "\n", + "**Note:** The fluxes in the `Object` catalog are already forced photometry measurements, and re-measuring them is not necessary for scientific analysis. The `Object` table is just used to get a set of coordinates for this tutorial.\n", + "\n", + "### 2.1. Load an image and catalog\n", + "\n", + "Identify coadded images overlapping a particular sky position and select one to run forced photometry on. Also load the corresponding `Object` table using the TAP service so that forced photometry can be run at the positions of already-detected objects.\n", "\n", - "Identify coadded images overlapping a particular sky position and select one to run forced photometry on. Also load the corresponding `object` table using the TAP service so that forced photometry can be run at the positions of already-detected objects." + "Define the right ascension and declination as the center of the ECDFS field, and use the $r$ band." ] }, { @@ -200,6 +222,14 @@ "band = 'r'" ] }, + { + "cell_type": "markdown", + "id": "35debe43-cd77-431b-a294-9e50a037738b", + "metadata": {}, + "source": [ + "Execute a butler query for deep coadd images." + ] + }, { "cell_type": "code", "execution_count": null, @@ -209,20 +239,29 @@ "source": [ "query = f\"band='{band}' AND \\\n", " patch.region OVERLAPS POINT({ra}, {dec})\"\n", - "\n", - "coadd_img_refs = butler.query_datasets('deep_coadd', where=query)\n", - "\n", - "ref = coadd_img_refs[0]" + "coadd_img_refs = butler.query_datasets('deep_coadd', where=query,\n", + " order_by='patch')" + ] + }, + { + "cell_type": "markdown", + "id": "cc662eff-84ed-4fa9-9baa-7901cda08038", + "metadata": {}, + "source": [ + "Assert that at least one matching image was returned.\n", + "Retrieve the first in the sorted list." ] }, { "cell_type": "code", "execution_count": null, - "id": "b8e8b42d-eda2-4c96-a412-22e99da51a5d", + "id": "7b1aa5b0-b59e-49c2-8891-f7d0d469851b", "metadata": {}, "outputs": [], "source": [ - "exposure = butler.get(ref)" + "assert len(coadd_img_refs) > 0\n", + "ref = coadd_img_refs[0]\n", + "deep_coadd = butler.get(ref)" ] }, { @@ -230,7 +269,7 @@ "id": "3a57d04c-00ae-4278-b637-9e312084419d", "metadata": {}, "source": [ - "Extract a list of objects from the `Object` table in a small radius near the position of interest. These will be used both as an input list for forced measurement, and as a comparison for the results." + "Query for point-like objects from the `Object` table in a 0.3 degree radius near the position of interest that overlap the deep coadd image's tract and patch. " ] }, { @@ -242,8 +281,7 @@ "source": [ "query = f\"\"\"SELECT TOP 50 objectId, coord_ra, coord_dec, {band}_psfMag, {band}_psfFlux\n", " FROM dp1.Object\n", - " WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),\n", - " CIRCLE('ICRS', {ra}, {dec}, 0.1)) = 1\n", + " WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), CIRCLE('ICRS', {ra}, {dec}, 0.3)) = 1\n", " AND refExtendedness < 0.5\n", " AND {band}_psfMag < 26 AND {band}_psfMag > 17\n", " AND tract = {ref.dataId['tract']}\n", @@ -254,17 +292,53 @@ "job.wait(phases=['COMPLETED', 'ERROR'])\n", "print('Job phase is', job.phase)\n", "if job.phase == 'ERROR':\n", - " job.raise_if_error()\n", + " job.raise_if_error()" + ] + }, + { + "cell_type": "markdown", + "id": "f0939bcb-a517-4ec3-b2c0-3fb3b8700681", + "metadata": {}, + "source": [ + "Retrieve the search results as `objtable`.\n", + "This table will be used both as an input list for forced measurement, and as a comparison for the results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1803ae74-76d1-42f1-9e29-39231188f026", + "metadata": {}, + "outputs": [], + "source": [ "assert job.phase == 'COMPLETED'\n", "objtable = job.fetch_result().to_table()" ] }, + { + "cell_type": "markdown", + "id": "930f2a06-8918-43f5-9a05-bca95da38b03", + "metadata": {}, + "source": [ + "Assert that the `objtable` is not empty." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f3234934-8c70-4a75-87b8-2b59e5c1eb1a", + "metadata": {}, + "outputs": [], + "source": [ + "assert len(objtable) > 0" + ] + }, { "cell_type": "markdown", "id": "e4243773-e324-4415-884c-857cdb497688", "metadata": {}, "source": [ - "Option to examine the `objtable`:" + "Option to display the `objtable`." ] }, { @@ -282,7 +356,7 @@ "id": "0abd209b-d6b2-4ec5-95d5-7171bb1465fc", "metadata": {}, "source": [ - "## 3. Run the forced measurement driver" + "### 2.2. Run the forced measurement driver" ] }, { @@ -321,7 +395,7 @@ "id": "df670f7b-6b6d-47f4-af2c-c9f419a5699d", "metadata": {}, "source": [ - "Option: examine the configuration." + "Option to examine the configuration." ] }, { @@ -357,7 +431,7 @@ "id": "5bd50b2d-613c-4d23-9a72-177f18faa6f0", "metadata": {}, "source": [ - "Run the task using the input table of source positions and IDs, the `deep_coadd` exposure image, and required parameters: column names for ID, RA, and Dec, plus a PSF footprint scaling factor (used to create synthetic footprints since detection is skipped in forced photometry)." + "Run the task using the input table of source positions and IDs, the `deep_coadd` image, and required parameters: column names for ID, RA, and Dec, plus a PSF footprint scaling factor (used to create synthetic footprints since detection is skipped in forced photometry)." ] }, { @@ -369,7 +443,7 @@ "source": [ "result = driver.runFromAstropy(\n", " objtable,\n", - " exposure,\n", + " deep_coadd,\n", " id_column_name=\"objectId\",\n", " ra_column_name=\"coord_ra\",\n", " dec_column_name=\"coord_dec\",\n", @@ -382,7 +456,31 @@ "id": "0c3176bd-4642-4b71-9e18-fd29a8f00125", "metadata": {}, "source": [ - "Because the `ForcedMeasurementDriverTask` relies on code that is typically used in early stages of data processing, it assumes that fluxes are instrumental fluxes and appends the suffixes `_instFlux` and `_instFluxErr` to the output columns. It also sets the units for these columns to \"ct\" (counts). Because we _know_ that we are measuring on calibrated images whose pixel values are in nJy (to confirm this, execute `exposure.metadata['BUNIT']` in a code cell to see what the units are), we rename the columns and set their units to nJy. The code that causes this behavior will be updated in the future to use keywords from the input image to set the units and suffixes." + "Because the `ForcedMeasurementDriverTask` relies on code that is typically used in early stages of data processing, it assumes that fluxes are instrumental fluxes and appends the suffixes `_instFlux` and `_instFluxErr` to the output columns.\n", + "It also sets the units for these columns to \"ct\" (counts).\n", + "\n", + "However, these assumptions about the types of fluxes are incorrect when running this task on calibrated visit and deep coadd images, whose pixel values are in nJy.\n", + "\n", + "Confirm that the pixel values are in units of nJy." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d52f2ecb-3212-4677-a3b3-656133244fcb", + "metadata": {}, + "outputs": [], + "source": [ + "deep_coadd.metadata['BUNIT']" + ] + }, + { + "cell_type": "markdown", + "id": "03781d81-4975-468b-9c1d-b4e6956053eb", + "metadata": {}, + "source": [ + "Use the function defined in Section 1.2 to fix the column names.\n", + "This step is temporary as the column names will be fixed in the future." ] }, { @@ -424,7 +522,10 @@ "\n", "Each flux measurement has an associated flux error and one or more flag columns to allow filtering of bad measurements. Additionally, there are many pixel-based flags (denoted `base_PixelFlags_*`) highlighting possible issues with the data for each source.\n", "\n", - "Columns that begin with `slot_` are duplicates of other columns in the table. These exist as a way of telling algorithms in the Science Pipelines which is the \"preferred\" measurement for things like positions, shapes, and fluxes. Because there are multiple ways of measuring shapes (for example), this lets processing algorithms know which of these to select. Note that the specific measurements populating (some of) the slots were configured above, for example when we specified `config.measurement.slots.psfFlux = \"base_PsfFlux\"`. " + "Columns that begin with `slot_` are duplicates of other columns in the table.\n", + "These exist as a way of telling algorithms in the Science Pipelines which is the \"preferred\" measurement for things like positions, shapes, and fluxes.\n", + "Because there are multiple ways of measuring shapes (for example), this lets processing algorithms know which of these to select.\n", + "Note that the specific measurements populating (some of) the slots were configured above, for example `config.measurement.slots.psfFlux = \"base_PsfFlux\"` was specified. " ] }, { @@ -434,7 +535,7 @@ "source": [ "The forced measurement functionality is limited to a handful of measurement plugins (for example, model-fitting algorithms will not work with this task because they require more ancillary information than it is set up to receive). Nonetheless, it offers a quick and easy way to extract measurements at arbitrary positions.\n", "\n", - "### 3.1 Compare the forced measurements to the Object table\n", + "### 2.3. Compare forced measurements to the Object table\n", "\n", "Make a plot showing the ratio of the measured PSF fluxes from the forced measurement task to the `psfFlux` for the same objects in the `Object` table. Plot these as a function of magnitude, using Astropy units to convert `slot_PsfFlux_flux` from the `result` table to AB magnitudes." ] @@ -449,51 +550,78 @@ "meas_mag = (result['slot_PsfFlux_flux'].value*u.nJy).to(u.ABmag)\n", "\n", "fig, ax = plt.subplots(figsize=(6, 4))\n", - "\n", - "ax.plot(meas_mag, result['slot_PsfFlux_flux']/objtable['r_psfFlux'], 'k.')\n", + "ax.plot(meas_mag, result['slot_PsfFlux_flux']/objtable['r_psfFlux'],\n", + " 'o', alpha=0.5, mew=1, mec='black', color='black')\n", "ax.hlines(1.0, 17.3, 26.3, linestyle=':', color='gray')\n", - "ax.minorticks_on()\n", "ax.set_xlim(17.3, 26.3)\n", - "ax.set_xlabel('measured PSF mag')\n", - "ax.set_ylabel('flux ratio (forced/Object)')\n", - "\n", - "plt.tight_layout()\n", + "ax.set_xlabel('force-measured PSF magnitude')\n", + "ax.set_ylabel('flux ratio (forced / Object)')\n", "plt.show()" ] }, { "cell_type": "markdown", - "id": "159d91c6-d600-4fa2-8dea-1a07e733059b", + "id": "4cd8bf7d-53d1-43a9-a299-7339c196242c", "metadata": {}, "source": [ - "## 4. Forced measurement on a visit image\n", + "> **Figure 1:** The flux ratio of the force-measured PSF flux divided by the PSF flux from the `Object` table, vs. the force-measured PSF magnitude (all for $r$-band).\n", "\n", - "The forced measurement task can be run on any `ExposureF` object. To demonstrate this, retrieve a `visit_image` and execute the driver to extract forced photometry on that image. The appropriate input catalog for `visit_image` measurements is the `Source` table." + "The plot above demonstrates that the **forced fluxes can sometimes be higher (brighter) than fluxes in the `Object` table**.\n", + "This is usually **due to contaminating flux** from nearby objects, as deblending is not performed as part of forced measurement, but it is performed when the fluxes in the `Object` table are derived." + ] + }, + { + "cell_type": "markdown", + "id": "75f941f2-74cf-44b8-8c98-5274714cc0b8", + "metadata": {}, + "source": [ + "Delete the image and catalogs from memory." ] }, { "cell_type": "code", "execution_count": null, - "id": "9f682830-1497-4909-99f5-43c280fd88e8", + "id": "0a774cc1-867d-4271-aeff-f6d4f9306c64", "metadata": {}, "outputs": [], "source": [ - "query = f\"band='{band}' AND \\\n", - " visit.region OVERLAPS POINT({ra}, {dec})\"\n", + "del deep_coadd, objtable, result" + ] + }, + { + "cell_type": "markdown", + "id": "159d91c6-d600-4fa2-8dea-1a07e733059b", + "metadata": {}, + "source": [ + "## 3. Measure sources in a visit image\n", "\n", - "visit_img_refs = butler.query_datasets('visit_image', where=query)\n", + "The above process of obtaining forced photometry for objects in a deep coadd image can be applied to sources in a visit image, using the same configuration parameters.\n", "\n", - "ref = visit_img_refs[0]" + "Retrieve the first $r$-band visit image that overlaps the same coordinates as defined above." ] }, { "cell_type": "code", "execution_count": null, - "id": "171e7626-85a7-4974-b1e6-d4964f8b814b", + "id": "9f682830-1497-4909-99f5-43c280fd88e8", "metadata": {}, "outputs": [], "source": [ - "exposure = butler.get(ref)" + "query = f\"band='{band}' AND \\\n", + " visit.region OVERLAPS POINT({ra}, {dec})\"\n", + "visit_img_refs = butler.query_datasets('visit_image', where=query, limit=1,\n", + " order_by='visit.timespan.begin')\n", + "assert len(visit_img_refs) == 1\n", + "ref = visit_img_refs[0]\n", + "visit_image = butler.get(ref)" + ] + }, + { + "cell_type": "markdown", + "id": "c0760743-4c48-4db1-95fb-b2d2dd05da30", + "metadata": {}, + "source": [ + "Query the TAP service for point-like sources detected in the visit image." ] }, { @@ -505,8 +633,7 @@ "source": [ "query = f\"\"\"SELECT TOP 50 sourceId, ra, dec, psfFlux\n", " FROM dp1.Source\n", - " WHERE CONTAINS(POINT('ICRS', ra, dec),\n", - " CIRCLE('ICRS', {ra}, {dec}, 0.1)) = 1\n", + " WHERE CONTAINS(POINT('ICRS', ra, dec), CIRCLE('ICRS', {ra}, {dec}, 0.3)) = 1\n", " AND extendedness < 0.5\n", " AND visit = {ref.dataId['visit']}\n", " AND detector = {ref.dataId['detector']}\"\"\"\n", @@ -517,8 +644,20 @@ "print('Job phase is', job.phase)\n", "if job.phase == 'ERROR':\n", " job.raise_if_error()\n", + "\n", "assert job.phase == 'COMPLETED'\n", - "srctable = job.fetch_result().to_table()" + "srctable = job.fetch_result().to_table()\n", + "assert len(srctable) > 0" + ] + }, + { + "cell_type": "markdown", + "id": "acc39c47-eb40-46a9-b130-c2873d14a41f", + "metadata": {}, + "source": [ + "Configure and run forced photometry.\n", + "\n", + "The only difference from Section 2 is that when the `driver` is called, `srctable` and `visit_image` are passed, and the `id_column_name`, `ra_column_name`, and `dec_column_name` are changed, as appropriate for the `Source` table column names: `sourceId`, `ra`, and `dec`." ] }, { @@ -533,7 +672,28 @@ "config.measurement.slots.centroid = \"base_TransformedCentroidFromCoord\"\n", "config.measurement.slots.shape = None\n", "config.measurement.doReplaceWithNoise = False\n", - "config.doApCorr = True" + "config.doApCorr = True\n", + "\n", + "driver = ForcedMeasurementDriverTask(config=config)\n", + "\n", + "src_result = driver.runFromAstropy(\n", + " srctable,\n", + " visit_image,\n", + " id_column_name=\"sourceId\",\n", + " ra_column_name=\"ra\",\n", + " dec_column_name=\"dec\",\n", + " psf_footprint_scaling=3.0,\n", + ")\n", + "\n", + "src_result = fix_flux_columns(src_result)" + ] + }, + { + "cell_type": "markdown", + "id": "1eb934e8-df87-4cc0-a096-8d5a29b42b17", + "metadata": {}, + "source": [ + "Option to view 5 rows of the results table." ] }, { @@ -543,91 +703,195 @@ "metadata": {}, "outputs": [], "source": [ - "driver = ForcedMeasurementDriverTask(config=config)" + "# src_result[:5]" ] }, { "cell_type": "code", "execution_count": null, - "id": "ca15977d-db4b-4f06-8f49-92678f92909d", + "id": "ad5ea790-39e1-4350-bdb9-361a5392ad10", "metadata": {}, "outputs": [], "source": [ - "src_result = driver.runFromAstropy(\n", - " srctable,\n", - " exposure,\n", - " id_column_name=\"sourceId\",\n", - " ra_column_name=\"ra\",\n", - " dec_column_name=\"dec\",\n", - " psf_footprint_scaling=3.0,\n", - ")" + "del visit_image, srctable, src_result" + ] + }, + { + "cell_type": "markdown", + "id": "03c59928-4b16-40b3-a656-3f840e74184d", + "metadata": {}, + "source": [ + "## 4. Measurements for user-specified coordinates\n", + "\n", + "Often, the motivation for doing forced photometry is to make forced measurements at a set of coordinates where no object or source has been detected with SNR $>5$ in the deep coadd or visit image (i.e., at coordinates that are not in the `Object` or `Source` tables).\n", + "\n", + "These coordinates might be for subthreshold detections (e.g., SNR $<5$) or for detections from other (non-Rubin) surveys, e.g., infrared, space-based.\n", + "\n", + "For this demonstration, a prepared input file of coordinates is loaded, and then forced photometry is run on a visit image.\n", + "\n", + "The input file was generated by re-running source detection with `config.thresholdValue = 4` on a specific $r$-band visit image (`visit` = 2024110800263, `detector` = 3), and saving the coordinates for sources with $4.0 < SNR < 4.5$ in a small region of the visit image to the file.\n", + "The sub-threshold detection magnitudes were also saved in the file.\n", + "\n", + "Load the file as an astropy table and print the column names." ] }, { "cell_type": "code", "execution_count": null, - "id": "2bdf49bc-b5ff-4d0a-936f-ccd77017aa7c", + "id": "476ad58a-851b-4de7-bd15-e600b15842f2", "metadata": {}, "outputs": [], "source": [ - "src_result = fix_flux_columns(src_result)" + "filepath = '/rubin/cst_repos/tutorial-notebooks-data/data/'\n", + "filename = 'dp1_105_3_subthreshold.dat'\n", + "subthresh = ascii.read(filepath+filename)\n", + "print('file columns: ', subthresh.colnames)\n", + "print('file length: ', len(subthresh))" + ] + }, + { + "cell_type": "markdown", + "id": "3e54b0d9-d25c-436c-ae2d-94095c636f32", + "metadata": {}, + "source": [ + "Option to view the loaded table." ] }, { "cell_type": "code", "execution_count": null, - "id": "15678271-af3e-4bea-a147-b989c970aed8", + "id": "72206982-8f19-45cf-96b9-d5e9f5ed6fd7", "metadata": {}, "outputs": [], "source": [ - "src_result[:5]" + "# subthresh" ] }, { "cell_type": "markdown", - "id": "af970fba-e2b5-42db-95a5-b24e6ec11caa", + "id": "0e6a872a-3b02-4d7f-86ce-7366f049da69", "metadata": {}, "source": [ - "### 4.1 Compare the forced measurements to the Source table\n", - "\n", - "Make a plot showing the ratio of the measured PSF fluxes from the forced measurement task to the `psfFlux` for the same objects in the `Source` table. Plot these as a function of magnitude, using Astropy units to convert `slot_PsfFlux_flux` from the `result` table to AB magnitudes." + "Retrieve the visit image from the Butler." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52f7fbda-0fde-4e72-93de-d9707dbdb1e8", + "metadata": {}, + "outputs": [], + "source": [ + "visit_image = butler.get('visit_image', visit=2024110800263, detector=3)" + ] + }, + { + "cell_type": "markdown", + "id": "16e78644-48a0-4cef-ac81-cbcaded56381", + "metadata": {}, + "source": [ + "Run the forced measurement driver with the same configuration as used in Section 3, but pass the `subthresh` table and use column names that match the input file." ] }, { "cell_type": "code", "execution_count": null, - "id": "b561fa35-fc6f-4b39-b448-f57f302a3084", + "id": "bbe2c9fb-dd25-46c2-9be1-906c9d152d96", "metadata": {}, "outputs": [], "source": [ - "meas_mag = (src_result['slot_PsfFlux_flux'].value*u.nJy).to(u.ABmag)\n", + "config = ForcedMeasurementDriverConfig()\n", + "config.measurement.slots.psfFlux = \"base_PsfFlux\"\n", + "config.measurement.slots.centroid = \"base_TransformedCentroidFromCoord\"\n", + "config.measurement.slots.shape = None\n", + "config.measurement.doReplaceWithNoise = False\n", + "config.doApCorr = True\n", "\n", - "plt.plot(meas_mag, src_result['slot_PsfFlux_flux']/srctable['psfFlux'], 'k.')\n", - "plt.hlines(1.0, 17.3, 26.3, linestyle=':', color='gray')\n", - "plt.minorticks_on()\n", - "plt.xlim(17.3, 26.3)\n", - "plt.xlabel('measured PSF mag')\n", - "plt.ylabel('flux ratio (forced/sourceTable)')\n", + "driver = ForcedMeasurementDriverTask(config=config)\n", "\n", - "plt.tight_layout()\n", - "plt.show()" + "src_result = driver.runFromAstropy(\n", + " subthresh,\n", + " visit_image,\n", + " id_column_name=\"id\",\n", + " ra_column_name=\"ra\",\n", + " dec_column_name=\"dec\",\n", + " psf_footprint_scaling=3.0,\n", + ")\n", + "\n", + "src_result = fix_flux_columns(src_result)" + ] + }, + { + "cell_type": "markdown", + "id": "a6af247b-b7c9-4768-af39-4c6c49eead79", + "metadata": {}, + "source": [ + "Assert that the input and output tables have the same length." ] }, { "cell_type": "code", "execution_count": null, - "id": "90e0d3b6-4761-4466-922e-bcb6be1188c1", + "id": "26f12ca6-9624-44c8-973a-8da08d5e5cdf", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "assert len(subthresh) == len(src_result)" + ] + }, + { + "cell_type": "markdown", + "id": "323a4504-ed40-493f-9eba-fd7e3305e784", + "metadata": {}, + "source": [ + "Option to display the first five lines of the output table." + ] }, { "cell_type": "code", "execution_count": null, - "id": "52f7fbda-0fde-4e72-93de-d9707dbdb1e8", + "id": "22037910-72e2-42a7-8568-2398449db47d", + "metadata": {}, + "outputs": [], + "source": [ + "# src_result[:5]" + ] + }, + { + "cell_type": "markdown", + "id": "16b266a9-d47f-4693-aa03-a0b899b97181", + "metadata": {}, + "source": [ + "Convert fluxes to magnitudes (and vice versa) and plot the flux ratio between the forced and sub-threshold detection fluxes as a function of the forced photometry magnitudes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0eb910c6-5c6b-4e75-a82e-29f7b2df2290", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "src_result['mag'] = (src_result['slot_PsfFlux_flux'].value*u.nJy).to(u.ABmag)\n", + "subthresh['r_flux'] = (subthresh['r_mag'].value*u.ABmag).to(u.nJy)\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 4))\n", + "ax.plot(subthresh['r_mag'], src_result['slot_PsfFlux_flux']/subthresh['r_flux'],\n", + " 'o', alpha=0.5, mew=1, mec='black', color='black')\n", + "ax.hlines(1.0, 17.3, 26.3, linestyle=':', color='gray')\n", + "ax.set_xlim(24.3, 24.5)\n", + "ax.set_xlabel('input file magnitude')\n", + "ax.set_ylabel('flux ratio (forced / input)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "a05b2fbc-ace9-4ca1-a37a-f8ff1af67e17", + "metadata": {}, + "source": [ + "> **Figure 2:** The flux ratio of the force-measured PSF flux divided by the PSF flux from the original sub-threshold detection (from the input file) vs. the input file's sub-threshold detection magnitude." + ] } ], "metadata": {