Posts with the Tag pyVO:

• Towards Data Discovery in pyVO

  2022-01-10 Markus Demleitner

  When I struggled with ways to properly integrate TAP services – which may have hundreds or thousands of different resources in one service – into the VO Registry without breaking what we already had, I realised that there are really two fundamentally different modes of using the VO Registry. In Discovering Data Collections's abstract I wrote:

  the Registry must support both VO-wide discovery of services by type ("service enumeration") and discovery by data collection ("data discovery").

  To illustrate the difference in a non-TAP case, suppose I have archived images of lensed quasars from Telescopes A, B, and C. All these image collections are resources in their own right and should be separately findable when people look for “resources with data from Telescope A“ or perhaps “images obtained between 2011-01-01 and 2011-12-31”.

  However, when a machine wants to find all images at a certain position, publishing the three resources through three different services would mean that that machine has to do three requests where one would work just as well. That is very relevant when you think about how the VO will evolve: At this point there are 342 SIAP services in the VO, and when you read this, that number may have grown further. Adding one service per collection will simply not scale when we want to keep the possibility of all-VO searches. Since I claim that is a very desirably thing, we need to enable collective services covering multiple subordinate resources.

  So, while in the first (“data discovery”) case one wants to query (or at least discover) the three resources separately, in the second case they should be ignored, and only a collective “images of lensed quasars” service should be queried.

  The technical solution to this requirement was creating “auxliary capabilities” as discribed in the endorsed note on discoving data collections cited above. But these of course need client support; VO clients up to now by and large do service enumeration, as that has been what we started with in the VO Registry. Client support would, roughly, mean that clients would present their users with data collections, and then offer the various ways to to access them.

  There are quite a number of technicalities involved in why that's not terribly straightforward for the “big” clients like TOPCAT and Aladin (though Aladin's discovery tree already comes rather close).

  Now that quite a number of people use pyVO interactively in jupyter notebooks, extending pyVO's registry interface to do data discovery in addition to the conventional service enumeration becomes an attractive target to have data discovery in practice.

  I have hence created pyVO PR #289. I think some the rough edges will need to be smoothed out before it can be merged, but meanwhile I'd be grateful if you could try it out already. To facilitate that, I have prepared a jupyter notebook that shows the basic ideas.

  Followup (2023-12-15)

  I have just prepared a slightly updated version of the notebook.

  To run it while the PR is not merged, you need to install the forked pyVO. In order to not clobber your main installation, you can install astropy using your package manager and then do the following (assuming your shell is bash or something suitably similar):

  virtualenv --system-site-packages try-discoverdata
  . try-discoverdata/bin/activate
  cd try-discoverdata
  git clone https://github.com/msdemlei/pyvo
  cd pyvo
  git checkout add-discoverdata
  python3 setup.py develop
  ipython3 notebook
  

  That should open a browser window in which you can open the notebook (you probably want to download it into the pyvo checkout in order to make the notebook selector see it). Enjoy!

