From 51f0fb9cfc4a9a858e922d60480fa4eb919973f1 Mon Sep 17 00:00:00 2001 From: Abdu Zoghbi Date: Mon, 20 Feb 2023 09:40:18 -0500 Subject: [PATCH] converted ipynb to md and added jupytext as a dependency; changes to ci --- .github/workflows/ci_tests_and_publish.yml | 3 +- CS_Catalog_Queries.ipynb | 632 ------------------ CS_Catalog_Queries.md | 389 +++++++++++ CS_Image_Access.ipynb | 350 ---------- CS_Image_Access.md | 216 +++++++ CS_Spectral_Access.ipynb | 211 ------ CS_Spectral_Access.md | 135 ++++ CS_UCDs.ipynb | 275 -------- CS_UCDs.md | 170 +++++ CS_VO_Tables.ipynb | 145 ----- CS_VO_Tables.md | 105 +++ Exercise_I.ipynb | 434 ------------- Exercise_I.md | 231 +++++++ Exercise_II.ipynb | 323 ---------- Exercise_II.md | 179 ++++++ Exercise_III.ipynb | 486 -------------- Exercise_III.md | 295 +++++++++ QuickReference.ipynb | 423 ------------ QuickReference.md | 273 ++++++++ UseCase_I.ipynb | 709 --------------------- UseCase_I.md | 452 +++++++++++++ UseCase_II.ipynb | 434 ------------- UseCase_II.md | 260 ++++++++ UseCase_III.ipynb | 661 ------------------- UseCase_III.md | 430 +++++++++++++ check_env.py | 3 +- conf.py | 3 + doc-requirements.txt | 1 + environment.yml | 3 +- tox.ini | 1 + 30 files changed, 3146 insertions(+), 5086 deletions(-) delete mode 100644 CS_Catalog_Queries.ipynb create mode 100644 CS_Catalog_Queries.md delete mode 100644 CS_Image_Access.ipynb create mode 100644 CS_Image_Access.md delete mode 100644 CS_Spectral_Access.ipynb create mode 100644 CS_Spectral_Access.md delete mode 100644 CS_UCDs.ipynb create mode 100644 CS_UCDs.md delete mode 100644 CS_VO_Tables.ipynb create mode 100644 CS_VO_Tables.md delete mode 100644 Exercise_I.ipynb create mode 100644 Exercise_I.md delete mode 100644 Exercise_II.ipynb create mode 100644 Exercise_II.md delete mode 100644 Exercise_III.ipynb create mode 100644 Exercise_III.md delete mode 100644 QuickReference.ipynb create mode 100644 QuickReference.md delete mode 100644 UseCase_I.ipynb create mode 100644 UseCase_I.md delete mode 100644 UseCase_II.ipynb create mode 100644 UseCase_II.md delete mode 100644 UseCase_III.ipynb create mode 100644 UseCase_III.md diff --git a/.github/workflows/ci_tests_and_publish.yml b/.github/workflows/ci_tests_and_publish.yml index 2e7da58..63a5337 100644 --- a/.github/workflows/ci_tests_and_publish.yml +++ b/.github/workflows/ci_tests_and_publish.yml @@ -74,6 +74,7 @@ jobs: - run: pip install -r doc-requirements.txt - run: | python check_env.py + jupytext --to notebook *.md pytest --nbval publish: @@ -99,4 +100,4 @@ jobs: with: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./_build/html/ - commit_message: ${{ github.event.head_commit.message }} \ No newline at end of file + commit_message: ${{ github.event.head_commit.message }} diff --git a/CS_Catalog_Queries.ipynb b/CS_Catalog_Queries.ipynb deleted file mode 100644 index 5c827fc..0000000 --- a/CS_Catalog_Queries.ipynb +++ /dev/null @@ -1,632 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Accessing astronomical catalogs\n", - "\n", - "There are two ways to access astronomical data catalogs that are provided as table data with a VO API.\n", - "\n", - "First, there is a __[Simple Cone Search (SCS) protocol](http://www.ivoa.net/documents/latest/ConeSearch.html)__ used to search a given table with a given position and radius, getting back a table of results. The interface requires only a position and search radius. \n", - "\n", - "For more complicated searches, the __[Table Access Protocol](http://www.ivoa.net/documents/TAP/)__ (TAP) protocol is a powerful tool to search any VO table. Here, we expand on its usage and that of the __[Astronomical Data Query Language](http://www.ivoa.net/documents/latest/ADQL.html)__ (ADQL) that it uses. \n", - "\n", - "\n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# suppress some specific warnings that are not important\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")\n", - "\n", - "import io\n", - "import numpy as np\n", - "\n", - "# Astropy imports\n", - "import astropy.units as u\n", - "import astropy.constants as const\n", - "from astropy.coordinates import SkyCoord\n", - "from astropy.cosmology import Planck15\n", - "from astropy.io import votable as apvot\n", - "import scipy.integrate\n", - "\n", - "## Generic VO access routines\n", - "import pyvo as vo\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "## 1. Simple cone search\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Starting with a single simple source first: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "coord = SkyCoord.from_name(\"m51\")\n", - "print(coord)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Below, we go through the exercise of how we can figure out the most relevant table. But for now, let's assume that we know that we want the CFA redshift catalog refered to as 'zcat'. VO services are listed in a central Registry that can be searched through a [web interface](http://vao.stsci.edu/keyword-search/) or using PyVO's `regsearch`. We use the registry to find the corresponding cone service and then submit our cone search. \n", - "\n", - "Registry services are of the following type: \n", - "* simple cone search: \"scs\"\n", - "* table access protocol: \"tap\" or \"table\"\n", - "* simple image search: \"sia\" or \"image\"\n", - "* simple spectral access: \"ssa\"\n", - "\n", - "There are a number of things in the registry related to 'zcat' so we find the specific one that we want, which is the CFA version:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='scs', keywords=['zcat'])\n", - "services.to_table()['ivoid', 'short_name', 'res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Supposing that we want the table with the short_name CFAZ, and we want to retrieve the data for all sources within an arcminute of our specified location:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Use the one that's CFAZ. \n", - "## Use list comprehension to check each service's short_name attribute and use the first.\n", - "cfaz_cone_service = [s for s in services if 'CFAZ' in s.short_name][0] \n", - "\n", - "## We are searching for sources within 10 arcminutes of M51. \n", - "results = cfaz_cone_service.search(pos=coord, radius=10*u.arcmin)\n", - "results.to_table()\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The SCS is quite straightforward and returns all of the columns of the given table (which can be anything) for the sources in the region queried." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "## 2. Table Access Protocol queries\n", - "\n", - "A TAP query is the most powerful way to search a catalog. A Simple Cone Search only allows you to ask for a position and radius, but TAP allows you to do much more, since the available tables contain much more information. \n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.1 TAP services\n", - "\n", - "Many services list a single TAP service in the Registry that can access many catalogs, boosting your efficiency, and letting you add constraints based on any column. This is the power of the TAP! \n", - "\n", - "Suppose for our example, we want to select bright galaxy candidates but don't know the coordinates. Therefore, we start from figuring out the best table to query. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As before, we use the `vo.regsearch()` for a servicetype 'tap'. There are a lot of TAP services in the registry, but they are listed slightly differently than cone services. The metadata on each catalog is usually published in the registry with its cone service, and then the full TAP service is listed as an \"auxiliary\" service. So to find a TAP service for a given catalog, we need to add the option *includeaux=True*. Alternatively, you can start with a single TAP service and then ask it specifically which tables it serves, but for this use case, that is less efficient. \n", - "\n", - "We'll first do a registry search for all auxiliary TAP services related to the \"CfA\" and \"redshift\". " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "tap_services = vo.regsearch(servicetype='tap', keywords=['cfa redshift'], includeaux=True)\n", - "tap_services.to_table()['ivoid', 'short_name', 'res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "There are many tables that mention these keywords. Pick some likely looking ones and look at the descriptions:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "for t in tap_services:\n", - " if \"cfa\" in t.res_title.lower() and \"magnitude\" in t.res_description.lower():\n", - " print(f\"{t.ivoid}: {t.res_description}\\n\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "From the above information, you can choose the table you want and then use the specified TAP service to query it as described below. \n", - "\n", - "But first we'll look at the other way of finding tables to query with TAP: by starting with the TAP services listed individually in the Registry. We see above that the HEASARC has the ZCAT, so what else does it have?" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can find out which tables a TAP serves and then look at the tables descriptions. The last line here sends a query directly to the service to ask it for a list of tables. (This can take a minute, since there may be a lot of tables.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Here, we're looking for a specific service, and we don't need the includeaux option:\n", - "tap_services = vo.regsearch(servicetype='tap',keywords=['heasarc'])\n", - "heasarc = tap_services[0]\n", - "heasarc_tables=heasarc.service.tables " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Then let's look for tables matching the terms we're interested in as above. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for tablename in heasarc_tables.keys():\n", - " if \"redshift\" in heasarc_tables[tablename].description.lower(): \n", - " heasarc_tables[tablename].describe()\n", - " print(\"Columns={}\".format(sorted([k.name for k in heasarc_tables[tablename].columns ])))\n", - " print(\"----\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "There are a number of tables that appear to be useful table for our goal, including the ZCAT, which contains columns with the information that we need to select a sample of the brightest nearby spiral galaxy candidates. \n", - "\n", - "Now that we know all the possible column information in the zcat catalog, we can do more than query on position (as in a cone search) but also on any other column (e.g., redshift, bmag, morph_type). The query has to be expressed in a language called __[ADQL](http://www.ivoa.net/documents/latest/ADQL.html)__. \n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "### 2.2 Expressing queries in ADQL" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "The basics of ADQL:\n", - "\n", - "* *SELECT * FROM my.interesting.catalog as cat...* \n", - "\n", - "says you want all (\"*\") columns from the catalog called \"my.interesting.catalog\", which you will refer to in the rest of the query by the more compact name of \"cat\". \n", - "\n", - "Instead of returning all columns, you can \n", - "\n", - "* *SELECT cat.RA, cat.DEC, cat.bmag from catalog as cat...* \n", - "\n", - "to only return the columns you're interested in. To use multiple catalogs, your query could start, e.g.,\n", - "\n", - "* *SELECT c1.RA,c1.DEC,c2.BMAG FROM catalog1 as c1 natural join catalog2 as c2...* \n", - "\n", - "says that you want to query two catalogs zipped together the \"natural\" way, i.e., by looking for a common column.\n", - "\n", - "To select only some rows of the catalog based on the value in a column, you can add: \n", - "\n", - "* *WHERE cat.bmag < 14* \n", - "\n", - "says that you want to retrieve only those entries in the catalog whose bmag column has a value less than 14.\n", - "\n", - "You can also append \n", - "\n", - "* *ORDER by cat.bmag* \n", - "\n", - "to return the result sorted ascending by one of the columns, adding *DESC* to the end for descending. \n", - "\n", - "A few special functions in the ADQL allow you to query regions:\n", - "\n", - "* *WHERE contains( point('ICRS', cat.ra, cat.dec), circle('ICRS', 210.5, -6.5, 0.5))=1*\n", - "\n", - "is how you would ask for any catalog entries whose RA,DEC lie within a circular region defined by RA,DEC 210.5,-6.5 and a radius of 0.5 (all in degrees). The 'ICRS' specifies the coordinate system. \n", - "\n", - "See the ADQL documentation for more. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "\n", - "### 2.3 A use case \n", - "\n", - "Here is a simple ADQL query where we print out the relevant columns for the bright (Bmag <14) sources found within 1 degree of M51 (we will discuss how to define the table and column names below):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Inside the format call, the {} are replaced by the given variables in order.\n", - "## So this asks for \n", - "## rows of public.zcat where that row's ra and dec (cat.ra and cat.dec from the catalog) \n", - "## are within radius 1deg of the given RA and DEC we got above for M51 \n", - "## (coord.ra.deg and coord.dec.deg from our variables defined above), and where \n", - "## the bmag column is less than 14. \n", - "query = \"\"\"SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where \n", - " contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{},{},1.0))=1 and\n", - " cat.bmag < 14\n", - " order by cat.radial_velocity_error \n", - " \"\"\".format(coord.ra.deg, coord.dec.deg)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "results=heasarc.service.run_async(query) \n", - "#results = heasarc.search(query)\n", - "results.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "See the __[information on the zcat](https://heasarc.gsfc.nasa.gov/W3Browse/galaxy-catalog/zcat.html)__ for column information. (We will use the 'radial_velocity' column rather than the 'redshift' column.) We note that spiral galaxies have morph_type between 1 - 9. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Therefore, we can generalize the query above to complete our exercise and select the brightest (bmag < 14), nearby (radial velocity < 3000), spiral ( morph_type = 1 - 9) galaxies as follows: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "query = \"\"\"SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where \n", - " cat.bmag < 14 and cat.morph_type between 1 and 9 and cat.Radial_Velocity < 3000 \n", - " order by cat.Radial_velocity \n", - " \"\"\".format(coord.ra.deg, coord.dec.deg)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "results = heasarc.service.run_async(query)\n", - "#results = heasarc.search(query)\n", - "results.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "\n", - "\n", - "### 2.4 TAP examples for a given service\n", - "\n", - "Each service may also provide some example queries. If they do, then you can see them with, e.g.:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "heasarc.service.examples[0]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Above, these examples look like a list of dictionaries. But they are actually a list of objects that can be executed:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for example in heasarc.service.examples:\n", - " print(example['QUERY'])\n", - " result=example.execute()\n", - " # Stop at one\n", - " break\n", - "result.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "## 3. Using the TAP to cross-correlate and combine" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 3.1 Cross-correlating to combine catalogs\n", - "\n", - "TAP can also be a powerful way to collect a lot of useful information from existing catalogs in one quick step. For this exercise, we will start with a list of sources, uploaded from our own table, and do a 'cross-correlation' with the *zcat* table. \n", - "\n", - "For more on creating and working with VO tables, see that [notebook](CS_VO_Tables.ipynb). Here, we just read one in that's already prepared: \n", - "\n", - "First, check that this service can handle uploaded tables. Not all do. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "heasarc.service.upload_methods" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The inline method is what PyVO will use. These take a while, i.e. half a minute." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "query=\"\"\"\n", - " SELECT cat.ra, cat.dec, Radial_Velocity, bmag, morph_type\n", - " FROM zcat cat, tap_upload.mysources mt \n", - " WHERE\n", - " contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',mt.ra,mt.dec,0.01))=1\n", - " and Radial_Velocity > 0\n", - " ORDER by cat.ra\"\"\"\n", - "zcattable = heasarc.service.run_async(query, uploads={'mysources': 'data/my_sources.xml'})\n", - "#zcattable = heasarc.search(query, uploads={'mysources': 'data/my_sources.xml'})\n", - "mytable = zcattable.to_table()\n", - "mytable" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Therefore we now have the Bmag, morphological type and radial velocities for all the sources in our list with a single TAP query. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "### 3.2 Cross-correlating with user-defined columns\n", - "\n", - "Our input list of sources contains galaxy pair candidates that may be interacting with each other. Therefore it would be interesting to know what the morphological type and the Bmagnitude are for the potential companions. \n", - "\n", - "In this advanced example, we want our search to be physically motivated since the criterion for galaxy interaction depends on the physical separation of the galaxies. Unlike the previous case, the search radius is not a constant, but varies for each candidate by the distance to the source. Specifically, we want to search for companions that are within 50 kpc of the candidate and therefore first need to find the angular diameter distance that corresponds to galaxy's distance (in our case the radial velocity).\n", - "\n", - "Therefore, we begin by taking our table of objects and adding an angDdeg column:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## The column 'radial_velocity' is c*z but doesn't include the unit; it is km/s\n", - "## Get the speed of light from astropy.constants and express in km/s\n", - "c = const.c.to(u.km/u.s).value \n", - "redshifts = mytable['radial_velocity']/c \n", - "mytable['redshift'] = redshifts\n", - "physdist = 0.05*u.Mpc # 50 kpc physical distance\n", - "\n", - "angDdist = Planck15.angular_diameter_distance(mytable['redshift'].data)\n", - "angDrad = np.arctan(physdist/angDdist)\n", - "mytable['angDdeg'] = angDrad.to(u.deg)\n", - "mytable" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we construct and run a query that uses the new angDdeg column in every row search. Note, we also don't want to list the original candidates since we know these are in the catalog and we want rather to find any companions. Therefore, we exclude the match if the radial velocities match exactly.\n", - "\n", - "This time, rather than write the table to disk, we'll keep it in memory and give Tap.query() a \"file-like\" object using io.BytesIO(). This can take half a minute: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## In memory only, use an IO stream. \n", - "vot_obj=io.BytesIO()\n", - "apvot.writeto(apvot.from_table(mytable),vot_obj)\n", - "## (Reset the \"file-like\" object to the beginning.)\n", - "vot_obj.seek(0)\n", - "query=\"\"\"SELECT mt.ra, mt.dec, cat.ra, cat.dec, cat.Radial_Velocity, cat.morph_type, cat.bmag \n", - " FROM zcat cat, tap_upload.mytable mt \n", - " WHERE\n", - " contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',mt.ra,mt.dec,mt.angDdeg))=1\n", - " and cat.Radial_Velocity > 0 and cat.radial_velocity != mt.radial_velocity\n", - " ORDER by cat.ra\"\"\"\n", - "# Currently broken due to a bug.\n", - "#mytable2 = heasarc.service.run_async(query, uploads={'mytable':vot_obj})\n", - "mytable2 = heasarc.search(query, uploads={'mytable':vot_obj})\n", - "vot_obj.close()\n", - "mytable2.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Therefore, by adding new information to our original data table, we could cross-correlate with the TAP. We find that, in our candidate list, there is one true pair of galaxies. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "\n", - "## 4 Synchronous versus asynchronous queries\n", - "\n", - "There is one technical detail about TAP queries that you will need to know. In the code cells above, there are two commands for sending the query, one of which is commented out. This is because, with the TAP, there are two ways to send such queries. The default when you use the `search()` method is to us a synchronous query, which means that the query is sent and the client waits for the response. For large and complicated queries, this may time out, or you may want to run several in parallel. So there are other options. \n", - "\n", - "The method `service.run_async()` uses an asynchronous query, which means that the query is sent, and then (under the hood without you needing to do anything), the method checks for a response. From your point of view, these methods look the same; PyVO is doing different things under the hood, but the method will not return until it has your result. \n", - "\n", - "You need to know about these two methods for a couple of reasons. First, some services will limit synchronous queries, i.e. they will not necessarily return *all* the results if there are too many of them. An asynchronous query should have no such restrictions. In the case of the HEASARC service that we use above, it does not matter, but you should be aware of this and be in the habit of using the asynchronous queries for complete results after an initial interactive exploration. \n", - "\n", - "The second reason to be aware of this is that asynchronous queries may be queued by the service, and they can take a lot longer if the service is very busy or the job is very large. (The synchronous option in this case may either time out, or it may return quickly but with incomplete results.)\n", - "\n", - "For very large queries, or for submitting queries in parallel, you may wish to use the `submit_job()`, `wait()`, and `fetch_results()` methods to avoid locking up your Python session. This is described in the [pyvo documentation](https://pyvo.readthedocs.io/en/latest/dal/index.html#jobs).\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.9" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": { - "height": "calc(100% - 180px)", - "left": "10px", - "top": "150px", - "width": "345.6px" - }, - "toc_section_display": true, - "toc_window_display": true - }, - "widgets": { - "state": {}, - "version": "1.1.1" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/CS_Catalog_Queries.md b/CS_Catalog_Queries.md new file mode 100644 index 0000000..b27f5d6 --- /dev/null +++ b/CS_Catalog_Queries.md @@ -0,0 +1,389 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.10.11 +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: + height: calc(100% - 180px) + left: 10px + top: 150px + width: 345.6px + toc_section_display: true + toc_window_display: true +widgets: + state: {} + version: 1.1.1 +--- + +# Accessing astronomical catalogs + +There are two ways to access astronomical data catalogs that are provided as table data with a VO API. + +First, there is a __[Simple Cone Search (SCS) protocol](http://www.ivoa.net/documents/latest/ConeSearch.html)__ used to search a given table with a given position and radius, getting back a table of results. The interface requires only a position and search radius. + +For more complicated searches, the __[Table Access Protocol](http://www.ivoa.net/documents/TAP/)__ (TAP) protocol is a powerful tool to search any VO table. Here, we expand on its usage and that of the __[Astronomical Data Query Language](http://www.ivoa.net/documents/latest/ADQL.html)__ (ADQL) that it uses. + +- [Accessing astronomical catalogs](#accessing-astronomical-catalogs) + - [1. Simple cone search](#1-simple-cone-search) + - [2. Table Access Protocol queries](#2-table-access-protocol-queries) + - [2.1 TAP services](#21-tap-services) + - [2.2 Expressing queries in ADQL](#22-expressing-queries-in-adql) + - [2.3 A use case](#23-a-use-case) + - [2.4 TAP examples for a given service](#24-tap-examples-for-a-given-service) + - [3. Using the TAP to cross-correlate and combine](#3-using-the-tap-to-cross-correlate-and-combine) + - [3.1 Cross-correlating to combine catalogs](#31-cross-correlating-to-combine-catalogs) + - [3.2 Cross-correlating with user-defined columns](#32-cross-correlating-with-user-defined-columns) + - [4. Synchronous versus asynchronous queries](#4-synchronous-versus-asynchronous-queries) + +```{code-cell} ipython3 +# suppress some specific warnings that are not important +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") + +import io +import numpy as np + +# Astropy imports +import astropy.units as u +import astropy.constants as const +from astropy.coordinates import SkyCoord +from astropy.cosmology import Planck15 +from astropy.io import votable as apvot +import scipy.integrate + +## Generic VO access routines +import pyvo as vo +``` + +## 1. Simple cone search + ++++ + +Starting with a single simple source first: + +```{code-cell} ipython3 +coord = SkyCoord.from_name("m51") +print(coord) +``` + +Below, we go through the exercise of how we can figure out the most relevant table. But for now, let's assume that we know that we want the CFA redshift catalog refered to as 'zcat'. VO services are listed in a central Registry that can be searched through a [web interface](http://vao.stsci.edu/keyword-search/) or using PyVO's `regsearch`. We use the registry to find the corresponding cone service and then submit our cone search. + +Registry services are of the following type: +* simple cone search: "scs" +* table access protocol: "tap" or "table" +* simple image search: "sia" or "image" +* simple spectral access: "ssa" + +There are a number of things in the registry related to 'zcat' so we find the specific one that we want, which is the CFA version: + +```{code-cell} ipython3 +services = vo.regsearch(servicetype='scs', keywords=['zcat']) +services.to_table()['ivoid', 'short_name', 'res_title'] +``` + +Supposing that we want the table with the short_name CFAZ, and we want to retrieve the data for all sources within an arcminute of our specified location: + +```{code-cell} ipython3 +## Use the one that's CFAZ. +## Use list comprehension to check each service's short_name attribute and use the first. +cfaz_cone_service = [s for s in services if 'CFAZ' in s.short_name][0] + +## We are searching for sources within 10 arcminutes of M51. +results = cfaz_cone_service.search(pos=coord, radius=10*u.arcmin) +results.to_table() +``` + +The SCS is quite straightforward and returns all of the columns of the given table (which can be anything) for the sources in the region queried. + ++++ + +## 2. Table Access Protocol queries + +A TAP query is the most powerful way to search a catalog. A Simple Cone Search only allows you to ask for a position and radius, but TAP allows you to do much more, since the available tables contain much more information. + ++++ + +### 2.1 TAP services + +Many services list a single TAP service in the Registry that can access many catalogs, boosting your efficiency, and letting you add constraints based on any column. This is the power of the TAP! + +Suppose for our example, we want to select bright galaxy candidates but don't know the coordinates. Therefore, we start from figuring out the best table to query. + ++++ + +As before, we use the `vo.regsearch()` for a servicetype 'tap'. There are a lot of TAP services in the registry, but they are listed slightly differently than cone services. The metadata on each catalog is usually published in the registry with its cone service, and then the full TAP service is listed as an "auxiliary" service. So to find a TAP service for a given catalog, we need to add the option *includeaux=True*. Alternatively, you can start with a single TAP service and then ask it specifically which tables it serves, but for this use case, that is less efficient. + +We'll first do a registry search for all auxiliary TAP services related to the "CfA" and "redshift". + +```{code-cell} ipython3 +tap_services = vo.regsearch(servicetype='tap', keywords=['cfa redshift'], includeaux=True) +tap_services.to_table()['ivoid', 'short_name', 'res_title'] +``` + +There are many tables that mention these keywords. Pick some likely looking ones and look at the descriptions: + +```{code-cell} ipython3 +for t in tap_services: + if "cfa" in t.res_title.lower() and "magnitude" in t.res_description.lower(): + print(f"{t.ivoid}: {t.res_description}\n") +``` + +From the above information, you can choose the table you want and then use the specified TAP service to query it as described below. + +But first we'll look at the other way of finding tables to query with TAP: by starting with the TAP services listed individually in the Registry. We see above that the HEASARC has the ZCAT, so what else does it have? + ++++ + +You can find out which tables a TAP serves and then look at the tables descriptions. The last line here sends a query directly to the service to ask it for a list of tables. (This can take a minute, since there may be a lot of tables.) + +```{code-cell} ipython3 +# Here, we're looking for a specific service, and we don't need the includeaux option: +tap_services = vo.regsearch(servicetype='tap',keywords=['heasarc']) +heasarc = tap_services[0] +heasarc_tables=heasarc.service.tables +``` + +Then let's look for tables matching the terms we're interested in as above. + +```{code-cell} ipython3 +for tablename in heasarc_tables.keys(): + if "redshift" in heasarc_tables[tablename].description.lower(): + heasarc_tables[tablename].describe() + print("Columns={}".format(sorted([k.name for k in heasarc_tables[tablename].columns ]))) + print("----") +``` + +There are a number of tables that appear to be useful table for our goal, including the ZCAT, which contains columns with the information that we need to select a sample of the brightest nearby spiral galaxy candidates. + +Now that we know all the possible column information in the zcat catalog, we can do more than query on position (as in a cone search) but also on any other column (e.g., redshift, bmag, morph_type). The query has to be expressed in a language called __[ADQL](http://www.ivoa.net/documents/latest/ADQL.html)__. + ++++ + +### 2.2 Expressing queries in ADQL + ++++ + +The basics of ADQL: + +* *SELECT * FROM my.interesting.catalog as cat...* + +says you want all ("*") columns from the catalog called "my.interesting.catalog", which you will refer to in the rest of the query by the more compact name of "cat". + +Instead of returning all columns, you can + +* *SELECT cat.RA, cat.DEC, cat.bmag from catalog as cat...* + +to only return the columns you're interested in. To use multiple catalogs, your query could start, e.g., + +* *SELECT c1.RA,c1.DEC,c2.BMAG FROM catalog1 as c1 natural join catalog2 as c2...* + +says that you want to query two catalogs zipped together the "natural" way, i.e., by looking for a common column. + +To select only some rows of the catalog based on the value in a column, you can add: + +* *WHERE cat.bmag < 14* + +says that you want to retrieve only those entries in the catalog whose bmag column has a value less than 14. + +You can also append + +* *ORDER by cat.bmag* + +to return the result sorted ascending by one of the columns, adding *DESC* to the end for descending. + +A few special functions in the ADQL allow you to query regions: + +* *WHERE contains( point('ICRS', cat.ra, cat.dec), circle('ICRS', 210.5, -6.5, 0.5))=1* + +is how you would ask for any catalog entries whose RA,DEC lie within a circular region defined by RA,DEC 210.5,-6.5 and a radius of 0.5 (all in degrees). The 'ICRS' specifies the coordinate system. + +See the ADQL documentation for more. + ++++ + +### 2.3 A use case + +Here is a simple ADQL query where we print out the relevant columns for the bright (Bmag <14) sources found within 1 degree of M51 (we will discuss how to define the table and column names below): + +```{code-cell} ipython3 +## Inside the format call, the {} are replaced by the given variables in order. +## So this asks for +## rows of public.zcat where that row's ra and dec (cat.ra and cat.dec from the catalog) +## are within radius 1deg of the given RA and DEC we got above for M51 +## (coord.ra.deg and coord.dec.deg from our variables defined above), and where +## the bmag column is less than 14. +query = """SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where + contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{},{},1.0))=1 and + cat.bmag < 14 + order by cat.radial_velocity_error + """.format(coord.ra.deg, coord.dec.deg) +``` + +```{code-cell} ipython3 +results=heasarc.service.run_async(query) +#results = heasarc.search(query) +results.to_table() +``` + +See the __[information on the zcat](https://heasarc.gsfc.nasa.gov/W3Browse/galaxy-catalog/zcat.html)__ for column information. (We will use the 'radial_velocity' column rather than the 'redshift' column.) We note that spiral galaxies have morph_type between 1 - 9. + ++++ + +Therefore, we can generalize the query above to complete our exercise and select the brightest (bmag < 14), nearby (radial velocity < 3000), spiral ( morph_type = 1 - 9) galaxies as follows: + +```{code-cell} ipython3 +query = """SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where + cat.bmag < 14 and cat.morph_type between 1 and 9 and cat.Radial_Velocity < 3000 + order by cat.Radial_velocity + """.format(coord.ra.deg, coord.dec.deg) +``` + +```{code-cell} ipython3 +results = heasarc.service.run_async(query) +#results = heasarc.search(query) +results.to_table() +``` + +### 2.4 TAP examples for a given service + +Each service may also provide some example queries. If they do, then you can see them with, e.g.: + +```{code-cell} ipython3 +heasarc.service.examples[0] +``` + +Above, these examples look like a list of dictionaries. But they are actually a list of objects that can be executed: + +```{code-cell} ipython3 +for example in heasarc.service.examples: + print(example['QUERY']) + result=example.execute() + # Stop at one + break +result.to_table() +``` + +## 3. Using the TAP to cross-correlate and combine + ++++ + +### 3.1 Cross-correlating to combine catalogs + +TAP can also be a powerful way to collect a lot of useful information from existing catalogs in one quick step. For this exercise, we will start with a list of sources, uploaded from our own table, and do a 'cross-correlation' with the *zcat* table. + +For more on creating and working with VO tables, see that [notebook](CS_VO_Tables.md). Here, we just read one in that's already prepared: + +First, check that this service can handle uploaded tables. Not all do. + +```{code-cell} ipython3 +heasarc.service.upload_methods +``` + +The inline method is what PyVO will use. These take a while, i.e. half a minute. + +```{code-cell} ipython3 +query=""" + SELECT cat.ra, cat.dec, Radial_Velocity, bmag, morph_type + FROM zcat cat, tap_upload.mysources mt + WHERE + contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',mt.ra,mt.dec,0.01))=1 + and Radial_Velocity > 0 + ORDER by cat.ra""" +zcattable = heasarc.service.run_async(query, uploads={'mysources': 'data/my_sources.xml'}) +#zcattable = heasarc.search(query, uploads={'mysources': 'data/my_sources.xml'}) +mytable = zcattable.to_table() +mytable +``` + +Therefore we now have the Bmag, morphological type and radial velocities for all the sources in our list with a single TAP query. + ++++ + +### 3.2 Cross-correlating with user-defined columns + +Our input list of sources contains galaxy pair candidates that may be interacting with each other. Therefore it would be interesting to know what the morphological type and the Bmagnitude are for the potential companions. + +In this advanced example, we want our search to be physically motivated since the criterion for galaxy interaction depends on the physical separation of the galaxies. Unlike the previous case, the search radius is not a constant, but varies for each candidate by the distance to the source. Specifically, we want to search for companions that are within 50 kpc of the candidate and therefore first need to find the angular diameter distance that corresponds to galaxy's distance (in our case the radial velocity). + +Therefore, we begin by taking our table of objects and adding an angDdeg column: + +```{code-cell} ipython3 +## The column 'radial_velocity' is c*z but doesn't include the unit; it is km/s +## Get the speed of light from astropy.constants and express in km/s +c = const.c.to(u.km/u.s).value +redshifts = mytable['radial_velocity']/c +mytable['redshift'] = redshifts +physdist = 0.05*u.Mpc # 50 kpc physical distance + +angDdist = Planck15.angular_diameter_distance(mytable['redshift'].data) +angDrad = np.arctan(physdist/angDdist) +mytable['angDdeg'] = angDrad.to(u.deg) +mytable +``` + +Now we construct and run a query that uses the new angDdeg column in every row search. Note, we also don't want to list the original candidates since we know these are in the catalog and we want rather to find any companions. Therefore, we exclude the match if the radial velocities match exactly. + +This time, rather than write the table to disk, we'll keep it in memory and give Tap.query() a "file-like" object using io.BytesIO(). This can take half a minute: + +```{code-cell} ipython3 +## In memory only, use an IO stream. +vot_obj=io.BytesIO() +apvot.writeto(apvot.from_table(mytable),vot_obj) +## (Reset the "file-like" object to the beginning.) +vot_obj.seek(0) +query="""SELECT mt.ra, mt.dec, cat.ra, cat.dec, cat.Radial_Velocity, cat.morph_type, cat.bmag + FROM zcat cat, tap_upload.mytable mt + WHERE + contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',mt.ra,mt.dec,mt.angDdeg))=1 + and cat.Radial_Velocity > 0 and cat.radial_velocity != mt.radial_velocity + ORDER by cat.ra""" +# Currently broken due to a bug. +#mytable2 = heasarc.service.run_async(query, uploads={'mytable':vot_obj}) +mytable2 = heasarc.search(query, uploads={'mytable':vot_obj}) +vot_obj.close() +mytable2.to_table() +``` + +Therefore, by adding new information to our original data table, we could cross-correlate with the TAP. We find that, in our candidate list, there is one true pair of galaxies. + ++++ + +## 4. Synchronous versus asynchronous queries + +There is one technical detail about TAP queries that you will need to know. In the code cells above, there are two commands for sending the query, one of which is commented out. This is because, with the TAP, there are two ways to send such queries. The default when you use the `search()` method is to us a synchronous query, which means that the query is sent and the client waits for the response. For large and complicated queries, this may time out, or you may want to run several in parallel. So there are other options. + +The method `service.run_async()` uses an asynchronous query, which means that the query is sent, and then (under the hood without you needing to do anything), the method checks for a response. From your point of view, these methods look the same; PyVO is doing different things under the hood, but the method will not return until it has your result. + +You need to know about these two methods for a couple of reasons. First, some services will limit synchronous queries, i.e. they will not necessarily return *all* the results if there are too many of them. An asynchronous query should have no such restrictions. In the case of the HEASARC service that we use above, it does not matter, but you should be aware of this and be in the habit of using the asynchronous queries for complete results after an initial interactive exploration. + +The second reason to be aware of this is that asynchronous queries may be queued by the service, and they can take a lot longer if the service is very busy or the job is very large. (The synchronous option in this case may either time out, or it may return quickly but with incomplete results.) + +For very large queries, or for submitting queries in parallel, you may wish to use the `submit_job()`, `wait()`, and `fetch_results()` methods to avoid locking up your Python session. This is described in the [pyvo documentation](https://pyvo.readthedocs.io/en/latest/dal/index.html#jobs). diff --git a/CS_Image_Access.ipynb b/CS_Image_Access.ipynb deleted file mode 100644 index 7c6b249..0000000 --- a/CS_Image_Access.ipynb +++ /dev/null @@ -1,350 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Searching for and retrieving images\n", - "\n", - "In this notebook, we show how to search for and retrieve images from VO services using the Registry and the __[Simple Image Access](http://www.ivoa.net/documents/SIA/)__ (SIA) protocol.\n", - "\n", - "\n", - "
\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**\\*Note:** for all of these notebooks, the results depend on real-time queries. Sometimes there are problems, either because a given service has changed, is undergoing maintenance, or the internet connectivity is having problems, etc. Always retry a couple of times, come back later and try again, and only then send us the problem report to investigate." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "\n", - "import numpy as np\n", - "\n", - "import matplotlib\n", - "import matplotlib.pyplot as plt\n", - "%matplotlib inline \n", - "\n", - "import pyvo as vo\n", - "\n", - "from astropy.io import fits\n", - "import astropy.coordinates as coord\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "\n", - "from IPython.display import Image as ipImage, display\n", - "\n", - "# There are a number of relatively unimportant warnings that show up, so for now, suppress them:\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "## 1. Finding SIA resources from the Registry\n", - "\n", - "First, how do we find out what services are available? These are listed in a registry at STScI (__[see here](http://www.ivoa.net/documents/RegTAP/)__). Our Registry function gives a simple interface for how to search for services. \n", - "\n", - "Let's search for services providing images in the ultraviolet bands:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "uv_services = vo.regsearch(servicetype='image',waveband='uv')\n", - "uv_services.to_table()['ivoid','short_name','res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This returns an astropy table containing information about the services available. We can then specify the service we want by using the corresponding row. We'll repeat the search with additional qualifiers to isolate the row we want (note that in the keyword search the \"%\" character can be used as a wild card):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "uvot_services = vo.regsearch(servicetype='image',waveband='uv',keywords=['swift'])\n", - "uvot_services.to_table()['ivoid','short_name','res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This shows us that the data we are interested in comes from the HEASARC's SkyView service, but the point of these VO tools is that you don't need to know that ahead of time or indeed to care where it comes from." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "\n", - "## 2. Using SIA to retrieve an image\n", - "\n", - "Now we look for images of our favorite source. See __[the SIA definition](http://www.ivoa.net/documents/WD/SIA/sia-20040524.html)__ for usage. In short, you can specify the central position and the size (degrees as one or two floats for the RA, DEC directions). It is up to the service to determine how to provide this. Optionally, you can limit it to the format you want, e.g., \"image/fits\" or \"image/png\" etc. \n", - "\n", - "What is returned to you is not the image itself but a list of images available and how to access them. This is easiest shown by example: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "coords = coord.SkyCoord.from_name(\"m51\")\n", - "\n", - "im_table = uvot_services[0].search(pos=coords,size=0.2,format='image/jpeg')\n", - "im_table.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Extract the fields you're interested in, e.g., the URLs of the images made by skyview. Note that specifying as we did SwiftUVOT, we get a number of different images, e.g., UVOT U, V, B, W1, W2, etc. For each survey, there are two URLs, first the FITS IMAGE and second the JPEG. \n", - "\n", - "Note that different services will return different column names, but all will have a column giving the URL to access the image. Though it has different column names in different services, it can always be accessed through the `getdataurl` function." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "url = im_table[0].getdataurl()\n", - "print(url)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 3. Viewing the resulting image\n", - "\n", - "\n", - "### JPG images\n", - "\n", - "Since we have asked for JPEG images, we can display an image in python easily by using its URL. Each row of the result has a getdataurl() method, and you can then hand the URL to an image displayer such as IPython.display:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "img = ipImage(url=im_table[0].getdataurl())\n", - "display(img)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Fits files\n", - "\n", - "Or download the FITS image and display it with imshow, or aplpy.\n", - "\n", - "(This often errors off with a time out message. Just try it again, possibly a couple of times.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Do the search again asking for FITS\n", - "im_table = uvot_services[0].search(pos=coords,size=0.2,format='image/fits')\n", - "\n", - "# Hand the url of the first result to fits.open()\n", - "hdu_list = fits.open(im_table[0].getdataurl())\n", - "hdu_list.info()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Using imshow" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(hdu_list[0].data, cmap='gray', origin='lower',vmax=0.1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 4. Example of data available through multiple services\n", - "\n", - "Suppose we want Sloan DSS data. A generic query finds us a number of possibilities (note that this doesn't work for keywords=['sdss']; be flexible and try several search terms):" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='image', keywords=['sloan'], waveband='optical')\n", - "services.to_table()[np.where(np.isin(services.to_table()['short_name'], 'SDSSDR7'))]['ivoid', 'short_name']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "So one of these is served by SDSS's SkyServer and the other by HEASARC's SkyView. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Using HEASARC" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "heasarc_dr7_service = [s for s in services if 'SDSSDR7' in s.short_name and 'heasarc' in s.ivoid][0]\n", - "\n", - "sdss_table_heasarc = heasarc_dr7_service.search(pos=coords,size=0.2,format='image/fits')\n", - "sdss_table_heasarc.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## If you only run this once, you can do it in memory in one line:\n", - "## This fetches the FITS as an astropy.io.fits object in memory\n", - "# hdu_list = sdss_table_heasarc[0].getdataobj()\n", - "## But if you might run this notebook repeatedly with limited bandwidth, \n", - "## download it once and cache it. \n", - "\n", - "# Get the filter g version\n", - "file_name=download_file(sdss_table_heasarc[0].getdataurl(), cache=True, timeout=600)\n", - "hdu_list = fits.open(file_name)\n", - "\n", - "plt.imshow(hdu_list[0].data, cmap='gray', origin='lower', vmax=1200, vmin=1010)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Using SDSS SkyServer" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "jhu_dr7_service = [s for s in services if 'SDSSDR7' in s.short_name and 'jhu' in s.ivoid][0]\n", - "\n", - "# Note: jhu_dr7_service access url has hard-wired \"format=image/fits\". \n", - "# If you specify anythign else, it errors. If you specify nothing, \n", - "# then the search() method puts \"format=all\", which errors. So specify \"format=None\" for now.\n", - "sdss_table_jhu=jhu_dr7_service.search(pos=coords,size=0.2, format=None)\n", - "sdss_table_jhu.to_table().show_in_notebook(display_length = 5)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get the filter g version\n", - "file_name=download_file(sdss_table_jhu[1].getdataurl(), cache=True, timeout=600)\n", - "hdu_list = fits.open(file_name)\n", - "\n", - "plt.imshow(hdu_list[0].data, cmap='gray', origin='lower',vmax=1200,vmin=1010)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "It turns out that SkyView is just getting images by using the SIAP internally to get the data from the SDSS service. The point of the VO protocols is that you don't need to know where the data are coming from. But they can be processed differently. " - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.13" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - }, - "widgets": { - "state": {}, - "version": "1.1.1" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/CS_Image_Access.md b/CS_Image_Access.md new file mode 100644 index 0000000..7cc6742 --- /dev/null +++ b/CS_Image_Access.md @@ -0,0 +1,216 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.9.13 +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: {} + toc_section_display: true + toc_window_display: true +widgets: + state: {} + version: 1.1.1 +--- + +# Searching for and retrieving images + +In this notebook, we show how to search for and retrieve images from VO services using the Registry and the __[Simple Image Access](http://www.ivoa.net/documents/SIA/)__ (SIA) protocol. + +- [Searching for and retrieving images](#searching-for-and-retrieving-images) + - [1. Finding SIA resources from the Registry](#1-finding-sia-resources-from-the-registry) + - [2. Using SIA to retrieve an image](#2-using-sia-to-retrieve-an-image) + - [3. Viewing the resulting image](#3-viewing-the-resulting-image) + - [JPG images](#jpg-images) + - [Fits files](#fits-files) + - [4. Example of data available through multiple services](#4-example-of-data-available-through-multiple-services) + - [Using HEASARC](#using-heasarc) + - [Using SDSS SkyServer](#using-sdss-skyserver) + ++++ + +**\*Note:** for all of these notebooks, the results depend on real-time queries. Sometimes there are problems, either because a given service has changed, is undergoing maintenance, or the internet connectivity is having problems, etc. Always retry a couple of times, come back later and try again, and only then send us the problem report to investigate. + +```{code-cell} ipython3 +import warnings + +import numpy as np + +import matplotlib +import matplotlib.pyplot as plt +%matplotlib inline + +import pyvo as vo + +from astropy.io import fits +import astropy.coordinates as coord +# For downloading files +from astropy.utils.data import download_file + +from IPython.display import Image as ipImage, display + +# There are a number of relatively unimportant warnings that show up, so for now, suppress them: +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +``` + +## 1. Finding SIA resources from the Registry + +First, how do we find out what services are available? These are listed in a registry at STScI (__[see here](http://www.ivoa.net/documents/RegTAP/)__). Our Registry function gives a simple interface for how to search for services. + +Let's search for services providing images in the ultraviolet bands: + +```{code-cell} ipython3 +uv_services = vo.regsearch(servicetype='image',waveband='uv') +uv_services.to_table()['ivoid','short_name','res_title'] +``` + +This returns an astropy table containing information about the services available. We can then specify the service we want by using the corresponding row. We'll repeat the search with additional qualifiers to isolate the row we want (note that in the keyword search the "%" character can be used as a wild card): + +```{code-cell} ipython3 +uvot_services = vo.regsearch(servicetype='image',waveband='uv',keywords=['swift']) +uvot_services.to_table()['ivoid','short_name','res_title'] +``` + +This shows us that the data we are interested in comes from the HEASARC's SkyView service, but the point of these VO tools is that you don't need to know that ahead of time or indeed to care where it comes from. + ++++ + +## 2. Using SIA to retrieve an image + +Now we look for images of our favorite source. See __[the SIA definition](http://www.ivoa.net/documents/WD/SIA/sia-20040524.html)__ for usage. In short, you can specify the central position and the size (degrees as one or two floats for the RA, DEC directions). It is up to the service to determine how to provide this. Optionally, you can limit it to the format you want, e.g., "image/fits" or "image/png" etc. + +What is returned to you is not the image itself but a list of images available and how to access them. This is easiest shown by example: + +```{code-cell} ipython3 +coords = coord.SkyCoord.from_name("m51") + +im_table = uvot_services[0].search(pos=coords,size=0.2,format='image/jpeg') +im_table.to_table() +``` + +Extract the fields you're interested in, e.g., the URLs of the images made by skyview. Note that specifying as we did SwiftUVOT, we get a number of different images, e.g., UVOT U, V, B, W1, W2, etc. For each survey, there are two URLs, first the FITS IMAGE and second the JPEG. + +Note that different services will return different column names, but all will have a column giving the URL to access the image. Though it has different column names in different services, it can always be accessed through the `getdataurl` function. + +```{code-cell} ipython3 +url = im_table[0].getdataurl() +print(url) +``` + +## 3. Viewing the resulting image + ++++ + +### JPG images + +Since we have asked for JPEG images, we can display an image in python easily by using its URL. Each row of the result has a getdataurl() method, and you can then hand the URL to an image displayer such as IPython.display: + +```{code-cell} ipython3 +img = ipImage(url=im_table[0].getdataurl()) +display(img) +``` + +### Fits files + +Or download the FITS image and display it with imshow, or aplpy. + +(This often errors off with a time out message. Just try it again, possibly a couple of times.) + +```{code-cell} ipython3 +# Do the search again asking for FITS +im_table = uvot_services[0].search(pos=coords,size=0.2,format='image/fits') + +# Hand the url of the first result to fits.open() +hdu_list = fits.open(im_table[0].getdataurl()) +hdu_list.info() +``` + +#### Using imshow + +```{code-cell} ipython3 +plt.imshow(hdu_list[0].data, cmap='gray', origin='lower',vmax=0.1) +``` + +## 4. Example of data available through multiple services + +Suppose we want Sloan DSS data. A generic query finds us a number of possibilities (note that this doesn't work for keywords=['sdss']; be flexible and try several search terms): + +```{code-cell} ipython3 +services = vo.regsearch(servicetype='image', keywords=['sloan'], waveband='optical') +services.to_table()[np.where(np.isin(services.to_table()['short_name'], 'SDSSDR7'))]['ivoid', 'short_name'] +``` + +So one of these is served by SDSS's SkyServer and the other by HEASARC's SkyView. + ++++ + +### Using HEASARC + +```{code-cell} ipython3 +heasarc_dr7_service = [s for s in services if 'SDSSDR7' in s.short_name and 'heasarc' in s.ivoid][0] + +sdss_table_heasarc = heasarc_dr7_service.search(pos=coords,size=0.2,format='image/fits') +sdss_table_heasarc.to_table() +``` + +```{code-cell} ipython3 +## If you only run this once, you can do it in memory in one line: +## This fetches the FITS as an astropy.io.fits object in memory +# hdu_list = sdss_table_heasarc[0].getdataobj() +## But if you might run this notebook repeatedly with limited bandwidth, +## download it once and cache it. + +# Get the filter g version +file_name=download_file(sdss_table_heasarc[0].getdataurl(), cache=True, timeout=600) +hdu_list = fits.open(file_name) + +plt.imshow(hdu_list[0].data, cmap='gray', origin='lower', vmax=1200, vmin=1010) +``` + +### Using SDSS SkyServer + +```{code-cell} ipython3 +jhu_dr7_service = [s for s in services if 'SDSSDR7' in s.short_name and 'jhu' in s.ivoid][0] + +# Note: jhu_dr7_service access url has hard-wired "format=image/fits". +# If you specify anythign else, it errors. If you specify nothing, +# then the search() method puts "format=all", which errors. So specify "format=None" for now. +sdss_table_jhu=jhu_dr7_service.search(pos=coords,size=0.2, format=None) +sdss_table_jhu.to_table().show_in_notebook(display_length = 5) +``` + +```{code-cell} ipython3 +# Get the filter g version +file_name=download_file(sdss_table_jhu[1].getdataurl(), cache=True, timeout=600) +hdu_list = fits.open(file_name) + +plt.imshow(hdu_list[0].data, cmap='gray', origin='lower',vmax=1200,vmin=1010) +``` + +It turns out that SkyView is just getting images by using the SIAP internally to get the data from the SDSS service. The point of the VO protocols is that you don't need to know where the data are coming from. But they can be processed differently. diff --git a/CS_Spectral_Access.ipynb b/CS_Spectral_Access.ipynb deleted file mode 100644 index 614350f..0000000 --- a/CS_Spectral_Access.ipynb +++ /dev/null @@ -1,211 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Retrieve spectra using Simple Spectral Access protocol\n", - "\n", - "This notebook is one of a set produced by NAVO to demonstrate data access with python tools. \n", - "\n", - "In this notebook, we show how to search for and retrieve spectra from VO services using the Registry and the __[Simple Spectral Access](http://www.ivoa.net/documents/SSA/)__ (SSA) protocol.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "\n", - "import matplotlib\n", - "import matplotlib.pyplot as plt\n", - "%matplotlib inline \n", - "\n", - "import requests, io\n", - "\n", - "from astropy.table import Table\n", - "import astropy.io.fits as fits\n", - "from astropy.coordinates import SkyCoord\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "\n", - "import pyvo as vo\n", - "\n", - "# There are a number of relatively unimportant warnings that show up, so for now, suppress them:\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Finding available Spectral Access Services\n", - "\n", - "First, we find out what spectral access services ('ssa') are available in the Registry offering x-ray data." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='ssa',waveband='x-ray')\n", - "services.to_table()['ivoid','short_name']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can look at only the Chandra entry:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "chandra_service = [s for s in services if 'Chandra' in s.short_name][0] \n", - "chandra_service.access_url" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Chandra Spectrum of Delta Ori\n", - "\n", - "Getting the list of spectra." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "delori = SkyCoord.from_name(\"Delta Ori\")\n", - "\n", - "spec_tables = chandra_service.search(pos=delori,diameter=0.1)\n", - "spec_tables.to_table().show_in_notebook()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Accessing one of the spectra." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## If you only run this once, you can do it in memory in one line:\n", - "## This fetches the FITS as an astropy.io.fits object in memory\n", - "# hdu_list = spec_tables[0].getdataobj()\n", - "## But if you might run this notebook repeatedly with limited bandwidth, \n", - "## download it once and cache it. \n", - "file_name = download_file(spec_tables[0].getdataurl(),cache=True)\n", - "hdu_list = fits.open(file_name)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Simple example of plotting a spectrum" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spec_table = Table(hdu_list[1].data)\n", - "spec_table" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "matplotlib.rcParams['figure.figsize'] = (12, 10)\n", - "\n", - "for i in range(len(spec_table)): \n", - " \n", - " ax = plt.subplot(6,2,i+1)\n", - " pha = plt.plot( spec_table['CHANNEL'][i],spec_table['COUNTS'][i])\n", - " ax.set_yscale('log')\n", - " \n", - " if spec_table['TG_PART'][i] == 1:\n", - " instr='HEG'\n", - " if spec_table['TG_PART'][i] == 2:\n", - " instr='MEG'\n", - " if spec_table['TG_PART'][i] == 3:\n", - " instr='LEG'\n", - " \n", - " ax.set_title(\"{grating}{order:+d}\".format(grating=instr, order=spec_table['TG_M'][i]))\n", - " \n", - " plt.tight_layout()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html). (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.ipynb) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.10" - }, - "nav_menu": {}, - "toc": { - "navigate_menu": true, - "number_sections": true, - "sideBar": true, - "threshold": 6, - "toc_cell": false, - "toc_section_display": "block", - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/CS_Spectral_Access.md b/CS_Spectral_Access.md new file mode 100644 index 0000000..8d323d1 --- /dev/null +++ b/CS_Spectral_Access.md @@ -0,0 +1,135 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.6.10 +nav_menu: {} +toc: + navigate_menu: true + number_sections: true + sideBar: true + threshold: 6 + toc_cell: false + toc_section_display: block + toc_window_display: true +--- + +# Retrieve spectra using Simple Spectral Access protocol + +This notebook is one of a set produced by NAVO to demonstrate data access with python tools. + +In this notebook, we show how to search for and retrieve spectra from VO services using the Registry and the __[Simple Spectral Access](http://www.ivoa.net/documents/SSA/)__ (SSA) protocol. + +```{code-cell} ipython3 +import numpy as np + +import matplotlib +import matplotlib.pyplot as plt +%matplotlib inline + +import requests, io + +from astropy.table import Table +import astropy.io.fits as fits +from astropy.coordinates import SkyCoord +# For downloading files +from astropy.utils.data import download_file + +import pyvo as vo + +# There are a number of relatively unimportant warnings that show up, so for now, suppress them: +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +``` + +## Finding available Spectral Access Services + +First, we find out what spectral access services ('ssa') are available in the Registry offering x-ray data. + +```{code-cell} ipython3 +services = vo.regsearch(servicetype='ssa',waveband='x-ray') +services.to_table()['ivoid','short_name'] +``` + +We can look at only the Chandra entry: + +```{code-cell} ipython3 +chandra_service = [s for s in services if 'Chandra' in s.short_name][0] +chandra_service.access_url +``` + +## Chandra Spectrum of Delta Ori + +Getting the list of spectra. + +```{code-cell} ipython3 +delori = SkyCoord.from_name("Delta Ori") + +spec_tables = chandra_service.search(pos=delori,diameter=0.1) +spec_tables.to_table().show_in_notebook() +``` + +Accessing one of the spectra. + +```{code-cell} ipython3 +## If you only run this once, you can do it in memory in one line: +## This fetches the FITS as an astropy.io.fits object in memory +# hdu_list = spec_tables[0].getdataobj() +## But if you might run this notebook repeatedly with limited bandwidth, +## download it once and cache it. +file_name = download_file(spec_tables[0].getdataurl(),cache=True) +hdu_list = fits.open(file_name) +``` + +## Simple example of plotting a spectrum + +```{code-cell} ipython3 +spec_table = Table(hdu_list[1].data) +spec_table +``` + +```{code-cell} ipython3 +matplotlib.rcParams['figure.figsize'] = (12, 10) + +for i in range(len(spec_table)): + + ax = plt.subplot(6,2,i+1) + pha = plt.plot( spec_table['CHANNEL'][i],spec_table['COUNTS'][i]) + ax.set_yscale('log') + + if spec_table['TG_PART'][i] == 1: + instr='HEG' + if spec_table['TG_PART'][i] == 2: + instr='MEG' + if spec_table['TG_PART'][i] == 3: + instr='LEG' + + ax.set_title("{grating}{order:+d}".format(grating=instr, order=spec_table['TG_M'][i])) + + plt.tight_layout() +``` + +This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html). (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.md) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) + +```{code-cell} ipython3 + +``` diff --git a/CS_UCDs.ipynb b/CS_UCDs.ipynb deleted file mode 100644 index e992f63..0000000 --- a/CS_UCDs.ipynb +++ /dev/null @@ -1,275 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# UCDs: working with heterogeneous tables\n", - "\n", - "Suppose you want to do something using a column that you expect to find in a bunch of different tables, like coordinates and time. It's a good bet that many if not most of the tables have coordinate columns, but there's no rule about what they have to be named. \n", - "\n", - "When doing detailed catalog queries with the TAP, you can obviously examine the columns of every table you're interested in to find the columns you want. Then you can hard-code the correct ones into each query for each table and service. \n", - "\n", - "Or, you can also search for keywords like \"ra\" or \"ascension\" in the columns and their descriptions to get the columns you want automatically that way. \n", - "\n", - "But is there are more generic way? [Unified Content Descriptors (UCDs)](http://www.ivoa.net/documents/latest/UCD.html) are a VO standard that allows table publishers to name their columns whatever they (or their contributors) want but to identify those that contain standard sorts of data. For example, the RA column could be called \"RA\", \"ra\", \"Right_Ascension\", etc. But in all cases, a VO service can label the column with its UCD, which is \"pos.eq.ra\". This information is not part of the table but part of the meta-data that the service may provide with that data. Though not required of all VO services, UCDs are commonly provided precisely to make such tasks as identifying the columns of interest easier to automate. \n", - "\n", - "This is easiest to show by example." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Generic VO access routines\n", - "import pyvo as vo\n", - "\n", - "# Ignore unimportant warnings\n", - "import warnings\n", - "warnings.filterwarnings('ignore', '.*Unknown element .*', vo.utils.xml.elements.UnknownElementWarning)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's look at some tables in a little more detail. Let's find the Hubble Source Catalog version 3 (HSCv3), assuming there's only one at MAST." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='tap', keywords=['mast'])\n", - "hsc=[s for s in services if 'HSCv3' in s.res_title][0]\n", - "\n", - "print(f'Title: {hsc.res_title}')\n", - "print(f'{hsc.res_description}')\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's see what tables are provided by this service for HSCv3. Note that this is another query to the service:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "tables = hsc.service.tables # Queries for details of the service's tables\n", - "print(f'{len(tables)} tables:')\n", - "for t in tables:\n", - " print(f'{t.name:30s} - {t.description}\\n----') # A more succinct option than t.describe()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's look at the first 10 columns of the DetailedCatalog table. Again, note that calling the columns attribute sends another query to the service to ask for the columns." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "columns=tables['dbo.DetailedCatalog'].columns\n", - "for c in columns:\n", - " print(f'{f\"{c.name} [{c.ucd}]\":30s} - {c.description}')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The PyVO method to get the columns will automatically fetch all the meta-data about those columns. It's up to the service provider to set them correctly, of course, but in this case, we see that the column named \"MatchRA\" is identified with the UCD \"pos.eq.ra\". \n", - "\n", - "So if we did not know the exact name used in HSCv3 for the RA, we could do something like this looking for the string \"RA\":" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ra_name=[c.name for c in columns if 'RA' in c.name or \"ascension\" in c.name.lower()]\n", - "print(ra_name)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "But a more general approach is to check for the correct UCD. It also has the further advantage that it can be used to label columns that should be used for certain purposes when there are multiple possibilities. For instance, this table has MatchRA and SourceRA. Let's check the UCD: \n", - "\n", - "(Note that the UCD is not required. If it isn't there, you get a None type, so code the check carefully)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ra_name=[c.name for c in columns if c.ucd and 'pos.eq.ra' in c.ucd][0]\n", - "dec_name=[c.name for c in columns if c.ucd and 'pos.eq.dec' in c.ucd][0]\n", - "ra_name,dec_name" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "What that shows you is that though there are two columns in this table that give RA information, only one has the 'pos.eq.ra' UCD. The documentation for this ought to explain the usage of these columns, and the UCD should not be used as a substitute for understanding the table. But it can be a useful tool. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In particular, you can use the UCDs to look for catalogs that might have the information you're interested in. Then you can code the same query to work for different tables (with different column names) in a loop. This sends a bunch of queries but doesn't take too long, a minute maybe. (One is particularly slow.) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Look for all TAP services with x-ray and optical data\n", - "collection={}\n", - "for s in vo.regsearch(servicetype='tap',keywords=['x-ray','optical']):\n", - " if \"wfau\" in s.ivoid: continue # These sometimes have issues\n", - " print(f\"Looking at service from {s.ivoid}\")\n", - " try:\n", - " tables=s.service.tables\n", - " except:\n", - " print(\"Problem with this service's tables endpoint. Continuing to next.\")\n", - " continue\n", - " # Find all the tables that have an RA,DEC and a start and end time\n", - " for t in tables:\n", - " names={}\n", - " for ucd in ['pos.eq.ra','pos.eq.dec','time.start','time.end']:\n", - " cols=[c.name for c in t.columns if c.ucd and ucd in c.ucd]\n", - " if len(cols) > 0: \n", - " names[ucd]=cols[0] # use the first that matches\n", - " if len(names.keys()) == 4: \n", - " print(f\" Table {t.name} has the right columns. Counting rows matching my time.\")\n", - " # For a first look, a very simple query counting rows in a \n", - " # time range of interest:\n", - " query=f\"select count({names['time.start']}) from {t.name}\" \\\n", - " f\" where {names['time.start']} > 52000 \" \n", - " try:\n", - " results=s.search(query)\n", - " except:\n", - " print(\"Problem executing query. Continuing to next.\")\n", - " continue\n", - " # For this simple query, the result is a single number, the count. \n", - " # But different services might name the result differently, so \n", - " # don't assume you know the column name. \n", - " print(\" Found {} results from {}\\n\".format(results.to_table()[0][0],t.name))\n", - " # If the query above asked for the matching data rather than the\n", - " # count, you might want to collect the results. \n", - " # Careful: here we're assuming the table names are unique\n", - " collection[t.name]=results\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can also use UCDs to look at the results. Above, we collected just the first 10 rows of the four columns we're interested in from every catalog that had them. But these tables still have their original column names. So the UCDs will still be useful, and PyVO provides a simple routine to convert from UCD to column (field) name. \n", - "\n", - "Note, however, that returning the UCDs as part of the result is not mandatory, and some services do not do it. So you'll have to check.\n", - "\n", - "Now we have a collection of rows from different tables with different columns. In the results object, we have access to a fieldname_with_ucd() function to get the column you want. Supposing we hadn't already looked for this in the above loop, let's now find out which of these tables has a magnitude column:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#ucd='pos.eq.ra'\n", - "ucd='phot.mag'\n", - "for tname,results in collection.items():\n", - " #print(f\"On table {tname}\")\n", - " # Sometimes this doesn't work well, so use a try:\n", - " try:\n", - " name=results.fieldname_with_ucd(ucd)\n", - " except:\n", - " pass\n", - " if name:\n", - " print(f\" Table {tname} has the {ucd} column named {name}\")\n", - " else:\n", - " print(f\" (Table {tname} didn't find the UCD.)\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Lastly, if you have a table of results from a TAP query (and if that service includes the UCDs), then you can get data based on UCDs with the getbyucd() method, which simply gets the corresponding element using fieldname_with_ucd():" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "results=hsc.service.search(\"select top 10 * from dbo.DetailedCatalog\")\n", - "[r.getbyucd('phot.mag') for r in results]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that we can see earlier in this notebook, when we looked at this table's contents, that there are two phot.mag fields in this table, MagAper2 and MagAuto. The getbyucd() and fieldname_with_ucd() routines do not currently allow you to handle multiple columns with the same UCD. The code can help you find what you want, but it depends on the meta data the service defines, and you still must look at the detailed information for each catalog you use to understand what it contains. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.10" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} diff --git a/CS_UCDs.md b/CS_UCDs.md new file mode 100644 index 0000000..f1fc249 --- /dev/null +++ b/CS_UCDs.md @@ -0,0 +1,170 @@ +--- +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.7.10 +--- + +# UCDs: working with heterogeneous tables + +Suppose you want to do something using a column that you expect to find in a bunch of different tables, like coordinates and time. It's a good bet that many if not most of the tables have coordinate columns, but there's no rule about what they have to be named. + +When doing detailed catalog queries with the TAP, you can obviously examine the columns of every table you're interested in to find the columns you want. Then you can hard-code the correct ones into each query for each table and service. + +Or, you can also search for keywords like "ra" or "ascension" in the columns and their descriptions to get the columns you want automatically that way. + +But is there are more generic way? [Unified Content Descriptors (UCDs)](http://www.ivoa.net/documents/latest/UCD.html) are a VO standard that allows table publishers to name their columns whatever they (or their contributors) want but to identify those that contain standard sorts of data. For example, the RA column could be called "RA", "ra", "Right_Ascension", etc. But in all cases, a VO service can label the column with its UCD, which is "pos.eq.ra". This information is not part of the table but part of the meta-data that the service may provide with that data. Though not required of all VO services, UCDs are commonly provided precisely to make such tasks as identifying the columns of interest easier to automate. + +This is easiest to show by example. + +```{code-cell} ipython3 +# Generic VO access routines +import pyvo as vo + +# Ignore unimportant warnings +import warnings +warnings.filterwarnings('ignore', '.*Unknown element .*', vo.utils.xml.elements.UnknownElementWarning) +``` + +Let's look at some tables in a little more detail. Let's find the Hubble Source Catalog version 3 (HSCv3), assuming there's only one at MAST. + +```{code-cell} ipython3 +services = vo.regsearch(servicetype='tap', keywords=['mast']) +hsc=[s for s in services if 'HSCv3' in s.res_title][0] + +print(f'Title: {hsc.res_title}') +print(f'{hsc.res_description}') +``` + +Now let's see what tables are provided by this service for HSCv3. Note that this is another query to the service: + +```{code-cell} ipython3 +tables = hsc.service.tables # Queries for details of the service's tables +print(f'{len(tables)} tables:') +for t in tables: + print(f'{t.name:30s} - {t.description}\n----') # A more succinct option than t.describe() +``` + +Let's look at the first 10 columns of the DetailedCatalog table. Again, note that calling the columns attribute sends another query to the service to ask for the columns. + +```{code-cell} ipython3 +columns=tables['dbo.DetailedCatalog'].columns +for c in columns: + print(f'{f"{c.name} [{c.ucd}]":30s} - {c.description}') +``` + +The PyVO method to get the columns will automatically fetch all the meta-data about those columns. It's up to the service provider to set them correctly, of course, but in this case, we see that the column named "MatchRA" is identified with the UCD "pos.eq.ra". + +So if we did not know the exact name used in HSCv3 for the RA, we could do something like this looking for the string "RA": + +```{code-cell} ipython3 +ra_name=[c.name for c in columns if 'RA' in c.name or "ascension" in c.name.lower()] +print(ra_name) +``` + +But a more general approach is to check for the correct UCD. It also has the further advantage that it can be used to label columns that should be used for certain purposes when there are multiple possibilities. For instance, this table has MatchRA and SourceRA. Let's check the UCD: + +(Note that the UCD is not required. If it isn't there, you get a None type, so code the check carefully) + +```{code-cell} ipython3 +ra_name=[c.name for c in columns if c.ucd and 'pos.eq.ra' in c.ucd][0] +dec_name=[c.name for c in columns if c.ucd and 'pos.eq.dec' in c.ucd][0] +ra_name,dec_name +``` + +What that shows you is that though there are two columns in this table that give RA information, only one has the 'pos.eq.ra' UCD. The documentation for this ought to explain the usage of these columns, and the UCD should not be used as a substitute for understanding the table. But it can be a useful tool. + ++++ + +In particular, you can use the UCDs to look for catalogs that might have the information you're interested in. Then you can code the same query to work for different tables (with different column names) in a loop. This sends a bunch of queries but doesn't take too long, a minute maybe. (One is particularly slow.) + +```{code-cell} ipython3 +# Look for all TAP services with x-ray and optical data +collection={} +for s in vo.regsearch(servicetype='tap',keywords=['x-ray','optical']): + if "wfau" in s.ivoid: continue # These sometimes have issues + print(f"Looking at service from {s.ivoid}") + try: + tables=s.service.tables + except: + print("Problem with this service's tables endpoint. Continuing to next.") + continue + # Find all the tables that have an RA,DEC and a start and end time + for t in tables: + names={} + for ucd in ['pos.eq.ra','pos.eq.dec','time.start','time.end']: + cols=[c.name for c in t.columns if c.ucd and ucd in c.ucd] + if len(cols) > 0: + names[ucd]=cols[0] # use the first that matches + if len(names.keys()) == 4: + print(f" Table {t.name} has the right columns. Counting rows matching my time.") + # For a first look, a very simple query counting rows in a + # time range of interest: + query=f"select count({names['time.start']}) from {t.name}" \ + f" where {names['time.start']} > 52000 " + try: + results=s.search(query) + except: + print("Problem executing query. Continuing to next.") + continue + # For this simple query, the result is a single number, the count. + # But different services might name the result differently, so + # don't assume you know the column name. + print(" Found {} results from {}\n".format(results.to_table()[0][0],t.name)) + # If the query above asked for the matching data rather than the + # count, you might want to collect the results. + # Careful: here we're assuming the table names are unique + collection[t.name]=results +``` + +You can also use UCDs to look at the results. Above, we collected just the first 10 rows of the four columns we're interested in from every catalog that had them. But these tables still have their original column names. So the UCDs will still be useful, and PyVO provides a simple routine to convert from UCD to column (field) name. + +Note, however, that returning the UCDs as part of the result is not mandatory, and some services do not do it. So you'll have to check. + +Now we have a collection of rows from different tables with different columns. In the results object, we have access to a fieldname_with_ucd() function to get the column you want. Supposing we hadn't already looked for this in the above loop, let's now find out which of these tables has a magnitude column: + +```{code-cell} ipython3 +#ucd='pos.eq.ra' +ucd='phot.mag' +for tname,results in collection.items(): + #print(f"On table {tname}") + # Sometimes this doesn't work well, so use a try: + try: + name=results.fieldname_with_ucd(ucd) + except: + pass + if name: + print(f" Table {tname} has the {ucd} column named {name}") + else: + print(f" (Table {tname} didn't find the UCD.)") +``` + +Lastly, if you have a table of results from a TAP query (and if that service includes the UCDs), then you can get data based on UCDs with the getbyucd() method, which simply gets the corresponding element using fieldname_with_ucd(): + +```{code-cell} ipython3 +results=hsc.service.search("select top 10 * from dbo.DetailedCatalog") +[r.getbyucd('phot.mag') for r in results] +``` + +Note that we can see earlier in this notebook, when we looked at this table's contents, that there are two phot.mag fields in this table, MagAper2 and MagAuto. The getbyucd() and fieldname_with_ucd() routines do not currently allow you to handle multiple columns with the same UCD. The code can help you find what you want, but it depends on the meta data the service defines, and you still must look at the detailed information for each catalog you use to understand what it contains. + +```{code-cell} ipython3 + +``` diff --git a/CS_VO_Tables.ipynb b/CS_VO_Tables.ipynb deleted file mode 100644 index a14412f..0000000 --- a/CS_VO_Tables.ipynb +++ /dev/null @@ -1,145 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Creating a VO Table from a CSV file " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "There are several ways of doing this, and there are a few object layers here, which can be confusing:\n", - "- Standard [astropy Table](https://docs.astropy.org/en/stable/table/) objects\n", - "- [Votable Table](https://docs.astropy.org/en/stable/api/astropy.io.votable.tree.Table.html#astropy.io.votable.tree.Table) objects\n", - "- [Votable VOTableFile](https://docs.astropy.org/en/stable/api/astropy.io.votable.tree.VOTableFile.html#astropy.io.votable.tree.VOTableFile) objects (may contain multiple votable Tables)\n", - "\n", - "Although some things can be done with generic astropy Tables, other VO operations can only be done with VO Tables or VOTableFile objects. \n", - "\n", - "This is easiest to see with an example. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import io, numpy\n", - "\n", - "from astropy import table as aptable\n", - "from astropy.io import votable as apvot" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Create a table with only two columns starting from an astropy Table" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "myaptable=aptable.Table(\n", - " numpy.array([\n", - " [19.0186, 46.7304],\n", - " [20.2887, 40.4703],\n", - " [125.886, 21.3377],\n", - " [136.002, 21.9679],\n", - " [141.057, 40.6372],\n", - " [146.700, 22.0116],\n", - " [148.785, 14.2922],\n", - " [149.751, 17.8168],\n", - " [175.039, 15.3270],\n", - " [191.542, 30.7317],\n", - " [194.913, 28.8959],\n", - " [199.026, 41.5011],\n", - " [206.577, 43.8511],\n", - " [209.963, 38.1821],\n", - " [213.556, 15.6214],\n", - " [219.967, 42.7421],\n", - " [226.693, 12.8502],\n", - " [237.489, 20.8057],\n", - " [241.519, 20.8014],\n", - " [317.088, 18.2002],\n", - " [329.235, 6.64845],\n", - " [333.830, 37.3012] ]), \n", - " names=[\"RA\",\"DEC\"])\n", - "\n", - "print(type(myaptable))\n", - "print(myaptable)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Then convert this to a VOTableFile object which contains a nested set of *resources* and *tables* (in this case, only one of each)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "myvotablefile = apvot.from_table(myaptable)\n", - "print(type(myvotablefile))\n", - "\n", - "for r in myvotablefile.resources:\n", - " mytable=r\n", - " for t in r.tables:\n", - " print(t)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.1" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 1 -} \ No newline at end of file diff --git a/CS_VO_Tables.md b/CS_VO_Tables.md new file mode 100644 index 0000000..4c858c5 --- /dev/null +++ b/CS_VO_Tables.md @@ -0,0 +1,105 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.9.1 +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: {} + toc_section_display: true + toc_window_display: true +--- + +# Creating a VO Table from a CSV file + ++++ + +There are several ways of doing this, and there are a few object layers here, which can be confusing: +- Standard [astropy Table](https://docs.astropy.org/en/stable/table/) objects +- [Votable Table](https://docs.astropy.org/en/stable/api/astropy.io.votable.tree.Table.html#astropy.io.votable.tree.Table) objects +- [Votable VOTableFile](https://docs.astropy.org/en/stable/api/astropy.io.votable.tree.VOTableFile.html#astropy.io.votable.tree.VOTableFile) objects (may contain multiple votable Tables) + +Although some things can be done with generic astropy Tables, other VO operations can only be done with VO Tables or VOTableFile objects. + +This is easiest to see with an example. + +```{code-cell} ipython3 +import io, numpy + +from astropy import table as aptable +from astropy.io import votable as apvot +``` + +## Create a table with only two columns starting from an astropy Table + +```{code-cell} ipython3 +myaptable=aptable.Table( + numpy.array([ + [19.0186, 46.7304], + [20.2887, 40.4703], + [125.886, 21.3377], + [136.002, 21.9679], + [141.057, 40.6372], + [146.700, 22.0116], + [148.785, 14.2922], + [149.751, 17.8168], + [175.039, 15.3270], + [191.542, 30.7317], + [194.913, 28.8959], + [199.026, 41.5011], + [206.577, 43.8511], + [209.963, 38.1821], + [213.556, 15.6214], + [219.967, 42.7421], + [226.693, 12.8502], + [237.489, 20.8057], + [241.519, 20.8014], + [317.088, 18.2002], + [329.235, 6.64845], + [333.830, 37.3012] ]), + names=["RA","DEC"]) + +print(type(myaptable)) +print(myaptable) +``` + +## Then convert this to a VOTableFile object which contains a nested set of *resources* and *tables* (in this case, only one of each) + +```{code-cell} ipython3 +myvotablefile = apvot.from_table(myaptable) +print(type(myvotablefile)) + +for r in myvotablefile.resources: + mytable=r + for t in r.tables: + print(t) +``` + +```{code-cell} ipython3 + +``` diff --git a/Exercise_I.ipynb b/Exercise_I.ipynb deleted file mode 100644 index 3012731..0000000 --- a/Exercise_I.ipynb +++ /dev/null @@ -1,434 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Science User Case - Inspecting a Candidate List\n", - "\n", - "Ogle et al. (2016) mined the NASA/IPAC Extragalactic Database (NED) to identify a new type of galaxy: Superluminous Spiral Galaxies.\n", - "\n", - "Here's the paper: https://ui.adsabs.harvard.edu/abs/2016ApJ...817..109O/abstract\n", - "\n", - "Table 1 lists the positions of these Super Spirals. Based on those positions, let's create multiwavelength cutouts for each super spiral to see what is unique about this new class of objects." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 1. Import the Python modules we'll be using." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Suppress unimportant warnings.\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")\n", - "warnings.filterwarnings('ignore', '.*RADECSYS=*', append=True)\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "\n", - "from astropy.coordinates import SkyCoord\n", - "from astropy.io import fits\n", - "from astropy.nddata import Cutout2D\n", - "import astropy.visualization as vis\n", - "from astropy.wcs import WCS\n", - "from astroquery.ipac.ned import Ned\n", - "\n", - "import pyvo as vo" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The next cell prepares the notebook to display our visualizations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2. Search NED for objects in this paper.\n", - "\n", - "Insert a Code Cell below by clicking on the \"Insert\" Menu and choosing \"Insert Cell Below\". Then consult QuickReference.md to figure out how to use astroquery to search NED for all objects in a paper, based on the refcode of the paper. Inspect the resulting astropy table." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: the QuickReference has this example:\n", - "#objects_in_paper = Ned.query_refcode('2018ApJ...858...62K')\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 3. Filter the NED results.\n", - "\n", - "The results from NED will include galaxies, but also other kinds of objects (e.g. galaxy clusters, galaxy groups). Print the 'Type' column to see the full range of classifications and filter the results so that we only keep the galaxies in the list." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 4. Search the NAVO Registry for image resources.\n", - "\n", - "The paper selected super spirals using WISE, SDSS, and GALEX images. Search the NAVO registry for all image resources, using the 'service_type' search parameter. How many image resources are currently available?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: the QuickReference has this example:\n", - "# services = vo.regsearch(servicetype='conesearch', keywords=['swift'])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 5. Search the NAVO Registry for image resources that will allow you to search for AllWISE images.\n", - "\n", - "There are hundreds of image resources...too many to quickly read through. Try adding the 'keywords' search parameter to your registry search, and find the image resource you would need to search the AllWISE images. Remember from the Known Issues that 'keywords' must be a list." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: the QuickReference has this example:\n", - "# services = vo.regsearch(servicetype='conesearch', keywords=['swift'])" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 6. Choose the AllWISE image service that you are interested in." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 7. Choose one of the galaxies in the NED list.\n", - "What is the position of this galaxy?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: the QuickReference has this example:\n", - "# m83_pos = SkyCoord('13h37m00.950s -29d51m55.51s')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 8. Search for a list of AllWISE images that cover this galaxy.\n", - "\n", - "How many images are returned? Which are you most interested in?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: the QuickReference has this example:\n", - "#results = services[1].search(pos=m83_pos, size=.2)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 9. Use the .to_table() method to view the results as an Astropy table." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 10. From the result in 8., select the first record for an image taken in WISE band W1 (3.6 micron)\n", - "\n", - "Hints:\n", - "* Loop over records and test on the `.bandpass_id` attribute of each record\n", - "* Print the `.title` and `.bandpass_id` of the record you find, to verify it is the right one." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 11. Visualize this AllWISE image.\n", - "Hint: Locate the galaxy in the image by overplotting a ring centered on the galaxy on the image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: the QuickReference has this example:\n", - "# file_name = download_file(results[0].getdataurl()) \n", - "# Hint: when displaying the fits file play with the vmax value, a \n", - "# vmax value around 10 will display a nice image" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 12. Plot a cutout of the AllWISE image, centered on your position.\n", - "\n", - "Try a 60 arcsecond cutout." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: using Cutout2D imported above from from astropy.nddata\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 13. Try visualizing a cutout of a GALEX image that covers your position.\n", - "\n", - "Repeat steps 5, 6, 8 through 12 for GALEX." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# Steps 5 & 6: Search the services found in step 4 for GALEX Sources\n", - "# Hint: Choose the GALEX service on STSCI" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Steps 8 & 9: Find the images that contain your chosen galaxy\n", - "# and display using .to_table()\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Step 10: Select one of the images and print the title and relevent \n", - "# info. For example you can select an image with an 'enrValue' of 2.35e-07\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Step 11: Visualize the image and overplot a ring on the image centered \n", - "# on your galaxy. (A very small vmax value around 0.01 is recomended)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Step 12: Use Cutout2D to get a 60 arcsecond cutout\n", - "# centered on the galaxy" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 14. Try visualizing a cutout of an SDSS image that covers your position.\n", - "\n", - "Hints:\n", - "* Search the registry using `keywords=['sloan']\n", - "* Find the service with a `short_name` of `'SDSS SIAP'`\n", - "* From Known Issues, recall that an empty string must be specified to the `format` parameter dues to a bug in the service.\n", - "* After obtaining your search results, select r-band images using the `.title` attribute of the records that are returned, since `.bandpass_id` is not populated." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "(As a workaround to a bug in the SDSS service, pass `format=''` as an argument to the search() function when using the SDSS service.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Search the registery to find a service" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Search the SDSS SIAP service for images containing your galaxy" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Find an image taken in the r filter" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Visualize the image and overplot a ring centered on the galaxy\n", - "# Hint: a vmax around 0.01 will produce a nice image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# use Cutout2D to get a 60 arcsecond cutout centered on the galaxy" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 15. Try looping over all positions and plotting multiwavelength cutouts.\n", - "Hint: Gather the data in ALLWISE, GALEX, and SDSS for each galaxy and plot as seperate cutouts centered on the galaxy position" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Warning: this takes a long time to run! You can limit it to the first three galaxies only, for example, for testing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.2" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/Exercise_I.md b/Exercise_I.md new file mode 100644 index 0000000..f35033d --- /dev/null +++ b/Exercise_I.md @@ -0,0 +1,231 @@ +--- +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.8.2 +--- + +# Science User Case - Inspecting a Candidate List + +Ogle et al. (2016) mined the NASA/IPAC Extragalactic Database (NED) to identify a new type of galaxy: Superluminous Spiral Galaxies. + +Here's the paper: https://ui.adsabs.harvard.edu/abs/2016ApJ...817..109O/abstract + +Table 1 lists the positions of these Super Spirals. Based on those positions, let's create multiwavelength cutouts for each super spiral to see what is unique about this new class of objects. + ++++ + +## 1. Import the Python modules we'll be using. + +```{code-cell} ipython3 +# Suppress unimportant warnings. +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +warnings.filterwarnings('ignore', '.*RADECSYS=*', append=True) + +import matplotlib.pyplot as plt +import numpy as np + +# For downloading files +from astropy.utils.data import download_file + +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.nddata import Cutout2D +import astropy.visualization as vis +from astropy.wcs import WCS +from astroquery.ipac.ned import Ned + +import pyvo as vo +``` + +The next cell prepares the notebook to display our visualizations. + +```{code-cell} ipython3 +%matplotlib inline +``` + +## 2. Search NED for objects in this paper. + +Insert a Code Cell below by clicking on the "Insert" Menu and choosing "Insert Cell Below". Then consult QuickReference.md to figure out how to use astroquery to search NED for all objects in a paper, based on the refcode of the paper. Inspect the resulting astropy table. + +```{code-cell} ipython3 +# Hint: the QuickReference has this example: +#objects_in_paper = Ned.query_refcode('2018ApJ...858...62K') +``` + +## 3. Filter the NED results. + +The results from NED will include galaxies, but also other kinds of objects (e.g. galaxy clusters, galaxy groups). Print the 'Type' column to see the full range of classifications and filter the results so that we only keep the galaxies in the list. + +```{code-cell} ipython3 + +``` + +## 4. Search the NAVO Registry for image resources. + +The paper selected super spirals using WISE, SDSS, and GALEX images. Search the NAVO registry for all image resources, using the 'service_type' search parameter. How many image resources are currently available? + +```{code-cell} ipython3 +# Hint: the QuickReference has this example: +# services = vo.regsearch(servicetype='conesearch', keywords=['swift']) +``` + +## 5. Search the NAVO Registry for image resources that will allow you to search for AllWISE images. + +There are hundreds of image resources...too many to quickly read through. Try adding the 'keywords' search parameter to your registry search, and find the image resource you would need to search the AllWISE images. Remember from the Known Issues that 'keywords' must be a list. + +```{code-cell} ipython3 +# Hint: the QuickReference has this example: +# services = vo.regsearch(servicetype='conesearch', keywords=['swift']) +``` + +## 6. Choose the AllWISE image service that you are interested in. + +```{code-cell} ipython3 + +``` + +## 7. Choose one of the galaxies in the NED list. +What is the position of this galaxy? + +```{code-cell} ipython3 +# Hint: the QuickReference has this example: +# m83_pos = SkyCoord('13h37m00.950s -29d51m55.51s') +``` + +## 8. Search for a list of AllWISE images that cover this galaxy. + +How many images are returned? Which are you most interested in? + +```{code-cell} ipython3 +# Hint: the QuickReference has this example: +#results = services[1].search(pos=m83_pos, size=.2) +``` + +## 9. Use the .to_table() method to view the results as an Astropy table. + +```{code-cell} ipython3 + +``` + +## 10. From the result in 8., select the first record for an image taken in WISE band W1 (3.6 micron) + +Hints: +* Loop over records and test on the `.bandpass_id` attribute of each record +* Print the `.title` and `.bandpass_id` of the record you find, to verify it is the right one. + +```{code-cell} ipython3 + +``` + +## 11. Visualize this AllWISE image. +Hint: Locate the galaxy in the image by overplotting a ring centered on the galaxy on the image + +```{code-cell} ipython3 +# Hint: the QuickReference has this example: +# file_name = download_file(results[0].getdataurl()) +# Hint: when displaying the fits file play with the vmax value, a +# vmax value around 10 will display a nice image +``` + +## 12. Plot a cutout of the AllWISE image, centered on your position. + +Try a 60 arcsecond cutout. + +```{code-cell} ipython3 +# Hint: using Cutout2D imported above from from astropy.nddata +``` + +## 13. Try visualizing a cutout of a GALEX image that covers your position. + +Repeat steps 5, 6, 8 through 12 for GALEX. + +```{code-cell} ipython3 +:tags: [output_scroll] + +# Steps 5 & 6: Search the services found in step 4 for GALEX Sources +# Hint: Choose the GALEX service on STSCI +``` + +```{code-cell} ipython3 +# Steps 8 & 9: Find the images that contain your chosen galaxy +# and display using .to_table() +``` + +```{code-cell} ipython3 +# Step 10: Select one of the images and print the title and relevent +# info. For example you can select an image with an 'enrValue' of 2.35e-07 +``` + +```{code-cell} ipython3 +# Step 11: Visualize the image and overplot a ring on the image centered +# on your galaxy. (A very small vmax value around 0.01 is recomended) +``` + +```{code-cell} ipython3 +# Step 12: Use Cutout2D to get a 60 arcsecond cutout +# centered on the galaxy +``` + +## 14. Try visualizing a cutout of an SDSS image that covers your position. + +Hints: +* Search the registry using `keywords=['sloan'] +* Find the service with a `short_name` of `'SDSS SIAP'` +* From Known Issues, recall that an empty string must be specified to the `format` parameter dues to a bug in the service. +* After obtaining your search results, select r-band images using the `.title` attribute of the records that are returned, since `.bandpass_id` is not populated. + ++++ + +(As a workaround to a bug in the SDSS service, pass `format=''` as an argument to the search() function when using the SDSS service.) + +```{code-cell} ipython3 +# Search the registery to find a service +``` + +```{code-cell} ipython3 +# Search the SDSS SIAP service for images containing your galaxy +``` + +```{code-cell} ipython3 +# Find an image taken in the r filter +``` + +```{code-cell} ipython3 +# Visualize the image and overplot a ring centered on the galaxy +# Hint: a vmax around 0.01 will produce a nice image +``` + +```{code-cell} ipython3 +# use Cutout2D to get a 60 arcsecond cutout centered on the galaxy +``` + +## 15. Try looping over all positions and plotting multiwavelength cutouts. +Hint: Gather the data in ALLWISE, GALEX, and SDSS for each galaxy and plot as seperate cutouts centered on the galaxy position + ++++ + +Warning: this takes a long time to run! You can limit it to the first three galaxies only, for example, for testing. + +```{code-cell} ipython3 + +``` diff --git a/Exercise_II.ipynb b/Exercise_II.ipynb deleted file mode 100644 index cfd9f6a..0000000 --- a/Exercise_II.ipynb +++ /dev/null @@ -1,323 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Preparing a proposal\n", - "\n", - "The Story: Suppose that you are preparing to write a proposal on NGC1365, aiming to investigate the intriguing black hole spin this galaxy with Chandra grating observations (see: https://www.space.com/19980-monster-black-hole-spin-discovery.html ) \n", - "\n", - "In writing proposals, there are often the same tasks that are required: including finding and analyzing previous observations of the proposal, and creating figures that include, e.g., multiwavelength images and spectrum for the source. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# As a hint, we include the code block for Python modules that you will likely need to import: \n", - "import matplotlib\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "%matplotlib inline \n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "from astropy.io import fits\n", - "\n", - "import pyvo as vo\n", - "\n", - "## There are a number of relatively unimportant warnings that \n", - "## show up, so for now, suppress them:\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 1: Find out what the previously quoted Chandra 2-10 keV flux of the central source is for NGC 1365. \n", - "\n", - "Hint: Do a Registry search for tables served by the HEASARC (where high energy data are archived) to find potential table with this information" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Start with a Registry query to find table services for the HEASARC.\n", - "# Then get the list of tables that this service serves.\n", - "# \n", - "# Hint: the QuickReference has this example:\n", - "#services = vo.regsearch(servicetype='tap', keywords=['heasarc'])\n", - "#tables = services[0].service.tables \n", - "#\n", - "# Hint2: the QuickReference also has this example:\n", - "#for c in tables['zcat'].columns:\n", - "# print(f'{c.name:30s} - {c.description}')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint: The Chansngcat ( https://heasarc.gsfc.nasa.gov/W3Browse/chandra/chansngcat.html ) table is likely the best table. Create a table with ra, dec, exposure time, and flux (and flux errors) from the public.chansngcat catalog for Chandra observations matched within 0.1 degree." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get the coordinate for NGC 1365 with astropy. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Construct a query that will get the ra, dec, exposure time, flux, and flux errors \n", - "# from this catalog in the region around this source and submit the query. \n", - "# (See the CS_Catalog_queries.ipynb )\n", - "# Hint: the QuickReference has this example:\n", - "#coord = SkyCoord.from_name(\"m83\")\n", - "#query = f'''\n", - "#SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where \n", - "#contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{coord.ra.deg},{coord.dec.deg},1.0))=1\n", - "#'''\n", - "#results = services[0].service.run_async(query)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "## Step 2: Make Images \n", - "\n", - "### Create ultraviolet and X-ray images\n", - "Hint: Start by checking what UV image services exist (e.g., GALEX?)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: start with a Registry search for relevant image services" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The keyword search for 'galex' returned a bunch of things that may have mentioned it, but let's just use the ones that have GALEX as their short name:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Filter on the short_name attribute of the list of services" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Though using the result as an Astropy Table makes it easier to look at the contents, to call the service itself, we cannot use the row of that table. You have to use the entry in the service result list itself. So use the table to browse, but select the list of services itself using the properties that have been defined as attributes such as short_name and ivoid:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# You may find more than one service. Look at both. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint: Next create a UV image for the source " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Do an image search for NGC 1365 in the UV services found above\n", - "#\n", - "# Hint: the QuickReference has this example:\n", - "#results = services[1].search(pos=m83_pos, size=.2)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get a FITS file and visualize the image \n", - "#\n", - "# Hint: the QuickReference has this example:\n", - "#file_name = download_file(results[0].getdataurl()) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint: Repeat steps for X-ray image. (Note: Ideally, we would find an image in the Chandra 'cxc' catalog) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "## Step 3: Make a spectrum \n", - "\n", - "### Find what Chandra spectral observations exist already for this source. \n", - "Hint: try searching for X-ray spectral data tables using the registry query" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Search the Registry to list services that contain X-ray spectral data" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint 2: Take a look at what data exist for our candidate, NGC 1365." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint 3: Download the data to make a spectrum. Note: you might end here and use Xspec to plot and model the spectrum. Or ... you can also try to take a quick look at the spectrum. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get it and look at it:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Or write it to disk" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Extension: Making a \"quick look\" spectrum. For our purposes, the 1st order of the HEG grating data would be sufficient." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Hint: You'll have to look into the details of the spectra. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html). (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.ipynb) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Congratulations! You have completed this notebook exercise." - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.1" - }, - "nav_menu": {}, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": "block", - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/Exercise_II.md b/Exercise_II.md new file mode 100644 index 0000000..a19ef2a --- /dev/null +++ b/Exercise_II.md @@ -0,0 +1,179 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.7.1 +nav_menu: {} +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: {} + toc_section_display: block + toc_window_display: true +--- + +# Preparing a proposal + +The Story: Suppose that you are preparing to write a proposal on NGC1365, aiming to investigate the intriguing black hole spin this galaxy with Chandra grating observations (see: https://www.space.com/19980-monster-black-hole-spin-discovery.html ) + +In writing proposals, there are often the same tasks that are required: including finding and analyzing previous observations of the proposal, and creating figures that include, e.g., multiwavelength images and spectrum for the source. + +```{code-cell} ipython3 +# As a hint, we include the code block for Python modules that you will likely need to import: +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +%matplotlib inline + +# For downloading files +from astropy.utils.data import download_file +from astropy.io import fits + +import pyvo as vo + +## There are a number of relatively unimportant warnings that +## show up, so for now, suppress them: +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +``` + +## Step 1: Find out what the previously quoted Chandra 2-10 keV flux of the central source is for NGC 1365. + +Hint: Do a Registry search for tables served by the HEASARC (where high energy data are archived) to find potential table with this information + +```{code-cell} ipython3 +# Start with a Registry query to find table services for the HEASARC. +# Then get the list of tables that this service serves. +# +# Hint: the QuickReference has this example: +#services = vo.regsearch(servicetype='tap', keywords=['heasarc']) +#tables = services[0].service.tables +# +# Hint2: the QuickReference also has this example: +#for c in tables['zcat'].columns: +# print(f'{c.name:30s} - {c.description}') +``` + +Hint: The Chansngcat ( https://heasarc.gsfc.nasa.gov/W3Browse/chandra/chansngcat.html ) table is likely the best table. Create a table with ra, dec, exposure time, and flux (and flux errors) from the public.chansngcat catalog for Chandra observations matched within 0.1 degree. + +```{code-cell} ipython3 +# Get the coordinate for NGC 1365 with astropy. +``` + +```{code-cell} ipython3 +# Construct a query that will get the ra, dec, exposure time, flux, and flux errors +# from this catalog in the region around this source and submit the query. +# (See the CS_Catalog_queries.md ) +# Hint: the QuickReference has this example: +#coord = SkyCoord.from_name("m83") +#query = f''' +#SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where +#contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{coord.ra.deg},{coord.dec.deg},1.0))=1 +#''' +#results = services[0].service.run_async(query) +``` + +## Step 2: Make Images + +### Create ultraviolet and X-ray images +Hint: Start by checking what UV image services exist (e.g., GALEX?) + +```{code-cell} ipython3 +# Hint: start with a Registry search for relevant image services +``` + +The keyword search for 'galex' returned a bunch of things that may have mentioned it, but let's just use the ones that have GALEX as their short name: + +```{code-cell} ipython3 +# Filter on the short_name attribute of the list of services +``` + +Though using the result as an Astropy Table makes it easier to look at the contents, to call the service itself, we cannot use the row of that table. You have to use the entry in the service result list itself. So use the table to browse, but select the list of services itself using the properties that have been defined as attributes such as short_name and ivoid: + +```{code-cell} ipython3 +# You may find more than one service. Look at both. +``` + +Hint: Next create a UV image for the source + +```{code-cell} ipython3 +# Do an image search for NGC 1365 in the UV services found above +# +# Hint: the QuickReference has this example: +#results = services[1].search(pos=m83_pos, size=.2) +``` + +```{code-cell} ipython3 +# Get a FITS file and visualize the image +# +# Hint: the QuickReference has this example: +#file_name = download_file(results[0].getdataurl()) +``` + +Hint: Repeat steps for X-ray image. (Note: Ideally, we would find an image in the Chandra 'cxc' catalog) + +```{code-cell} ipython3 + +``` + +## Step 3: Make a spectrum + +### Find what Chandra spectral observations exist already for this source. +Hint: try searching for X-ray spectral data tables using the registry query + +```{code-cell} ipython3 +# Search the Registry to list services that contain X-ray spectral data +``` + +Hint 2: Take a look at what data exist for our candidate, NGC 1365. + +```{code-cell} ipython3 + +``` + +Hint 3: Download the data to make a spectrum. Note: you might end here and use Xspec to plot and model the spectrum. Or ... you can also try to take a quick look at the spectrum. + +```{code-cell} ipython3 +# Get it and look at it: +``` + +```{code-cell} ipython3 +## Or write it to disk +``` + +Extension: Making a "quick look" spectrum. For our purposes, the 1st order of the HEG grating data would be sufficient. + +```{code-cell} ipython3 +# Hint: You'll have to look into the details of the spectra. +``` + +This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html). (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.md) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) + ++++ + +Congratulations! You have completed this notebook exercise. diff --git a/Exercise_III.ipynb b/Exercise_III.ipynb deleted file mode 100644 index fe1d0c7..0000000 --- a/Exercise_III.ipynb +++ /dev/null @@ -1,486 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Creating a stellar color-magnitude (or Hertzsprung-Russell) diagram\n", - "\n", - "The [Hertzsprung-Russell diagram](https://en.wikipedia.org/wiki/Hertzsprung–Russell_diagram) is a fundamental diagram in astronomy that displays important relationships between the stellar color (or temperature) and absolute brightness (or luminosity). \n", - "\n", - "In this exercise, we will use existing stellar catalogs to produce the H-R diagram. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# As a hint, we include the code block for Python modules that you will likely need to import: \n", - "import matplotlib\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "%matplotlib inline \n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "from astropy.io import fits\n", - "\n", - "import pyvo as vo\n", - "from pyvo import registry\n", - "\n", - "## There are a number of relatively unimportant warnings that \n", - "## show up, so for now, suppress them:\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 1: Find appropriate catalogs\n", - "\n", - "We want to find a star catalog that has the available data to produce the H-R diagram, i.e., the absolute magnitudes (or both apparent magnitudes AND distances, so we can calculate the absolute magnitudes) in two optical bands (e.g., B and V). This would give us color. Or we need B- OR V- band magnitude and the stellar temperature. \n", - "\n", - "To simplify this problem, we want to find a catalog of an open cluster of stars, where all the stars were born around the same time and are located in one cluster. This simplifies the issue of getting accurate distances to the stars. One famous cluster is the Pleiades in the constellation of Taurus. So first we start by searching for an existing catalog with data on Pleiades that will provide the necessary information about the stars: magnitudes in two bands (e.g., B and V), which can be used to measure color, or temperature of the star plus one magnitude. \n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### DATA DISCOVERY STEPS: \n", - "\n", - "Here is useful link for [how the pyvo registry search works](https://pyvo.readthedocs.io/en/latest/registry/index.html)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# Write some code to perform a registry search using \n", - "# appropriate keywords and specify the service type. \n", - "# Also apply the includeaux=True option to return \n", - "# the maximum number of services: \n", - "\n", - "\n", - "# Print the size of the output table and get a summary of the sources:\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note: The includeaux=True includes auxiliary services. \n", - "\n", - "Because we specified the service type, this returns only the TAP services for each of the avalible resorces matching our search criteria. So, the 'interfaces' column will only show TAP service as an avalible service, even if other services are avalible from the same resource. We know we need the service to be a TAP service since we know that we want to eventually access the search using additional column information (i.e., beyond RA and Dec, which is all that Simple Cone Search returns). " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Next, we need to find which of these has the columns of interest, i.e. magnitudes in two bands to create the color-magnitude diagram. \n", - "\n", - "We can re-run the registry search, but further restrict the results by column UCD. We want tables that have magnitude columns; the most basic UCD to describe a magnitude column is phot.mag" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Perform another registry search, this time searching\n", - "# tables that have at least 1 magnitude column. \n", - "# Note: '%' serves as a wild card when searching columns\n", - "\n", - "# How many tables do you get? " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note: the '%' serves as a wild card when searching by UCD\n", - "\n", - "The IOVA standard enables resources to be as sepecific as they would like when defining the UCD of columns. For example, 'phot.mag' and 'phot.mag;em.opt.V' can both be used to describe a column containing the V magnitudes of objects. If a resource uses the latter to describe a column, a search using 'phot.mag' will not return that resource. A wild card would need to be used or the exact column UCD. The UCD search requires an exact match for a resource to be returned, so using the wild card will make it easier to discover a wider variety of resources. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "So using this we can reduce the matched tables to ones that are a bit more catered to our experiment. Note, that there is redundancy in some resources since these are available via multiple services and/or publishers. Therefore a bit more cleaning can be done to provide only the unique matches. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Print the 'ivoid' column for these tables \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Write a function that searches the result and only\n", - "# outputs unique tables\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can read more information about the results we found. For each resource element (i.e. row in the table above), there are useful attributes, which are [described here]( https://pyvo.readthedocs.io/en/latest/api/pyvo.registry.regtap.RegistryResource.html#pyvo.registry.regtap.RegistryResource)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# Output the descriptions of the resulting matches:\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - " RESULT: Based on these, the second one (by Eichhorn et al) looks like a good start. \n", - "\n", - "### At this point, you can proceed to Step 2. \n", - "\n", - "-- OR --\n", - "\n", - "### Try a different data discovery method! " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - "### Alternative Method: Use ADS to search for appropriate paper and access data via NED.\n", - "\n", - "There are multiple paths for the data discovery. So it may also be that you know the paper that has the data you are interested in and want to access via the bibcode or authors, etc. \n", - "\n", - "In this case, let's assume that we have the information that the Eichhorn+1970 paper has the data that we need to create the H-R diagram: https://ui.adsabs.harvard.edu/abs/1970MmRAS..73..125E/abstract\n", - "\n", - "We can either search by bibcode (1970MmRAS..73..125E) or \"Eichhorn\" to get the access_urls that will allow us to work with the data. \n", - "\n", - "Before this step, if may help to see the names of the fields available to use. Notice the following fields: \n", - "\n", - "\"source_value\" contains the bibcode information that we want; \"creator_seq\" lists the authors; \n", - "\n", - "and \n", - "\n", - "\"access_url\" provides the url from where the data can be accessed. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Print the fieldnames for the tap services. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "First, Try using bibcode: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bibcode = '1970MmRAS..73..125E' # Eichhorn\n", - "\n", - "# \n", - "# print \"short_name\", \"source_value\" and \"access_url\" information \n", - "# for the table entry that matches this bibcode. \n", - "#\n", - "# Note that using the to_table() lets you search the result \n", - "# easily using all columns. But in the end, you want to get\n", - "# back not an astropy table row, which you cannot use, but the\n", - "# original RegistryResult that has the callable TAP service. \n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that the URL is a generic TAP url for Vizier. All of its tables can be accessed by that same TAP services. It'll be in the ADQL query itself that you specify the table name. We'll see this below." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, try using Author name: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "author = 'Eichhorn'\n", - "\n", - "# \n", - "# print \"short_name\", \"reference_url\" and \"access_url\" information \n", - "# for the table entry that matches this author name. \n", - "#\n", - "# Hint: The pyVO Registry resource page may help:\n", - "# https://pyvo.readthedocs.io/en/latest/api/pyvo.registry.regtap.RegistryResource.html#pyvo.registry.regtap.RegistryResource.search\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "These examples provide a few ways to access the information of interest. \n", - "\n", - "Below are a few other ways to see what the tap_service table contains. \n", - "\n", - "1. To view the column information: tap_services.to_table().columns() shows the metadata contained in the tap service. We will reference some of this columns below as we try to find the appropriate table. \n", - "\n", - "2. tap_services[index].describe(): The table with the tap_services output has, in our case, 83 tables listed and each includes metadata containing some human readable description. You can get the description for one case or for all the records by iterating through the resource. In the former case, we show the description for the Eichhorn data. The latter case also follows. \n", - " " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# Print the full description for the Eichhorn+1970 example. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 2: Acquire the relevant data and make a plot!\n", - "In order to query the table, we need the table name, note this is NOT the same as the short name we found above:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get the table name for the Eichhorn+1970 data\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Query the data to output the table information (i.e. the data in the \n", - "# columns in the table)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can access the column data as array using the .getcolumn(colname) attribute, where the colname is given in the table above. In particular the \"CI\" is the color index and \"Ptm\" is the photovisual magnitude. See [here](https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/90) for details about the columns. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# get color and magnitude column information for\n", - "# stars in this table. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plotting... \n", - "Note: The magnitudes here are apparent and therefore in plotting, the color-magnitude diagram is typically brightness increasing upwards (higher y-axis) so we will flip the y-axis here. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.ylim(15, 0)\n", - "plt.ylabel(\"V [apparent mag]\")\n", - "plt.xlabel(\"B-V\")\n", - "#\n", - "# plot color and magnitude data: \n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 3. Compare with other color-magnitude diagrams for Pleiades:\n", - "\n", - "There is nice discussion here: http://www.southastrodel.com/Page03009a.htm about the color-magnitude diagram. Their Fig 4 looks slightly cleaner because part of this investigation was to select the 270 stars that are vetted members and restricted to stellar types more massive than K0. \n", - "\n", - "The dataset is from Raboud+1998 (1998A&A...329..101R)\n", - "\n", - "Therefore in this next step, we will use the bibcode to select this data and overplot with the previous data to compare. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bibcode = '1998A&A...329..101R' # Raboud\n", - "\n", - "# Repeat steps above with this other dataset:\n", - "# First find the table usign the bib code" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Then query the table using the table name and get the\n", - "# appropriate columns" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Then plot!\n", - "\n", - "plt.ylim(15, 0)\n", - "plt.ylabel(\"V [apparent mag]\")\n", - "plt.xlabel(\"B-V\")\n", - "\n", - "# Plot the Eichhorn data as black circles and \n", - "# Raboud data as red squares. \n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## BONUS: Step 4: The CMD as a distance indicator! \n", - "\n", - "Since the y-axis above is apparent magnitude, we can use the obvious features (e.g., main sequence curve) to translate the apparent magnitudes to absolute magnitudes (by comparing to published H-R diagrams given in absolute magnitudes) and measure the distance to Pleiades! \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "sun_color = 0.65 # from http://www.astro.ucla.edu/~wright/magcolor.htm\n", - "\n", - "# Based on the color-magnitude diagram, what would be\n", - "# the Sun's apparent magnitude at the distance of the Pleiades?\n", - "# i.e. if you overplot the Sun on the CMD above\n", - "# what's the Sun's magnitude, corresponding to the sun's color\n", - "# given as B-V=0.65? \n", - "\n", - "# Overplot the sun in the plot above for reference. \n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Measure the distance to the Pleaides, using the Sun\n", - "# as a reference. We know the Sun's absolute magnitude: \n", - "Vabs = 4.8 ## Sun @ B-V = 0.65 (taken from Wikipedia)\n", - "\n", - "# Using the Sun's apparent magnitude from your estimate above\n", - "# what's the distance to Pleiades in pc? " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "True distance to Pleaides is 136.2 pc (https://en.wikipedia.org/wiki/Pleiades ). Not bad!" - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.2" - }, - "nav_menu": {}, - "toc": { - "navigate_menu": true, - "number_sections": true, - "sideBar": true, - "threshold": 6, - "toc_cell": false, - "toc_section_display": "block", - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/Exercise_III.md b/Exercise_III.md new file mode 100644 index 0000000..08c5d84 --- /dev/null +++ b/Exercise_III.md @@ -0,0 +1,295 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.8.2 +nav_menu: {} +toc: + navigate_menu: true + number_sections: true + sideBar: true + threshold: 6 + toc_cell: false + toc_section_display: block + toc_window_display: true +--- + +# Creating a stellar color-magnitude (or Hertzsprung-Russell) diagram + +The [Hertzsprung-Russell diagram](https://en.wikipedia.org/wiki/Hertzsprung–Russell_diagram) is a fundamental diagram in astronomy that displays important relationships between the stellar color (or temperature) and absolute brightness (or luminosity). + +In this exercise, we will use existing stellar catalogs to produce the H-R diagram. + +```{code-cell} ipython3 +# As a hint, we include the code block for Python modules that you will likely need to import: +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +%matplotlib inline + +# For downloading files +from astropy.utils.data import download_file +from astropy.io import fits + +import pyvo as vo +from pyvo import registry + +## There are a number of relatively unimportant warnings that +## show up, so for now, suppress them: +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +``` + +## Step 1: Find appropriate catalogs + +We want to find a star catalog that has the available data to produce the H-R diagram, i.e., the absolute magnitudes (or both apparent magnitudes AND distances, so we can calculate the absolute magnitudes) in two optical bands (e.g., B and V). This would give us color. Or we need B- OR V- band magnitude and the stellar temperature. + +To simplify this problem, we want to find a catalog of an open cluster of stars, where all the stars were born around the same time and are located in one cluster. This simplifies the issue of getting accurate distances to the stars. One famous cluster is the Pleiades in the constellation of Taurus. So first we start by searching for an existing catalog with data on Pleiades that will provide the necessary information about the stars: magnitudes in two bands (e.g., B and V), which can be used to measure color, or temperature of the star plus one magnitude. + ++++ + +### DATA DISCOVERY STEPS: + +Here is useful link for [how the pyvo registry search works](https://pyvo.readthedocs.io/en/latest/registry/index.html). + +```{code-cell} ipython3 +:tags: [output_scroll] + +# Write some code to perform a registry search using +# appropriate keywords and specify the service type. +# Also apply the includeaux=True option to return +# the maximum number of services: + + +# Print the size of the output table and get a summary of the sources: +``` + +Note: The includeaux=True includes auxiliary services. + +Because we specified the service type, this returns only the TAP services for each of the avalible resorces matching our search criteria. So, the 'interfaces' column will only show TAP service as an avalible service, even if other services are avalible from the same resource. We know we need the service to be a TAP service since we know that we want to eventually access the search using additional column information (i.e., beyond RA and Dec, which is all that Simple Cone Search returns). + ++++ + +#### Next, we need to find which of these has the columns of interest, i.e. magnitudes in two bands to create the color-magnitude diagram. + +We can re-run the registry search, but further restrict the results by column UCD. We want tables that have magnitude columns; the most basic UCD to describe a magnitude column is phot.mag + +```{code-cell} ipython3 +# Perform another registry search, this time searching +# tables that have at least 1 magnitude column. +# Note: '%' serves as a wild card when searching columns + +# How many tables do you get? +``` + +Note: the '%' serves as a wild card when searching by UCD + +The IOVA standard enables resources to be as sepecific as they would like when defining the UCD of columns. For example, 'phot.mag' and 'phot.mag;em.opt.V' can both be used to describe a column containing the V magnitudes of objects. If a resource uses the latter to describe a column, a search using 'phot.mag' will not return that resource. A wild card would need to be used or the exact column UCD. The UCD search requires an exact match for a resource to be returned, so using the wild card will make it easier to discover a wider variety of resources. + ++++ + +So using this we can reduce the matched tables to ones that are a bit more catered to our experiment. Note, that there is redundancy in some resources since these are available via multiple services and/or publishers. Therefore a bit more cleaning can be done to provide only the unique matches. + +```{code-cell} ipython3 +# Print the 'ivoid' column for these tables +``` + +```{code-cell} ipython3 +# Write a function that searches the result and only +# outputs unique tables +``` + +We can read more information about the results we found. For each resource element (i.e. row in the table above), there are useful attributes, which are [described here]( https://pyvo.readthedocs.io/en/latest/api/pyvo.registry.regtap.RegistryResource.html#pyvo.registry.regtap.RegistryResource) + +```{code-cell} ipython3 +:tags: [output_scroll] + +# Output the descriptions of the resulting matches: +``` + + RESULT: Based on these, the second one (by Eichhorn et al) looks like a good start. + +### At this point, you can proceed to Step 2. + +-- OR -- + +### Try a different data discovery method! + ++++ {"tags": ["output_scroll"]} + +### Alternative Method: Use ADS to search for appropriate paper and access data via NED. + +There are multiple paths for the data discovery. So it may also be that you know the paper that has the data you are interested in and want to access via the bibcode or authors, etc. + +In this case, let's assume that we have the information that the Eichhorn+1970 paper has the data that we need to create the H-R diagram: https://ui.adsabs.harvard.edu/abs/1970MmRAS..73..125E/abstract + +We can either search by bibcode (1970MmRAS..73..125E) or "Eichhorn" to get the access_urls that will allow us to work with the data. + +Before this step, if may help to see the names of the fields available to use. Notice the following fields: + +"source_value" contains the bibcode information that we want; "creator_seq" lists the authors; + +and + +"access_url" provides the url from where the data can be accessed. + +```{code-cell} ipython3 +# Print the fieldnames for the tap services. +``` + +First, Try using bibcode: + +```{code-cell} ipython3 +bibcode = '1970MmRAS..73..125E' # Eichhorn + +# +# print "short_name", "source_value" and "access_url" information +# for the table entry that matches this bibcode. +# +# Note that using the to_table() lets you search the result +# easily using all columns. But in the end, you want to get +# back not an astropy table row, which you cannot use, but the +# original RegistryResult that has the callable TAP service. +``` + +Note that the URL is a generic TAP url for Vizier. All of its tables can be accessed by that same TAP services. It'll be in the ADQL query itself that you specify the table name. We'll see this below. + ++++ + +Next, try using Author name: + +```{code-cell} ipython3 +author = 'Eichhorn' + +# +# print "short_name", "reference_url" and "access_url" information +# for the table entry that matches this author name. +# +# Hint: The pyVO Registry resource page may help: +# https://pyvo.readthedocs.io/en/latest/api/pyvo.registry.regtap.RegistryResource.html#pyvo.registry.regtap.RegistryResource.search +``` + +These examples provide a few ways to access the information of interest. + +Below are a few other ways to see what the tap_service table contains. + +1. To view the column information: tap_services.to_table().columns() shows the metadata contained in the tap service. We will reference some of this columns below as we try to find the appropriate table. + +2. tap_services[index].describe(): The table with the tap_services output has, in our case, 83 tables listed and each includes metadata containing some human readable description. You can get the description for one case or for all the records by iterating through the resource. In the former case, we show the description for the Eichhorn data. The latter case also follows. + + +```{code-cell} ipython3 +:tags: [output_scroll] + +# Print the full description for the Eichhorn+1970 example. +``` + +## Step 2: Acquire the relevant data and make a plot! +In order to query the table, we need the table name, note this is NOT the same as the short name we found above: + +```{code-cell} ipython3 +# Get the table name for the Eichhorn+1970 data +``` + +```{code-cell} ipython3 +# Query the data to output the table information (i.e. the data in the +# columns in the table) +``` + +We can access the column data as array using the .getcolumn(colname) attribute, where the colname is given in the table above. In particular the "CI" is the color index and "Ptm" is the photovisual magnitude. See [here](https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/90) for details about the columns. + +```{code-cell} ipython3 +# get color and magnitude column information for +# stars in this table. +``` + +### Plotting... +Note: The magnitudes here are apparent and therefore in plotting, the color-magnitude diagram is typically brightness increasing upwards (higher y-axis) so we will flip the y-axis here. + +```{code-cell} ipython3 +plt.ylim(15, 0) +plt.ylabel("V [apparent mag]") +plt.xlabel("B-V") +# +# plot color and magnitude data: + +``` + +## Step 3. Compare with other color-magnitude diagrams for Pleiades: + +There is nice discussion here: http://www.southastrodel.com/Page03009a.htm about the color-magnitude diagram. Their Fig 4 looks slightly cleaner because part of this investigation was to select the 270 stars that are vetted members and restricted to stellar types more massive than K0. + +The dataset is from Raboud+1998 (1998A&A...329..101R) + +Therefore in this next step, we will use the bibcode to select this data and overplot with the previous data to compare. + +```{code-cell} ipython3 +bibcode = '1998A&A...329..101R' # Raboud + +# Repeat steps above with this other dataset: +# First find the table usign the bib code +``` + +```{code-cell} ipython3 +# Then query the table using the table name and get the +# appropriate columns +``` + +```{code-cell} ipython3 +# Then plot! + +plt.ylim(15, 0) +plt.ylabel("V [apparent mag]") +plt.xlabel("B-V") + +# Plot the Eichhorn data as black circles and +# Raboud data as red squares. +``` + +## BONUS: Step 4: The CMD as a distance indicator! + +Since the y-axis above is apparent magnitude, we can use the obvious features (e.g., main sequence curve) to translate the apparent magnitudes to absolute magnitudes (by comparing to published H-R diagrams given in absolute magnitudes) and measure the distance to Pleiades! + +```{code-cell} ipython3 + +sun_color = 0.65 # from http://www.astro.ucla.edu/~wright/magcolor.htm + +# Based on the color-magnitude diagram, what would be +# the Sun's apparent magnitude at the distance of the Pleiades? +# i.e. if you overplot the Sun on the CMD above +# what's the Sun's magnitude, corresponding to the sun's color +# given as B-V=0.65? + +# Overplot the sun in the plot above for reference. + +``` + +```{code-cell} ipython3 +# Measure the distance to the Pleaides, using the Sun +# as a reference. We know the Sun's absolute magnitude: +Vabs = 4.8 ## Sun @ B-V = 0.65 (taken from Wikipedia) + +# Using the Sun's apparent magnitude from your estimate above +# what's the distance to Pleiades in pc? +``` + +True distance to Pleaides is 136.2 pc (https://en.wikipedia.org/wiki/Pleiades ). Not bad! diff --git a/QuickReference.ipynb b/QuickReference.ipynb deleted file mode 100644 index 149da9a..0000000 --- a/QuickReference.ipynb +++ /dev/null @@ -1,423 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Quick Reference" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 0. Setup\n", - "Please make sure your environment is set up according to the [instructions here](https://nasa-navo.github.io/navo-workshop/00_SETUP.html).\n", - "\n", - "Ensure you have the latest version of the workshop material by updating your environment:\n", - "TBD" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 1. Overview\n", - "NASA services can be queried from Python in multiple ways.\n", - "* Generic Virtual Observatory (VO) queries.\n", - " * Call sequence is consistent, including for non-NASA resources.\n", - " * Use the [`pyvo` package](https://pyvo.readthedocs.io/en/latest/).\n", - " * [Known issues/caveats](https://github.com/NASA-NAVO/navo-workshop/blob/main/KNOWN_ISSUES.md).\n", - "* [Astroquery](https://astroquery.readthedocs.io/en/latest/) interfaces:\n", - " * Call sequences not quite as consistent, but follow similar patterns.\n", - "* Ad hoc archive-specific interfaces\n", - "\n", - "## 2. VO Services\n", - "This workshop will introduce 4 types of VO queries:\n", - "* **VO Registry** - Discover what services are available worldwide\n", - "* **Simple Cone Search** - Search for catalog object within a specified cone region\n", - "* **Simple Image Access** - Search for image products within a spatial region\n", - "* **Simple Spectral Access** - Search for spectral products within a spatial region\n", - "* **Table Access** - SQL-like queries to databases\n", - "\n", - "### 2.0 Import Necessary Packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Generic VO access routines\n", - "import pyvo as vo\n", - "\n", - "# For specifying coordinates and angles\n", - "from astropy.coordinates import SkyCoord\n", - "from astropy.coordinates import Angle\n", - "from astropy import units as u\n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "\n", - "# Ignore unimportant warnings\n", - "import warnings\n", - "warnings.filterwarnings('ignore', '.*Unknown element mirrorURL.*', vo.utils.xml.elements.UnknownElementWarning)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.1 Look Up Services in VO Registry\n", - "Simple example: Find Simple Cone Search (conesearch) services related to SWIFT." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='conesearch', keywords=['swift'])\n", - "services" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### 2.1.1 Use different arguments/values to modify the simple example\n", - "| Argument | Description | Examples |\n", - "| :-----: | :----------- | :-------- |\n", - "| **servicetype** | Type of service | `conesearch` or `scs` for **Simple Cone Search**
`image` or `sia` for **Simple Image Access**
`spectrum` or `ssa` for **Simple Spectral Access**
`table` or `tap` for **Table Access Protocol**|\n", - "| **keyword** | List of one or more keyword(s) to match service's metadata. Both ORs and ANDs may be specified.
| `['galex', 'swift']` matches 'galex' or 'swift'
`['hst survey']` matches services mentioning both 'hst' and 'survey' |\n", - "| **waveband** | Resulting services have data in the specified waveband(s) | ‘radio’, ‘millimeter’, ‘infrared’, ‘optical’, ‘uv’, ‘euv’, ‘x-ray’ ‘gamma-ray’ |\n", - "\n", - "#### 2.1.2 Inspect the results.\n", - "##### Using pyvo\n", - "Although not lists, `pyvo` results can be iterated over to see each individual result. The results are specialized based on the type of query, providing access to the important properties of the results. Some useful accessors with registry results are:\n", - "* `short_name` - A short name\n", - "* `res_title` - A more descriptive title\n", - "* `res_description` - A more verbose description\n", - "* `reference_url` - A link for more information\n", - "* `ivoid` - A unique identifier for the service. Gives some indication of what organization is serving the data." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Print the number of results and the 1st 4 short names and titles.\n", - "print(f'Number of results: {len(services)}\\n')\n", - "for s in list(services)[:4]: # (Treat services as list to get the subset of rows)\n", - " print(f'{s.short_name} - {s.res_title}')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "##### Filtering results\n", - "Of the services we found, which one(s) have 'stsci.edu' in their unique identifier?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "stsci_services = [s for s in services if 'stsci.edu' in s.ivoid]\n", - "for s in stsci_services:\n", - " print (f'(STScI): {s.short_name} - {s.res_title}')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "##### Using astropy\n", - "With the `to_table()` method, `pyvo` results can also be converted to Astropy `Table` objects which offer a variety of addional features. See [astropy.table](http://docs.astropy.org/en/stable/table/) for more on working with Astropy Tables." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Convert to an Astropy Table\n", - "services_table = services.to_table()\n", - "\n", - "# Print the column names and display 1st 3 rows with a subset of columns\n", - "print(f'\\nColumn Names:\\n{services_table.colnames}\\n') \n", - "services_table['short_name', 'res_title', 'res_description'][:3] " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.2 Cone search\n", - "Example: Find a cone search service for the USNO-B catalog and search it around M51 with a .1 degree radius. (More inspection could be done on the service list instead of blindly choosing the first service.) \n", - "\n", - "The position (`pos`) is best specified with [SkyCoord](http://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html) objects. \n", - "\n", - "The size of the region is specified with the `radius` keyword and may be decimal degrees or an Astropy [Angle](http://docs.astropy.org/en/stable/api/astropy.coordinates.Angle.html#astropy.coordinates.Angle). " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m51_pos = SkyCoord.from_name(\"m51\")\n", - "services = vo.regsearch(servicetype='conesearch', keywords='usno-b')\n", - "results = services[0].search(pos=m51_pos, radius=0.1)\n", - "# Astropy Table is useful for displaying cone search results.\n", - "results.to_table() " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.3 Image search\n", - "Example: Find an image search service for GALEX, and search it around coordinates 13:37:00.950,-29:51:55.51 (M83) with a radius of .2 degrees. Download the first file in the results.\n", - "#### Find an image service" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='image', keywords=['galex'])\n", - "services.to_table()['ivoid', 'short_name', 'res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Search one of the services\n", - "The first service looks good. Search it!\n", - "\n", - "For more details on using `SkyCoord`, see [its documentation](http://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord)\n", - "\n", - "**NOTE**: For image searches, the size of the region is defined by the `size` keyword which is more like a diameter than a radius." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "m83_pos = SkyCoord('13h37m00.950s -29d51m55.51s')\n", - "results = services[1].search(pos=m83_pos, size=.2)\n", - "\n", - "# We can look at the results.\n", - "results.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Download an image\n", - "For the first result, print the file format and download the file. If repeatedly executing this code, add `cache=True` to `download_file()` to prevent repeated downloads.\n", - "\n", - "See [`download_file()` documentation here.](https://docs.astropy.org/en/stable/api/astropy.utils.data.download_file.html#astropy.utils.data.download_file)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(results[0].format)\n", - "file_name = download_file(results[0].getdataurl()) \n", - "file_name" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.4 Spectral search\n", - "Example: Find a spectral service for x-ray data. Query it around Delta Ori with a search **diameter** of 10 arc minutes, and download the first data product. Note that the results table can be inspected for potentially useful columns.\n", - "\n", - "Spectral search is very similar to image search. In this example, note:\n", - "* **`diameter`** defines the size of the search region\n", - "* `waveband` used in `regsearch()`\n", - "* Astropy `Angle` used to specify radius units other than degrees." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Search for a spectrum search service that has x-ray data.\n", - "services = vo.regsearch(servicetype='spectrum', waveband='x-ray')\n", - "\n", - "# Assuming there are services and the first one is OK...\n", - "results = services[0].search(pos=SkyCoord.from_name(\"Delta Ori\"), \n", - " diameter=Angle(10 * u.arcmin))\n", - "\n", - "# Assuming there are results, download the first file.\n", - "print(f'Title: {results[0].title}, Format: {results[0].format}')\n", - "file_name = download_file(results[0].getdataurl()) \n", - "file_name" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2.5 Table search\n", - "Example: Find the HEASARC Table Access Protocol (TAP) service, get some information about the available tables." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "services = vo.regsearch(servicetype='tap', keywords=['heasarc'])\n", - "print(f'{len(services)} service(s) found.')\n", - "# We found only one service. Print some info about the service and its tables.\n", - "print(f'{services[0].describe()}')\n", - "tables = services[0].service.tables # Queries for details of the service's tables\n", - "print(f'{len(tables)} tables:')\n", - "for t in tables:\n", - " print(f'{t.name:30s} - {t.description}') # A more succinct option than t.describe()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Column Information\n", - "For any table, we can list the column names and descriptions." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for c in tables['zcat'].columns:\n", - " print(f'{c.name:30s} - {c.description}')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Perform a Query\n", - "Example: Perform a cone search on the ZCAT catalog at M83 with a 1.0 degree radius." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "coord = SkyCoord.from_name(\"m83\")\n", - "query = f'''\n", - "SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where \n", - "contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{coord.ra.deg},{coord.dec.deg},1.0))=1\n", - "'''\n", - "results = services[0].service.run_async(query)\n", - "\n", - "results.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 3. Astroquery \n", - "Many archives have Astroquery or Astroquery-compliant modules for data access, including:\n", - "\n", - "* Astroquery\n", - " * [HEASARC Queries (astroquery.heasarc)](https://astroquery.readthedocs.io/en/latest/heasarc/heasarc.html)\n", - " * [HITRAN Queries (astroquery.hitran)](https://astroquery.readthedocs.io/en/latest/hitran/hitran.html)\n", - " * [IRSA Image Server program interface (IBE) Queries (astroquery.ibe)](https://astroquery.readthedocs.io/en/latest/ibe/ibe.html)\n", - " * [IRSA Queries (astroquery.ipac.irsa)](https://astroquery.readthedocs.io/en/latest/ipac/irsa/irsa.html)\n", - " * [IRSA Dust Extinction Service Queries (astroquery.ipac.irsa.irsa_dust)](https://astroquery.readthedocs.io/en/latest/ipac/irsa/irsa_dust.html)\n", - " * [JPL Spectroscopy Queries (astroquery.jplspec)](https://astroquery.readthedocs.io/en/latest/jplspec/jplspec.html)\n", - " * [MAST Queries (astroquery.mast)](https://astroquery.readthedocs.io/en/latest/mast/mast.html)\n", - " * [NASA ADS Queries (astroquery.nasa_ads)](https://astroquery.readthedocs.io/en/latest/nasa_ads/nasa_ads.html)\n", - " * [NED Queries (astroquery.ipac.ned)](https://astroquery.readthedocs.io/en/latest/ipac/ned/ned.html)\n", - "* Astroquery-compliant\n", - " * [KOA Queries (pykoa.koa)](https://koa.ipac.caltech.edu/UserGuide/PyKOA/PyKOA.html)\n", - "\n", - "For more, see https://astroquery.readthedocs.io/en/latest/\n", - "\n", - "### 3.1 NED\n", - "Example: Get an Astropy Table containing the objects from paper 2018ApJ...858...62K. For more on the API, see [astroquery](https://astroquery.readthedocs.io/en/latest/ipac/ned/ned.html)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from astroquery.ipac.ned import Ned\n", - "objects_in_paper = Ned.query_refcode('2018ApJ...858...62K')\n", - "objects_in_paper" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.9" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} \ No newline at end of file diff --git a/QuickReference.md b/QuickReference.md new file mode 100644 index 0000000..87fb41e --- /dev/null +++ b/QuickReference.md @@ -0,0 +1,273 @@ +--- +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.9.9 +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: {} + toc_section_display: true + toc_window_display: true +--- + +# Quick Reference + ++++ + +## 0. Setup +Please make sure your environment is set up according to the [instructions here](https://nasa-navo.github.io/navo-workshop/00_SETUP.html). + +Ensure you have the latest version of the workshop material by updating your environment: +TBD + ++++ + +## 1. Overview +NASA services can be queried from Python in multiple ways. +* Generic Virtual Observatory (VO) queries. + * Call sequence is consistent, including for non-NASA resources. + * Use the [`pyvo` package](https://pyvo.readthedocs.io/en/latest/). + * [Known issues/caveats](https://github.com/NASA-NAVO/navo-workshop/blob/main/KNOWN_ISSUES.md). +* [Astroquery](https://astroquery.readthedocs.io/en/latest/) interfaces: + * Call sequences not quite as consistent, but follow similar patterns. +* Ad hoc archive-specific interfaces + +## 2. VO Services +This workshop will introduce 4 types of VO queries: +* **VO Registry** - Discover what services are available worldwide +* **Simple Cone Search** - Search for catalog object within a specified cone region +* **Simple Image Access** - Search for image products within a spatial region +* **Simple Spectral Access** - Search for spectral products within a spatial region +* **Table Access** - SQL-like queries to databases + +### 2.0 Import Necessary Packages + +```{code-cell} ipython3 +# Generic VO access routines +import pyvo as vo + +# For specifying coordinates and angles +from astropy.coordinates import SkyCoord +from astropy.coordinates import Angle +from astropy import units as u + +# For downloading files +from astropy.utils.data import download_file + +# Ignore unimportant warnings +import warnings +warnings.filterwarnings('ignore', '.*Unknown element mirrorURL.*', vo.utils.xml.elements.UnknownElementWarning) +``` + +### 2.1 Look Up Services in VO Registry +Simple example: Find Simple Cone Search (conesearch) services related to SWIFT. + +```{code-cell} ipython3 +services = vo.regsearch(servicetype='conesearch', keywords=['swift']) +services +``` + +#### 2.1.1 Use different arguments/values to modify the simple example +| Argument | Description | Examples | +| :-----: | :----------- | :-------- | +| **servicetype** | Type of service | `conesearch` or `scs` for **Simple Cone Search**
`image` or `sia` for **Simple Image Access**
`spectrum` or `ssa` for **Simple Spectral Access**
`table` or `tap` for **Table Access Protocol**| +| **keyword** | List of one or more keyword(s) to match service's metadata. Both ORs and ANDs may be specified.
| `['galex', 'swift']` matches 'galex' or 'swift'
`['hst survey']` matches services mentioning both 'hst' and 'survey' | +| **waveband** | Resulting services have data in the specified waveband(s) | ‘radio’, ‘millimeter’, ‘infrared’, ‘optical’, ‘uv’, ‘euv’, ‘x-ray’ ‘gamma-ray’ | + +#### 2.1.2 Inspect the results. +##### Using pyvo +Although not lists, `pyvo` results can be iterated over to see each individual result. The results are specialized based on the type of query, providing access to the important properties of the results. Some useful accessors with registry results are: +* `short_name` - A short name +* `res_title` - A more descriptive title +* `res_description` - A more verbose description +* `reference_url` - A link for more information +* `ivoid` - A unique identifier for the service. Gives some indication of what organization is serving the data. + +```{code-cell} ipython3 +# Print the number of results and the 1st 4 short names and titles. +print(f'Number of results: {len(services)}\n') +for s in list(services)[:4]: # (Treat services as list to get the subset of rows) + print(f'{s.short_name} - {s.res_title}') +``` + +##### Filtering results +Of the services we found, which one(s) have 'stsci.edu' in their unique identifier? + +```{code-cell} ipython3 +stsci_services = [s for s in services if 'stsci.edu' in s.ivoid] +for s in stsci_services: + print (f'(STScI): {s.short_name} - {s.res_title}') +``` + +##### Using astropy +With the `to_table()` method, `pyvo` results can also be converted to Astropy `Table` objects which offer a variety of addional features. See [astropy.table](http://docs.astropy.org/en/stable/table/) for more on working with Astropy Tables. + +```{code-cell} ipython3 +# Convert to an Astropy Table +services_table = services.to_table() + +# Print the column names and display 1st 3 rows with a subset of columns +print(f'\nColumn Names:\n{services_table.colnames}\n') +services_table['short_name', 'res_title', 'res_description'][:3] +``` + +### 2.2 Cone search +Example: Find a cone search service for the USNO-B catalog and search it around M51 with a .1 degree radius. (More inspection could be done on the service list instead of blindly choosing the first service.) + +The position (`pos`) is best specified with [SkyCoord](http://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html) objects. + +The size of the region is specified with the `radius` keyword and may be decimal degrees or an Astropy [Angle](http://docs.astropy.org/en/stable/api/astropy.coordinates.Angle.html#astropy.coordinates.Angle). + +```{code-cell} ipython3 +m51_pos = SkyCoord.from_name("m51") +services = vo.regsearch(servicetype='conesearch', keywords='usno-b') +results = services[0].search(pos=m51_pos, radius=0.1) +# Astropy Table is useful for displaying cone search results. +results.to_table() +``` + +### 2.3 Image search +Example: Find an image search service for GALEX, and search it around coordinates 13:37:00.950,-29:51:55.51 (M83) with a radius of .2 degrees. Download the first file in the results. +#### Find an image service + +```{code-cell} ipython3 +services = vo.regsearch(servicetype='image', keywords=['galex']) +services.to_table()['ivoid', 'short_name', 'res_title'] +``` + +#### Search one of the services +The first service looks good. Search it! + +For more details on using `SkyCoord`, see [its documentation](http://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html#astropy.coordinates.SkyCoord) + +**NOTE**: For image searches, the size of the region is defined by the `size` keyword which is more like a diameter than a radius. + +```{code-cell} ipython3 +m83_pos = SkyCoord('13h37m00.950s -29d51m55.51s') +results = services[1].search(pos=m83_pos, size=.2) + +# We can look at the results. +results.to_table() +``` + +#### Download an image +For the first result, print the file format and download the file. If repeatedly executing this code, add `cache=True` to `download_file()` to prevent repeated downloads. + +See [`download_file()` documentation here.](https://docs.astropy.org/en/stable/api/astropy.utils.data.download_file.html#astropy.utils.data.download_file) + +```{code-cell} ipython3 +print(results[0].format) +file_name = download_file(results[0].getdataurl()) +file_name +``` + +### 2.4 Spectral search +Example: Find a spectral service for x-ray data. Query it around Delta Ori with a search **diameter** of 10 arc minutes, and download the first data product. Note that the results table can be inspected for potentially useful columns. + +Spectral search is very similar to image search. In this example, note: +* **`diameter`** defines the size of the search region +* `waveband` used in `regsearch()` +* Astropy `Angle` used to specify radius units other than degrees. + +```{code-cell} ipython3 +# Search for a spectrum search service that has x-ray data. +services = vo.regsearch(servicetype='spectrum', waveband='x-ray') + +# Assuming there are services and the first one is OK... +results = services[0].search(pos=SkyCoord.from_name("Delta Ori"), + diameter=Angle(10 * u.arcmin)) + +# Assuming there are results, download the first file. +print(f'Title: {results[0].title}, Format: {results[0].format}') +file_name = download_file(results[0].getdataurl()) +file_name +``` + +### 2.5 Table search +Example: Find the HEASARC Table Access Protocol (TAP) service, get some information about the available tables. + +```{code-cell} ipython3 +:tags: [output_scroll] + +services = vo.regsearch(servicetype='tap', keywords=['heasarc']) +print(f'{len(services)} service(s) found.') +# We found only one service. Print some info about the service and its tables. +print(f'{services[0].describe()}') +tables = services[0].service.tables # Queries for details of the service's tables +print(f'{len(tables)} tables:') +for t in tables: + print(f'{t.name:30s} - {t.description}') # A more succinct option than t.describe() +``` + +#### Column Information +For any table, we can list the column names and descriptions. + +```{code-cell} ipython3 +for c in tables['zcat'].columns: + print(f'{c.name:30s} - {c.description}') +``` + +#### Perform a Query +Example: Perform a cone search on the ZCAT catalog at M83 with a 1.0 degree radius. + +```{code-cell} ipython3 +coord = SkyCoord.from_name("m83") +query = f''' +SELECT ra, dec, Radial_Velocity, radial_velocity_error, bmag, morph_type FROM public.zcat as cat where +contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{coord.ra.deg},{coord.dec.deg},1.0))=1 +''' +results = services[0].service.run_async(query) + +results.to_table() +``` + +## 3. Astroquery +Many archives have Astroquery or Astroquery-compliant modules for data access, including: + +* Astroquery + * [HEASARC Queries (astroquery.heasarc)](https://astroquery.readthedocs.io/en/latest/heasarc/heasarc.html) + * [HITRAN Queries (astroquery.hitran)](https://astroquery.readthedocs.io/en/latest/hitran/hitran.html) + * [IRSA Image Server program interface (IBE) Queries (astroquery.ibe)](https://astroquery.readthedocs.io/en/latest/ibe/ibe.html) + * [IRSA Queries (astroquery.ipac.irsa)](https://astroquery.readthedocs.io/en/latest/ipac/irsa/irsa.html) + * [IRSA Dust Extinction Service Queries (astroquery.ipac.irsa.irsa_dust)](https://astroquery.readthedocs.io/en/latest/ipac/irsa/irsa_dust.html) + * [JPL Spectroscopy Queries (astroquery.jplspec)](https://astroquery.readthedocs.io/en/latest/jplspec/jplspec.html) + * [MAST Queries (astroquery.mast)](https://astroquery.readthedocs.io/en/latest/mast/mast.html) + * [NASA ADS Queries (astroquery.nasa_ads)](https://astroquery.readthedocs.io/en/latest/nasa_ads/nasa_ads.html) + * [NED Queries (astroquery.ipac.ned)](https://astroquery.readthedocs.io/en/latest/ipac/ned/ned.html) +* Astroquery-compliant + * [KOA Queries (pykoa.koa)](https://koa.ipac.caltech.edu/UserGuide/PyKOA/PyKOA.html) + +For more, see https://astroquery.readthedocs.io/en/latest/ + +### 3.1 NED +Example: Get an Astropy Table containing the objects from paper 2018ApJ...858...62K. For more on the API, see [astroquery](https://astroquery.readthedocs.io/en/latest/ipac/ned/ned.html) + +```{code-cell} ipython3 +from astroquery.ipac.ned import Ned +objects_in_paper = Ned.query_refcode('2018ApJ...858...62K') +objects_in_paper +``` diff --git a/UseCase_I.ipynb b/UseCase_I.ipynb deleted file mode 100644 index 4c9cdab..0000000 --- a/UseCase_I.ipynb +++ /dev/null @@ -1,709 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Science User Case - Inspecting a Candidate List\n", - "\n", - "[Ogle et al. (2016)](https://ui.adsabs.harvard.edu/abs/2016ApJ...817..109O/abstract) mined the NASA/IPAC Extragalactic Database (NED) to identify a new type of galaxy: Superluminous Spiral Galaxies.\n", - "\n", - "Table 1 lists the positions of these Super Spirals. Based on those positions, let's create multiwavelength cutouts for each super spiral to see what is unique about this new class of objects." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 1. Import the Python modules we'll be using." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Suppress unimportant warnings.\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")\n", - "warnings.filterwarnings('ignore', '.*RADECSYS=*', append=True)\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "\n", - "from astropy.coordinates import SkyCoord\n", - "from astropy.io import fits\n", - "from astropy.nddata import Cutout2D\n", - "import astropy.visualization as vis\n", - "from astropy.wcs import WCS\n", - "from astroquery.ipac.ned import Ned\n", - "\n", - "import pyvo as vo" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The next cell prepares the notebook to display our visualizations." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 2. Search NED for objects in this paper.\n", - "\n", - "Insert a Code Cell below by clicking on the \"Insert\" Menu and choosing \"Insert Cell Below\". Then consult QuickReference.md to figure out how to use astroquery to search NED for all objects in a paper, based on the refcode of the paper. Inspect the resulting astropy table." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "objects_in_paper = Ned.query_refcode('2016ApJ...817..109O')\n", - "objects_in_paper.show_in_notebook()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 3. Filter the NED results.\n", - "\n", - "The results from NED will include galaxies, but also other kinds of objects (e.g. galaxy clusters, galaxy groups). Print the 'Type' column to see the full range of classifications and filter the results so that we only keep the galaxies in the list." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "objects_in_paper['Type']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# Keep only the galaxies from the list\n", - "galaxies = objects_in_paper[np.array(objects_in_paper['Type']) == 'G']\n", - "\n", - "galaxies.show_in_notebook()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 4. Search the NAVO Registry for image resources.\n", - "\n", - "The paper selected super spirals using WISE, SDSS, and GALEX images. Search the NAVO registry for all image resources, using the 'service_type' search parameter. How many image resources are currently available?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_services = vo.regsearch(servicetype='image')\n", - "\n", - "print(f'{len(image_services)} result(s) found.')\n", - "\n", - "image_services.to_table()['ivoid', 'short_name', 'res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 5. Search the NAVO Registry for image resources that will allow you to search for AllWISE images.\n", - "\n", - "There are hundreds of image resources...too many to quickly read through. Try adding the 'keywords' search parameter to your registry search, and find the image resource you would need to search the AllWISE images. Remember from the Known Issues that 'keywords' must be a list." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "allwise_image_services = vo.regsearch(servicetype='image', keywords=['allwise'])\n", - "\n", - "print(f'{len(allwise_image_services)} result(s) found.')\n", - "\n", - "allwise_image_services.to_table()['ivoid', 'short_name', 'res_title']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 6. Choose the AllWISE image service that you are interested in." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "allwise_image_service = allwise_image_services[0]\n", - "allwise_image_service.service" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 7. Choose one of the galaxies in the NED list.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ra = galaxies['RA'][0]\n", - "dec = galaxies['DEC'][0]\n", - "pos = SkyCoord(ra, dec, unit = 'deg')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ra,dec" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 8. Search for a list of AllWISE images that cover this galaxy.\n", - "\n", - "How many images are returned? Which are you most interested in?" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "allwise_image_table = allwise_image_service.search(pos=pos, size=0)\n", - "allwise_image_table" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 9. Use the .to_table() method to view the results as an Astropy table." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "allwise_images = allwise_image_table.to_table()\n", - "allwise_images" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 10. From the result in 8., select the first record for an image taken in WISE band W1 (3.6 micron)\n", - "\n", - "Hints:\n", - "* Loop over records and test on the `.bandpass_id` attribute of each record\n", - "* Print the `.title` and `.bandpass_id` of the record you find, to verify it is the right one." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for allwise_image_record in allwise_image_table:\n", - " if 'W1' in allwise_image_record.bandpass_id:\n", - " break\n", - "print(allwise_image_record.title, allwise_image_record.bandpass_id)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 11. Visualize this AllWISE image." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## If you only run this once, you can do it in memory in one line:\n", - "## This fetches the FITS as an astropy.io.fits object in memory\n", - "#allwise_w1_image = allwise_image_record.getdataobj()\n", - "## But if you might run this notebook repeatedly with limited bandwidth, \n", - "## download it once and cache it. \n", - "file_name = download_file(allwise_image_record.getdataurl(), cache=True) \n", - "allwise_w1_image = fits.open(file_name)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig = plt.figure()\n", - "\n", - "wcs = WCS(allwise_w1_image[0].header)\n", - "ax = fig.add_subplot(1, 1, 1, projection=wcs)\n", - "ax.imshow(allwise_w1_image[0].data, cmap='gray_r', origin='lower', vmax = 10)\n", - "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 12. Plot a cutout of the AllWISE image, centered on your position.\n", - "\n", - "Try a 60 arcsecond cutout." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "size = 60\n", - "cutout = Cutout2D(allwise_w1_image[0].data, pos, (size, size), wcs=wcs)\n", - "wcs = cutout.wcs\n", - "\n", - "fig = plt.figure()\n", - "\n", - "ax = fig.add_subplot(1, 1, 1, projection=wcs)\n", - "ax.imshow(cutout.data, cmap='gray_r', origin='lower', vmax = 10)\n", - "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 13. Try visualizing a cutout of a GALEX image that covers your position.\n", - "\n", - "Repeat steps 4, 5, 6, 8 through 12 for GALEX." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "galex_image_services = vo.regsearch(keywords=['galex'], servicetype='image')\n", - "print(f'{len(galex_image_services)} result(s) found.')\n", - "galex_image_services.to_table()['ivoid', 'short_name', 'res_title']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "galex_image_service = galex_image_services[0]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "galex_image_table = galex_image_service.search(pos=pos, size=0.0, intersect='covers')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for i in range(len(galex_image_table)):\n", - " if (('image/fits' in galex_image_table[i].format) and\n", - " (galex_image_table['energy_bounds_center'][i]==2.35e-07) and\n", - " (galex_image_table['productType'][i] == 'SCIENCE')):\n", - " break\n", - "galex_image_record = galex_image_table[i]\n", - "print(galex_image_record.title, galex_image_record.bandpass_id)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## See above regarding two ways to do this: \n", - "#galex_nuv_image = fits.open(galex_image_record.getdataurl())\n", - "file_name = download_file(galex_image_record.getdataurl(), cache=True) \n", - "galex_nuv_image=fits.open(file_name)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "image_data = galex_nuv_image[0].data\n", - "print('Min:', np.min(image_data))\n", - "print('Max:', np.max(image_data))\n", - "print('Mean:', np.mean(image_data))\n", - "print('Stdev:', np.std(image_data))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig = plt.figure()\n", - "ax = fig.add_subplot(1, 1, 1, projection=WCS(galex_nuv_image[0].header))\n", - "ax.imshow(galex_nuv_image[0].data, cmap='gray_r', origin='lower', vmin=0.0, vmax=0.01)\n", - "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cutout = Cutout2D(galex_nuv_image[0].data, pos, size, wcs=WCS(galex_nuv_image[0].header))\n", - "\n", - "fig = plt.figure()\n", - "\n", - "ax = fig.add_subplot(1, 1, 1, projection=cutout.wcs)\n", - "ax.imshow(cutout.data, cmap='gray_r', origin='lower', vmin = 0.0, vmax = 0.01)\n", - "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 14. Try visualizing a cutout of an SDSS image that covers your position.\n", - "\n", - "Hints:\n", - "* Search the registry using `keywords=['sloan']\n", - "* Find the service with a `short_name` of `'SDSS SIAP'`\n", - "* From Known Issues, recall that an empty string must be specified to the `format` parameter dues to a bug in the service.\n", - "* After obtaining your search results, select r-band images using the `.title` attribute of the records that are returned, since `.bandpass_id` is not populated." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sdss_image_services = vo.regsearch(keywords=['sloan'], servicetype='image')\n", - "sdss_image_services.to_table()['ivoid', 'short_name', 'res_title', 'source_value']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Use list comprehension to check each service's short_name attribute.\n", - "# Given the above, we know the first match is the right one. \n", - "sdss_image_service = [s for s in sdss_image_services if 'SDSS SIAP' in s.short_name ][0]\n", - "sdss_image_service.short_name" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As a workaround to a bug in the SDSS service, pass `format=None` when searching." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "sdss_image_table = sdss_image_service.search(pos=pos, size=0.0, format=None, intersect='covers')\n", - "len(sdss_image_table['Title'])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for sdss_rband_record in sdss_image_table:\n", - " if 'Sloan Digital Sky Survey - Filter r' in sdss_rband_record.title:\n", - " break\n", - "print(sdss_rband_record.title, sdss_rband_record.bandpass_id)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## See above regarding two ways to do this\n", - "# sdss_rband_image = fits.open(sdss_rband_record.getdataurl())\n", - "file_name = download_file(sdss_rband_record.getdataurl(), cache=True) \n", - "sdss_rband_image=fits.open(file_name)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig = plt.figure()\n", - "\n", - "ax = fig.add_subplot(1, 1, 1, projection=WCS(sdss_rband_image[0].header))\n", - "\n", - "interval = vis.PercentileInterval(99.9)\n", - "vmin,vmax = interval.get_limits(sdss_rband_image[0].data)\n", - "norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))\n", - "ax.imshow(sdss_rband_image[0].data, cmap = 'gray_r', norm = norm, origin = 'lower') \n", - "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cutout = Cutout2D(sdss_rband_image[0].data, pos, size, wcs=WCS(sdss_rband_image[0].header))\n", - "fig = plt.figure()\n", - "ax = fig.add_subplot(1, 1, 1, projection=cutout.wcs)\n", - "vmin,vmax = interval.get_limits(sdss_rband_image[0].data)\n", - "norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))\n", - "ax.imshow(cutout.data, cmap = 'gray_r', norm = norm, origin = 'lower') \n", - "ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## 15. Try looping over all positions and plotting multiwavelength cutouts." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Warning: this cell takes a long time to run! We limit it to the first three galaxies only." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Pick the first 3 galaxies.\n", - "galaxy_subset = galaxies[0:3]\n", - "\n", - "# For each galaxy,\n", - "for galaxy in galaxy_subset:\n", - "\n", - " # Establish the position.\n", - " ra = galaxy['RA']\n", - " dec = galaxy['DEC']\n", - " pos = SkyCoord(ra, dec, unit = 'deg') \n", - " \n", - " # Set up the plot for this position.\n", - " fig = plt.figure(figsize=(20,6))\n", - " plt.suptitle('POSITION = ' + str(ra) + ', ' + str(dec), fontsize=16)\n", - "\n", - " # GALEX\n", - " \n", - " # Find the GALEX images that overlap the position.\n", - " galex_image_table = galex_image_service.search(pos=pos, size=0.25)\n", - " \n", - " # Find the GALEX All-Sky Image Survey (AIS) Near-UV FITS coadd. \n", - " galex_image_record = None\n", - " for record in galex_image_table:\n", - " if (('image/fits' in record.format) and\n", - " (record['energy_bounds_center'] == 2.35e-07) and\n", - " (record['productType'] == 'SCIENCE')):\n", - " galex_image_record = record\n", - " break\n", - " \n", - " if galex_image_record is not None:\n", - " # Create a cutout.\n", - " file_name = download_file(galex_image_record.getdataurl(), cache=True) \n", - " gimage = fits.open(file_name)\n", - " galex_cutout = Cutout2D(gimage[0].data, pos, size, wcs=WCS(gimage[0].header))\n", - "\n", - " # Plot the cutout in the first position of a 1x3 (rowsxcols) grid.\n", - " ax = fig.add_subplot(1, 3, 1, projection=galex_cutout.wcs)\n", - " ax.set_title(galex_image_record.title)\n", - " ax.imshow(galex_cutout.data, cmap='gray_r', origin='lower', vmin = 0.0, vmax = 0.01)\n", - " ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n", - " else:\n", - " # We didn't find a suitable image, so leave that subplot blank.\n", - " ax = fig.add_subplot(1, 3, 1, projection=galex_cutout.wcs)\n", - " ax.set_title('GALEX image not found')\n", - " \n", - " # SDSS\n", - " \n", - " # Find the SDSS images that overlap the position.\n", - " sdss_image_table = sdss_image_service.search(pos=pos, size=0, format=None)\n", - " \n", - " # Find the first SDSS r-band image.\n", - " sdss_rband_record = None\n", - " for record in sdss_image_table:\n", - " if 'Sloan Digital Sky Survey - Filter r' in record.title:\n", - " sdss_rband_record = record\n", - " break\n", - " \n", - " if sdss_rband_record is not None:\n", - " # Create a cutout.\n", - " file_name = download_file(sdss_rband_record.getdataurl(), cache=True) \n", - " sdss_rband_image=fits.open(file_name)\n", - "\n", - " sdss_cutout = Cutout2D(sdss_rband_image[0].data, pos, size,\n", - " wcs=WCS(sdss_rband_image[0].header))\n", - "\n", - " # Plot the cutout in the second position of a 1x3 grid.\n", - " vmin,vmax = interval.get_limits(sdss_cutout.data)\n", - " norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))\n", - " ax = fig.add_subplot(1, 3, 2, projection=sdss_cutout.wcs)\n", - " ax.imshow(sdss_cutout.data, cmap = 'gray_r', norm = norm, origin = 'lower') \n", - " ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n", - " ax.set_title(sdss_rband_record.title)\n", - " else:\n", - " # We didn't find a suitable image, so leave that subplot blank.\n", - " ax = fig.add_subplot(1, 3, 2, projection=galex_cutout.wcs)\n", - " ax.set_title('SDSS rband image not found')\n", - " \n", - " # AllWISE\n", - " \n", - " # Find the AllWISE images that overlap the position.\n", - " allwise_image_table = allwise_image_service.search(pos=pos, size=0)\n", - " \n", - " # Find the first AllWISE W1 channel image.\n", - " allwise_image_record = None\n", - " for record in allwise_image_table:\n", - " if 'W1' in record.bandpass_id:\n", - " allwise_image_record = record\n", - " break\n", - " \n", - " if allwise_image_record is not None:\n", - " # Create a cutout.\n", - " file_name = download_file(allwise_image_record.getdataurl(), cache=True) \n", - " allwise_w1_image=fits.open(file_name)\n", - "\n", - " allwise_cutout = Cutout2D(allwise_w1_image[0].data, pos, (size, size),\n", - " wcs=WCS(allwise_w1_image[0].header))\n", - "\n", - " # Plot the cutout in the third position of a 1x3 grid.\n", - " ax = fig.add_subplot(1, 3, 3, projection=allwise_cutout.wcs)\n", - " ax.imshow(allwise_cutout.data, cmap='gray_r', origin='lower', vmax = 10)\n", - " ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none')\n", - " ax.set_title(allwise_image_record.title)\n", - " else:\n", - " # We didn't find a suitable image, so leave that subplot blank.\n", - " ax = fig.add_subplot(1, 3, 3, projection=galex_cutout.wcs)\n", - " ax.set_title('AllWISE W1 image not found')\n" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.13" - }, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": true, - "toc_window_display": true - }, - "widgets": { - "state": {}, - "version": "1.1.1" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/UseCase_I.md b/UseCase_I.md new file mode 100644 index 0000000..732d4a5 --- /dev/null +++ b/UseCase_I.md @@ -0,0 +1,452 @@ +--- +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.9.13 +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: {} + toc_section_display: true + toc_window_display: true +widgets: + state: {} + version: 1.1.1 +--- + +# Science User Case - Inspecting a Candidate List + +[Ogle et al. (2016)](https://ui.adsabs.harvard.edu/abs/2016ApJ...817..109O/abstract) mined the NASA/IPAC Extragalactic Database (NED) to identify a new type of galaxy: Superluminous Spiral Galaxies. + +Table 1 lists the positions of these Super Spirals. Based on those positions, let's create multiwavelength cutouts for each super spiral to see what is unique about this new class of objects. + ++++ + +## 1. Import the Python modules we'll be using. + +```{code-cell} ipython3 +# Suppress unimportant warnings. +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +warnings.filterwarnings('ignore', '.*RADECSYS=*', append=True) + +import matplotlib.pyplot as plt +import numpy as np + +# For downloading files +from astropy.utils.data import download_file + +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.nddata import Cutout2D +import astropy.visualization as vis +from astropy.wcs import WCS +from astroquery.ipac.ned import Ned + +import pyvo as vo +``` + +The next cell prepares the notebook to display our visualizations. + +```{code-cell} ipython3 +%matplotlib inline +``` + +## 2. Search NED for objects in this paper. + +Insert a Code Cell below by clicking on the "Insert" Menu and choosing "Insert Cell Below". Then consult QuickReference.md to figure out how to use astroquery to search NED for all objects in a paper, based on the refcode of the paper. Inspect the resulting astropy table. + +```{code-cell} ipython3 +:tags: [output_scroll] + +objects_in_paper = Ned.query_refcode('2016ApJ...817..109O') +objects_in_paper.show_in_notebook() +``` + +## 3. Filter the NED results. + +The results from NED will include galaxies, but also other kinds of objects (e.g. galaxy clusters, galaxy groups). Print the 'Type' column to see the full range of classifications and filter the results so that we only keep the galaxies in the list. + +```{code-cell} ipython3 +objects_in_paper['Type'] +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +# Keep only the galaxies from the list +galaxies = objects_in_paper[np.array(objects_in_paper['Type']) == 'G'] + +galaxies.show_in_notebook() +``` + +## 4. Search the NAVO Registry for image resources. + +The paper selected super spirals using WISE, SDSS, and GALEX images. Search the NAVO registry for all image resources, using the 'service_type' search parameter. How many image resources are currently available? + +```{code-cell} ipython3 +image_services = vo.regsearch(servicetype='image') + +print(f'{len(image_services)} result(s) found.') + +image_services.to_table()['ivoid', 'short_name', 'res_title'] +``` + +## 5. Search the NAVO Registry for image resources that will allow you to search for AllWISE images. + +There are hundreds of image resources...too many to quickly read through. Try adding the 'keywords' search parameter to your registry search, and find the image resource you would need to search the AllWISE images. Remember from the Known Issues that 'keywords' must be a list. + +```{code-cell} ipython3 +allwise_image_services = vo.regsearch(servicetype='image', keywords=['allwise']) + +print(f'{len(allwise_image_services)} result(s) found.') + +allwise_image_services.to_table()['ivoid', 'short_name', 'res_title'] +``` + +## 6. Choose the AllWISE image service that you are interested in. + +```{code-cell} ipython3 +allwise_image_service = allwise_image_services[0] +allwise_image_service.service +``` + +## 7. Choose one of the galaxies in the NED list. + +```{code-cell} ipython3 +ra = galaxies['RA'][0] +dec = galaxies['DEC'][0] +pos = SkyCoord(ra, dec, unit = 'deg') +``` + +```{code-cell} ipython3 +ra,dec +``` + +## 8. Search for a list of AllWISE images that cover this galaxy. + +How many images are returned? Which are you most interested in? + +```{code-cell} ipython3 +allwise_image_table = allwise_image_service.search(pos=pos, size=0) +allwise_image_table +``` + +## 9. Use the .to_table() method to view the results as an Astropy table. + +```{code-cell} ipython3 +allwise_images = allwise_image_table.to_table() +allwise_images +``` + +## 10. From the result in 8., select the first record for an image taken in WISE band W1 (3.6 micron) + +Hints: +* Loop over records and test on the `.bandpass_id` attribute of each record +* Print the `.title` and `.bandpass_id` of the record you find, to verify it is the right one. + +```{code-cell} ipython3 +for allwise_image_record in allwise_image_table: + if 'W1' in allwise_image_record.bandpass_id: + break +print(allwise_image_record.title, allwise_image_record.bandpass_id) +``` + +## 11. Visualize this AllWISE image. + +```{code-cell} ipython3 +## If you only run this once, you can do it in memory in one line: +## This fetches the FITS as an astropy.io.fits object in memory +#allwise_w1_image = allwise_image_record.getdataobj() +## But if you might run this notebook repeatedly with limited bandwidth, +## download it once and cache it. +file_name = download_file(allwise_image_record.getdataurl(), cache=True) +allwise_w1_image = fits.open(file_name) +``` + +```{code-cell} ipython3 +fig = plt.figure() + +wcs = WCS(allwise_w1_image[0].header) +ax = fig.add_subplot(1, 1, 1, projection=wcs) +ax.imshow(allwise_w1_image[0].data, cmap='gray_r', origin='lower', vmax = 10) +ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') +``` + +## 12. Plot a cutout of the AllWISE image, centered on your position. + +Try a 60 arcsecond cutout. + +```{code-cell} ipython3 +size = 60 +cutout = Cutout2D(allwise_w1_image[0].data, pos, (size, size), wcs=wcs) +wcs = cutout.wcs + +fig = plt.figure() + +ax = fig.add_subplot(1, 1, 1, projection=wcs) +ax.imshow(cutout.data, cmap='gray_r', origin='lower', vmax = 10) +ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') +``` + +## 13. Try visualizing a cutout of a GALEX image that covers your position. + +Repeat steps 4, 5, 6, 8 through 12 for GALEX. + +```{code-cell} ipython3 +galex_image_services = vo.regsearch(keywords=['galex'], servicetype='image') +print(f'{len(galex_image_services)} result(s) found.') +galex_image_services.to_table()['ivoid', 'short_name', 'res_title'] +``` + +```{code-cell} ipython3 +galex_image_service = galex_image_services[0] +``` + +```{code-cell} ipython3 +galex_image_table = galex_image_service.search(pos=pos, size=0.0, intersect='covers') +``` + +```{code-cell} ipython3 +for i in range(len(galex_image_table)): + if (('image/fits' in galex_image_table[i].format) and + (galex_image_table['energy_bounds_center'][i]==2.35e-07) and + (galex_image_table[i]['productType'] == 'SCIENCE')): + break +galex_image_record = galex_image_table[i] +print(galex_image_record.title, galex_image_record.bandpass_id) +``` + +```{code-cell} ipython3 +## See above regarding two ways to do this: +#galex_nuv_image = fits.open(galex_image_record.getdataurl()) +file_name = download_file(galex_image_record.getdataurl(), cache=True) +galex_nuv_image=fits.open(file_name) +``` + +```{code-cell} ipython3 +image_data = galex_nuv_image[0].data +print('Min:', np.min(image_data)) +print('Max:', np.max(image_data)) +print('Mean:', np.mean(image_data)) +print('Stdev:', np.std(image_data)) +``` + +```{code-cell} ipython3 +fig = plt.figure() +ax = fig.add_subplot(1, 1, 1, projection=WCS(galex_nuv_image[0].header)) +ax.imshow(galex_nuv_image[0].data, cmap='gray_r', origin='lower', vmin=0.0, vmax=0.01) +ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') +``` + +```{code-cell} ipython3 +cutout = Cutout2D(galex_nuv_image[0].data, pos, size, wcs=WCS(galex_nuv_image[0].header)) + +fig = plt.figure() + +ax = fig.add_subplot(1, 1, 1, projection=cutout.wcs) +ax.imshow(cutout.data, cmap='gray_r', origin='lower', vmin = 0.0, vmax = 0.01) +ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') +``` + +## 14. Try visualizing a cutout of an SDSS image that covers your position. + +Hints: +* Search the registry using `keywords=['sloan'] +* Find the service with a `short_name` of `'SDSS SIAP'` +* From Known Issues, recall that an empty string must be specified to the `format` parameter dues to a bug in the service. +* After obtaining your search results, select r-band images using the `.title` attribute of the records that are returned, since `.bandpass_id` is not populated. + +```{code-cell} ipython3 +sdss_image_services = vo.regsearch(keywords=['sloan'], servicetype='image') +sdss_image_services.to_table()['ivoid', 'short_name', 'res_title', 'source_value'] +``` + +```{code-cell} ipython3 +# Use list comprehension to check each service's short_name attribute. +# Given the above, we know the first match is the right one. +sdss_image_service = [s for s in sdss_image_services if 'SDSS SIAP' in s.short_name ][0] +sdss_image_service.short_name +``` + +As a workaround to a bug in the SDSS service, pass `format=None` when searching. + +```{code-cell} ipython3 +sdss_image_table = sdss_image_service.search(pos=pos, size=0.0, format=None, intersect='covers') +len(sdss_image_table['Title']) +``` + +```{code-cell} ipython3 +for sdss_rband_record in sdss_image_table: + if 'Sloan Digital Sky Survey - Filter r' in sdss_rband_record.title: + break +print(sdss_rband_record.title, sdss_rband_record.bandpass_id) +``` + +```{code-cell} ipython3 +## See above regarding two ways to do this +# sdss_rband_image = fits.open(sdss_rband_record.getdataurl()) +file_name = download_file(sdss_rband_record.getdataurl(), cache=True) +sdss_rband_image=fits.open(file_name) +``` + +```{code-cell} ipython3 +fig = plt.figure() + +ax = fig.add_subplot(1, 1, 1, projection=WCS(sdss_rband_image[0].header)) + +interval = vis.PercentileInterval(99.9) +vmin,vmax = interval.get_limits(sdss_rband_image[0].data) +norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000)) +ax.imshow(sdss_rband_image[0].data, cmap = 'gray_r', norm = norm, origin = 'lower') +ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') +``` + +```{code-cell} ipython3 +cutout = Cutout2D(sdss_rband_image[0].data, pos, size, wcs=WCS(sdss_rband_image[0].header)) +fig = plt.figure() +ax = fig.add_subplot(1, 1, 1, projection=cutout.wcs) +vmin,vmax = interval.get_limits(sdss_rband_image[0].data) +norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000)) +ax.imshow(cutout.data, cmap = 'gray_r', norm = norm, origin = 'lower') +ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') +``` + +## 15. Try looping over all positions and plotting multiwavelength cutouts. + ++++ + +Warning: this cell takes a long time to run! We limit it to the first three galaxies only. + +```{code-cell} ipython3 +# Pick the first 3 galaxies. +galaxy_subset = galaxies[0:3] + +# For each galaxy, +for galaxy in galaxy_subset: + + # Establish the position. + ra = galaxy['RA'] + dec = galaxy['DEC'] + pos = SkyCoord(ra, dec, unit = 'deg') + + # Set up the plot for this position. + fig = plt.figure(figsize=(20,6)) + plt.suptitle('POSITION = ' + str(ra) + ', ' + str(dec), fontsize=16) + + # GALEX + + # Find the GALEX images that overlap the position. + galex_image_table = galex_image_service.search(pos=pos, size=0.25) + + # Find the GALEX All-Sky Image Survey (AIS) Near-UV FITS coadd. + galex_image_record = None + for record in galex_image_table: + if (('image/fits' in record.format) and + (record['energy_bounds_center'] == 2.35e-07) and + (record['productType'] == 'SCIENCE')): + galex_image_record = record + break + + if galex_image_record is not None: + # Create a cutout. + file_name = download_file(galex_image_record.getdataurl(), cache=True) + gimage = fits.open(file_name) + galex_cutout = Cutout2D(gimage[0].data, pos, size, wcs=WCS(gimage[0].header)) + + # Plot the cutout in the first position of a 1x3 (rowsxcols) grid. + ax = fig.add_subplot(1, 3, 1, projection=galex_cutout.wcs) + ax.set_title(galex_image_record.title) + ax.imshow(galex_cutout.data, cmap='gray_r', origin='lower', vmin = 0.0, vmax = 0.01) + ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') + else: + # We didn't find a suitable image, so leave that subplot blank. + ax = fig.add_subplot(1, 3, 1, projection=galex_cutout.wcs) + ax.set_title('GALEX image not found') + + # SDSS + + # Find the SDSS images that overlap the position. + sdss_image_table = sdss_image_service.search(pos=pos, size=0, format=None) + + # Find the first SDSS r-band image. + sdss_rband_record = None + for record in sdss_image_table: + if 'Sloan Digital Sky Survey - Filter r' in record.title: + sdss_rband_record = record + break + + if sdss_rband_record is not None: + # Create a cutout. + file_name = download_file(sdss_rband_record.getdataurl(), cache=True) + sdss_rband_image=fits.open(file_name) + + sdss_cutout = Cutout2D(sdss_rband_image[0].data, pos, size, + wcs=WCS(sdss_rband_image[0].header)) + + # Plot the cutout in the second position of a 1x3 grid. + vmin,vmax = interval.get_limits(sdss_cutout.data) + norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000)) + ax = fig.add_subplot(1, 3, 2, projection=sdss_cutout.wcs) + ax.imshow(sdss_cutout.data, cmap = 'gray_r', norm = norm, origin = 'lower') + ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') + ax.set_title(sdss_rband_record.title) + else: + # We didn't find a suitable image, so leave that subplot blank. + ax = fig.add_subplot(1, 3, 2, projection=galex_cutout.wcs) + ax.set_title('SDSS rband image not found') + + # AllWISE + + # Find the AllWISE images that overlap the position. + allwise_image_table = allwise_image_service.search(pos=pos, size=0) + + # Find the first AllWISE W1 channel image. + allwise_image_record = None + for record in allwise_image_table: + if 'W1' in record.bandpass_id: + allwise_image_record = record + break + + if allwise_image_record is not None: + # Create a cutout. + file_name = download_file(allwise_image_record.getdataurl(), cache=True) + allwise_w1_image=fits.open(file_name) + + allwise_cutout = Cutout2D(allwise_w1_image[0].data, pos, (size, size), + wcs=WCS(allwise_w1_image[0].header)) + + # Plot the cutout in the third position of a 1x3 grid. + ax = fig.add_subplot(1, 3, 3, projection=allwise_cutout.wcs) + ax.imshow(allwise_cutout.data, cmap='gray_r', origin='lower', vmax = 10) + ax.scatter(ra, dec, transform=ax.get_transform('fk5'), s=500, edgecolor='red', facecolor='none') + ax.set_title(allwise_image_record.title) + else: + # We didn't find a suitable image, so leave that subplot blank. + ax = fig.add_subplot(1, 3, 3, projection=galex_cutout.wcs) + ax.set_title('AllWISE W1 image not found') +``` diff --git a/UseCase_II.ipynb b/UseCase_II.ipynb deleted file mode 100644 index bb03952..0000000 --- a/UseCase_II.ipynb +++ /dev/null @@ -1,434 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Preparing a proposal\n", - "\n", - "The Story: Suppose that you are preparing to write a proposal on NGC1365, aiming to investigate the intriguing black hole spin this galaxy with Chandra grating observations (see: [Monster Blackhole Spin Revealed](https://www.space.com/19980-monster-black-hole-spin-discovery.html)) \n", - "\n", - "In writing proposals, there are often the same tasks that are required: including finding and analyzing previous observations of the proposal, and creating figures that include, e.g., multiwavelength images and spectrum for the source. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# As a hint, we include the code block for Python modules that you will likely need to import: \n", - "import matplotlib\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "%matplotlib inline \n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "from astropy.io import fits\n", - "\n", - "import pyvo as vo\n", - "\n", - "## There are a number of relatively unimportant warnings that \n", - "## show up, so for now, suppress them:\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 1: Find out what the previously quoted Chandra 2-10 keV flux of the central source is for NGC 1365. \n", - "\n", - "Hint: Do a Registry search for tables served by the HEASARC (where high energy data are archived) to find potential table with this information" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# This gets all services matching, which should be a list of only one:\n", - "tap_services=vo.regsearch(servicetype='table',keywords=['heasarc'])\n", - "# This fetches the list of tables that this service serves:\n", - "heasarc_tables=tap_services[0].service.tables" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint: The [Chansngcat](https://heasarc.gsfc.nasa.gov/W3Browse/chandra/chansngcat.html) table is likely the best table. Create a table with ra, dec, exposure time, and flux (and flux errors) from the public.chansngcat catalog for Chandra observations matched within 0.1 degree." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "for tablename in heasarc_tables.keys():\n", - " if \"chansng\" in tablename: \n", - " print(\"Table {} has columns={}\\n\".format(\n", - " tablename,\n", - " sorted([k.name for k in heasarc_tables[tablename].columns ])))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get the coordinate for NGC 1365\n", - "import astropy.coordinates as coord\n", - "pos=coord.SkyCoord.from_name(\"ngc1365\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Construct a query that will get the ra, dec, exposure time, flux, and flux errors \n", - "# from this catalog in the region around this source:\n", - "query=\"\"\"SELECT ra, dec, exposure, flux, flux_lower, flux_upper FROM public.chansngcat as cat \n", - " where contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{},{},0.1))=1 \n", - " and cat.exposure > 0 order by cat.exposure\"\"\".format(pos.ra.deg, pos.dec.deg)\n", - "# Submit the query. (See the CS_Catalog_queries.ipynb for\n", - "# information about these two search options.)\n", - "results=tap_services[0].service.run_async(query)\n", - "#results=tap_services[0].search(query)\n", - "# Look at the results\n", - "results.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "## Step 2: Make Images: \n", - "\n", - "### Create ultraviolet and X-ray images\n", - "Hint: Start by checking what UV image services exist (e.g., GALEX?)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Note that to browse the columns, use the .to_table() method\n", - "uv_services=vo.regsearch(servicetype='image',keywords='galex', waveband='uv')\n", - "uv_services.to_table()['ivoid','short_name']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The keyword search for 'galex' returned a bunch of things that may have mentioned it, but let's just use the ones that have GALEX as their short name:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "uv_services.to_table()[\n", - " np.array(['GALEX' in u.short_name for u in uv_services])\n", - " ]['ivoid', 'short_name']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Though using the result as an Astropy Table makes it easier to look at the contents, to call the service itself, we cannot use the row of that table. You have to use the entry in the service result list itself. So use the table to browse, but select the list of services itself using the properties that have been defined as attributes such as short_name and ivoid:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "galex_stsci=[s for s in uv_services if 'GALEX' in s.short_name and 'stsci' in s.ivoid][0]\n", - "galex_heasarc=[s for s in uv_services if 'GALEX' in s.short_name and 'heasarc' in s.ivoid][0]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint: Next create a UV image for the source " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Do an image search for NGC 1365 in the UV service found above\n", - "im_table_stsci=galex_stsci.search(pos=pos,size=0.1)\n", - "im_table_stsci.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Let's see what HEASARC offers, and this time limit it to FITS \n", - "# this option doesn't currently work for STScI's service)\n", - "im_table_heasarc=galex_heasarc.search(pos=pos,size=0.1,format='image/fits')\n", - "im_table_heasarc.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## If you only run this once, you can do it in memory in one line:\n", - "## This fetches the FITS as an astropy.io.fits object in memory\n", - "#dataobj=im_table_heasarc[0].getdataobj()\n", - "## But if you might run this notebook repeatedly with limited bandwidth, \n", - "## download it once and cache it. \n", - "file_name = download_file(im_table_heasarc[0].getdataurl(), cache=True, timeout=600)\n", - "dataobj=fits.open(file_name)\n", - "print(type(dataobj))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get the FITS file (which is index 0 for the NUV image or index=2 for the FUV image)\n", - "from pylab import figure, cm\n", - "from matplotlib.colors import LogNorm\n", - "plt.matshow(dataobj[0].data, origin='lower', cmap=cm.gray_r, norm=LogNorm(vmin=0.005, vmax=0.3))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint: Repeat steps for X-ray image. (Note: Ideally, we would find an image in the Chandra 'cxc' catalog) " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "x_services=vo.regsearch(servicetype='image',keywords=['chandra'], waveband='x-ray')\n", - "print(x_services.to_table()['short_name','ivoid'])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Do an image search for NGC 1365 in the X-ray CDA service found above\n", - "xim_table=x_services[0].search(pos=pos,size=0.2)\n", - "## Some of these are FITS and some JPEG. Look at the columns:\n", - "print( xim_table.to_table().columns )\n", - "first_fits_image_row = [x for x in xim_table if 'image/fits' in x.format][0] " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Create an image from the first FITS file (index=1) by downloading:\n", - "## See above for options\n", - "#xhdu_list=first_fits_image_row.getdataobj()\n", - "file_name = download_file(first_fits_image_row.getdataurl(), cache=True, timeout=600)\n", - "xhdu_list=fits.open(file_name)\n", - "\n", - "\n", - "plt.imshow(xhdu_list[0].data, origin='lower', cmap='cool', norm=LogNorm(vmin=0.1, vmax=500.))\n", - "plt.xlim(460, 560)\n", - "plt.ylim(460, 560)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "## Step 3: Make a spectrum: \n", - "\n", - "### Find what Chandra spectral observations exist already for this source. \n", - "Hint: try searching for X-ray spectral data tables using the registry query" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Use the TAP protocol to list services that contain X-ray spectral data\n", - "xsp_services=vo.regsearch(servicetype='ssa',waveband='x-ray')\n", - "xsp_services.to_table()['short_name','ivoid','waveband']" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint 2: Take a look at what data exist for our candidate, NGC 1365." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spec_tables=xsp_services[0].search(pos=pos,radius=0.2,verbose=True)\n", - "spec_tables.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Hint 3: Download the data to make a spectrum. Note: you might end here and use Xspec to plot and model the spectrum. Or ... you can also try to take a quick look at the spectrum. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Get it and look at it:\n", - "#hdu_list=spec_tables[0].getdataobj()\n", - "file_name = download_file(spec_tables[0].getdataurl(), cache=True, timeout=600)\n", - "hdu_list=fits.open(file_name)\n", - "\n", - "spectra=hdu_list[1].data\n", - "print(spectra.columns)\n", - "print(len(spectra))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## Or write it to disk\n", - "import os\n", - "if not os.path.isdir('downloads'):\n", - " os.makedirs(\"downloads\")\n", - "fname=spec_tables[0].make_dataset_filename()\n", - "# Known issue where the suffix is incorrect:\n", - "fname=fname.replace('None','fits')\n", - "with open('downloads/{}'.format(fname),'wb') as outfile:\n", - " outfile.write(spec_tables[0].getdataset().read())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Extension: Making a \"quick look\" spectrum. For our purposes, the 1st order of the HEG grating data would be sufficient." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "j=1\n", - "for i in range(len(spectra)):\n", - " matplotlib.rcParams['figure.figsize'] = (8, 3)\n", - " if abs(spectra['TG_M'][i]) == 1 and (spectra['TG_PART'][i]) == 1:\n", - " ax=plt.subplot(1,2,j)\n", - " pha = plt.plot( spectra['CHANNEL'][i],spectra['COUNTS'][i])\n", - " ax.set_yscale('log')\n", - " if spectra['TG_PART'][i] == 1:\n", - " instr='HEG'\n", - " ax.set_title(\"{grating}{order:+d}\".format(grating=instr, order=spectra['TG_M'][i]))\n", - " plt.tight_layout()\n", - " j=j+1\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html). (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.ipynb) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Congratulations! You have completed this notebook exercise." - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.7" - }, - "nav_menu": {}, - "toc": { - "base_numbering": 1, - "nav_menu": {}, - "number_sections": false, - "sideBar": true, - "skip_h1_title": false, - "title_cell": "Table of Contents", - "title_sidebar": "Contents", - "toc_cell": false, - "toc_position": {}, - "toc_section_display": "block", - "toc_window_display": true - }, - "widgets": { - "state": {}, - "version": "1.1.1" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/UseCase_II.md b/UseCase_II.md new file mode 100644 index 0000000..da865f9 --- /dev/null +++ b/UseCase_II.md @@ -0,0 +1,260 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.9.7 +nav_menu: {} +toc: + base_numbering: 1 + nav_menu: {} + number_sections: false + sideBar: true + skip_h1_title: false + title_cell: Table of Contents + title_sidebar: Contents + toc_cell: false + toc_position: {} + toc_section_display: block + toc_window_display: true +widgets: + state: {} + version: 1.1.1 +--- + +# Preparing a proposal + +The Story: Suppose that you are preparing to write a proposal on NGC1365, aiming to investigate the intriguing black hole spin this galaxy with Chandra grating observations (see: [Monster Blackhole Spin Revealed](https://www.space.com/19980-monster-black-hole-spin-discovery.html)) + +In writing proposals, there are often the same tasks that are required: including finding and analyzing previous observations of the proposal, and creating figures that include, e.g., multiwavelength images and spectrum for the source. + +```{code-cell} ipython3 +# As a hint, we include the code block for Python modules that you will likely need to import: +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +%matplotlib inline + +# For downloading files +from astropy.utils.data import download_file +from astropy.io import fits + +import pyvo as vo + +## There are a number of relatively unimportant warnings that +## show up, so for now, suppress them: +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +``` + +## Step 1: Find out what the previously quoted Chandra 2-10 keV flux of the central source is for NGC 1365. + +Hint: Do a Registry search for tables served by the HEASARC (where high energy data are archived) to find potential table with this information + +```{code-cell} ipython3 +# This gets all services matching, which should be a list of only one: +tap_services=vo.regsearch(servicetype='table',keywords=['heasarc']) +# This fetches the list of tables that this service serves: +heasarc_tables=tap_services[0].service.tables +``` + +Hint: The [Chansngcat](https://heasarc.gsfc.nasa.gov/W3Browse/chandra/chansngcat.html) table is likely the best table. Create a table with ra, dec, exposure time, and flux (and flux errors) from the public.chansngcat catalog for Chandra observations matched within 0.1 degree. + +```{code-cell} ipython3 +for tablename in heasarc_tables.keys(): + if "chansng" in tablename: + print("Table {} has columns={}\n".format( + tablename, + sorted([k.name for k in heasarc_tables[tablename].columns ]))) +``` + +```{code-cell} ipython3 +# Get the coordinate for NGC 1365 +import astropy.coordinates as coord +pos=coord.SkyCoord.from_name("ngc1365") +``` + +```{code-cell} ipython3 +# Construct a query that will get the ra, dec, exposure time, flux, and flux errors +# from this catalog in the region around this source: +query="""SELECT ra, dec, exposure, flux, flux_lower, flux_upper FROM public.chansngcat as cat + where contains(point('ICRS',cat.ra,cat.dec),circle('ICRS',{},{},0.1))=1 + and cat.exposure > 0 order by cat.exposure""".format(pos.ra.deg, pos.dec.deg) +# Submit the query. (See the CS_Catalog_queries.md for +# information about these two search options.) +results=tap_services[0].service.run_async(query) +#results=tap_services[0].search(query) +# Look at the results +results.to_table() +``` + +## Step 2: Make Images: + +### Create ultraviolet and X-ray images +Hint: Start by checking what UV image services exist (e.g., GALEX?) + +```{code-cell} ipython3 +## Note that to browse the columns, use the .to_table() method +uv_services=vo.regsearch(servicetype='image',keywords='galex', waveband='uv') +uv_services.to_table()['ivoid','short_name'] +``` + +The keyword search for 'galex' returned a bunch of things that may have mentioned it, but let's just use the ones that have GALEX as their short name: + +```{code-cell} ipython3 +uv_services.to_table()[ + np.array(['GALEX' in u.short_name for u in uv_services]) + ]['ivoid', 'short_name'] +``` + +Though using the result as an Astropy Table makes it easier to look at the contents, to call the service itself, we cannot use the row of that table. You have to use the entry in the service result list itself. So use the table to browse, but select the list of services itself using the properties that have been defined as attributes such as short_name and ivoid: + +```{code-cell} ipython3 +galex_stsci=[s for s in uv_services if 'GALEX' in s.short_name and 'stsci' in s.ivoid][0] +galex_heasarc=[s for s in uv_services if 'GALEX' in s.short_name and 'heasarc' in s.ivoid][0] +``` + +Hint: Next create a UV image for the source + +```{code-cell} ipython3 +# Do an image search for NGC 1365 in the UV service found above +im_table_stsci=galex_stsci.search(pos=pos,size=0.1) +im_table_stsci.to_table() +``` + +```{code-cell} ipython3 +# Let's see what HEASARC offers, and this time limit it to FITS +# this option doesn't currently work for STScI's service) +im_table_heasarc=galex_heasarc.search(pos=pos,size=0.1,format='image/fits') +im_table_heasarc.to_table() +``` + +```{code-cell} ipython3 +## If you only run this once, you can do it in memory in one line: +## This fetches the FITS as an astropy.io.fits object in memory +#dataobj=im_table_heasarc[0].getdataobj() +## But if you might run this notebook repeatedly with limited bandwidth, +## download it once and cache it. +file_name = download_file(im_table_heasarc[0].getdataurl(), cache=True, timeout=600) +dataobj=fits.open(file_name) +print(type(dataobj)) +``` + +```{code-cell} ipython3 +# Get the FITS file (which is index 0 for the NUV image or index=2 for the FUV image) +from pylab import figure, cm +from matplotlib.colors import LogNorm +plt.matshow(dataobj[0].data, origin='lower', cmap=cm.gray_r, norm=LogNorm(vmin=0.005, vmax=0.3)) +``` + +Hint: Repeat steps for X-ray image. (Note: Ideally, we would find an image in the Chandra 'cxc' catalog) + +```{code-cell} ipython3 +x_services=vo.regsearch(servicetype='image',keywords=['chandra'], waveband='x-ray') +print(x_services.to_table()['short_name','ivoid']) +``` + +```{code-cell} ipython3 +## Do an image search for NGC 1365 in the X-ray CDA service found above +xim_table=x_services[0].search(pos=pos,size=0.2) +## Some of these are FITS and some JPEG. Look at the columns: +print( xim_table.to_table().columns ) +first_fits_image_row = [x for x in xim_table if 'image/fits' in x.format][0] +``` + +```{code-cell} ipython3 +## Create an image from the first FITS file (index=1) by downloading: +## See above for options +#xhdu_list=first_fits_image_row.getdataobj() +file_name = download_file(first_fits_image_row.getdataurl(), cache=True, timeout=600) +xhdu_list=fits.open(file_name) + + +plt.imshow(xhdu_list[0].data, origin='lower', cmap='cool', norm=LogNorm(vmin=0.1, vmax=500.)) +plt.xlim(460, 560) +plt.ylim(460, 560) +``` + +## Step 3: Make a spectrum: + +### Find what Chandra spectral observations exist already for this source. +Hint: try searching for X-ray spectral data tables using the registry query + +```{code-cell} ipython3 +# Use the TAP protocol to list services that contain X-ray spectral data +xsp_services=vo.regsearch(servicetype='ssa',waveband='x-ray') +xsp_services.to_table()['short_name','ivoid','waveband'] +``` + +Hint 2: Take a look at what data exist for our candidate, NGC 1365. + +```{code-cell} ipython3 +spec_tables=xsp_services[0].search(pos=pos,radius=0.2,verbose=True) +spec_tables.to_table() +``` + +Hint 3: Download the data to make a spectrum. Note: you might end here and use Xspec to plot and model the spectrum. Or ... you can also try to take a quick look at the spectrum. + +```{code-cell} ipython3 +# Get it and look at it: +#hdu_list=spec_tables[0].getdataobj() +file_name = download_file(spec_tables[0].getdataurl(), cache=True, timeout=600) +hdu_list=fits.open(file_name) + +spectra=hdu_list[1].data +print(spectra.columns) +print(len(spectra)) +``` + +```{code-cell} ipython3 +## Or write it to disk +import os +if not os.path.isdir('downloads'): + os.makedirs("downloads") +fname=spec_tables[0].make_dataset_filename() +# Known issue where the suffix is incorrect: +fname=fname.replace('None','fits') +with open('downloads/{}'.format(fname),'wb') as outfile: + outfile.write(spec_tables[0].getdataset().read()) +``` + +Extension: Making a "quick look" spectrum. For our purposes, the 1st order of the HEG grating data would be sufficient. + +```{code-cell} ipython3 +j=1 +for i in range(len(spectra)): + matplotlib.rcParams['figure.figsize'] = (8, 3) + if abs(spectra['TG_M'][i]) == 1 and (spectra['TG_PART'][i]) == 1: + ax=plt.subplot(1,2,j) + pha = plt.plot( spectra['CHANNEL'][i],spectra['COUNTS'][i]) + ax.set_yscale('log') + if spectra['TG_PART'][i] == 1: + instr='HEG' + ax.set_title("{grating}{order:+d}".format(grating=instr, order=spectra['TG_M'][i])) + plt.tight_layout() + j=j+1 +``` + +This can then be analyzed in your favorite spectral analysis tool, e.g., [pyXspec](https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html). (For the winter 2018 AAS workshop, we demonstrated this in a [notebook](https://github.com/NASA-NAVO/aas_workshop_2018/blob/master/heasarc/heasarc_Spectral_Access.md) that you can consult for how to use pyXspec, but the pyXspec documentation will have more information.) + ++++ + +Congratulations! You have completed this notebook exercise. diff --git a/UseCase_III.ipynb b/UseCase_III.ipynb deleted file mode 100644 index de029dc..0000000 --- a/UseCase_III.ipynb +++ /dev/null @@ -1,661 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Creating a stellar color-magnitude (or Hertzsprung-Russell) diagram\n", - "\n", - "The [Hertzsprung-Russell diagram](https://en.wikipedia.org/wiki/Hertzsprung–Russell_diagram) is a fundamental diagram in astronomy that displays important relationships between the stellar color (or temperature) and absolute brightness (or luminosity). \n", - "\n", - "In this exercise, we will use existing stellar catalogs to produce the H-R diagram. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# As a hint, we include the code block for Python modules that you will likely need to import: \n", - "import matplotlib\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "%matplotlib inline \n", - "\n", - "# For downloading files\n", - "from astropy.utils.data import download_file\n", - "from astropy.io import fits\n", - "\n", - "import pyvo as vo\n", - "from pyvo import registry\n", - "\n", - "## There are a number of relatively unimportant warnings that \n", - "## show up, so for now, suppress them:\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n", - "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.*\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 1: Find appropriate catalogs\n", - "\n", - "We want to find a star catalog that has the available data to produce the H-R diagram, i.e., the absolute magnitudes (or both apparent magnitudes AND distances, so we can calculate the absolute magnitudes) in two optical bands (e.g., B and V). This would give us color. Or we need B- OR V- band magnitude and the stellar temperature. \n", - "\n", - "To simplify this problem, we want to find a catalog of an open cluster of stars, where all the stars were born around the same time and are located in one cluster. This simplifies the issue of getting accurate distances to the stars. One famous cluster is the Pleiades in the constellation of Taurus. So first we start by searching for an existing catalog with data on Pleiades that will provide the necessary information about the stars: magnitudes in two bands (e.g., B and V), which can be used to measure color, or temperature of the star plus one magnitude. \n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### DATA DISCOVERY STEPS: \n", - "\n", - "Here is useful link for [how the pyvo registry search works](https://pyvo.readthedocs.io/en/latest/registry/index.html)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "tap_services = registry.search(servicetype = 'tap', keywords=['star pleiades'], includeaux=True)\n", - "print(len(tap_services))\n", - "tap_services.get_summary()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note: The includeaux=True includes auxiliary services. \n", - "\n", - "Because we specified the service type, this returns only the TAP services for each of the avalible resorces matching our search criteria. So, the 'interfaces' column will only show TAP service as an avalible service, even if other services are avalible from the same resource. We know we need the service to be a TAP service since we know that we want to eventually access the search using additional column information (i.e., beyond RA and Dec, which is all that Simple Cone Search returns). " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Next, we need to find which of these has the columns of interest, i.e. magnitudes in two bands to create the color-magnitude diagram. \n", - "\n", - "We can re-run the registry search, but further restrict the results by column UCD. We want tables that have magnitude columns; the most basic UCD to describe a magnitude column is phot.mag" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "tap_services = registry.search(servicetype='tap', keywords=['star pleiades'], \n", - " ucd = ['phot.mag%'], includeaux=True)\n", - "print(len(tap_services))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note: the '%' serves as a wild card when searching by UCD\n", - "\n", - "The IOVA standard enables resources to be as sepecific as they would like when defining the UCD of columns. For example, 'phot.mag' and 'phot.mag;em.opt.V' can both be used to describe a column containing the V magnitudes of objects. If a resource uses the latter to describe a column, a search using 'phot.mag' will not return that resource. A wild card would need to be used or the exact column UCD. The UCD search requires an exact match for a resource to be returned, so using the wild card will make it easier to discover a wider variety of resources. " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - "So using this we can reduce the matched tables to ones that are a bit more catered to our experiment. Note, that there can be redundancy in some resources since these are available via multiple services and/or publishers. Therefore a bit more cleaning can be done to provide only the unique matches. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "tap_services.to_table()['ivoid']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "def getunique( result ):\n", - " short_name = []\n", - " unique_ind = []\n", - " for i in range(len(result)):\n", - " short = result[i].short_name \n", - " if short not in short_name: \n", - " short_name.append(short)\n", - " unique_ind.append(i)\n", - " else:\n", - " print(i)\n", - " \n", - " return(unique_ind)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "uniq_ind=getunique(tap_services)\n", - "print(len(uniq_ind))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "tap_services.to_table()[uniq_ind]['ivoid']" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - "This shows that in this case, all of our TAP results are unique. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can read more information about the results we found. For each resource element (i.e. row in the table above), there are useful attributes, which are [described here]( https://pyvo.readthedocs.io/en/latest/api/pyvo.registry.regtap.RegistryResource.html#pyvo.registry.regtap.RegistryResource)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# To read the descriptions of the resulting matches:\n", - "\n", - "for i in uniq_ind: \n", - " print(\" *** \\n\")\n", - " print(tap_services[i].creators)\n", - " print(tap_services[i].res_description)\n", - " \n", - " " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - " RESULT: Based on these, the second one (by Eichhorn et al) looks like a good start. \n", - "\n", - "### At this point, you can proceed to Step 2. \n", - "\n", - "-- OR --\n", - "\n", - "### Try a different data discovery method! " - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - "### Alternative Method: Use ADS to search for appropriate paper and access data via NED.\n", - "\n", - "There are multiple paths for the data discovery. So it may also be that you know the paper that has the data you are interested in and want to access via the bibcode or authors, etc. \n", - "\n", - "In this case, let's assume that we have the information that the Eichhorn+1970 paper has the data that we need to create the H-R diagram: https://ui.adsabs.harvard.edu/abs/1970MmRAS..73..125E/abstract\n", - "\n", - "We can either search by bibcode (1970MmRAS..73..125E) or \"Eichhorn\" to get the access_urls that will allow us to work with the data. \n", - "\n", - "Before this step, if may help to see the names of the fields available to use. Notice the following fields: \n", - "\n", - "\"source_value\" contains the bibcode information that we want; \"creator_seq\" lists the authors; \n", - "\n", - "and \n", - "\n", - "\"access_url\" provides the url from where the data can be accessed. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "## You already have this from above:\n", - "# tap_services = registry.search(servicetype='tap', keywords=['star pleiades'], \n", - "# ucd = ['phot.mag%'], includeaux=True)\n", - "print(tap_services.fieldnames)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "First, Try using bibcode: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bibcode = '1970MmRAS..73..125E' # Eichhorn\n", - "idx=-1\n", - "for s in tap_services:\n", - " idx+=1\n", - " if bibcode in s['source_value']:\n", - " myidx=idx\n", - " print(f\"{s.short_name}, {s.source_value}, {s.access_url}\")\n", - " break \n", - "# Note that using the to_table() lets you search the result \n", - "# easily using all columns. But in the end, you want to get\n", - "# back not an astropy table row, which you cannot use, but the\n", - "# original RegistryResult that has the callable TAP service. \n", - "myTAP=tap_services[myidx]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that the URL is a generic TAP url for Vizier. All of its tables can be accessed by that same TAP services. It'll be in the ADQL query itself that you specify the table name. We'll see this below." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, try using Author name: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "author = 'Eichhorn'\n", - "\n", - "for record in tap_services: \n", - " names=record.creators\n", - " if 'Eichhorn' in names[0]: \n", - " print(\"For %s: \" %record.short_name)\n", - " print(\" Access URL: %s\" %record['access_urls'][0]) \n", - " print(\" Reference URL: %s\" %record.reference_url)\n", - "\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "In the code above, the record is a Registry Resource. You can access the attribute, \"creators\", from the resource, which is relevant for our example here since this is a direct way to get the author names. The other attributes, \"access_url\" and \"reference_url\", provides two types of URLs. The former can be used to access the service resource (as described above) and the latter points to a human-readable document describing this resource. \n", - "\n", - "If you click on the Reference URL links, you can find the abstract, Readme file, Vizier table, and much more information associated with this catalog.\n", - "\n", - "Additional information about PyVO Registry Resource and its attributes is available from this page. \n", - "\n", - "\n", - "

***

" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "These examples provide a few ways to access the information of interest. \n", - "\n", - "Below are a few other ways to see what the tap_service table contains. \n", - "\n", - "1. To view the column information: tap_services.to_table().columns() shows the metadata contained in the tap service. We will reference some of this columns below as we try to find the appropriate table. \n", - "\n", - "2. tap_services[index].describe(): The table with the tap_services output has, in our case, 83 tables listed and each includes metadata containing some human readable description. You can get the description for one case or for all the records by iterating through the resource. In the former case, we show the description for the Eichhorn data, whose index is uniq_ind[1]. The latter case also follows. \n", - " " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "print( tap_services.to_table().columns )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "tap_services[uniq_ind[1]].describe() # For Eichhorm+1970 example. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true, - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# To iterate over all the tables: \n", - "for tapsvc in tap_services: \n", - " print(\"--------------------------------------------- \\n\")\n", - " tapsvc.describe()" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - "## Step 2: Acquire the relevant data and make a plot!\n", - "In order to query the table, we need the table name, note this is NOT the same as the short name we found above:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "tables = tap_services[uniq_ind[1]].service.tables\n", - "\n", - "short_name = \"I/90\"\n", - "# find table name: \n", - "for name in tables.keys():\n", - " if short_name in name: \n", - " print(name)" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": ["output_scroll"] - }, - "source": [ - "We can write code to eliminate the other cases (e.g., VI or VIII...) but we wanted to keep this cell to illustrate that the table name (which is required for the query) will likely include the short_name appended to \"/catalog\" (or \"/table\"). \n", - "\n", - "But the other roman numeral catalogs are obviously different catalogs. Therefore try the below for a better match: " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# find (more restricted) table name: \n", - "for name in tables.keys():\n", - " if name.startswith(short_name): \n", - " print(name)\n", - " tablename=name" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "query = 'SELECT * FROM \"%s\"' %tablename\n", - "print(query)\n", - "results = tap_services[short_name].search(query)\n", - "results.to_table()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can access the column data as array using the .getcolumn(colname) attribute, where the colname is given in the table above. In particular the \"CI\" is the color index and \"Ptm\" is the photovisual magnitude. See [here](https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/90) for details about the columns. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "color = results.getcolumn('CI')\n", - "mag = results.getcolumn('Ptm') " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Plotting... \n", - "Note: The magnitudes here are apparent and therefore in plotting, the color-magnitude diagram is typically brightness increasing upwards (higher y-axis) so we will flip the y-axis here. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.ylim(15, 0)\n", - "plt.ylabel(\"V [apparent mag]\")\n", - "plt.xlabel(\"B-V\")\n", - "\n", - "plt.plot(color, mag, 'o', color='black')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Step 3. Compare with other color-magnitude diagrams for Pleiades:\n", - "\n", - "There is nice discussion here: http://www.southastrodel.com/Page03009a.htm about the color-magnitude diagram. Their Fig 4 looks slightly cleaner because part of this investigation was to select the 270 stars that are vetted members and restricted to stellar types more massive than K0. \n", - "\n", - "The dataset is from Raboud+1998 (1998A&A...329..101R)\n", - "\n", - "Therefore in this next step, we will use the bibcode to select this data and overplot with the previous data to compare. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bibcode = '1998A&A...329..101R' # Raboud\n", - "all_bibcodes = tap_services.getcolumn('source_value')\n", - "all_shortnames = tap_services.getcolumn('short_name')\n", - "\n", - "match = np.where(all_bibcodes == bibcode)\n", - "\n", - "# Show relevant short_name (for Raboud paper): \n", - "short_name = all_shortnames[match][0]\n", - "print(short_name)\n", - "print(\"----------\")\n", - "ind = int(match[0])\n", - "\n", - "tap_services[ind].describe()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": ["output_scroll"] - }, - "outputs": [], - "source": [ - "# Doing steps above to view table from Raboud+1998\n", - "\n", - "# This 'tables' will return the same service tables as the 'tables' defined \n", - "# ealier in Step 2, since both the Eichhorn+1970 and Raboud+1998 come from the \n", - "# same service. Therfore, we did not need to redefine 'tables', but we kept it\n", - "# for completion, as different tables don't always use the same service. \n", - "tables = tap_services[ind].service.tables \n", - "\n", - "# find table name: \n", - "for name in tables.keys():\n", - " if short_name in name: \n", - " tablename = name\n", - " \n", - "# Another way to find the table name:\n", - "# Raboud_table = tap_services[short_name].get_tables()\n", - "# tablename = [name for name in Raboud_table.keys()][0]\n", - " \n", - "query = 'SELECT * FROM \"%s\"' %tablename\n", - "print(query)\n", - "results = tap_services[ind].search(query)\n", - "results.to_table()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "R98_color = results.getcolumn('B-V')\n", - "R98_mag = results.getcolumn('Vmag') \n", - "\n", - "plt.ylim(15, 0)\n", - "plt.ylabel(\"V [apparent mag]\")\n", - "plt.xlabel(\"B-V\")\n", - "plt.plot(color, mag, 'o', markersize=4.0, color='black') ## This is Eichhorn data\n", - "plt.plot(R98_color, R98_mag, 's', markersize=5.0, color='red') ## This is new data from Raboud+98\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## BONUS: Step 4: The CMD as a distance indicator! \n", - "\n", - "Since the y-axis above is apparent magnitude, we can use the obvious features (e.g., main sequence curve) to translate the apparent magnitudes to absolute magnitudes (by comparing to published H-R diagrams given in absolute magnitudes) and measure the distance to Pleiades! \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "R98_color = results.getcolumn('B-V')\n", - "R98_mag = results.getcolumn('Vmag') \n", - "\n", - "sun_color = 0.65 # from http://www.astro.ucla.edu/~wright/magcolor.htm\n", - "sun_mag = 10.4 # Played with this value until it looked centered in relation at the B-V color above (yellow star!) \n", - "\n", - "plt.ylim(15, 0)\n", - "plt.ylabel(\"V [apparent mag]\")\n", - "plt.xlabel(\"B-V\")\n", - "plt.plot(color, mag, 'o', markersize=4.0, color='black') ## This is Eichhorn data\n", - "plt.plot(R98_color, R98_mag, 's', markersize=5.0, color='red') ## This is new data from Raboud+98\n", - "plt.plot(sun_color, sun_mag, '*', markersize=15.0, color='yellow') ## This is our estimated center point\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Another measure... use the Sun: \n", - "Vabs = 4.8 ## Sun @ B-V = 0.65 (taken from Wikipedia)\n", - "Vapp = 10.4 ## Based on rough reading of plot above at B-V = 0.65\n", - "\n", - "dm= Vapp - Vabs # distance module = 5log d / 10pc. \n", - "dist = 10. ** (dm / 5. + 1.)\n", - "print(\"%10.1f pc \" %dist)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "True distance to Pleaides is 136.2 pc ( https://en.wikipedia.org/wiki/Pleiades ). Not bad! " - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.8.2" - }, - "nav_menu": {}, - "toc": { - "navigate_menu": true, - "number_sections": true, - "sideBar": true, - "threshold": 6, - "toc_cell": false, - "toc_section_display": "block", - "toc_window_display": true - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/UseCase_III.md b/UseCase_III.md new file mode 100644 index 0000000..c5f4c1d --- /dev/null +++ b/UseCase_III.md @@ -0,0 +1,430 @@ +--- +anaconda-cloud: {} +jupytext: + notebook_metadata_filter: all + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.14.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +language_info: + codemirror_mode: + name: ipython + version: 3 + file_extension: .py + mimetype: text/x-python + name: python + nbconvert_exporter: python + pygments_lexer: ipython3 + version: 3.8.2 +nav_menu: {} +toc: + navigate_menu: true + number_sections: true + sideBar: true + threshold: 6 + toc_cell: false + toc_section_display: block + toc_window_display: true +--- + +# Creating a stellar color-magnitude (or Hertzsprung-Russell) diagram + +The [Hertzsprung-Russell diagram](https://en.wikipedia.org/wiki/Hertzsprung–Russell_diagram) is a fundamental diagram in astronomy that displays important relationships between the stellar color (or temperature) and absolute brightness (or luminosity). + +In this exercise, we will use existing stellar catalogs to produce the H-R diagram. + +```{code-cell} ipython3 +# As a hint, we include the code block for Python modules that you will likely need to import: +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +%matplotlib inline + +# For downloading files +from astropy.utils.data import download_file +from astropy.io import fits + +import pyvo as vo +from pyvo import registry + +## There are a number of relatively unimportant warnings that +## show up, so for now, suppress them: +import warnings +warnings.filterwarnings("ignore", module="astropy.io.votable.*") +warnings.filterwarnings("ignore", module="pyvo.utils.xml.*") +``` + +## Step 1: Find appropriate catalogs + +We want to find a star catalog that has the available data to produce the H-R diagram, i.e., the absolute magnitudes (or both apparent magnitudes AND distances, so we can calculate the absolute magnitudes) in two optical bands (e.g., B and V). This would give us color. Or we need B- OR V- band magnitude and the stellar temperature. + +To simplify this problem, we want to find a catalog of an open cluster of stars, where all the stars were born around the same time and are located in one cluster. This simplifies the issue of getting accurate distances to the stars. One famous cluster is the Pleiades in the constellation of Taurus. So first we start by searching for an existing catalog with data on Pleiades that will provide the necessary information about the stars: magnitudes in two bands (e.g., B and V), which can be used to measure color, or temperature of the star plus one magnitude. + ++++ + +### DATA DISCOVERY STEPS: + +Here is useful link for [how the pyvo registry search works](https://pyvo.readthedocs.io/en/latest/registry/index.html). + +```{code-cell} ipython3 +:tags: [output_scroll] + +tap_services = registry.search(servicetype = 'tap', keywords=['star pleiades'], includeaux=True) +print(len(tap_services)) +tap_services.get_summary() +``` + +Note: The includeaux=True includes auxiliary services. + +Because we specified the service type, this returns only the TAP services for each of the avalible resorces matching our search criteria. So, the 'interfaces' column will only show TAP service as an avalible service, even if other services are avalible from the same resource. We know we need the service to be a TAP service since we know that we want to eventually access the search using additional column information (i.e., beyond RA and Dec, which is all that Simple Cone Search returns). + ++++ + +#### Next, we need to find which of these has the columns of interest, i.e. magnitudes in two bands to create the color-magnitude diagram. + +We can re-run the registry search, but further restrict the results by column UCD. We want tables that have magnitude columns; the most basic UCD to describe a magnitude column is phot.mag + +```{code-cell} ipython3 +:tags: [output_scroll] + +tap_services = registry.search(servicetype='tap', keywords=['star pleiades'], + ucd = ['phot.mag%'], includeaux=True) +print(len(tap_services)) +``` + +Note: the '%' serves as a wild card when searching by UCD + +The IOVA standard enables resources to be as sepecific as they would like when defining the UCD of columns. For example, 'phot.mag' and 'phot.mag;em.opt.V' can both be used to describe a column containing the V magnitudes of objects. If a resource uses the latter to describe a column, a search using 'phot.mag' will not return that resource. A wild card would need to be used or the exact column UCD. The UCD search requires an exact match for a resource to be returned, so using the wild card will make it easier to discover a wider variety of resources. + ++++ {"tags": ["output_scroll"]} + +So using this we can reduce the matched tables to ones that are a bit more catered to our experiment. Note, that there can be redundancy in some resources since these are available via multiple services and/or publishers. Therefore a bit more cleaning can be done to provide only the unique matches. + +```{code-cell} ipython3 +:tags: [output_scroll] + +tap_services.to_table()['ivoid'] +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +def getunique( result ): + short_name = [] + unique_ind = [] + for i in range(len(result)): + short = result[i].short_name + if short not in short_name: + short_name.append(short) + unique_ind.append(i) + else: + print(i) + + return(unique_ind) +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +uniq_ind=getunique(tap_services) +print(len(uniq_ind)) +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +tap_services.to_table()[uniq_ind]['ivoid'] +``` + ++++ {"tags": ["output_scroll"]} + +This shows that in this case, all of our TAP results are unique. + ++++ + +We can read more information about the results we found. For each resource element (i.e. row in the table above), there are useful attributes, which are [described here]( https://pyvo.readthedocs.io/en/latest/api/pyvo.registry.regtap.RegistryResource.html#pyvo.registry.regtap.RegistryResource) + +```{code-cell} ipython3 +:tags: [output_scroll] + +# To read the descriptions of the resulting matches: + +for i in uniq_ind: + print(" *** \n") + print(tap_services[i].creators) + print(tap_services[i].res_description) + + +``` + ++++ {"tags": ["output_scroll"]} + + RESULT: Based on these, the second one (by Eichhorn et al) looks like a good start. + +### At this point, you can proceed to Step 2. + +-- OR -- + +### Try a different data discovery method! + ++++ {"tags": ["output_scroll"]} + +### Alternative Method: Use ADS to search for appropriate paper and access data via NED. + +There are multiple paths for the data discovery. So it may also be that you know the paper that has the data you are interested in and want to access via the bibcode or authors, etc. + +In this case, let's assume that we have the information that the Eichhorn+1970 paper has the data that we need to create the H-R diagram: https://ui.adsabs.harvard.edu/abs/1970MmRAS..73..125E/abstract + +We can either search by bibcode (1970MmRAS..73..125E) or "Eichhorn" to get the access_urls that will allow us to work with the data. + +Before this step, if may help to see the names of the fields available to use. Notice the following fields: + +"source_value" contains the bibcode information that we want; "creator_seq" lists the authors; + +and + +"access_url" provides the url from where the data can be accessed. + +```{code-cell} ipython3 +## You already have this from above: +# tap_services = registry.search(servicetype='tap', keywords=['star pleiades'], +# ucd = ['phot.mag%'], includeaux=True) +print(tap_services.fieldnames) +``` + +First, Try using bibcode: + +```{code-cell} ipython3 +bibcode = '1970MmRAS..73..125E' # Eichhorn +idx=-1 +for s in tap_services: + idx+=1 + if bibcode in s['source_value']: + myidx=idx + print(f"{s.short_name}, {s.source_value}, {s.access_url}") + break +# Note that using the to_table() lets you search the result +# easily using all columns. But in the end, you want to get +# back not an astropy table row, which you cannot use, but the +# original RegistryResult that has the callable TAP service. +myTAP=tap_services[myidx] +``` + +Note that the URL is a generic TAP url for Vizier. All of its tables can be accessed by that same TAP services. It'll be in the ADQL query itself that you specify the table name. We'll see this below. + ++++ + +Next, try using Author name: + +```{code-cell} ipython3 +:tags: [output_scroll] + +author = 'Eichhorn' + +for record in tap_services: + names=record.creators + if 'Eichhorn' in names[0]: + print("For %s: " %record.short_name) + print(" Access URL: %s" %record['access_urls'][0]) + print(" Reference URL: %s" %record.reference_url) + + +``` + +In the code above, the record is a Registry Resource. You can access the attribute, "creators", from the resource, which is relevant for our example here since this is a direct way to get the author names. The other attributes, "access_url" and "reference_url", provides two types of URLs. The former can be used to access the service resource (as described above) and the latter points to a human-readable document describing this resource. + +If you click on the Reference URL links, you can find the abstract, Readme file, Vizier table, and much more information associated with this catalog. + +Additional information about PyVO Registry Resource and its attributes is available from this page. + + +

***

+ ++++ + +These examples provide a few ways to access the information of interest. + +Below are a few other ways to see what the tap_service table contains. + +1. To view the column information: tap_services.to_table().columns() shows the metadata contained in the tap service. We will reference some of this columns below as we try to find the appropriate table. + +2. tap_services[index].describe(): The table with the tap_services output has, in our case, 83 tables listed and each includes metadata containing some human readable description. You can get the description for one case or for all the records by iterating through the resource. In the former case, we show the description for the Eichhorn data, whose index is uniq_ind[1]. The latter case also follows. + + +```{code-cell} ipython3 +:tags: [output_scroll] + +print( tap_services.to_table().columns ) +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +tap_services[uniq_ind[1]].describe() # For Eichhorm+1970 example. +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +# To iterate over all the tables: +for tapsvc in tap_services: + print("--------------------------------------------- \n") + tapsvc.describe() +``` + ++++ {"tags": ["output_scroll"]} + +## Step 2: Acquire the relevant data and make a plot! +In order to query the table, we need the table name, note this is NOT the same as the short name we found above: + +```{code-cell} ipython3 +:tags: [output_scroll] + +tables = tap_services[uniq_ind[1]].service.tables + +short_name = "I/90" +# find table name: +for name in tables.keys(): + if short_name in name: + print(name) +``` + ++++ {"tags": ["output_scroll"]} + +We can write code to eliminate the other cases (e.g., VI or VIII...) but we wanted to keep this cell to illustrate that the table name (which is required for the query) will likely include the short_name appended to "/catalog" (or "/table"). + +But the other roman numeral catalogs are obviously different catalogs. Therefore try the below for a better match: + +```{code-cell} ipython3 +:tags: [output_scroll] + +# find (more restricted) table name: +for name in tables.keys(): + if name.startswith(short_name): + print(name) + tablename=name +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +query = 'SELECT * FROM "%s"' %tablename +print(query) +results = tap_services[short_name].search(query) +results.to_table() +``` + +We can access the column data as array using the .getcolumn(colname) attribute, where the colname is given in the table above. In particular the "CI" is the color index and "Ptm" is the photovisual magnitude. See [here](https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/90) for details about the columns. + +```{code-cell} ipython3 +color = results.getcolumn('CI') +mag = results.getcolumn('Ptm') +``` + +### Plotting... +Note: The magnitudes here are apparent and therefore in plotting, the color-magnitude diagram is typically brightness increasing upwards (higher y-axis) so we will flip the y-axis here. + +```{code-cell} ipython3 +plt.ylim(15, 0) +plt.ylabel("V [apparent mag]") +plt.xlabel("B-V") + +plt.plot(color, mag, 'o', color='black') +``` + +## Step 3. Compare with other color-magnitude diagrams for Pleiades: + +There is nice discussion here: http://www.southastrodel.com/Page03009a.htm about the color-magnitude diagram. Their Fig 4 looks slightly cleaner because part of this investigation was to select the 270 stars that are vetted members and restricted to stellar types more massive than K0. + +The dataset is from Raboud+1998 (1998A&A...329..101R) + +Therefore in this next step, we will use the bibcode to select this data and overplot with the previous data to compare. + +```{code-cell} ipython3 +bibcode = '1998A&A...329..101R' # Raboud +all_bibcodes = tap_services.getcolumn('source_value') +all_shortnames = tap_services.getcolumn('short_name') + +match = np.where(all_bibcodes == bibcode) + +# Show relevant short_name (for Raboud paper): +short_name = all_shortnames[match][0] +print(short_name) +print("----------") +ind = int(match[0]) + +tap_services[ind].describe() +``` + +```{code-cell} ipython3 +:tags: [output_scroll] + +# Doing steps above to view table from Raboud+1998 + +# This 'tables' will return the same service tables as the 'tables' defined +# ealier in Step 2, since both the Eichhorn+1970 and Raboud+1998 come from the +# same service. Therfore, we did not need to redefine 'tables', but we kept it +# for completion, as different tables don't always use the same service. +tables = tap_services[ind].service.tables + +# find table name: +for name in tables.keys(): + if short_name in name: + tablename = name + +# Another way to find the table name: +# Raboud_table = tap_services[short_name].get_tables() +# tablename = [name for name in Raboud_table.keys()][0] + +query = 'SELECT * FROM "%s"' %tablename +print(query) +results = tap_services[ind].search(query) +results.to_table() +``` + +```{code-cell} ipython3 +R98_color = results.getcolumn('B-V') +R98_mag = results.getcolumn('Vmag') + +plt.ylim(15, 0) +plt.ylabel("V [apparent mag]") +plt.xlabel("B-V") +plt.plot(color, mag, 'o', markersize=4.0, color='black') ## This is Eichhorn data +plt.plot(R98_color, R98_mag, 's', markersize=5.0, color='red') ## This is new data from Raboud+98 +``` + +## BONUS: Step 4: The CMD as a distance indicator! + +Since the y-axis above is apparent magnitude, we can use the obvious features (e.g., main sequence curve) to translate the apparent magnitudes to absolute magnitudes (by comparing to published H-R diagrams given in absolute magnitudes) and measure the distance to Pleiades! + +```{code-cell} ipython3 +R98_color = results.getcolumn('B-V') +R98_mag = results.getcolumn('Vmag') + +sun_color = 0.65 # from http://www.astro.ucla.edu/~wright/magcolor.htm +sun_mag = 10.4 # Played with this value until it looked centered in relation at the B-V color above (yellow star!) + +plt.ylim(15, 0) +plt.ylabel("V [apparent mag]") +plt.xlabel("B-V") +plt.plot(color, mag, 'o', markersize=4.0, color='black') ## This is Eichhorn data +plt.plot(R98_color, R98_mag, 's', markersize=5.0, color='red') ## This is new data from Raboud+98 +plt.plot(sun_color, sun_mag, '*', markersize=15.0, color='yellow') ## This is our estimated center point +``` + +```{code-cell} ipython3 +# Another measure... use the Sun: +Vabs = 4.8 ## Sun @ B-V = 0.65 (taken from Wikipedia) +Vapp = 10.4 ## Based on rough reading of plot above at B-V = 0.65 + +dm= Vapp - Vabs # distance module = 5log d / 10pc. +dist = 10. ** (dm / 5. + 1.) +print("%10.1f pc " %dist) +``` + +True distance to Pleaides is 136.2 pc ( https://en.wikipedia.org/wiki/Pleiades ). Not bad! diff --git a/check_env.py b/check_env.py index 2fba41e..06b37b7 100644 --- a/check_env.py +++ b/check_env.py @@ -18,7 +18,8 @@ 'astropy': ('5.2', None), 'scipy': ('1.0', None), 'pyvo': ('1.4', None), - 'astroquery': ('0.4.7dev', None) + 'astroquery': ('0.4.7dev', None), + 'jupytext': ('1.14', None) } diff --git a/conf.py b/conf.py index 315500d..609bca3 100644 --- a/conf.py +++ b/conf.py @@ -70,3 +70,6 @@ # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] + +# myst configurations +myst_heading_anchors = 4 diff --git a/doc-requirements.txt b/doc-requirements.txt index 878da0e..ec3dd6c 100644 --- a/doc-requirements.txt +++ b/doc-requirements.txt @@ -5,3 +5,4 @@ sphinx-book-theme sphinx-copybutton nbval pytest +jupytext diff --git a/environment.yml b/environment.yml index a524499..b8b3d50 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,8 @@ dependencies: - jupyterlab>=3.1 - pyvo>=1.4 - scipy>=1.0 + - jupytext>=1.14 - pip - pip: - astroquery - - --pre \ No newline at end of file + - --pre diff --git a/tox.ini b/tox.ini index aa77669..a342b6d 100644 --- a/tox.ini +++ b/tox.ini @@ -31,6 +31,7 @@ deps = commands = pip freeze + !buildhtml: jupytext --from myst --to notebook *.md !buildhtml: pytest --nbval buildhtml: sphinx-build -b html . _build/html -D nb_execution_mode=auto -nWT --keep-going