switch to visit image

Author: MelissaGraham

DP1/100_How_to_Use_RSP_Tools/103_Image_access_and_display/103_4_Small_Image_Cutout.ipynb

@@ -30,13 +30,12 @@
    "source": [
     "**Learning objective:** How to use the Rubin image cutout service with DP1.\n",
     "\n",
-    "**LSST data products:** `deep_coadd`\n",
+    "**LSST data products:** `visit_image`\n",
     "\n",
     "**Packages:** `lsst.rsp.utils`, `pyvo`, `lsst.rsp.service`, `lsst.afw.display`\n",
     "\n",
     "**Credit:** Originally developed by the Rubin Community Science team.\n",
     "Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.\n",
-    "This notebook builds on DP0.2 science demonstration of image cutouts (tutorial notebook 13a), and technical exploration of the cutout service developed by Leanne Guy (in preparation), and builds on an earlier notebook written by Alex Drlica-Wagner <a href=\"https://github.com/rubin-dp0/cst-dev/blob/main/ADW_gravelpit/TAP_Image_Access.ipynb\">linked here</a>.\n",
     "\n",
     "**Get Support:**\n",
     "Everyone is encouraged to ask questions or raise issues in the \n",
@@ -109,7 +108,8 @@
     "from pyvo.dal.adhoc import DatalinkResults, SodaQuery\n",
     "\n",
     "from astropy import units as u\n",
-    "from astropy.coordinates import Angle"
+    "from astropy.coordinates import Angle\n",
+    "from astropy.time import Time"
    ]
   },
   {
@@ -160,12 +160,16 @@
    "id": "3d4d9571-2f48-440e-9c40-90f14f6eac8c",
    "metadata": {},
    "source": [
-    "## 2. Locating the image for cutout\n",
+    "## 2. Find the visit image\n",
     "\n",
-    "First, define a point on the sky as the center of the image cutout. Once the RA and Dec are defined, create a SpherePoint object to define the location on the sky. Use these in a SIA query to identify the overlapping deep_coadd image in Section 2.4.  \n",
+    "The cutout services needs the `access_url` for the image from which a cutout is desired.\n",
+    "The SIA or TAP services can be used to find the desired image and retrieve its `access_url`.\n",
     "\n",
-    "Define search coordinates right ascension (`target_ra`) and declination (`target_dec`), in degrees, to use in all queries.\n",
-    "These coordinates are the location of the Hubble Ultra Deep Field, near the center of the Extended Chandra Deep Field South (ECDFS) DP1 map."
+    "For this example, make an r-band (effective wavelength 622.1 nm) cutout centered on a set of coordinates in the ECDFS field,\n",
+    "for a visit image that was obtained between MJD 60623.256 and 60623.259.\n",
+    "\n",
+    "Define the coordinates right ascension (`target_ra`) and declination (`target_dec`) in degrees.\n",
+    "Define the band and the start and end times."
    ]
   },
   {
@@ -175,10 +179,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "target_ra = 53.1567053\n",
-    "target_dec = -27.7815854\n",
-    "\n",
-    "spherePoint = geom.SpherePoint(target_ra*geom.degrees, target_dec*geom.degrees)"
+    "target_ra = 53.1246023\n",
+    "target_dec = -27.7404715\n",
+    "eff_wl = 622.1e-09\n",
+    "time1 = Time(60623.256, format=\"mjd\", scale=\"tai\")\n",
+    "time2 = Time(60623.259, format=\"mjd\", scale=\"tai\")"
    ]
   },
   {
@@ -214,9 +219,7 @@
    "id": "08b71293-d3ab-49ca-9a95-2566488eec33",
    "metadata": {},
    "source": [
-    "Query for `deep_coadd` images by setting `calib_level=3` and `dpsubtype='lsst.deep_coadd'`.\n",
-    "\n",
-    "This query will return 6 `deep_coadd` images (one for each filter)."
+    "This query will return 1 `visit_image` (by design)."
    ]
   },
   {
@@ -226,7 +229,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = service.search(pos=circle, calib_level=3, dpsubtype='lsst.deep_coadd')\n",
+    "results = service.search(pos=circle, calib_level=2, dpsubtype='lsst.visit_image',\n",
+    "                         band=eff_wl, time=(time1, time2))\n",
     "print(len(results))"
    ]
   },