• LAMOST5 meets Datalink

  2019-12-11 Markus Demleitner

  One of the busiest spectral survey instruments operated right now is the Large Sky Area Multi-Object Fiber Spectrograph Telescope (LAMOST). And its data in the VO, more or less: DR2 and DR3 have been brought into the VO by our Czech colleagues, but since they currently lack resources to update their services to the latest releases, they have kindly given me their DaCHS resource descriptor, and so I had a head start for publishing DR5 in Heidelberg.

  With some minor updates, here it is now: Over nine million medium-resolution spectra covering large parts of the northen sky – the spatial coverage is like this:

  

  There's lots of fun to be had with this; of course, there's an SSA service, so when you point Aladin or Splat at some part of the covered sky and look for spectra, chances are you'll see LAMOST spectra, and when working on some of our tutorials (this one, for example), it happened that LAMOST actually had what I was looking for when writing them.

  But I'd like to use the opportunity to mention two other modes of accessing the data.

  

  Tablesample and TOPCAT's Plot Table activation action

  Say you'd like to look at spectra of M stars and would like to have some sample from across the sky, fire up TOPCAT, point its TAP client the GAVO DC TAP service (http://dc.g-vo.org/tap) and run something like:

  select
    ssa_pubDID, accref, raj2000, dej2000, ssa_targsubclass
  from lamost5.data tablesample(1)
  where
    ssa_targsubclass like 'M%'
  

  This is using the TABLESAMPLE modifier in the from clause, which isn't standard ADQL yet. As mentioned in the DaCHS 1.4 announcement, DaCHS has a prototype implementation of what's been discussed on the IVOA's DAL mailing list: pick a part of a table rather than the full one. It takes a percentage as an argument, and tells the server to choose about this percentage of the table's records using a reasonable and fast heuristic. Note that this won't give you perfect statistical sampling, but if it's not “good enough” for some purpose, I'd like to learn about that purpose.

  Drawing a proper statistical sample, on the other hand, would take minutes on the GAVO database server – with tablesample, I had the roughly 6000 spectra the above query returns essentially instantaneously, and from eyeballing a sky plot of them, I'd say their distribution is close enough to that of the full DR5. So: tablesample is your friend.

  For a quick look at the spectra themselves, in TOPCAT click Views/Activation Actions, check “Plot Table” and make sure TOPCAT proposes the accref column as “Table Location” (if you don't see these items, update your TOPCAT – it's worth it). Now click on a row or perhaps a dot on a plot and behold an M spectrum.

  Cutouts via Datalink

  LAMOST releases spectra in FITS format pretty much like the ones you may know from SDSS. The trick above works because we instead hand out proper, IVOA Spectral Data Model-compliant spectra through SSA and TAP. However, if you need to go back to the original files, you can, using Datalink. If you're unsure what this Datalink thing is: call me vain, but I still like my 2015 ADASS poster explaining that. In TOPCAT, you'd be using the “Invoke Service” activation action to get to the datalinks.

  If you have actual work to do, offloading repetetive work to the computer is what you want, and fortunately, pyVO knows about datalink, too. I give you this is hard to discover so far, and the interface is... a tiny bit clunky. Until some kind soul cleans up the pyVO datalink act, a poster Stefan and I showed at the 2017 ADASS might give you an idea which buttons to press. Or read on and see how things work for LAMOST5.

  The shortest way to datalinks is a TAP query that at least retrieves the ssa_pubdid column (that's a must; Datalink can't work without it) and, on the result, run the iter_datalinks method. This returns an object in which you can find the associated data items (in this case, a preview and the original FITS with the #progenitor semantics), plus the cutout service.

  Hence, a minimal example for pulling the legacy FITS links out of the first three items in lamost5.data would look like this:

  import pyvo
  
  svc = pyvo.dal.TAPService("http://dc.g-vo.org/tap")
  for dl in svc.run_sync("select top 3 ssa_pubdid"
          " from lamost5.data").iter_datalinks():
      print(next(dl.bysemantics("#progenitor")
          )["access_url"].decode("ascii"))
  

  This is a bit different from listing 2 in the poster linked above because it's python3, so getting the first element from iterator an iterator looks a bit different, and (curse astropy.votable for returning VOTable chars as bytes rather than strings!) you'll want to turn the URL into a proper string manually.

  Another, actually more interesting, thing you can do with Datalink is cut out regions of interest. The LAMOST spectra are fairly long (though of course still small by image standards), so if you're only interested in a single line, you can save a bit of storage and bandwidth over blindly pulling the whole thing.

  For instance, if you wanted to pull the vicinity of the H and K Fraunhofer lines from the matches in the loop in the snippet above, you could say:

  from astropy import units as u
  proc = next(dl.iter_procs())
  cutout = proc.processed(band=(392*u.nm,398*u.nm))
  

  And this is what I've done for the decorative left border above: it's the H and K line profiles for 0.1% of the stars LAMOST has classified as G8. Building the image didn't take more than a few seconds (where I'd like the cutouts to be faster by a factor of 10; I guess that's about an afternoon of work for me, so if it'd save you more than that afternoon, poke me to do it).

  What's coming back is tables. By the time python has digested these, they're numpy record arrays. Thus, you can immediately bring in your beloved scipy (or whatever). For instance, if for some reason you're convinced that the H and K lines should be fit by identical Gaussians in the boring case and would like find objects for which that's patently untrue and that hence could be un-boring, here's how you could do that:

  def spectral_model(wl, c1, c2, depth, width):
      return (1
          -depth*numpy.exp(-numpy.square(wl-c1)
              /numpy.square(width))
          -depth*numpy.exp(-numpy.square(wl-c2)
              /numpy.square(width)))
  
  for pubdid, prof in get_profiles(
          "G8", (392*u.nm,398*u.nm), 0.01, 4):
      prof["flux"] /= max(prof["flux"])
      popt, pcov = curve_fit(
          spectral_model, prof["spectral"], prof["flux"],
          sigma=prof["flux_error"],
          p0=[3968, 3934, 1, 1])
      if pcov[3][3]>1:
          break
  

  – where get_profiles is essentially doing the TAP plus datalink routine above, except I'm swallowing spectra with too much noise and I have the function transform the spectral coordinate into the objects' rest frames. If you're curious how I'm doing this just based on the IVOA Spectral Data Model, check the source linked at the bottom of this post.

  I've just run this, and the first spectrum that the machinery flagged as suspicious was this:

  

  – which doesn't look like I've made a discovery just yet. But that doesn't mean there's not a lot to find within LAMOST5's lines...

  To get you up to speed quickly: here's the actual python3 code I ran for the “analysis” and the plot.

