add an example showing a vizier overlay on a skyview image
@@ -476,3 +476,63 @@ Retrieve data from a particular co-I or PI from the ESO archive
import shutil
for fn in files:
shutil.move(fn+'.TAR','.')
+
+Example 9
++++++++++
+
+Retrieve an image from skyview and overlay a Vizier catalog on it.
+This example approximately reproduces Figure 1 of 2016ApJ...826...16E,
+except with a different background.::
+
+
+ from astropy import coordinates, units as u, wcs
+ from astroquery.skyview import SkyView
+ from astroquery.vizier import Vizier
+ import pylab as pl
+
+ center = coordinates.SkyCoord.from_name('Orion KL')
+
+ # Grab an image from SkyView of the Orion KL nebula region
+ imglist = SkyView.get_images(position=center, survey='2MASS-J')
+
+ # the returned value is a list of images, but there is only one
+ img = imglist[0]
+
+ # 'img' is now a fits.HDUList object; the 0th entry is the image
+ mywcs = wcs.WCS(img[0].header)
+
+ fig = pl.figure(1)
+ fig.clf() # just in case one was open before
+ # use astropy's wcsaxes tool to create an RA/Dec image
+ ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=mywcs)
+ ax.set_xlabel("RA")
+ ax.set_ylabel("Dec")
+
+ ax.imshow(img[0].data, cmap='gray_r', interpolation='none', origin='lower',
+ norm=pl.matplotlib.colors.LogNorm())
+
+
+ # retrieve a specific table from Vizier to overplot
+ tablelist = Vizier.query_region(center, radius=5*u.arcmin, catalog='J/ApJ/826/16/table1')
+ # again, the result is a list of tables, so we'll get the first one
+ result = tablelist[0]
+
+ # convert the ra/dec entries in the table to astropy coordinates
+ tbl_crds = coordinates.SkyCoord(result['RAJ2000'], result['DEJ2000'],
+ unit=(u.hour, u.deg), frame='fk5')
+
+ # we want this table too:
+ tablelist2 = Vizier(row_limit=10000).query_region(center, radius=5*u.arcmin, catalog='J/ApJ/540/236')
+ result2 = tablelist2[0]
+ tbl_crds2 = coordinates.SkyCoord(result2['RAJ2000'], result2['DEJ2000'],
+ unit=(u.hour, u.deg), frame='fk5')
+
+
+ # overplot the data in the image
+ ax.plot(tbl_crds.ra, tbl_crds.dec, 'o', transform=ax.get_transform('fk5'),
+ mec='b', mfc='none')
+ ax.plot(tbl_crds2.ra, tbl_crds2.dec, 'o', transform=ax.get_transform('fk5'),
+ mec='r', mfc='none')
+ # zoom in on the relevant region
+ ax.axis([117,217,200,300])
+