@@ -235,7 +239,7 @@
    "id": "8f62bb47-7a7c-49cd-969a-6d60b5488abe",
    "metadata": {},
    "source": [
-    "Option to display the results as an Astropy table."
+    "Display the results as an Astropy table."
    ]
   },
   {
@@ -245,62 +249,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# results.to_table()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "2e641dc9-57ae-4dc5-9bee-6e54f48d88dc",
-   "metadata": {},
-   "source": [
-    "Select the r-band deep coadd image."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "82a36217-98b7-4633-8535-1a055d7d2ff9",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "tx = np.where(results['lsst_band'] == 'r')[0]\n",
-    "print(tx)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "98235dce-1afe-42c8-9d01-0a82d8b8ab2a",
-   "metadata": {},
-   "source": [
-    "Define the corresponding `index` for the `results`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "59c760db-7ecd-42e2-b8bf-1ae4910ab4a5",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "index = tx[0].item()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "977ed070-845c-4145-908a-6d68f2c913bf",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "results[index]"
+    "results.to_table()"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "5ddf1d6c-7acb-4467-8a32-4659bd45ab5b",
    "metadata": {},
    "source": [
-    "In the table, the `access_url` contains the web URL datalink for the r-band `deep_coadd`. This datalink will be needed to generate the image cutout."
+    "In the table, the `access_url` contains the web URL datalink for the image. This datalink will be needed to generate the image cutout."
    ]
   },
   {
@@ -315,13 +272,17 @@
     "An example TAP query using this method that retrieves the `access_url` is:\n",
     "\n",
     "```\n",
-    "SELECT dataproduct_type, dataproduct_subtype, calib_level, lsst_band,\n",
-    "       instrument_name, access_url, access_format \n",
+    "SELECT dataproduct_type,dataproduct_subtype,calib_level,lsst_band,em_min,em_max,lsst_tract,lsst_patch,\n",
+    "       lsst_filter,lsst_visit,lsst_detector,t_exptime,t_min,t_max,s_ra,s_dec,s_fov,obs_id,\n",
+    "       obs_collection,o_ucd,facility_name,instrument_name,obs_title,s_region,access_url,\n",
+    "       access_format \n",
     "FROM ivoa.ObsCore \n",
-    "WHERE calib_level = 3\n",
-    "AND dataproduct_type = 'image'\n",
-    "AND dataproduct_subtype = 'lsst.deep_coadd'\n",
-    "AND CONTAINS(POINT('ICRS', 53.1567053, -27.7815854), s_region)=1\n",
+    "WHERE CONTAINS(POINT('ICRS', 53.1567053, -27.7815854), s_region)=1\n",
+    "      AND obs_collection = 'LSST.DP1' AND calib_level = 2\n",
+    "      AND dataproduct_type = 'image' AND instrument_name = 'LSSTComCam'\n",
+    "      AND dataproduct_subtype = 'lsst.visit_image'\n",
+    "      AND ( t_min <= 60623.259 AND 60623.256 <= t_max )\n",
+    "      AND ( 622e-9 BETWEEN em_min AND em_max )\n",
     "```"
    ]
   },
@@ -346,7 +307,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "datalink_url = results[index].access_url\n",
+    "datalink_url = results[0].access_url\n",
     "dl_result = DatalinkResults.from_result_url(datalink_url, session=get_pyvo_auth())\n",
     "\n",
     "f\"Datalink status: {dl_result.status}. Datalink service url: {datalink_url}\""
@@ -377,7 +338,7 @@
    "id": "c6fa46da-88c8-4b5e-b86c-4b5878f3badf",
    "metadata": {},
    "source": [
-    "The variable `sq` now holds the result of the SODA query using the data link (which currently still points the full LSST patch and tract deep_coadd image, at its remote location in the database). The cell below will now demonstrate how to extract a cutout from `sq`.  \n",
+    "The variable `sq` now holds the result of the SODA query using the data link (which currently still points the full LSST `visit_image`, at its remote location in the database). The cell below will now demonstrate how to extract a cutout from `sq`.  \n",
     "\n",
     "### 3.1. Define cutout center and edge\n",
     "\n",
@@ -391,7 +352,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Radius = 0.06 * u.deg\n",
+    "spherePoint = geom.SpherePoint(target_ra*geom.degrees, target_dec*geom.degrees)\n",
+    "Radius = 0.01 * u.deg\n",
     "sq.circle = (spherePoint.getRa().asDegrees() * u.deg,\n",
     "             spherePoint.getDec().asDegrees() * u.deg,\n",
     "             Radius)"
@@ -447,7 +409,7 @@
    "id": "458a4d25-ca95-4f92-b76e-986b7a03b320",
    "metadata": {},
    "source": [
-    "> Figure 1: The `deep_coadd` cutout image, displayed in grayscale with a scale bar at right.\n",
+    "> Figure 1: The cutout image, displayed in grayscale with a scale bar at right.\n",
     "\n",
     "#### 3.2.1. Option to save cutout to disk\n",
     "\n",
@@ -502,7 +464,9 @@
     "\n",
     "POLYGON=12 34 14 34 14 36 12 36\n",
     "\n",
-    "Since the center of the cutout is already defined in RA and Dec in the cells above (`spherePoint`), this example will define each x,y set as RA+/-edge and Dec+/-edge.   "
+    "Since the center of the cutout is already defined in RA and Dec in the cells above (`spherePoint`), this example will define each x,y set as RA+/-edge and Dec+/-edge.\n",
+    "\n",
+    "> **Warning:** Visit images are rotated, and although it is the \"Declination edge\" that is defined to be smaller, that corresponds to the x-axis when the cutout is displayed below, due to image rotation."
    ]
   },
   {
@@ -516,15 +480,17 @@
     "                              dl_result.get_adhocservice_by_id(\"cutout-sync\"),\n",
     "                              session=get_pyvo_auth())\n",
     "\n",
-    "edge1 = 0.06 * u.deg\n",
-    "edge2 = 0.02 * u.deg\n",
+    "ra_edge = 0.02 * u.deg\n",
+    "de_edge = 0.005 * u.deg\n",
     "\n",
-    "sqp.polygon = (spherePoint.getRa().asDegrees() * u.deg - edge1,\n",
-    "               spherePoint.getDec().asDegrees() * u.deg - edge2,\n",
-    "               spherePoint.getRa().asDegrees() * u.deg - edge1,\n",
-    "               spherePoint.getDec().asDegrees() * u.deg + edge2,\n",
-    "               spherePoint.getRa().asDegrees() * u.deg + edge1,\n",
-    "               spherePoint.getDec().asDegrees() * u.deg + edge2)\n",
+    "sqp.polygon = (spherePoint.getRa().asDegrees() * u.deg - ra_edge,\n",
+    "               spherePoint.getDec().asDegrees() * u.deg - de_edge,\n",
+    "               spherePoint.getRa().asDegrees() * u.deg - ra_edge,\n",
+    "               spherePoint.getDec().asDegrees() * u.deg + de_edge,\n",
+    "               spherePoint.getRa().asDegrees() * u.deg + ra_edge,\n",
+    "               spherePoint.getDec().asDegrees() * u.deg + de_edge,\n",
+    "               spherePoint.getRa().asDegrees() * u.deg + ra_edge,\n",
+    "               spherePoint.getDec().asDegrees() * u.deg - de_edge)\n",
     "\n",
     "cutout_bytes = sqp.execute_stream().read()\n",
     "sqp.raise_if_error()\n",
@@ -551,8 +517,7 @@
    "id": "c83e5413-2edf-41fd-8542-b622ba8324b0",
    "metadata": {},
    "source": [
-    "> Figure 2: A rectangular polygon cutout from the same `deep_coadd` image used in Figure 1.\n",
-    "\n"
+    "> Figure 2: A rectangular polygon cutout from a rotated `visit_image`."
    ]
   },
   {
@@ -562,7 +527,7 @@
    "source": [
     "### 3.4. Correcting for cos(d)\n",
     "\n",
-    "There is an important difference to note between the circle and polygon shape definitions. The angular distance on the sky that defines the circular cutout size already accounts for the difference in angular distance in the RA direction is smaller by a factor of cos(declination), where declination is in units radians. The difference increases with higher declination. However, the polygon definition does not automatically account for this cosine factor. Thus, circle and polygon cutout definitions using the same cutout edge length will not match size in the RA direction (for deep_coadds, the x-axis). The 2 cells below demonstrate how to make this correction to the polygon cutout definition to create symmetric cutouts with polygon. Here, reset the edge sizes to be the same as `Radius` from the circle definition above.\n",
+    "There is an important difference to note between the circle and polygon shape definitions. The angular distance on the sky that defines the circular cutout size already accounts for the difference in angular distance in the RA direction is smaller by a factor of cos(declination), where declination is in units radians. The difference increases with higher declination. However, the polygon definition does not automatically account for this cosine factor. Thus, circle and polygon cutout definitions using the same cutout edge length will not match size in the RA direction. The 2 cells below demonstrate how to make this correction to the polygon cutout definition to create symmetric cutouts with polygon. Here, reset the edge sizes to be the same as `Radius` from the circle definition above.\n",
     "\n",
     "First, generate a polygon cutout without factoring in cos(dec)."
    ]