• Small Telescopes, Large Surveys

  2019-03-13 Markus Demleitner

  

  Plate technology at Bamberg observatory: a blink comparator with one plate mounted, and a survey camera that was once used at Boyden Station, an astronomer outpost in 60ies South Africa.

  I'm currently at the workshop “Large surveys with small telescopes: past, present, and future” (or Astroplate III for short) in Bamberg, where people are discussing using and re-using the rich heritage of historical observations (hence the “plate” part) as well growing that heritage in the age of large CCDs, fast computers and large disks.

  Using and re-using is of course what the Virtual Observatory is about, and we've been keeping fairly large plate collections in our data center for quite a while (among them the Archives of Landessternwarte Königstuhl or the Palomar-Leiden Trojan surveys, and there is the WFPDB TAP-accessibly). Therefore, people from GAVO Heidelberg have been to all past astroplate conferences.

  For this one, I brought a brand-new tutorial on plate scans in the VO, which, I hope, also works as a general introduction to image discovery in the VO using SIAP, Datalink, and Obscore. If you're doing image stuff now and then, please have a quick look at the thing – I am particularly grateful for hints on what to improve or perhaps particularly obvious use cases for the material discussed.

  Such VO proselytising aside, the conference is discussing the wide variety of creative, low-cost data collectors out there as well as computer-aided re-analysis extracting new knowledge from decades-old data. If I had to choose a single come-to-think-of-it moment, it would be Norbert Zacharias' observation that if you have a well-behaved object and you'd like to know where it was in 1900, it's now more accurate to extrapolate Gaia astrometry to the epoch of observation than to measure it on the plate itself. Which is saying a lot about the amazing feat of engineering that Gaia is.

  This is not, however, an argument for dumping the old data. Usually, it is exactly what is not so well-behaved (like those) that's interesting – both in terms of astrometry and in terms of photometry (for which there's a lot more unruly behaviour in the first place). To figure out how objects don't behave well, and, for objects disguising as well-behaved only on time scales of the (say) Gaia mission, which these are, the key is “old” data. The freshness of which we're discussing this week.

