Cluster analysis around the Herbig AeBe star HD 200775¶
This use case is based on work by Pérez Blanco+ 2019.
Follow the instructions here to install the pyESASky widget. More information can be found here.
Workflow:
- Query the Vioque+ 2018 catalogue of Herbig Ae/Be stars to obtain the stellar parameters, using the VizieR astroquery module.
- Query the Gaia archive for sources around HD 200775 and identify possible cluster candidates by filtering on proper motions and quality, using the Gaia astroquery module.
- Query the Herschel archive and download a PACS image of the region, using the ESASky astroquery module.
- Overlay the Gaia results on the Herschel image.
- Visualise the results in the pyESASky widget.
Modification history:¶
- Notebook updated to reflect backend changes as of August 2021.
In [3]:
# Import all the required python modules:
import astropy.units as u
from astropy.coordinates.sky_coordinate import SkyCoord
from astropy.units import Quantity
from astroquery.vizier import Vizier
from astroquery.gaia import Gaia
from astroquery.esasky import ESASky
from astropy.wcs import WCS
from astropy.visualization import (MinMaxInterval, SqrtStretch, ImageNormalize, ManualInterval)
from pyesasky.pyesasky import ESASkyWidget
from pyesasky.catalogue import Catalogue
from pyesasky.catalogueDescriptor import CatalogueDescriptor
from pyesasky.metadataDescriptor import MetadataDescriptor
from pyesasky.metadataType import MetadataType
from pyesasky.cooFrame import CooFrame
In [4]:
#%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
Step 1. Query VizieR¶
In [5]:
# Search for 'Gaia DR2 study of Herbig Ae/Be stars (Vioque+, 2018)':
catalog_list = Vizier.find_catalogs('Vioque+,2018')
print({k:v.description for k,v in catalog_list.items()})
In [6]:
# Get the above list of catalogues:
Vizier.ROW_LIMIT = -1
catalogs = Vizier.get_catalogs(catalog_list.keys())
print(catalogs)
In [7]:
# Access the first table:
Vioque = catalogs['J/A+A/620/A128/hqsample']
print(Vioque[54])
Step 2. Query the Gaia archive¶
The following performs an asynchronous query (asynchronous rather than synchronous queries should be performed when retrieving more than 2000 rows) using the Astronomical Data Query Language (ADQL) on the Gaia DR2 catalogue for sources within a search radius of 0.2 degrees around the Herbig Ae/Be star, firstly with no quality filtering on the catalogue.
See here for more information and examples of ADQL queries.
In [8]:
tables = Gaia.load_tables(only_names=True)
In [9]:
# Do the following to load and look at the available Gaia table names:
for table in (tables):
print (table.get_qualified_name())
In [10]:
tables = Gaia.load_tables(only_names=True)
job = Gaia.launch_job_async("SELECT * FROM gaiadr2.gaia_source \
WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),CIRCLE('ICRS',315.40383,68.16327,0.2))=1 \
;", dump_to_file=False)
In [11]:
g = job.get_results()
print (g['source_id', 'pmra', 'pmdec', 'parallax'])
In [12]:
# Plot the results in a proper motion plot of proper motion in RA (pmra) versus proper motion in DEC (pmdec)
# in the range pmra [-20,20] and pmdec [-20,20]
plt.figure(figsize=(10, 8))
graph = plt.scatter(g['pmra'], g['pmdec'], alpha=0.9, c=g['parallax'], cmap=plt.cm.coolwarm)
cb = plt.colorbar(graph)
cb.set_label('parallax',fontsize=14)
#plt.scatter(g['pmra'], g['pmdec'],color='r',alpha=0.3)
plt.xlim(-20,20)
plt.ylim(-20,20)
plt.title('Gaia DR2 sources proper motion plot around HD 200775',fontsize=16)
plt.xlabel(r'pmRA',fontsize=14)
plt.ylabel(r'pmDEC',fontsize=14)
plt.show()
In [13]:
# Retrieve the Gaia parameters parallax, pmra and pmdec for HD 200775. From Vioque+ 2018 we know that
# parallax = 2.770664
job2 = Gaia.launch_job_async("SELECT ra, dec, parallax, pmra, pmdec FROM gaiadr2.gaia_source \
WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec),CIRCLE('ICRS',315.40383,68.16327,0.2))=1 \
AND parallax BETWEEN 2.770660 AND 2.770670 ;", dump_to_file=False)
r = job2.get_results()
print (r['dec', 'parallax', 'pmra', 'pmdec'])
In [14]:
# Run a second asynchronous query on the Gaia archive, filtering the results by quality, proper motion and parallax:
job3 = Gaia.launch_job_async("SELECT * FROM gaiadr2.gaia_source as gaia \
INNER JOIN gaiadr2.ruwe \
ON gaia.source_id = gaiadr2.ruwe.source_id \
WHERE CONTAINS(POINT('ICRS',gaia.ra,gaia.dec),CIRCLE('ICRS',315.40383,68.16327,0.2))=1 \
AND gaia.phot_g_mean_mag < 19 \
AND gaia.parallax_over_error > 10 \
AND gaia.phot_g_mean_flux_over_error > 50 \
AND gaia.visibility_periods_used > 8 \
AND gaiadr2.ruwe.ruwe < 1.4 \
AND gaia.parallax BETWEEN 2 AND 3.5 \
AND gaia.pmra BETWEEN 5 AND 11 \
AND gaia.pmdec BETWEEN -4 AND 1;", dump_to_file=False)
candidates = job3.get_results()
len(candidates)
Out[14]:
In [15]:
plt.figure(figsize=(10, 8))
plt.scatter(candidates['pmra'], candidates['pmdec'], s=40, edgecolor='g', facecolor='none')
graph = plt.scatter(g['pmra'], g['pmdec'], alpha=0.9, c=g['parallax'], cmap=plt.cm.coolwarm)
cb = plt.colorbar(graph)
cb.set_label('parallax',fontsize=14)
plt.scatter(candidates['pmra'], candidates['pmdec'], s=40, edgecolor='g', facecolor='none')
plt.xlim(-20,20)
plt.ylim(-20,20)
plt.title('Gaia DR2 sources proper motion plot around HD 200775',fontsize=16)
plt.xlabel(r'pmRA',fontsize=14)
plt.ylabel(r'pmDEC',fontsize=14)
plt.legend(["HD 200775 cluster candidates"])
plt.show()
Step 3. Query the Herschel archive¶
In [16]:
# Query ESASky for the available Herschel maps
maps = ESASky.query_object_maps(position='HD 200775', missions=['HERSCHEL'])
print (maps)
In [17]:
# Inspect the table
maps['HERSCHEL'].info
Out[17]:
In [18]:
# Inspect the columns
maps['HERSCHEL']['observation_id', 'instrument', 'filter'].pprint()
In [19]:
# Download the images
maps_data = ESASky.get_maps(query_table_list=maps,missions='HERSCHEL')
In [20]:
# Inspection of the Herschel PACS 70um header
her_hdu = maps_data["HERSCHEL"][0]["70"]
her_hdu.info()
her_hdu[0].header
Out[20]:
Step 4. Overlay the Gaia sources on the Herschel image¶
In [21]:
# Overlaying the Gaia sources on the Herschel image
her_image = her_hdu['image'].data
norm = ImageNormalize(her_image,interval = ManualInterval(-0.05,0.3)) #play around with the ManualInterval numbers to change the image contrast
fig = plt.figure(figsize=(10,10),dpi=100)
wcs_h = WCS(her_hdu['image'].header)
ax = fig.add_subplot(111,projection=wcs_h)
ax.set_title("HD 200775 region in Herschel PACS 70 µm")
ax.imshow(her_image,cmap=plt.cm.copper,origin='lower',interpolation='nearest', norm=norm)
p1 = ax.scatter(g['ra'],g['dec'],transform=ax.get_transform('world'), \
s=30, edgecolor='#1E90FF', facecolor='none', label='Gaia DR2 sources')
p2 = ax.scatter(candidates['ra'],candidates['dec'],transform=ax.get_transform('world'), \
s=60, edgecolor='#00ff00', facecolor='none', label='HD 200775 cluster candidates')
ax.set_xlabel("RA")
ax.set_ylabel("Dec")
ax.legend(["Gaia DR2 sources","HD 200775 cluster candidates"])
Out[21]:
In [22]:
# Overlaying the Gaia sources, cluster candidates and their proper motions on the Herschel image
wcs_h = WCS(her_hdu['image'].header)
wcs_h
fig = plt.figure(figsize=(8,8),dpi=100)
ax = fig.add_subplot(111,projection=wcs_h)
ax.set_title("HD 200775 region in Herschel PACS 70 µm")
ax.imshow(her_image,cmap=plt.cm.gray,origin='lower',interpolation='nearest', norm=norm)
p1 = ax.scatter(g['ra'],g['dec'],transform=ax.get_transform('world'), \
s=30, edgecolor='#1E90FF', facecolor='none', label='Gaia DR2 sources')
p2 = ax.scatter(candidates['ra'],candidates['dec'],transform=ax.get_transform('world'), \
s=60, edgecolor='#00ff00', facecolor='none', label='HD 200775 cluster candidates')
Q = ax.quiver(np.array(candidates['ra']),np.array(candidates['dec']),-np.array(candidates['pmra']),np.array(candidates['pmdec']),
transform=ax.get_transform('world'),color='plum',alpha=0.55) # Define the X-Y axis (ra,dec) and the arrow direction (pmra,pmdec)
ax.set_xlabel("RA")
ax.set_ylabel("Dec")
ax.legend(["Gaia DR2 sources","HD 200775 cluster candidates"])
Out[22]:
Step 5. Visualise the results in pyESASky¶
In [23]:
esasky = ESASkyWidget()
esasky
In [29]:
# Go to HD 200775
esasky.goToTargetName('HD 200775')
In [31]:
# Zoom to 15 arcminutes
esasky.setFoV(0.5)
In [28]:
# Set the background HiPS to Herschel PACS colour:
#esasky.setHiPS("Herschel","https://cdn.skies.esac.esa.int/Herschel/PACS-color/")
# Set the background HiPS to 2MASS colour:
#esasky.setHiPS("2MASS","https://cdn.skies.esac.esa.int/2MASSColor/")
# Set the background HiPS to AllWISE colour:
#esasky.setHiPS("AllWISE","https://cdn.skies.esac.esa.int/AllWISEColor/")
# Set the background HiPS to Spitzer colour:
esasky.setHiPS("Spitzer","https://cdn.skies.esac.esa.int/Spitzer/IRAC134_RGB_bright/")
In [ ]:
#Overlay the Gaia DR2 cluster candidate sources
# overlayCatalogueFromAstropyTable('<catName>', '<cooFrame>', <HTMLcolor>, '<(astropy.table)>', '<RAColumn>', '<DecColumn>', '<nameColumn>')
# where:
# - <catName> : name of the catalogue that will be used in pyESASky as label
# - <cooFRAME> : coordinate frame
# - <HTMLcolor> : HTML color. It could be a "Color name", "Hex color code" or "RGB color code"
# - <lineWidth> : width used to draw sources. From 1 to 10
# - <idColumn> : name of the column containing a unique identifier for sources if any. None if not applicable
# - <nameColumn> : name of the column with the name of the source
# - <RAColumn> : name of the RA column in degrees
# - <DecColumn> : name of the Dec column in degrees
esasky.overlayCatalogueFromAstropyTable('Candidates', 'J2000', '#00ff00', 5, candidates, '','','')
All the normal ESASky functionality is available in the pyESASky widget. Select the catalogues icon and Gaia DR2 to plot all Gaia DR2 sources in the field of view. Change the shape and increase the size of the sources via selecting the colour icon in the results tab.