esasky.astroquery advanced use case¶
This notebook illustrates some more advanced usages of the ESASky implementation in astroquery
Author: Ivan Valtchanov¶
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.
#
# import some necessary packages
#
import matplotlib.pyplot as plt
from matplotlib.legend_handler import HandlerLine2D
%matplotlib inline
import sys
reload(sys)
sys.setdefaultencoding('utf8')
import numpy as np
from astropy.wcs import WCS
from astropy import visualization
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 the 3XMM point-source catalogue sources
- Overplot the Planck PSZ catalogue sources
# ESASky.list_maps()
# ESASky.list_catalogs()
maps = ESASky.query_object_maps('M51',["XMM-EPIC","HERSCHEL"])
print (maps)
print (maps["XMM-EPIC"].info)
maps["XMM-EPIC"]["observation_id","duration"].pprint()
We only need one of the 6 XMM-EPIC observation, let's take the longest one: 0303420101, it is with index 3.
# need to copy if you don't want to modify the original maps
nxmm = len(maps["XMM-EPIC"])
maps["XMM-EPIC"].remove_rows(np.delete(range(nxmm),3))
maps["XMM-EPIC"]["observation_id","duration"].pprint()
Now let's select the Herschel 250 µm map too, 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 7.
nher = len(maps["HERSCHEL"])
maps["HERSCHEL"].remove_rows(np.delete(range(nher),7))
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 deafult folder.
maps_data = ESASky.get_maps(maps)
Let's have a look at the maps_data
print (maps_data["HERSCHEL"])
print (maps_data["XMM-EPIC"])
xmm_hdu = maps_data["XMM-EPIC"][0]
xmm_hdu.info()
xmm_hdu[0].header
her_hdu = maps_data["HERSCHEL"][0]["250"]
her_hdu.info()
her_hdu[0].header
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"
# need to create the Gaussian kernel, it needs the st.dev. in pixels, so we need ot get the WCS of the image first
wcs = WCS(xmm_hdu[0].header)
pix = wcs.wcs.cdelt[1]*3600.0 # pixel size in arcsec
fwhm = 15/pix
conv = np.sqrt(8.0*np.log(2.0))
stdev = fwhm/conv
gauss = Gaussian2DKernel(stddev=stdev)
xmm_sm5 = convolve(xmm_hdu[0].data,gauss)
fig = plt.figure(figsize=(10,10),dpi=100)
pp = 98.0 # colour cut percentage
ax = fig.add_subplot(121,projection=wcs)
ax.set_title("XMM raw")
img_scaled = visualization.scale_image(xmm_hdu[0].data,percent=pp)
ax.imshow(img_scaled,cmap=plt.cm.gray,origin='lower',interpolation='nearest')
ay = fig.add_subplot(122,projection=wcs)
ay.set_title("XMM smoothed and contour")
img_scaled2 = visualization.scale_image(xmm_sm5,percent=pp)
ay.imshow(img_scaled2,cmap=plt.cm.gray,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(WCS(xmm_hdu[0].header)),
levels=np.linspace(3,100,20),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")
img_scaled = visualization.scale_image(her_hdu['image'].data,percent=pp)
ax.imshow(img_scaled,cmap=plt.cm.gray,origin='lower',interpolation='nearest')
cs = ax.contour(xmm_sm5, transform=ax.get_transform(wcs),
levels=np.linspace(3,50,10),colors='yellow', label="XMM [0.2-7] keV")
plt.xlabel = 'RA'
plt.ylabel = 'Dec'
ax.legend([q for q in cs.collections],["XMM EPIC"])
cats = ESASky.query_region_catalogs('M51','10 arcmin')
print (cats)
catx= cats['PLANCK-PCCS2-HFI']
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")
img_scaled = visualization.scale_image(her_hdu['image'].data,percent=pp)
ax.imshow(img_scaled,cmap=plt.cm.gray,origin='lower',interpolation='nearest')
cs = ax.contour(xmm_sm5, transform=ax.get_transform(wcs),
levels=np.linspace(6,50,10),colors='yellow', label="XMM [0.2-7] keV")
p1 = ax.scatter(src857['ra'],src857['dec'],transform=ax.get_transform('world'), \
s=300, edgecolor='red', facecolor='none', label='Planck PCCS2-HFI')
plt.xlabel = 'RA'
plt.ylabel = 'Dec'
ax.legend([[q for q in cs.collections][0],p1],["XMM EPIC","Planck PCCS2-HFI857"])
Final notes¶
This simple example illustrated a simple usage of the esasky module wihtin astroquery and showed how to use ESA's archival data on XMM-EPIC, Herschel-SPIRE and Planck point source catalogue.
The final figure in the example shows some potential quality issues with the Planck PCCS2-HFI as the sources seem to be offset from the centre of M51 (NGC 5194) and its sattelite M51B (NGC 5195).
There is a relatively good match between the SPIRE 250 µm map and the XMM-EPIC contours.