• Gaia DR2: A light version and light curves

  2018-05-07 Markus Demleitner

  

  Topcat is doing datalink, and our little python script has plotted a two-color time series of RMC 18 (or so I think).

  If anyone ever writes a history of the VO, the second data release of Gaia on April 25, 2018 will probably mark its coming-of-age – at least if you, like me, consider the Registry the central element of the VO. It was spectacular to view the spike of tens of Registry queries per second right around 12:00 CEST, the moment the various TAP services handing out the data made it public (with great aplomb, of course).

  In GAVO's Data Center we also carry Gaia DR2 data. Our host institute, the Zentrum für Astronomie in Heidelberg, also has a dedicated Gaia server. This gives relieves us from having to be a true mirror of the upstream data release. And since the source catalog has lots and lots of columns that most users will not be using most of the time, we figured a “light” version of the source catalog might fill an interesting ecological niche: Behold gaia.dr2light on the GAVO DC TAP service, containing essentially just the basic astrometric parameters and the diagonal of the covariance matrix.

  That has two advantages: Result sets with SELECT * are a lot less unwieldy (but: just don't do this with Gaia DR2), and, more importantly, a lighter table puts less load on the server. You see, conventional databases read entire rows when processing data, and having just 30% of the columns means we will be 3 times faster on I/O-bound tasks (assuming the same hardware, of course). Hence, and contrary to several other DR2-carrying sites, you can perform full sequential scans before timing out on our TAP service on gaia.dr2light. If, on the other hand, you need to do debugging or full-covariance-matrix error calculations: The full DR2 gaia_source table is available in many places in the VO. Just use the Registry.

  Photometry via TAP

  A piece of Gaia DR2 that's not available in this form anywhere else is the lightcurves; that's per-transit photometry in the G, BP, and RP band for about 0.5 million objects that the reduction system classified as variable. ESAC publishes these through datalink from within their gaia_source table, and what you get back is a VOTable that has the photometry in the three bands interleaved.

  I figured it might be useful if that data were available in a TAP-queriable table with lightcurves in the database. And that's how gaia.dr2epochflux came into being. In there, you have three triples of arrays: the epochs (g_transit_time, bp_obs_time, and rp_obs_time), the fluxes (g_transit_flux, bp_flux, and rp_flux), and their errors (you can probably guess their names). So, to retrieve G lightcurves where available together with a gaia_source query of your liking, you could write something like:

  SELECT g.*, g_transit_time, g_transit_flux
  FROM gaia.dr2light AS g
  LEFT OUTER JOIN gaia.dr2epochflux
  USING (source_id)
  WHERE ...whatever...
  

  – the LEFT OUTER JOIN arranges things such that the g_transit_time and g_transit_flux columns simply are NULL when there are no lightcurves; with a normal (“inner”) join, rows without lightcurves would not be returned in such a query.

  To give you an idea of what you can do with this, suppose you would like to discover new variable blue supergiants in the Gaia data (who knows – you might discover the precursor of the next nearby supernova!). You could start with establishing color cuts and train your favourite machine learning device on light curves of variable blue supergiants. Here's how to get (and, for simplicity, plot) time series of stars classified as blue supergiants by Simbad for which Gaia DR2 lightcurves are available, using pyvo and a little async trick:

  from matplotlib import pyplot as plt
  import pyvo
  
  def main():
    simbad = pyvo.dal.TAPService(
      "http://simbad.u-strasbg.fr:80/simbad/sim-tap")
    gavodc = pyvo.dal.TAPService("http://dc.g-vo.org/tap")
  
    # Get blue supergiants from Simbad
    simjob = simbad.submit_job("""
      select main_id, ra, dec
      from basic
      where otype='BlueSG*'""")
    simjob.run()
  
    # Get lightcurves from Gaia
    try:
      simjob.wait()
      time_series = gavodc.run_sync("""
        SELECT b.*, bp_obs_time, bp_flux, rp_obs_time, rp_flux
        FROM (SELECT
           main_id, source_id, g.ra, g.dec
           FROM
          gaia.dr2light as g
           JOIN TAP_UPLOAD.t1 AS tc
           ON (0.002>DISTANCE(tc.ra, tc.dec, g.ra, g.dec))
        OFFSET 0) AS b
        JOIN gaia.dr2epochflux
        USING (source_id)
        """,
        uploads={"t1": simjob.result_uri})
    finally:
      simjob.delete()
  
    # Now plot one after the other
    for row in time_series.table:
      plt.plot(row["bp_obs_time"], row["bp_flux"])
      plt.plot(row["rp_obs_time"], row["rp_flux"])
      plt.show(block=False)
      raw_input("{}; press return for next...".format(row["main_id"]))
      plt.cla()
  
  if __name__=="__main__":
    main()
  

  If you bother to read the code, you'll notice that we transfer the Simbad result directly to the GAVO data center without first downloading it. That's fairly boring in this case, where the table is small. But if you have a narrow pipe for one reason or another and some 10⁵ rows, passing around async result URLs is a useful trick.

  In this particular case the whole thing returns just four stars, so perhaps that's not a terribly useful target for your learning machine. But this piece of code should get you started to where there's more data.

  You should read the column descriptions and footnotes in the query results (or from the reference URL) – this tells you how to interpret the times and how to make magnitudes from the fluxes if you must. You probably can't hear it any more, but just in case: If you can, process fluxes rather than magnitudes from Gaia, because the errors are painful to interpret in magnitudes when the fluxes are small (try it!).

  Note how the photometry data is stored in arrays in the database, and that VOTables can just transport these. The bad news is that support for manipulating arrays in ADQL is pretty much zero at this point; this means that, when you have trained your ML device, you'll probably have to still download lots and lots of light curves rather than write some elegant ADQL to do the filtering server-side. However, I'd be highly interested to work out how some tastefully chosen user defined functions might enable offloading at least a good deal of that analysis to the database. So – if you know what you'd like to do, by all means let me know. Perhaps there's something I can do for you.

  Incidentally, I'll talk a bit more about ADQL arrays in a blog post coming up in a few weeks (I think). Don't miss it, subscribe to our feed).

  Datalink

  In the results from queries involving gaia.dr2epochflux, we also provide datalinks. These let you retrieve lightcurves that already have mags and that are more easily plotted. Perhaps more importantly, they link back to the full ESAC lightcurves that, in addition, give you a lot more debug information and are required if you want to reliably identify photometry points with the identifiers of the transits that generated them.

  Datalink support in clients still is not great, but it's growing nicely. Your ideas for workflows that should be supported are (again) most welcome – and have a good chance of being adopted. So, try things out, for instance by getting the most recent TOPCAT (as of this writing) and do the following:

  1. Open the VO/TAP dialog from the menu bar and double click the GAVO DC TAP service.

  2. Enter:

       SELECT source_id, ra, dec,
       phot_bp_mean_mag, phot_rp_mean_mag, phot_g_mean_mag,
       g_transit_time, g_transit_flux,
       rp_obs_time, rp_flux
       FROM gaia.dr2epochflux
       JOIN gaia.dr2light
       USING (source_id)
       WHERE parallax>50
     
     into “ADQL” text to retrieve lightcurves for the more nearby
     variables (in reality, you'd have to be a bit more careful with the
     distances, but you already knew that).
     

  3. plot something like phot_bp_mean_mag-phot_rp_mean_mag vs. phot_g_mean_mag (and adapt the plot to fit your viewing habits).

  4. Open the dialog for Views/Activation Actions (from the menu bar or the tool bar – same thing), check “Invoke Service”, choose “View Datalink Table”.

  5. Whenever you click on a a point in your CMD, a window will pop up in which you can choose between the time series in the various bands, and you can pull in the data from ESAC; to load a table, select “Load Table” from the actions near the foot of the datalink table and click “Invoke”.

  Yeah. It's clunky. Help us make it better with your fresh ideas for interfaces (and don't be cross with us if we have to marry them with what's technically feasible and readily generalised).

  SSAP and Obscore

  If you're fed up with bleeding-edge tech, the light curves are also available through good old SSAP and Obscore. To use that, just get Splat (or another SSA client, preferably with a bit of time series support). Look for a Gaia DR2 time series service (you may have to update the service list before you find it), enter (in keeping with our LBV theme) S Dor as position and hit “Lookup” followed by “Send Query”. Just click on any result to just view the time series – and then apply Splat's rich tool set to it.

  Update (8.5.2018): Clusters

  Here's another quick application – how about looking for variable stars in clusters? This piece of ADQL should get you started:

  SELECT TOP 100
    source_id, ra, dec, parallax, g.pmra, g.pmdec,
    m.name, m.pmra AS c_pmra, m.pmde AS c_pmde,
    m.e_pm AS c_e_pm,
    1/dist AS cluster_parallax
  FROM
    gaia.dr2epochflux
    JOIN gaia.dr2light AS g USING (source_id)
    JOIN mwsc.main AS m
    ON (1=CONTAINS(
      POINT(g.ra, g.dec),
      CIRCLE(m.raj2000, m.dej2000, rcluster)))
  WHERE IN_UNIT(pmdec, 'deg/yr') BETWEEN m.pmde-m.e_pm*3 AND m.pmde+m.e_pm*3
  

  – yes, you'll want to constrain pmra, too, and the distance, and properly deal with error and all. But you get simple lightcurves for free. Just add them in the SELECT clause!

« Page 2 / 2