DOC: adding docs and plot for Irsa.query_sia()

Author: bsipocz

docs/ipac/irsa/irsa.rst

@@ -220,6 +220,62 @@ returned as a `~pyvo.dal.TAPResults` object. Its ``to_table`` or ``to_qtable`` m
     202.809023  46.964558  15.874     0.081 ...  15.322     0.188     AAC    000
 
 
+Simple image access queries
+---------------------------
+
+`~astroquery.ipac.irsa.query_sia` provides a way to access IRSA's Simple
+Image Access VO service. In the following example we are looking for Spitzer
+Enhanced Imaging products in the centre of the COSMOS field as an `~astropy.table.Table`.
+
+.. doctest-remote-data::
+
+   >>> from astroquery.ipac.irsa import Irsa
+   >>> from astropy.coordinates import SkyCoord
+   >>> from astropy import units as u
+   >>>
+   >>> coord = SkyCoord('150.01d 2.2d', frame='icrs')
+   >>> spitzer_images = Irsa.query_sia(pos=(coord, 1 * u.arcmin), collection='spitzer_seip').to_table()
+
+Now open a cutout image for one of the science images. You could either use
+the the IRSA on-premise data or the cloud version of it using the
+``access_url`` or ``cloud_access`` columns:
+
+.. doctest-remote-data::
+
+   >>> from astropy.io import fits
+   >>> from astropy.nddata import Cutout2D
+   >>> from astropy.wcs import WCS
+   >>> science_image = spitzer_images[spitzer_images['dataproduct_subtype'] == 'science'][0]
+   >>> with fits.open(science_image['access_url'], use_fsspec=True) as hdul:
+   ...     cutout = Cutout2D(hdul[0].section, position=coord, size=2 * u.arcmin, wcs=WCS(hdul[0].header))
+
+Now plot the cutout.
+
+.. doctest-skip::
+
+   >>> import matplotlib.pyplot as plt
+   >>> plt.imshow(cutout.data, cmap='grey')
+   >>> plt.show()
+
+.. plot::
+
+   from astroquery.ipac.irsa import Irsa
+   from astropy.coordinates import SkyCoord
+   from astropy import units as u
+   from astropy.io import fits
+   from astropy.nddata import Cutout2D
+   from astropy.wcs import WCS
+   import matplotlib.pyplot as plt
+   coord = SkyCoord('150.01d 2.2d', frame='icrs')
+   spitzer_images = Irsa.query_sia(pos=(coord, 1 * u.arcmin), collection='spitzer_seip').to_table()
+   science_image = spitzer_images[spitzer_images['dataproduct_subtype'] == 'science'][0]
+   with fits.open(science_image['access_url'], use_fsspec=True) as hdul:
+        cutout = Cutout2D(hdul[0].section, position=coord, size=2 *
+   u.arcmin, wcs=WCS(hdul[0].header))
+   plt.imshow(cutout.data, cmap='grey')
+   plt.show()
+
+
 Other Configurations
 --------------------