This notebook illustrates some example usages of the ESASky implementation in astroquery.
First you need to install astroquery and esasky.
Astroquery can be installed with pip install astroquery, the latest version should come with esasky. Alternatively, you can grab the latest astroquery with esasky from here.
The documentation for astroquery.esasky is available here.
Use Case 1: Retrieve imaging data for a single object¶
In this use case, imaging data are retrieved for a single object, indicated by its name (resolved by Simbad) or coordinates.
We start by importing the ESASky astroquery module and other necessary packages:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from astroquery.esasky import ESASky
First, let's check the available maps:
ESASky.list_maps()
Let's do a search around the position of our object:
maps = ESASky.query_object_maps('M51')
print (maps)
The output is a TableList with the keys corresponding to the mission (or in some cases, mission and instrument) names for which there are images available covering our target position.
We can also do a search by coordinates:
maps = ESASky.query_object_maps('13 29 52.7 +47 11 43')
print (maps)
The list of accepted coordinate formats can be found here. The method has a tolerance of 5 arcsec to allow for positional errors.
Let's check the content of the 'XMM-OM-OPTICAL' table:
maps['XMM-OM-UV'].info
As we see, it contains a table with URLs for the maps, and some additional info, like map ra, dec, bandpass ('filter'), or the exposure time ('duration').
Let's now see what bands are available:
maps['XMM-OM-OPTICAL']['observation_id', 'instrument', 'filter', 'duration'].pprint()
The result is stored in memory. To get the actual images:
maps_data = ESASky.get_maps(maps,'XMM-OM-OPTICAL', download_dir = '/Users/belen/tmp')
maps_data
The output maps_data is a product containing all the FITS files in memory, which are also downloaded to a local folder. If no folder is specified, the current working directory is used by default.
It is possible to download all the available images at once by typing 'all' as the mission name in the ESASky.get_maps() method.
We are ready to work with these maps. For example, let's inspect the header of one of them:
hdu = maps_data["XMM-OM-OPTICAL"][3]
hdu[0].header
Let's now display it:
image = hdu[0].data
plt.imshow(image, cmap='gray', norm=LogNorm())
plt.colorbar()
Use Case 2: Retrieve catalogue data for a single object¶
We will now inspect and retrieve the catalogue data available for a given object, using either a name (resolved by Simbad) or coordinates.
As we did in Use Case 1 with the images, let's first get a list of the available catalogues:
ESASky.list_catalogs()
We can look for an object in these catalogues by name or coordinates:
cats = ESASky.query_object_catalogs('M51')
print (cats)
cats = ESASky.query_object_catalogs('13 29 52.7 +47 11 43')
print (cats)
As in the previous use case, the query results are stored in a TableList with the keys corresponding to the catalog name.
It is also possible to specify the catalogues to search:
cats = ESASky.query_object_catalogs('M51',["XMM-EPIC","XMM-SLEW", "HSC"])
print (cats)
In this example, we only get two tables, because there are no available sources in the XMM-SLEW catalogue.
Let's now visualise the XMM-EPIC results table:
hsc_table = cats["HSC"]
hsc_table.info()
print (hsc_table)
We can also choose the columns to display:
print (hsc_table['name', 'ra', 'dec'])
Let's plot these sources:
plt.xlabel('RA (deg)')
plt.ylabel('DEC (deg)')
plt.xlim(reversed(plt.xlim(202.467,202.472)))
plt.ylim(47.193,47.197)
plt.scatter(hsc_table['ra'], hsc_table['dec'])
plt.grid(True)
Use Case 3: Retrieve data in a sky region¶
It is also possible retrieve imaging and catalogue data in a circular sky region. The procedure is very similar to the case of a single object discussed above, but using the query_region_maps() and query_region_catalogs() methods instead of query_object_maps() and query_object_catalogs().
Let's first see the case where imaging data are retrieved. To query the region, we have to enter the central coordinates and the radius:
maps = ESASky.query_region_maps('13 29 52.7 +47 11 43', '25 arcmin')
print (maps)
The difference with the query_object_maps() module is that now we have to explicitly indicate the search radius; if not, an error message is returned.
To search in a region around a given object, we can also enter the object's name instead of its coordinates:
maps = ESASky.query_region_maps('M51', '25 arcmin')
print (maps)
If we compare the result of this query with that of Use Case 1, we will see that now we get more data, because we are using a larger search radius (recall that the search by object uses a radius of only 5 arcsec). In particular, we now have INTEGRAL data, and the number of observations from the other missions has increased.
The query for catalogues works similarly:
cats = ESASky.query_region_catalogs('M51',"25 arcmin",["XMM-EPIC", "HSC"])
print (cats)
We also note here how the number of rows in each table is larger than before, because of the larger search radius.
By default, the query returns a maximum of 10,000 rows per table. Note that this limit has been reached for the HSC. If we want to make sure that all the rows are retrieved, we can set the row_limit parameter to -1:
cats = ESASky.query_region_catalogs('M51',"25 arcmin",["XMM-EPIC", "HSC"], row_limit=-1)
print (cats)
Now the number of rows in the HSC table is 100,000, which is the absolute maximum that can be queried.
Let's now see how the XMM-EPIC table looks like:
xmm_epic_table = cats["XMM-EPIC"]
print (xmm_epic_table)
We are ready to work with the data. For example, let's plot a histogram of the EPIC fluxes:
NBINS = 100
plt.xlim(0, 3e-13)
plt.ylim(0,300)
flux = plt.hist(xmm_epic_table['ep_8_flux'], NBINS)
