Example ESASky use¶
This notebook illustrates some relatively more advanced usages of the ESASky implementation as a module of astroquery.
Authors: Ivan Valtchanov, Nuria Álvarez and Henrik Norman.
First you need to install astroquery which contains esasky.
Astroquery can be installed with pip install --pre astroquery, the latest version should come with esasky.
The documentation for astroquery.esasky is available here.
Modification history:¶
- Notebook updated to reflect backend changes as of July 2021.
#
# import some necessary packages
#
import matplotlib.pyplot as plt
from matplotlib.legend_handler import HandlerLine2D
from astropy.io import fits
%matplotlib inline
import numpy as np
from astropy.wcs import WCS
from astropy.visualization import PercentileInterval, ImageNormalize, AsinhStretch, LogStretch, LinearStretch
from astroquery.esasky import ESASky
from astropy.convolution import convolve, Kernel, Gaussian2DKernel
Workflow¶
This is a simple workflow.
- Load an XMM-Newton EPIC image of a target of choice.
- Smooth the image with a Gaussian with FWHM=5".
- Generate 10 contours.
- Plot the original photon image in [0.2-7] keV (ESASky data) and overplot the contours.
- Load a Herschel 250 µm image of the same target.
- Plot the Herschel image and overplot the XMM smoothed image contours.
- Overplot the Planck PSZ catalogue sources.
ESASky.list_maps()
maps = ESASky.query_object_maps(position='M51', missions=["XMM","HERSCHEL"])
print (maps)
ESASky.list_catalogs()
#
# get information on the content of the results table returned by the query
#
print(maps["XMM"].info)
maps_list = ESASky.list_maps()
maps_list
#
# and list the observations and their duration
#
maps["XMM"]["observation_id","instrument","duration"].pprint()
We only need one of the XMM observation, let's take the longest one: 0830191401, it is with index 9.
# need to copy if you don't want to modify the original maps
nxmm = len(maps["XMM"])
maps["XMM"].remove_rows(np.delete(range(nxmm),9))
maps["XMM"]["observation_id","duration"].pprint()
Now let's select the Herschel 250 µm map before we submit the download request to ESASky.
print (maps["HERSCHEL"].info)
maps["HERSCHEL"]["observation_id","filter","duration"].pprint()
Let's pick up SPIRE observation 1342188589, this is index 6.
nher = len(maps["HERSCHEL"])
maps["HERSCHEL"].remove_rows(np.delete(range(nher),6))
maps["HERSCHEL"]["observation_id","filter","duration"].pprint()
And we are ready to download the maps. They will be available in the memory but also saved to the default folder.
maps_data = ESASky.get_maps(query_table_list=maps)
Let's have a look at the Herschel maps_data. As we only kept one observation then there will be only one set of maps (index=0) and they will be stored in a dictionary with keys '250', '350' and '500'.
print (maps_data["HERSCHEL"][0].keys())
And then the XMM maps_data, which is also for one entry (index=0) but it's only one map in a Primary HDU.
print (maps_data["XMM"][0])
Let's check the HDUs.
#
# XMM
#
xmm_hdu = maps_data["XMM"][0]
xmm_hdu.info()
The XMM image is 3-D with different 3-D extensions. The content of the 3-color FITS files from the XMM pipeline is explained in the XMM Pipeline Specification Document. The primary extension is the background and exposure corrected image of MOS1+MOS2+pn; RAWIMAGE is the count-rate image, FWCIMAGE is the background image (Filter-Wheel-Closed data), OOTIMAGE is the Out-of-Time events image.
The PRIMARY extension contains images with three energy bands: 0 => [0.3,0.7] keV, 1 => [0.7,1.2] keV and 2 => [1.2,7.0] keV. These are produced with etruecolor XMM-SAS task.
For the tests we'll use the second slice (the "green" image, [0.7,1.2] keV) and the RAWIMAGE extension.
Note The FITS HDU info displays the dimensions of each extension as (648,648,3), but when we extract slices (i.e. per band) in python the band dimension is the first one, i.e. to get the image in band [0.3,0.7] keV use hdu['RAWIMAGE'].data[0,:,:].
#
# Herschel
#
her_hdu = maps_data["HERSCHEL"][0]["250"]
her_hdu.info()
The Herschel HDU contains a lot of extensions, with image, error and coverage.
Now we have the maps loaded in the session and we can start with some of the points in the workflow.
- Smooth the XMM-EPIC map with a Gaussian with FWHM=5"
Because the data is 3-D, the WCS is also 3-D and we have to slice it and extract the sliced WCS.
# need to crate the Gaussian kernel, it needs the st.dev. in pixels, so we need to get the WCS of the image first
wcs = WCS(xmm_hdu[0].header)
#
# slice and extract the WCS
#
w = wcs.dropaxis(2)
pix = w.wcs.cdelt[1]*3600.0 # pixel size in arcsec
#
g_image = xmm_hdu['RAWIMAGE'].data[1,:,:]
#
fwhm = 15.0/pix
conv = np.sqrt(8.0*np.log(2.0))
stdev = fwhm/conv
gauss = Gaussian2DKernel(x_stddev=stdev, y_stddev=stdev)
xmm_sm5 = convolve(g_image,gauss)
Let's visualise the "green" ([0.7,1.2] keV) XMM RAWIMAGE extension and add contours from the smoothed image.
fig = plt.figure(figsize=(10,10),dpi=100)
pp = 98.0 # colour cut percentage
ax = fig.add_subplot(121,projection=w)
ax.set_title("XMM primary")
norm = ImageNormalize(g_image,interval=PercentileInterval(pp), stretch=LinearStretch())
ax.imshow(g_image,cmap=plt.cm.gray,origin='lower',norm=norm,interpolation='nearest')
ay = fig.add_subplot(122,projection=w)
ay.set_title("XMM smoothed and contour")
norm = ImageNormalize(xmm_sm5,interval=PercentileInterval(pp), stretch=LinearStretch())
ay.imshow(xmm_sm5,cmap=plt.cm.gray,norm=norm,origin='lower',interpolation='nearest')
#print (np.min(xmm_sm5),np.max(xmm_sm5))
# now plot the contours
ay.contour(xmm_sm5, transform=ay.get_transform(w),
levels=np.logspace(-4,-1,10),colors='yellow')
Now, let's show the Herschel 250 µm image.
fig = plt.figure(figsize=(10,10),dpi=100)
pp = 98.0 # colour cut percentage
wcs_h = WCS(her_hdu['image'].header)
ax = fig.add_subplot(111,projection=wcs_h)
ax.set_title("Herschel SPIRE 250 µm")
norm_her = ImageNormalize(her_hdu['image'].data[~np.isnan(her_hdu['image'].data)],interval=PercentileInterval(pp), stretch=AsinhStretch())
ax.imshow(her_hdu['image'].data,cmap=plt.cm.gray,norm=norm_her,origin='lower',interpolation='nearest')
cs = ax.contour(xmm_sm5, transform=ax.get_transform(w),
levels=np.logspace(-4,-1,10),colors='yellow')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
#
ax.legend([q for q in cs.collections],["XMM [0.7-1.2] keV"]);
Now, let's load some catalogues.
cats = ESASky.query_region_catalogs(position='M51', radius='10 arcmin')
print (cats)
catx= cats['PLANCK-PCCS2-HFI']
print(catx.info)
Select only the 857 GHz sources.
ix857 = np.where(catx['frequency'] == 857)[0]
src857 = catx[ix857]
src857
Now let's combine everything into a final figure.
fig = plt.figure(figsize=(10,10),dpi=100)
pp = 95.0 # colour cut percentage
wcs_h = WCS(her_hdu['image'].header)
ax = fig.add_subplot(111,projection=wcs_h)
ax.set_title("Herschel SPIRE 250 µm")
ax.imshow(her_hdu['image'].data,cmap=plt.cm.gray,norm=norm_her,origin='lower',interpolation='nearest')
cs = ax.contour(xmm_sm5, transform=ax.get_transform(w),
levels=np.logspace(-4,-1,10),colors='yellow')
p1 = ax.scatter(src857['ra'],src857['dec'],transform=ax.get_transform('world'), \
s=300, edgecolor='red', facecolor='none', label='Planck PCCS2-HFI')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.legend([[q for q in cs.collections][0],p1],["XMM [0.2-7] keV","Planck PCCS2-HFI857"]);