{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "considered-spanking",
   "metadata": {},
   "source": [
    "# Data Discovery in using pyVO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "registered-mirror",
   "metadata": {},
   "source": [
    "This is a demo for how to use the proposed (as of Jan 2022; this may be merged when you read this) new features of pyvo.registry.  See https://blog.g-vo.org/towards-data-discovery-in-pyvo.html for some background.\n",
    "\n",
    "Until add-discoverdata is merged into pyvo, you can find HTML documentation for this at http://docs.g-vo.org/temporary-pyvo/registry/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "published-fountain",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up things; we're also ignoring over-zealous\n",
    "# astropy warnings against bleeding-edge VOTable.\n",
    "from pyvo import registry, dal\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore', module=\"astropy.io.votable.*\")\n",
    "warnings.filterwarnings('ignore', module=\"pyvo.io.vosi.vodataservice\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modified-mitchell",
   "metadata": {},
   "source": [
    "The most general way to run registry queries is by passing registry.search Constraints.  It is quite a bit more flexible than the conventional keywords, but admittedly somewhat more verbose.\n",
    "\n",
    "For instance, to find data giving redshifts on quasars, you could say:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "beginning-explanation",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs = registry.search(\n",
    "    registry.Freetext(\"quasar\"),\n",
    "    registry.UCD(\"src.redshift\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "smooth-electric",
   "metadata": {},
   "source": [
    "In many cases, you can also use an equivalent based on keyword arguments, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mineral-national",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs = registry.search(keywords=\"quasar\", \n",
    "                       ucd=\"src.redshift\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "substantial-emission",
   "metadata": {},
   "source": [
    "The list of constraints available (and explanations what they do) is available at http://docs.g-vo.org/temporary-pyvo/registry/#basic-interface.\n",
    "\n",
    "What's coming back here is a collection of “Resources” (in general corresponding to data collections).  The simplest way to have a look at them is through the ``to_table`` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rotary-brain",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs.get_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "skilled-carter",
   "metadata": {},
   "source": [
    "While this particular list is perhaps a bit unwieldy, this lets you relatively quickly browse what is available.  In particular, the last column tells you how you can talk to a service serving the data.\n",
    "\n",
    "Once you have found data you are interested in, you can pick it out of the list using the numeric index (which, however, is unstable between sessions and thus we don't do it here), using the short name (for which there *could* be clashes, though they should be rare) or through the ivoid (which is globally unique, but somewhat lengthy).  In this example, we are using the short name.\n",
    "\n",
    "Let's say we want to work with the resource III/175, “Optical spectroscopy of radio sources”.  By the last column, there is a cone search, TAP, and web service that provides access to it.\n",
    "\n",
    "The most immediate way to get to the data usually is the cone search, which gives something like a dump of a catalogue around a position (using 0,0,180 will give you the full catalogue most of the time).  To get to the service from your registry result, use the ``get_service`` method of the registry resource, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "touched-ratio",
   "metadata": {},
   "outputs": [],
   "source": [
    "svc = rscs[\"III/175\"].get_service(\"conesearch\")\n",
    "svc.search((126, -20), 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "responsible-bradley",
   "metadata": {},
   "source": [
    "A more powerful interface is TAP, which lets you send database queries to the service (forget about the “#aux” in the interface name for now).  To do something sensible in TAP, you need to know the name(s) of the table(s) making up the resource.  You can figure these out using the registry record's ``get_tables`` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brave-biotechnology",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs[\"III/175\"].get_tables()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "comprehensive-consolidation",
   "metadata": {},
   "source": [
    "Let's have a look at what columns one of these tables has – this is a normal astropy table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "young-hundred",
   "metadata": {},
   "outputs": [],
   "source": [
    "td = rscs[\"III/175\"].get_tables()['III/175/table1']\n",
    "td.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olympic-second",
   "metadata": {},
   "source": [
    "From here, you could inspect the various BaseParams for units, descriptions, and the like, but for this level of interactivity, you may want to use TOPCAT.  Give it the service's access URL:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sudden-jerusalem",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs[\"III/175\"].get_interface(\"tap\").access_url"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mounted-indianapolis",
   "metadata": {},
   "source": [
    "Oh, and while I was preparing this notebook, the metadata of this resource still had a bug, which showed itself as warnings of the type\n",
    "\n",
    "```\n",
    "WARNING: W02: ?:?:?: W02: '' is not a valid datatype according to the VOSI spec [pyvo.io.vosi.vodataservice]\n",
    "```\n",
    "\n",
    "While you might ignore warnings, at least with errors it is usually a good idea to notify the operators.  To see who to talk to, use the ``get_contact`` method of the record:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "floral-mountain",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs[\"III/175\"].get_contact()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "disabled-compilation",
   "metadata": {},
   "source": [
    "To actually run queries, get a TAP service and do queries based on the columns that you found.  For instance, here is how to find what types there are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "transsexual-firmware",
   "metadata": {},
   "outputs": [],
   "source": [
    "svc = rscs[\"III/175\"].get_service(\"tap\")\n",
    "svc.run_sync('SELECT DISTINCT type FROM \"III/175/table1\"').to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "structural-residence",
   "metadata": {},
   "source": [
    "To figure out the correlation between the 5 GHz flux and the optical magnitude for Quasars, you could do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ranking-today",
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_and_mag = svc.run_sync('SELECT m, S5Ghz FROM \"III/175/table1\"').to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "inner-award",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "stats.pearsonr(flux_and_mag[\"m\"], flux_and_mag[\"S5GHz\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "demographic-employee",
   "metadata": {},
   "source": [
    "That there's an anticorrelation (the first value returned) is not surprising (magnitudes grow as flux decrease).  The low statistics for correlatedness (the second value) might be more surprising.  My next step would be to do a plot.  But I've vowed to not have a plot in here, so let's go on to the last interface type on our sample service: web.\n",
    "\n",
    "That corresponds to something you can operate with your web browser, and hence there's just one thing pyVO can do: Open a browser.  That happens when you call the service's ``search`` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "speaking-latest",
   "metadata": {},
   "outputs": [],
   "source": [
    "rscs[\"III/175\"].get_service(\"web\").search()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "trying-bubble",
   "metadata": {},
   "source": [
    "By the way, this is *not* the way to look for a webpage *on* the service.  That web pages is available (provided the publishers did their homework) in a resources' reference_url attribute.  To get there, you could do: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nonprofit-record",
   "metadata": {},
   "outputs": [],
   "source": [
    "import webbrowser\n",
    "webbrowser.open(rscs[\"III/175\"].reference_url, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "challenging-discount",
   "metadata": {},
   "source": [
    "There are more constraints available than just free text and UCD.\n",
    "A particularly interesting one is the spatial coverage. For instance, you could look for data on flare stars around the Orion nebula like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "circular-express",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import SkyCoord\n",
    "flrscs = registry.search(\n",
    "    registry.Freetext(\"flare\"),\n",
    "    registry.Spatial((SkyCoord.from_name(\"M42\"), 2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "binding-brook",
   "metadata": {},
   "outputs": [],
   "source": [
    "flrscs.get_summary().show_in_notebook(display_length=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hydraulic-rating",
   "metadata": {},
   "source": [
    "The services here a bit more diverse than with our first example.  For instance, there are image services, as you will see when you skim the last column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sunrise-tournament",
   "metadata": {},
   "outputs": [],
   "source": [
    "matches = flrscs[\"flare_survey.dat\"].get_service(\"sia\").search(\n",
    "    pos=SkyCoord.from_name(\"M42\"),\n",
    "    size=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "applicable-inspection",
   "metadata": {},
   "source": [
    "In order to have at least a few images, let's use datalink to fetch a few previews of our matches (this doesn't work on all services; if it doesn't complain to the operators, demanding datalink support – see the thing with get_contact above)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "realistic-evans",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image, display\n",
    "for dl in matches.iter_datalinks():\n",
    "    for row in dl.bysemantics(\"#preview\"):\n",
    "        display(Image(url=row[\"access_url\"], width=200,\n",
    "                            embed=True, format=\"jpeg\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "committed-wheel",
   "metadata": {},
   "source": [
    "There are similar constraints for the Spectral and Time axes.  For instance, to look for resources talking about spectra and the Balmer break, you could say:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "substantial-nightmare",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "registry.search(\n",
    "    registry.Freetext(\"spectra\"),\n",
    "    registry.Spectral(364*u.nm)).get_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "renewable-single",
   "metadata": {},
   "source": [
    "Note that in particular for time and spectral coverage, as of 2022 many data providers in the VO have not updated their resource records to provide such information; hence, you will have to expect missing resources.  For spectral coverage, see also the ``Waveband`` constraint, which is older and therefore better supported."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continuous-telephone",
   "metadata": {},
   "source": [
    "Behind the scenes, all this just does ADQL queries via TAP.  So, whenever the pre-canned queries from the Registry module are not enough (e.g., because you want to do table uploads or need exotic constraints), you can simply switch to using TAP directly.  To help you with that, you can use the ``build_regtap_query`` function to get an ADQL query to start with.  For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "worth-catch",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(registry.get_RegTAP_query(\n",
    "    registry.Spatial((30, 40)),\n",
    "    registry.Servicetype('tap'),\n",
    "    registry.Datamodel(\"obscore\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "southwest-highway",
   "metadata": {},
   "source": [
    "This isn't overly pretty, but once you've had a look at the RegTAP documentation at https://ivoa.net/documents/RegTAP/, it should start to make sense.  By cutting and pasting, you could create a registry query using an uploaded object list, perhaps a bit like this (ignore the next code cells if you've not played with TAP uploads yet and/or feel uncomfortable near to large amounts of ADQL):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dedicated-snowboard",
   "metadata": {},
   "outputs": [],
   "source": [
    "objects = dal.TAPService(\"http://dc.g-vo.org/tap\").run_sync(\n",
    "    \"SELECT source_id, ra, dec FROM gaia.dr3lite TABLESAMPLE(0.000005)\"\n",
    ").to_table()\n",
    "objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "musical-council",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyvo.registry import regtap\n",
    "\n",
    "rt_query = \"\"\"\n",
    "SELECT DISTINCT\n",
    "ivoid, res_title, \n",
    "res_description, access_url FROM\n",
    "rr.resource\n",
    "NATURAL LEFT OUTER JOIN rr.capability\n",
    "NATURAL LEFT OUTER JOIN rr.interface\n",
    "NATURAL LEFT OUTER JOIN rr.res_detail\n",
    "NATURAL LEFT OUTER JOIN rr.stc_spatial\n",
    "JOIN TAP_UPLOAD.t1\n",
    "ON\n",
    "  (1 = CONTAINS(MOC(6, POINT(TAP_UPLOAD.t1.ra, TAP_UPLOAD.t1.dec)), coverage))\n",
    "WHERE\n",
    "  (detail_xpath = '/capability/dataModel/@ivo-id' AND 1 = ivo_nocasematch(detail_value, 'ivo://ivoa.net/std/obscore%'))\n",
    "  AND (standard_id IN ('ivo://ivoa.net/std/tap'))\n",
    "\"\"\"\n",
    "ocrscs = regtap.get_RegTAP_service(\n",
    "    ).run_sync(rt_query, uploads={\"t1\": objects}).to_table()\n",
    "ocrscs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interim-entry",
   "metadata": {},
   "source": [
    "Note, however, that in particular Obscore services are notoriously bad at properly defining their physical coverage, so this sort of query is probably more appropriate for TAP tables and perhaps image or spectral services."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "complete-lebanon",
   "metadata": {},
   "source": [
    "Finally, as soon as you constrain the service type, ``registry.search`` will work exactly as before add-discoverdata.  So, you can still do an all-VO query for, say, X-ray images as you always could (thanks to ``get_service``, it's now even a bit simpler):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cordless-george",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_matches = 0\n",
    "for res in registry.search(\n",
    "        keywords=\"rosat\", waveband=\"X-Ray\", servicetype=\"image\"):\n",
    "    try:\n",
    "        print(f\"Querying {res.short_name}...\")\n",
    "        mats = res.get_service().search(pos=(30, 20), size=0.3)\n",
    "        print(f\"...yielded {len(mats)}\")\n",
    "        total_matches += len(mats)\n",
    "    except Exception as msg:\n",
    "        print(f\"Service {res.short_name} failed: {msg}\")\n",
    "print(f\"Total found: {total_matches}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "corresponding-pharmacy",
   "metadata": {},
   "source": [
    "Comments, questions and ideas for improvement are very welcome. Contact:\n",
    "msdemlei@ari.uni-heidelberg.de (PGP key: 0x555FA86CC57AE128)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "welsh-fifth",
   "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}