diff --git a/DP2/300_Science_demos/301_DP2_observations/301_1_DP2_overview.ipynb b/DP2/300_Science_demos/301_DP2_observations/301_1_DP2_overview.ipynb index 0fa4bbad..3fbdddbf 100644 --- a/DP2/300_Science_demos/301_DP2_observations/301_1_DP2_overview.ipynb +++ b/DP2/300_Science_demos/301_DP2_observations/301_1_DP2_overview.ipynb @@ -21,8 +21,8 @@ "For the Rubin Science Platform at data.lsst.cloud. \\\n", "Data Release: [Data Preview 2](http://dp2.lsst.io) \\\n", "Container Size: Large \\\n", - "LSST Science Pipelines version: r30.0.10 \\\n", - "Last verified to run: 2026-07-24 \\\n", + "LSST Science Pipelines version: r30.0.11 \\\n", + "Last verified to run: 2026-08-17 \\\n", "Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)" ] }, @@ -53,7 +53,7 @@ "\n", "This notebook provides an overview of the observations that were obtained during commissioning and science validation with the LSST Camera ([lsstcam.lsst.io](https://lsstcam.lsst.io/)), processed with v30 of the Rubin Science Pipelines ([pipelines.lsst.io](https://pipelines.lsst.io/)), and released as Data Preview 2 (DP2; [dp2.lsst.io](https://dp2.lsst.io/)).\n", "\n", - "The code in this notebook was used to generated the [DP2 observations overview webpage](https://dp2.lsst.io/).\n", + "The code in this notebook was used to generate the [DP2 observations overview webpage](https://dp2.lsst.io/).\n", "\n", "**Key concepts for the DP2 dataset**\n", "\n", @@ -63,7 +63,7 @@ "**Skymap, tract, and patch:**\n", "One tract is one square region of LSST's all-sky tesselation, the \"skymap\".\n", "Tracts are $~1.66$ deg per side, identified by a unique four-digit number, and are subdivided into 100 overlapping patches.\n", - "One deep coadd image is created per patch, provided sufficient input visits exist.\n", + "One deep coadd image is created per patch in each band, provided sufficient input visits exist.\n", "The tract and patch numbers together with the filter uniquely identify a deep coadd image.\n", "\n", "**For DP2, the area covered by visits is greater than the area covered by deep coadd images.**\n", @@ -100,8 +100,10 @@ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.colors as mcolors\n", + "from astropy.coordinates import SkyCoord\n", "from astropy.table import Table, join\n", "from astropy.time import Time\n", + "import astropy.units as u\n", "import skyproj\n", "import healpy as hp\n", "import hpgeom as hpg\n", @@ -205,7 +207,7 @@ "id": "dd00d543-04cf-45a4-891e-e6f80313b59d", "metadata": {}, "source": [ - "Define dictionaries with the Data Preview 2 names and centers for Deep Drilling Fields (DDFs) and small survey fields (SFS). \n", + "Define dictionaries with the Data Preview 2 names and centers for Deep Drilling Fields (DDFs) and small field surveys (SFS). \n", "The coordinates are Right Ascension and Declination as International Celestial Reference System (ICRS) coordinates in decimal degrees.\n", "\n", "A radius in degrees is defined with the coordinates: it is *approximate*, and encompasses only the *boresight* (field-of-view center) coordinates for all individual visits of that field. It does not encompass every detector, or every skymap patch, or every detected source or object associated with that field." @@ -236,7 +238,7 @@ " \"Rubin_SV_225_-40\": [225.0, -39.5, 1.9],\n", " \"Rubin_SV_280_-48\": [280.1, -48.0, 2.0],\n", " \"Rubin_SV_300_-41\": [300.3, -41.0, 1.8],\n", - " \"Rubin_SV_320_-15\": [320.2, -15.1, 4.5],\n", + " \"Rubin_SV_320_-15\": [320.2, -15.1, 4.51],\n", "}\n", "region_names = list(regions.keys())" ] @@ -246,7 +248,7 @@ "id": "a500c7ca-134b-45e0-80ba-dec5e7397ccb", "metadata": {}, "source": [ - "Define a function that will remove a figure and all it's data, to help clear memory." + "Define a function that will remove a figure and all its data, to help clear memory." ] }, { @@ -376,7 +378,7 @@ "id": "ce49c560-766b-407c-b74e-a804b85a9d36", "metadata": {}, "source": [ - "Use the `skyproj` package with a McBryde skyprojection to visualize the distribution of DP2 visits on the sky in HEALPix that are approximately the same size as the LSSTCam FOV (`nside`=19)." + "Use the `skyproj` package with a McBryde sky projection to visualize the distribution of DP2 visits on the sky in HEALPix that are approximately the same size as the LSSTCam FOV (`nside`=19)." ] }, { @@ -387,12 +389,12 @@ "outputs": [], "source": [ "orig_cmap = plt.get_cmap('Blues')\n", - "new_cmap = truncate_colormap(orig_cmap, 0.10, 0.75)\n", + "new_cmap = truncate_colormap(orig_cmap, 0.20, 0.8)\n", "\n", "fig, ax = plt.subplots(figsize=(12, 8))\n", "sp = skyproj.McBrydeSkyproj(ax=ax)\n", - "vras = np.asarray(visit_table['ra'], dtype='float')\n", - "vdecs = np.asarray(visit_table['dec'], dtype='float')\n", + "vras = visit_table['ra'].value.data\n", + "vdecs = visit_table['dec'].value.data\n", "sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap=new_cmap)\n", "sp.draw_colorbar(label='Number of visits (any filter)', shrink=0.5, pad=0.01)\n", "for r, name in enumerate(region_names):\n", @@ -409,7 +411,7 @@ "sp.ax.legend(loc=\"lower center\", bbox_to_anchor=(0.5, 1.1), ncol=5, handletextpad=0)\n", "plt.tight_layout()\n", "plt.show()\n", - "del sp, vras, vdecs" + "del sp" ] }, { @@ -461,15 +463,16 @@ "metadata": {}, "outputs": [], "source": [ + "sc_all = SkyCoord(ra = vras*u.deg, dec = vdecs*u.deg, frame='icrs')\n", + "\n", "for i, name in enumerate(region_names):\n", " ra, dec, rad = regions[name]\n", - " cosdec = np.cos(np.deg2rad(dec))\n", - " offsets = np.sqrt((cosdec*(visit_table['ra']-ra))**2 + (visit_table['dec']-dec)**2)\n", - " tx = np.where(offsets < rad)[0]\n", + " sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_all).value <= rad)[0]\n", " for filt in filter_names:\n", " field_filt_dist[filt][i] = len(np.where(visit_table['band'][tx] == filt)[0])\n", " field_filt_dist['Total'][i] = len(tx)\n", - " del ra, dec, rad, offsets, tx\n", + " del ra, dec, rad, tx\n", "\n", "i = np.where(field_filt_dist['Field'] == \"All\")[0]\n", "for filt in filter_names:\n", @@ -520,14 +523,13 @@ "\n", "for i, name in enumerate(region_names):\n", " ra, dec, rad = regions[name]\n", - " cosdec = np.cos(np.deg2rad(dec))\n", - " offsets = np.sqrt((cosdec*(visit_table['ra']-ra))**2 + (visit_table['dec']-dec)**2)\n", - " tx = np.where(offsets < rad)[0]\n", + " sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_all).value <= rad)[0]\n", " field_epochs['visits'][i] = len(tx)\n", " unique_nights = np.unique(np.floor(visit_table['expMidptMJD'][tx]))\n", " field_epochs['nights'][i] = len(unique_nights)\n", " field_epochs['mean(visits/night)'][i] = int(np.round(len(tx)/len(unique_nights), 0))\n", - " del ra, dec, rad, offsets, tx, unique_nights\n", + " del ra, dec, rad, tx, unique_nights\n", "\n", "i = np.where(field_epochs['Field'] == \"All\")[0]\n", "field_epochs['visits'][i] = len(visit_table)\n", @@ -627,7 +629,7 @@ "id": "c863cf20-d032-4340-a135-42109ca22515", "metadata": {}, "source": [ - "> **Figure 3:** The cumulative distribution of visit airmass for all visits, for all filters combined. The 50% percentile for airmass is 1.2, and is marked with a dashed grey line." + "> **Figure 3:** The cumulative distribution of visit airmass for all visits, for all filters combined. The 50th percentile for airmass is 1.2, and is marked with a dashed grey line." ] }, { @@ -800,7 +802,7 @@ "source": [ "Create tables of the mean seeing and magnitude limit, averaged over all visits in a given filter, per region.\n", "\n", - "Because seeing and magnitude limit are filter-dependent, but an average over all filters is not calculated." + "Because seeing and magnitude limit are filter-dependent, an average over all filters is not calculated." ] }, { @@ -819,9 +821,8 @@ "\n", "for i, name in enumerate(region_names):\n", " ra, dec, rad = regions[name]\n", - " cosdec = np.cos(np.deg2rad(dec))\n", - " offsets = np.sqrt((cosdec*(joined_table['ra']-ra))**2 + (joined_table['dec']-dec)**2)\n", - " tx = np.where(offsets < rad)[0]\n", + " sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_all).value <= rad)[0]\n", " for filt in filter_names:\n", " fx = np.where(joined_table['band'][tx] == filt)[0]\n", " if len(fx) > 0:\n", @@ -831,7 +832,7 @@ " field_seeing[filt][i] = np.nan\n", " field_maglim[filt][i] = np.nan\n", " del fx\n", - " del ra, dec, rad, offsets, tx\n", + " del ra, dec, rad, tx\n", "\n", "i = np.where(field_seeing['Field'] == \"All\")[0]\n", "for filt in filter_names:\n", @@ -923,7 +924,7 @@ "id": "a3aace34-9b8c-4c80-a4a7-531885d9d05c", "metadata": {}, "source": [ - "Get the unique tract number for DP2, and print their total number." + "Get the unique tract numbers for DP2, and print their total number." ] }, { @@ -979,9 +980,9 @@ " vertex_list = info.getVertexList()\n", " ras = []\n", " decs = []\n", - " for vertice in vertex_list:\n", - " ras.append(vertice.getRa().asDegrees())\n", - " decs.append(vertice.getDec().asDegrees())\n", + " for vertex in vertex_list:\n", + " ras.append(vertex.getRa().asDegrees())\n", + " decs.append(vertex.getDec().asDegrees())\n", " sp.draw_polygon(ras, decs, edgecolor='darkgrey', alpha=1, linewidth=0.5, facecolor=None)\n", " del info, vertex_list, ras, decs\n", "for r, name in enumerate(region_names):\n", @@ -1015,7 +1016,7 @@ "source": [ "#### 3.1.2. Single-field map of patches\n", "\n", - "Draw one box per patch for the DDF ECDFS field (2 in Figure 6).\n", + "Draw one box per patch for the DDF ECDFS field (numbered 3 in Figure 6).\n", "\n", "Recall from Section 1.2. that the third component of `regions[]` is the radius within which the boresight center coordinates for all visits of the field are included.\n", "To capture all tracts that contain patches for which these visits might have contributed to a deep coadd image,\n", @@ -1034,10 +1035,9 @@ "name = \"DDF_ECDFS\"\n", "ra, dec, rad = regions[name]\n", "radius = rad + np.sqrt(2*1.66**2)\n", - "delta_ra = (ra - patches_table['s_ra'])*(np.cos(np.deg2rad(dec)))\n", - "delta_dec = dec - patches_table['s_dec']\n", - "offset = np.sqrt(delta_ra**2 + delta_dec**2)\n", - "tx = np.where(offset < radius)[0]\n", + "sc_patches = SkyCoord(ra = patches_table['s_ra'], dec = patches_table['s_dec'], frame='icrs')\n", + "sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + "tx = np.where(sc_region.separation(sc_patches).value <= radius)[0]\n", "tract_color_dict = {4848: '#1f77b4', 4849: '#ff7f0e', 5063: '#2ca02c'}\n", "\n", "fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n", @@ -1051,15 +1051,16 @@ " if patches_table['lsst_tract'][x] in tract_color_dict:\n", " color = tract_color_dict[patches_table['lsst_tract'][x]]\n", " ax[a].plot(ras, decs, color=color, linewidth=0.5, alpha=1)\n", - " ax[a].set_xlabel(\"Right Ascension\", fontsize=12)\n", - " ax[a].set_ylabel(\"Declination\", fontsize=12)\n", + " ax[a].set_xlabel(\"Right Ascension (deg)\", fontsize=12)\n", + " ax[a].set_ylabel(\"Declination (deg)\", fontsize=12)\n", " del s_region, temp, ras, decs\n", "ax[1].set_xlim([ra-0.5, ra+0.5])\n", "ax[1].set_ylim([dec-0.5, dec+0.5])\n", + "ax[0].invert_xaxis()\n", + "ax[1].invert_xaxis()\n", "plt.tight_layout()\n", "plt.show()\n", - "del name, ra, dec, rad, radius, tx\n", - "del delta_ra, delta_dec, offset" + "del name, ra, dec, rad, radius, tx, sc_patches, sc_region" ] }, { @@ -1067,7 +1068,7 @@ "id": "195299b2-42d1-4818-9115-53996e737f86", "metadata": {}, "source": [ - "> **Figure 7:** At left, all pre-defined patches for the EDCFS. At right, the zoom-in better shows that both tracts and individual patches within a tract overlap at their edges. The center three tracts are assigned colors just to help with visualizing the overlap." + "> **Figure 7:** At left, all pre-defined patches for the EDCFS. At right, the zoom-in better shows that both tracts and individual patches within a tract overlap at their edges. The center three tracts are assigned colors to help with visualizing the overlap." ] }, { @@ -1088,11 +1089,12 @@ "# name = \"DDF_ECDFS\"\n", "# ra, dec, rad = regions[name]\n", "# radius = rad + np.sqrt(2*1.66**2)\n", - "# tx = np.where(np.sqrt(((ra - patches_table['s_ra'])*(np.cos(np.deg2rad(dec))))**2 +\n", - "# (dec - patches_table['s_dec'])**2) < radius)[0]\n", + "# sc_patches = SkyCoord(ra = patches_table['s_ra'], dec = patches_table['s_dec'], frame='icrs')\n", + "# sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + "# tx = np.where(sc_region.separation(sc_patches).value <= rad)[0]\n", "# tract_list = np.unique(patches_table['lsst_tract'][tx])\n", "# print(name, tract_list)\n", - "# del name, ra, dec, rad, radius, tx, tract_list" + "# del name, ra, dec, rad, radius, tx, tract_list, sc_patches, sc_region" ] }, { @@ -1244,10 +1246,9 @@ "name = \"DDF_ECDFS\"\n", "coords = regions[name]\n", "radius = coords[2] + np.sqrt(2*1.66**2)\n", - "delta_ra = (coords[0] - patches_table['s_ra'])*(np.cos(np.deg2rad(coords[1])))\n", - "delta_dec = coords[1] - patches_table['s_dec']\n", - "offset = np.sqrt(delta_ra**2 + delta_dec**2)\n", - "tx = np.where(offset < radius)[0]\n", + "sc_patches = SkyCoord(ra = patches_table['s_ra'], dec = patches_table['s_dec'], frame='icrs')\n", + "sc_region = SkyCoord(ra = coords[0]*u.deg, dec = coords[1]*u.deg, frame='icrs')\n", + "tx = np.where(sc_region.separation(sc_patches).value <= radius)[0]\n", "pixels = np.unique(np.array(hpg.angle_to_pixel(default_nside_coverage,\n", " patches_table['s_ra'][tx],\n", " patches_table['s_dec'][tx])))\n", @@ -1403,6 +1404,14 @@ "source": [ "# field_coadd_size" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "df9921d9-905e-42d1-ad16-3cd92c34a80d", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/DP2/300_Science_demos/301_DP2_observations/301_2_Small_field_areas.ipynb b/DP2/300_Science_demos/301_DP2_observations/301_2_Small_field_areas.ipynb index 0476e6f9..0ac1be9b 100644 --- a/DP2/300_Science_demos/301_DP2_observations/301_2_Small_field_areas.ipynb +++ b/DP2/300_Science_demos/301_DP2_observations/301_2_Small_field_areas.ipynb @@ -21,8 +21,8 @@ "For the Rubin Science Platform at data.lsst.cloud. \\\n", "Data Release: [Data Preview 2](http://dp2.lsst.io) \\\n", "Container Size: Large \\\n", - "LSST Science Pipelines version: v30.0.10 \\\n", - "Last verified to run: 2026-07-24 \\\n", + "LSST Science Pipelines version: v30.0.11 \\\n", + "Last verified to run: 2026-08-17 \\\n", "Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)" ] }, @@ -303,7 +303,7 @@ "id": "a500c7ca-134b-45e0-80ba-dec5e7397ccb", "metadata": {}, "source": [ - "Define a function that will remove a figure and all it's data, to help clear memory." + "Define a function that will remove a figure and all its data, to help clear memory." ] }, { @@ -521,6 +521,7 @@ "metadata": {}, "outputs": [], "source": [ + "sc_all = SkyCoord(ra = visit_table['ra'], dec = visit_table['dec'], frame='icrs')\n", "data_start = np.floor(np.min(visit_table['expMidptMJD']))\n", "data_end = np.ceil(np.max(visit_table['expMidptMJD']))\n", "temp = []\n", @@ -528,8 +529,8 @@ "fig, ax = plt.subplots(len(regions), 2, figsize=(10, 20))\n", "for r, name in enumerate(region_names):\n", " ra, dec, rad = regions[name]\n", - " offsets = np.sqrt((visit_table['ra']-ra)**2 + (visit_table['dec']-dec)**2)\n", - " tx = np.where(offsets < rad)[0]\n", + " sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_all).value <= rad)[0]\n", " for f, filt in enumerate(filter_names):\n", " fx = np.where(visit_table['band'][tx] == filt)[0]\n", " if len(fx) > 0:\n", @@ -551,7 +552,7 @@ " tdiffs.append(comb[1] - comb[0])\n", " time_diffs = np.array(tdiffs)\n", " temp.append((len(tx), np.mean(time_diffs)))\n", - " del ra, dec, offsets, tx, dates, tdiffs, time_diffs\n", + " del ra, dec, sc_region, tx, dates, tdiffs, time_diffs\n", " ax[r, 0].set_xlabel('MJD - ' + str(data_start))\n", " ax[r, 0].set_ylabel('Cumulative visits')\n", " ax[r, 0].set_xlim([0, data_end - data_start])\n", @@ -678,12 +679,14 @@ "metadata": {}, "outputs": [], "source": [ + "sc_joined = SkyCoord(ra = joined_table['ra'], dec = joined_table['dec'], frame='icrs')\n", "temp = []\n", + "\n", "fig, ax = plt.subplots(len(regions), 2, figsize=(10, 20))\n", "for r, name in enumerate(region_names):\n", " ra, dec, rad = regions[name]\n", - " offsets = np.sqrt((joined_table['ra']-ra)**2 + (joined_table['dec']-dec)**2)\n", - " tx = np.where(offsets < rad)[0]\n", + " sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_joined).value <= rad)[0]\n", " for f, filt in enumerate(filter_names):\n", " fx = np.where(joined_table['band'][tx] == filt)[0]\n", " if len(fx) > 0:\n", @@ -699,7 +702,7 @@ " np.mean(joined_table['mean_maglim'][fx]),\n", " np.std(joined_table['mean_maglim'][fx])))\n", " del fx\n", - " del ra, dec, offsets, tx\n", + " del ra, dec, sc_region, tx\n", " ax[r, 0].set_xlabel('Mean seeing FWHM [arcsec]')\n", " ax[r, 0].set_ylabel('Number of visits')\n", " ax[r, 0].set_title(name)\n", @@ -955,7 +958,9 @@ "\n", "Use a relatively small search radius of just 0.2 degrees for faster queries.\n", "\n", - "Plot the magnitude distributions for PSF and cModel magnitudes per filter, for each field." + "### 4.1. Extendedness\n", + "\n", + "Plot the magnitude distributions for PSF and cModel magnitudes per filter, for each field, separating point sources and extended sources via the `refSizeExtendedness` parameter." ] }, { @@ -1022,7 +1027,6 @@ " ax[1].legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n", " fig.subplots_adjust(wspace=0)\n", " fig.suptitle(name)\n", - " # plt.tight_layout()\n", " plt.show()\n", "\n", " del region, query\n", @@ -1047,7 +1051,7 @@ "id": "46277fad-7109-402b-90e2-cdf68ec8a1e2", "metadata": {}, "source": [ - "Display, for each field, the number of point-like (PL) and potentially-extended (EX) objects in the search area, and estimate the 2D sky density in arcminutes per square degree. This is only an estimate, and not a robust measurement; the search radius was only 0.2 deg." + "Display, for each field, the number of point-like (PL) and potentially-extended (EX) objects in the search area, and estimate the 2D sky density in number of sources per square arcminute. This is only an estimate, and not a robust measurement; the search radius was only 0.2 deg." ] }, { @@ -1085,12 +1089,225 @@ "region_galactic[tx]" ] }, + { + "cell_type": "markdown", + "id": "c446c034-0713-4ca8-8015-273a4c0d81cb", + "metadata": {}, + "source": [ + "### 4.2. Color-magnitude and color-color diagrams\n", + "\n", + "Extract `Object` catalog data with a larger radius (0.5 deg.) for a high-latitude field (Rubin_SV_212_-7) and a relatively low-latitude field (Rubin_SV_225_-40). Compare their color-magnitude and color-color diagrams." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d816774-90dd-41b8-9933-f70f2bc8a5de", + "metadata": {}, + "outputs": [], + "source": [ + "fields = ['Rubin_SV_212_-7', 'Rubin_SV_225_-40']\n", + "search_radius = 0.5\n", + "objcats = {}\n", + "\n", + "for field in fields:\n", + " region = regions[field]\n", + " query = \"\"\"SELECT g_psfMag, r_psfMag, i_psfMag,\n", + " g_cModelMag, r_cModelMag, i_cModelMag,\n", + " refSizeExtendedness\n", + " FROM dp2.Object\n", + " WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),\n", + " CIRCLE('ICRS', {}, {}, {})) = 1\n", + " AND i_psfFlux_flag = 0 and i_cModel_flag = 0\n", + " \"\"\".format(region[0], region[1], search_radius)\n", + " job = tap_service.submit_job(query)\n", + " job.run()\n", + " job.wait(phases=['COMPLETED', 'ERROR'])\n", + " assert job.phase == 'COMPLETED'\n", + " object_table = job.fetch_result().to_table()\n", + " objcats[field] = object_table" + ] + }, + { + "cell_type": "markdown", + "id": "2206c2d4-8904-4ea5-9df6-310cf7c78025", + "metadata": {}, + "source": [ + "Set the `refSizeExtendedness` threshold to 0.25, to select a more pure sample of point sources." + ] + }, { "cell_type": "code", "execution_count": null, "id": "e27fc30c-3451-4b99-a118-e2343d638250", "metadata": {}, "outputs": [], + "source": [ + "cat_sv212_m7 = objcats['Rubin_SV_212_-7']\n", + "cat_sv225_m40 = objcats['Rubin_SV_225_-40']\n", + "\n", + "ptsource_sv212_m7 = np.where(cat_sv212_m7['refSizeExtendedness'] < 0.25)[0]\n", + "extsource_sv212_m7 = np.where(cat_sv212_m7['refSizeExtendedness'] >= 0.25)[0]\n", + "ptsource_sv225_m40 = np.where(cat_sv225_m40['refSizeExtendedness'] < 0.25)[0]\n", + "extsource_sv225_m40 = np.where(cat_sv225_m40['refSizeExtendedness'] >= 0.25)[0]" + ] + }, + { + "cell_type": "markdown", + "id": "f3678274-a62d-46cf-b15c-b35e3a76a726", + "metadata": {}, + "source": [ + "Plot the color-magnitude ($r$ vs. $g-r$) and color-color ($r-i$ vs. $g-r$) diagrams for each field, separating point and extended sources." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "143f1dd7-512d-42f5-bde4-4654dc9eb5c2", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(7, 6),\n", + " height_ratios=[1, 2.5])\n", + "plt.subplots_adjust(hspace=0, wspace=0)\n", + "\n", + "grmin, grmax = -0.9, 2.3\n", + "rimin, rimax = -1.3, 2.8\n", + "magmin, magmax = 15.8, 26.8\n", + "\n", + "ax1.hexbin(cat_sv212_m7[ptsource_sv212_m7]['g_psfMag']-cat_sv212_m7[ptsource_sv212_m7]['r_psfMag'],\n", + " cat_sv212_m7[ptsource_sv212_m7]['r_psfMag']-cat_sv212_m7[ptsource_sv212_m7]['i_psfMag'],\n", + " gridsize=(150, 150),\n", + " extent=(grmin, grmax, rimin, rimax), bins='log', cmap='Grays')\n", + "ax1.set_title('Point sources (PSF mags)')\n", + "ax1.set_xlim(grmin, grmax)\n", + "ax1.set_ylim(rimin, rimax)\n", + "ax1.set_ylabel(r'(r−i)')\n", + "ax1.set_xticklabels([])\n", + "ax1.minorticks_on()\n", + "\n", + "ax2.hexbin(cat_sv212_m7[extsource_sv212_m7]['g_cModelMag']-cat_sv212_m7[extsource_sv212_m7]['r_cModelMag'],\n", + " cat_sv212_m7[extsource_sv212_m7]['r_cModelMag']-cat_sv212_m7[extsource_sv212_m7]['i_cModelMag'],\n", + " gridsize=(150, 150),\n", + " extent=(grmin, grmax, rimin, rimax), bins='log', cmap='Grays')\n", + "ax2.set_title('Extended sources (cModel mags)')\n", + "ax2.set_xlim(grmin, grmax)\n", + "ax2.set_ylim(rimin, rimax)\n", + "ax2.set_xticklabels([])\n", + "ax2.set_yticklabels([])\n", + "ax2.minorticks_on()\n", + "\n", + "ax3.hexbin(cat_sv212_m7[ptsource_sv212_m7]['g_psfMag']-cat_sv212_m7[ptsource_sv212_m7]['r_psfMag'],\n", + " cat_sv212_m7[ptsource_sv212_m7]['r_psfMag'], gridsize=(150, 200),\n", + " extent=(grmin, grmax, magmin, magmax), bins='log', cmap='Grays')\n", + "ax3.set_xlim(grmin, grmax)\n", + "ax3.set_ylim(magmax, magmin)\n", + "ax3.set_xlabel(r'(g−r)')\n", + "ax3.set_ylabel(r'r magnitude')\n", + "ax3.minorticks_on()\n", + "\n", + "ax4.hexbin(cat_sv212_m7[extsource_sv212_m7]['g_cModelMag']-cat_sv212_m7[extsource_sv212_m7]['r_cModelMag'],\n", + " cat_sv212_m7[extsource_sv212_m7]['r_cModelMag'], gridsize=(150, 200),\n", + " extent=(grmin, grmax, magmin, magmax), bins='log', cmap='Grays')\n", + "ax4.set_xlim(grmin, grmax)\n", + "ax4.set_ylim(magmax, magmin)\n", + "ax4.set_xlabel(r'(g−r)')\n", + "ax4.set_yticklabels([])\n", + "ax4.minorticks_on()\n", + "\n", + "plt.suptitle('Rubin_SV_212_-7; Galactic latitude: 51.3 deg.')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "56ec7ebf-dea3-4dfc-a3af-aea22801c438", + "metadata": {}, + "source": [ + "> **Figure 8:** Color-magnitude (CMD; $r$ vs. $g-r$) and color-color ($r-i$ vs. $g-r$) diagrams for the SV_212_-7 field. The left column shows likely stars (`refSizeExtendedness` < 0.25) using `PSF` magnitudes, and the right column shows likely galaxies using `cModel` magnitudes. Top panels display color-color diagrams, and the bottom panels show the CMDs. This high-latitude field has few stellar features in the color-color diagrams, as it is mainly populated by galaxies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "108ebc3f-4c6e-4f7b-aab2-9feb447904da", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(7, 6),\n", + " height_ratios=[1, 2.5])\n", + "plt.subplots_adjust(hspace=0, wspace=0)\n", + "\n", + "grmin, grmax = -0.9, 2.3\n", + "rimin, rimax = -1.3, 2.8\n", + "magmin, magmax = 15.8, 26.8\n", + "\n", + "ax1.hexbin(cat_sv225_m40[ptsource_sv225_m40]['g_psfMag']-cat_sv225_m40[ptsource_sv225_m40]['r_psfMag'],\n", + " cat_sv225_m40[ptsource_sv225_m40]['r_psfMag']-cat_sv225_m40[ptsource_sv225_m40]['i_psfMag'],\n", + " gridsize=(150, 150),\n", + " extent=(grmin, grmax, rimin, rimax), bins='log', cmap='Grays')\n", + "ax1.set_title('Point sources (PSF mags)')\n", + "ax1.set_xlim(grmin, grmax)\n", + "ax1.set_ylim(rimin, rimax)\n", + "ax1.set_ylabel(r'(r−i)')\n", + "ax1.set_xticklabels([])\n", + "ax1.minorticks_on()\n", + "\n", + "ax2.hexbin(cat_sv225_m40[extsource_sv225_m40]['g_cModelMag']-cat_sv225_m40[extsource_sv225_m40]['r_cModelMag'],\n", + " cat_sv225_m40[extsource_sv225_m40]['r_cModelMag']-cat_sv225_m40[extsource_sv225_m40]['i_cModelMag'],\n", + " gridsize=(150, 150),\n", + " extent=(grmin, grmax, rimin, rimax), bins='log', cmap='Grays')\n", + "ax2.set_title('Extended sources (cModel mags)')\n", + "ax2.set_xlim(grmin, grmax)\n", + "ax2.set_ylim(rimin, rimax)\n", + "ax2.set_xticklabels([])\n", + "ax2.set_yticklabels([])\n", + "ax2.minorticks_on()\n", + "\n", + "ax3.hexbin(cat_sv225_m40[ptsource_sv225_m40]['g_psfMag']-cat_sv225_m40[ptsource_sv225_m40]['r_psfMag'],\n", + " cat_sv225_m40[ptsource_sv225_m40]['r_psfMag'], gridsize=(150, 200),\n", + " extent=(grmin, grmax, magmin, magmax), bins='log', cmap='Grays')\n", + "ax3.set_xlim(grmin, grmax)\n", + "ax3.set_ylim(magmax, magmin)\n", + "ax3.set_xlabel(r'(g−r)')\n", + "ax3.set_ylabel(r'r magnitude')\n", + "ax3.minorticks_on()\n", + "\n", + "ax4.hexbin(cat_sv225_m40[extsource_sv225_m40]['g_cModelMag']-cat_sv225_m40[extsource_sv225_m40]['r_cModelMag'],\n", + " cat_sv225_m40[extsource_sv225_m40]['r_cModelMag'], gridsize=(150, 200),\n", + " extent=(grmin, grmax, magmin, magmax), bins='log', cmap='Grays')\n", + "ax4.set_xlim(grmin, grmax)\n", + "ax4.set_ylim(magmax, magmin)\n", + "ax4.set_xlabel(r'(g−r)')\n", + "ax4.set_yticklabels([])\n", + "ax4.minorticks_on()\n", + "\n", + "plt.suptitle('Rubin_SV_225_-40; Galactic latitude: 16.9 deg.')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "05502c15-9823-4e22-a20f-749be500187e", + "metadata": {}, + "source": [ + "> **Figure 9:** Color-magnitude (CMD; $r$ vs. $g-r$) and color-color ($r-i$ vs. $g-r$) diagrams for the SV_225_-40 field. The left column shows likely stars (`refSizeExtendedness` < 0.25) using `PSF` magnitudes, and the right column shows likely galaxies using `cModel` magnitudes. Top panels display color-color diagrams, and the bottom panels show the CMDs. This low-latitude field has a clearly defined stellar locus in the color-color diagram." + ] + }, + { + "cell_type": "markdown", + "id": "c0ecee91-d091-45de-8642-d902976bb72b", + "metadata": {}, + "source": [ + "As expected, the lower latitude field has many more Milky Way stars and a very clearly defined stellar locus in the point-source color-color diagram. Some Milky Way stars are evident in the high-latitude field, but the majority of objects are galaxies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86af5a58-08c8-4da4-99e7-e7790ec82d07", + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/DP2/300_Science_demos/301_DP2_observations/301_3_Deep_drilling_fields.ipynb b/DP2/300_Science_demos/301_DP2_observations/301_3_Deep_drilling_fields.ipynb index a5b02b8e..26e42890 100644 --- a/DP2/300_Science_demos/301_DP2_observations/301_3_Deep_drilling_fields.ipynb +++ b/DP2/300_Science_demos/301_DP2_observations/301_3_Deep_drilling_fields.ipynb @@ -21,8 +21,8 @@ "For the Rubin Science Platform at data.lsst.cloud. \\\n", "Data Release: [Data Preview 2](http://dp2.lsst.io) \\\n", "Container Size: Large \\\n", - "LSST Science Pipelines version: v30.0.10 \\\n", - "Last verified to run: 2026-07-24 \\\n", + "LSST Science Pipelines version: v30.0.11 \\\n", + "Last verified to run: 2026-08-17 \\\n", "Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)" ] }, @@ -147,7 +147,9 @@ "import matplotlib.pyplot as plt\n", "import matplotlib.colors as mcolors\n", "import itertools\n", + "from astropy.coordinates import SkyCoord\n", "from astropy.table import Table, join\n", + "import astropy.units as u\n", "import skyproj\n", "import healpy as hp\n", "import hpgeom as hpg\n", @@ -296,7 +298,7 @@ "id": "a500c7ca-134b-45e0-80ba-dec5e7397ccb", "metadata": {}, "source": [ - "Define a function that will remove a figure and all it's data, to help clear memory." + "Define a function that will remove a figure and all its data, to help clear memory." ] }, { @@ -426,7 +428,7 @@ "id": "ce49c560-766b-407c-b74e-a804b85a9d36", "metadata": {}, "source": [ - "Use the `skyproj` package with a McBryde skyprojection to visualize the distribution of DP2 visits on the sky in HEALPix that are approximately the same size as the LSSTCam FOV (`nside`=19)." + "Use the `skyproj` package with a McBryde sky projection to visualize the distribution of DP2 visits on the sky in HEALPix that are approximately the same size as the LSSTCam FOV (`nside`=19)." ] }, { @@ -437,7 +439,7 @@ "outputs": [], "source": [ "orig_cmap = plt.get_cmap('Blues')\n", - "new_cmap = truncate_colormap(orig_cmap, 0.10, 0.75)\n", + "new_cmap = truncate_colormap(orig_cmap, 0.20, 0.8)\n", "\n", "fig, ax = plt.subplots(figsize=(12, 8))\n", "sp = skyproj.McBrydeSkyproj(ax=ax)\n", @@ -486,6 +488,7 @@ "metadata": {}, "outputs": [], "source": [ + "sc_all = SkyCoord(ra = visit_table['ra'], dec = visit_table['dec'], frame='icrs')\n", "data_start = np.floor(np.min(visit_table['expMidptMJD']))\n", "data_end = np.ceil(np.max(visit_table['expMidptMJD']))\n", "temp = []\n", @@ -493,8 +496,8 @@ "fig, ax = plt.subplots(len(regions), 2, figsize=(10, 20))\n", "for r, name in enumerate(region_names):\n", " ra, dec, rad = regions[name]\n", - " offsets = np.sqrt((visit_table['ra']-ra)**2 + (visit_table['dec']-dec)**2)\n", - " tx = np.where(offsets < rad)[0]\n", + " sc_region = SkyCoord(ra = ra*u.deg, dec = dec*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_all).value <= rad)[0]\n", " for f, filt in enumerate(filter_names):\n", " fx = np.where(visit_table['band'][tx] == filt)[0]\n", " if len(fx) > 0:\n", @@ -516,7 +519,7 @@ " tdiffs.append(comb[1] - comb[0])\n", " time_diffs = np.array(tdiffs)\n", " temp.append((len(tx), np.mean(time_diffs)))\n", - " del ra, dec, offsets, tx, dates, tdiffs, time_diffs\n", + " del ra, dec, sc_region, tx, dates, tdiffs, time_diffs\n", " ax[r, 0].set_xlabel('MJD - ' + str(data_start))\n", " ax[r, 0].set_ylabel('Cumulative visits')\n", " ax[r, 0].set_xlim([0, data_end - data_start])\n", @@ -845,14 +848,15 @@ "metadata": {}, "outputs": [], "source": [ + "sc_patches = SkyCoord(ra = patches_table['s_ra'], dec = patches_table['s_dec'], frame='icrs')\n", + "\n", "shortlist = [name for name in region_names if name != 'DDF_XMM_LSS']\n", + "\n", "for r, name in enumerate(shortlist):\n", " coords = regions[name]\n", " radius = coords[2] + np.sqrt(2*1.66**2)\n", - " delta_ra = (coords[0] - patches_table['s_ra'])*(np.cos(np.deg2rad(coords[1])))\n", - " delta_dec = coords[1] - patches_table['s_dec']\n", - " offset = np.sqrt(delta_ra**2 + delta_dec**2)\n", - " tx = np.where(offset < radius)[0]\n", + " sc_region = SkyCoord(ra = coords[0]*u.deg, dec = coords[1]*u.deg, frame='icrs')\n", + " tx = np.where(sc_region.separation(sc_patches).value <= radius)[0]\n", " pixels = np.unique(np.array(hpg.angle_to_pixel(32,\n", " patches_table['s_ra'][tx],\n", " patches_table['s_dec'][tx])))\n", @@ -893,7 +897,7 @@ "\n", " plt.tight_layout()\n", " plt.show()\n", - " del coords, radius, delta_ra, delta_dec, offset\n", + " del coords, radius\n", " del tx, pixels, hspmap\n", " del query, bind, refs, deep_coadd, show_tract, show_patch\n", " del sp\n", @@ -989,7 +993,6 @@ " ax[1].legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n", " fig.subplots_adjust(wspace=0)\n", " fig.suptitle(name)\n", - " # plt.tight_layout()\n", " plt.show()\n", "\n", " del region, query\n", @@ -1033,10 +1036,202 @@ "nobj_i_table" ] }, + { + "cell_type": "markdown", + "id": "d1c2f8b6-c8d0-4fe7-9a50-fec6abd62532", + "metadata": {}, + "source": [ + "# 5. Cluster red sequence in ELAIS_S1\n", + "\n", + "PSZ2 G309.43-72.86 ([Léget 2026](https://arxiv.org/abs/2606.09938)) is a known rich cluster at redshift=0.35 covered by the Data Preview 2 (DP2).\n", + "Its coordinates are within the DDF_ELAIS_S1 field.\n", + "\n", + "Query for objects near this position, and select red sequence galaxies, then confirm that they cluster on the sky as expected." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26e718d4-9c35-401a-8cb9-534e063b27c9", + "metadata": {}, + "outputs": [], + "source": [ + "ra_cl, dec_cl = 10.196, -44.134\n", + "\n", + "range_deg = 0.1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76dd03c6-d4c8-45fa-801b-8eb6666bcbef", + "metadata": {}, + "outputs": [], + "source": [ + "query = \"SELECT objectId, coord_ra, coord_dec, refSizeExtendedness, \" \\\n", + " \"deblend_peak_center_x as x, deblend_peak_center_y as y, \" \\\n", + " \"r_cModelMag, i_cModelMag, z_cModelMag, \" \\\n", + " \"r_cModelMagErr, i_cModelMagErr, z_cModelMagErr, \" \\\n", + " \"r_cModelFlux/r_cModelFluxErr AS r_cModel_SNR, \" \\\n", + " \"i_cModelFlux/i_cModelFluxErr AS i_cModel_SNR, \" \\\n", + " \"z_cModelFlux/z_cModelFluxErr AS z_cModel_SNR, \" \\\n", + " \"r_cModel_flag, i_cModel_flag, z_cModel_flag \" \\\n", + " \"FROM dp2.Object \" \\\n", + " \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \" \\\n", + " \"CIRCLE('ICRS', %f, %f, %f)) = 1 \"%(ra_cl, dec_cl, range_deg)\n", + "\n", + "job = tap_service.submit_job(query)\n", + "job.run()\n", + "job.wait(phases=['COMPLETED', 'ERROR'])\n", + "print('Job phase is', job.phase)\n", + "if job.phase == 'ERROR':\n", + " job.raise_if_error()\n", + "assert job.phase == 'COMPLETED'\n", + "obj_cat = job.fetch_result().to_table()\n", + "print('Number of objects: ', len(obj_cat))" + ] + }, + { + "cell_type": "markdown", + "id": "1a210736-0760-4a40-8fe8-430edd4f173b", + "metadata": {}, + "source": [ + "Select extended objects, and filter out flagged and low SNR objects." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c2ad33f5-a328-4e7a-b3bb-e3240239e424", + "metadata": {}, + "outputs": [], + "source": [ + "band_list = [\"r\", \"i\", \"z\"]\n", + "\n", + "sel = obj_cat[\"refSizeExtendedness\"] > 0.5\n", + "\n", + "for band in band_list:\n", + " sel &= (obj_cat[f\"{band}_cModel_flag\"] == False)\n", + " sel &= obj_cat[f'{band}_cModel_SNR'] > 3.\n", + "\n", + "obj_cat_sel = obj_cat[sel]" + ] + }, + { + "cell_type": "markdown", + "id": "742bc8d2-ecda-4668-9374-6b0ec26220dd", + "metadata": {}, + "source": [ + "Plot color-magnitude diagrams, then select the red sequence (the \"plume\" of stars at roughly the same color)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48e56519-567f-46a6-bca1-be6d44b93f21", + "metadata": {}, + "outputs": [], + "source": [ + "color_dict = {\n", + " \"ri\": 0.6,\n", + " \"iz\": 0.3,\n", + "}\n", + "\n", + "eps = 0.125\n", + "sel_rs = np.array([True] * len(obj_cat_sel))\n", + "mag_max = 23\n", + "\n", + "for ind in range(len(band_list)-1):\n", + "\n", + " band1 = band_list[ind]\n", + " band2 = band_list[ind+1]\n", + "\n", + " fig, ax = plt.subplots(figsize=(4, 3), layout=\"constrained\")\n", + "\n", + " ax.scatter(obj_cat[f\"{band1}_cModelMag\"],\n", + " obj_cat[f\"{band1}_cModelMag\"] - obj_cat[f\"{band2}_cModelMag\"],\n", + " marker=\".\", s=1, c='Gray', label='all objects')\n", + " ax.scatter(obj_cat_sel[f\"{band1}_cModelMag\"],\n", + " obj_cat_sel[f\"{band1}_cModelMag\"] - obj_cat_sel[f\"{band2}_cModelMag\"],\n", + " marker=\".\", s=1, c='Red', label='extended')\n", + "\n", + " ax.legend()\n", + " ax.set_xlim([18, 25])\n", + " ax.set_ylim([-1.5, 2.5])\n", + "\n", + " ax.set_xlabel(band1)\n", + " ax.set_ylabel(f\"{band1}-{band2}\")\n", + "\n", + " color = color_dict[f\"{band1}{band2}\"]\n", + " ax.axhline(color - eps, ls=\"--\", c=\"k\", alpha=0.3)\n", + " ax.axhline(color + eps, ls=\"--\", c=\"k\", alpha=0.3)\n", + "\n", + " sel_rs &= np.abs(obj_cat_sel[f\"{band1}_cModelMag\"] - obj_cat_sel[f\"{band2}_cModelMag\"] - color) < eps\n", + " sel_rs &= obj_cat_sel[f\"{band1}_cModelMag\"] < mag_max\n", + "\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "2802a612-5f3a-4bac-b79f-fd68a9aa619b", + "metadata": {}, + "source": [ + "> **Figure 8:** Color-magnitude diagrams for all objects (gray symbols) and potentially-extended (red) objects in a portion of the ELAIS-S1 field. Each panel contains a clear \"plume\" of extended objects at a roughly fixed color, spanning a large range in magnitude. This is the \"red sequence\" of the galaxy cluster." + ] + }, + { + "cell_type": "markdown", + "id": "5928dfe7-c12a-4bf6-90c1-09c5c52d8466", + "metadata": {}, + "source": [ + "Plot positions of all objects, and overlay the selected red sequence candidates." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "444074eb-d354-4bc9-aa8c-ae4efbb8c88f", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(4, 3), layout=\"constrained\")\n", + "ax.scatter(obj_cat[\"coord_ra\"],\n", + " obj_cat[\"coord_dec\"],\n", + " s=1., alpha=0.3, c='Gray',\n", + " label='all objects')\n", + "ax.scatter(obj_cat_sel[\"coord_ra\"][sel_rs],\n", + " obj_cat_sel[\"coord_dec\"][sel_rs],\n", + " s=3., c='Red',\n", + " label='extended')\n", + "ax.plot([ra_cl], [dec_cl], 'k*', label='cluster center')\n", + "ax.legend(loc='upper left')\n", + "ax.invert_xaxis()\n", + "ax.set_xlabel(\"RA [deg]\")\n", + "ax.set_ylabel(\"DEC [deg]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "14a43dab-afc3-424a-a408-3bb0a292beb2", + "metadata": {}, + "source": [ + "> **Figure 9:** Spatial distribution for all objects (gray points) and potentially-extended (red) objects. The extended objects are clumped near the position of the cluster center, which is shown as a black star." + ] + }, + { + "cell_type": "markdown", + "id": "c5fa38a8-043e-483e-82ef-d99d4006f527", + "metadata": {}, + "source": [ + "Notice that the red objects (selected to be galaxies in the red sequence) cluster near the center, suggesting that our cuts have successfully identified the galaxy cluster." + ] + }, { "cell_type": "code", "execution_count": null, - "id": "e27fc30c-3451-4b99-a118-e2343d638250", + "id": "afb82a98-fc07-4979-9ccf-ced7982db75b", "metadata": {}, "outputs": [], "source": []