This notebook shows how to access the Hubble Catalogs of Variables (HCV). The HCV is a large catalog of faint variable objects extracted from version 3 of the Hubble Source Catalog. The HCV project at the National Observatory of Athens was funded by the European Space Agency (PI: Alceste Bonanos). The data products for the HCV are available both at the ESA Hubble Archive at ESAC through the HCV Explorer interface and at STScI.
This notebook uses the MAST HSC catalog interface, which supports queries to the current and previous versions of the Hubble Source Catalog. It allows searches of several different tables including the HCV summary and detailed tables. It also has an associated API that is used for data access in this notebook.
For similar examples using the MAST CasJobs interface, a SQL database query interface that is more complex to use but more powerful than the API, see the HCV_CasJobs_demo notebook.
This notebook is available for download.
Running the notebook from top to bottom takes less than 1 minute (depending on the speed of your computer and network connection).
This notebook requires the use of Python 3.
This needs some special modules in addition to the common requirements of astropy
, numpy
and scipy
. For anaconda versions of Python the installation commands are:
conda install requests conda install pillow
Run the commands one at a time since conda may ask for confirmation.
%matplotlib inline
import astropy, pylab, time, sys, os, requests, json
import numpy as np
from PIL import Image
from io import BytesIO
from astropy.table import Table, join
from astropy.io import ascii
# Set page width to fill browser for longer output lines
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
# set width for pprint
astropy.conf.max_width = 150
hcvcone(ra,dec,radius [,keywords])
function searches the HCV catalog near a position.hcvsearch()
function performs general non-positional queries.hcvmetadata()
function gives information about the columns available in a table. resolve(name)
function uses the MAST Name Resolver (which relies on both SIMBAD and NED) to get the RA,Dec position for an object.hscapiurl = "https://catalogs.mast.stsci.edu/api/v0.1/hsc"
def hcvcone(ra,dec,radius,table="hcvsummary",release="v3",format="csv",magtype="magaper2",
columns=None, baseurl=hscapiurl, verbose=False,
**kw):
"""Do a cone search of the HSC catalog (including the HCV)
Parameters
----------
ra (float): (degrees) J2000 Right Ascension
dec (float): (degrees) J2000 Declination
radius (float): (degrees) Search radius (<= 0.5 degrees)
table (string): hcvsummary, hcv, summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
format: csv, votable, json
columns: list of column names to include (None means use defaults)
baseurl: base URL for the request
verbose: print info about request
**kw: other parameters (e.g., 'numimages.gte':2)
"""
data = kw.copy()
data['ra'] = ra
data['dec'] = dec
data['radius'] = radius
return hcvsearch(table=table,release=release,format=format,magtype=magtype,
columns=columns,baseurl=baseurl,verbose=verbose,**data)
def hcvsearch(table="hcvsummary",release="v3",magtype="magaper2",format="csv",
columns=None, baseurl=hscapiurl, verbose=False,
**kw):
"""Do a general search of the HSC catalog (possibly without ra/dec/radius)
Parameters
----------
table (string): hcvsummary, hcv, summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
format: csv, votable, json
columns: list of column names to include (None means use defaults)
baseurl: base URL for the request
verbose: print info about request
**kw: other parameters (e.g., 'numimages.gte':2). Note this is required!
"""
data = kw.copy()
if not data:
raise ValueError("You must specify some parameters for search")
if format not in ("csv","votable","json"):
raise ValueError("Bad value for format")
url = "{}.{}".format(cat2url(table,release,magtype,baseurl=baseurl),format)
if columns:
# check that column values are legal
# create a dictionary to speed this up
dcols = {}
for col in hcvmetadata(table,release,magtype)['name']:
dcols[col.lower()] = 1
badcols = []
for col in columns:
if col.lower().strip() not in dcols:
badcols.append(col)
if badcols:
raise ValueError('Some columns not found in table: {}'.format(', '.join(badcols)))
# two different ways to specify a list of column values in the API
# data['columns'] = columns
data['columns'] = '[{}]'.format(','.join(columns))
# either get or post works
# r = requests.post(url, data=data)
r = requests.get(url, params=data)
if verbose:
print(r.url)
r.raise_for_status()
if format == "json":
return r.json()
else:
return r.text
def hcvmetadata(table="hcvsummary",release="v3",magtype="magaper2",baseurl=hscapiurl):
"""Return metadata for the specified catalog and table
Parameters
----------
table (string): hcvsummary, hcv, summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
baseurl: base URL for the request
Returns an astropy table with columns name, type, description
"""
url = "{}/metadata".format(cat2url(table,release,magtype,baseurl=baseurl))
r = requests.get(url)
r.raise_for_status()
v = r.json()
# convert to astropy table
tab = Table(rows=[(x['name'],x['type'],x['description']) for x in v],
names=('name','type','description'))
return tab
def cat2url(table="hcvsummary",release="v3",magtype="magaper2",baseurl=hscapiurl):
"""Return URL for the specified catalog and table
Parameters
----------
table (string): hcvsummary, hcv, summary, detailed, propermotions, or sourcepositions
release (string): v3 or v2
magtype (string): magaper2 or magauto (only applies to summary table)
baseurl: base URL for the request
Returns a string with the base URL for this request
"""
checklegal(table,release,magtype)
if table == "summary":
url = "{baseurl}/{release}/{table}/{magtype}".format(**locals())
else:
url = "{baseurl}/{release}/{table}".format(**locals())
return url
def checklegal(table,release,magtype):
"""Checks if this combination of table, release and magtype is acceptable
Raises a ValueError exception if there is problem
"""
releaselist = ("v2", "v3")
if release not in releaselist:
raise ValueError("Bad value for release (must be one of {})".format(
', '.join(releaselist)))
if release=="v2":
tablelist = ("summary", "detailed")
else:
tablelist = ("summary", "detailed", "propermotions", "sourcepositions",
"hcvsummary", "hcv")
if table not in tablelist:
raise ValueError("Bad value for table (for {} must be one of {})".format(
release, ", ".join(tablelist)))
if table == "summary":
magtypelist = ("magaper2", "magauto")
if magtype not in magtypelist:
raise ValueError("Bad value for magtype (must be one of {})".format(
", ".join(magtypelist)))
def resolve(name):
"""Get the RA and Dec for an object using the MAST name resolver
Parameters
----------
name (str): Name of object
Returns RA, Dec tuple with position
"""
resolverRequest = {'service':'Mast.Name.Lookup',
'params':{'input':name,
'format':'json'
},
}
resolvedObjectString = mastQuery(resolverRequest)
resolvedObject = json.loads(resolvedObjectString)
# The resolver returns a variety of information about the resolved object,
# however for our purposes all we need are the RA and Dec
try:
objRa = resolvedObject['resolvedCoordinate'][0]['ra']
objDec = resolvedObject['resolvedCoordinate'][0]['decl']
except IndexError as e:
raise ValueError("Unknown object '{}'".format(name))
return (objRa, objDec)
def mastQuery(request, url='https://mast.stsci.edu/api/v0/invoke'):
"""Perform a MAST query.
Parameters
----------
request (dictionary): The MAST request json object
url (string): The service URL
Returns the returned data content
"""
# Encoding the request as a json string
requestString = json.dumps(request)
r = requests.post(url, data={'request': requestString})
r.raise_for_status()
return r.text
target = 'IC 1613'
ra, dec = resolve(target)
print(target,ra,dec)
This searches the HCV summary table for objects within 0.5 degrees of the galaxy center. Note that this returns both variable and non-variable objects.
radius = 0.5 # degrees
t0 = time.time()
tab = ascii.read(hcvcone(ra,dec,radius,table="hcvsummary"))
print("Completed in {:.1f} sec".format(time.time()-t0))
# clean up the output format
tab['MeanMag'].format = "{:.3f}"
tab['MeanCorrMag'].format = "{:.3f}"
tab['MAD'].format = "{:.4f}"
tab['Chi2'].format = "{:.4f}"
tab['RA'].format = "{:.6f}"
tab['Dec'].format = "{:.6f}"
# show some of the variable sources
tab[tab['AutoClass']>0]
Several of the table columns have information on the variability.
AutoClass
and ExpertClass
have summary information on the variability for a given MatchID
object.AutoClass
: Classification as provided by the system: 0=constant 1=single filter variable candidate (SFVC) 2=multi-filter variable candidate (MFVC)ExpertClass
: Classification as provided by expert: 0=not classified by expert, 1=high confidence variable, 2=probable variable, 4=possible artifactMAD
and Chi2
are variability indices using the median absolute deviation and the $\chi^2$ parameter for the given filter.VarQualFlag
is a variability quality flag (see Section 5 of the paper). The five letters correspond to CI, D, MagerrAper2, MagAper2-MagAuto, p2p; AAAAA corresponds to the highest quality flag.FilterDetFlag
is the filter detection flag: 1=source is variable in this filter, 0=source is not variable in this filter.See the HCV paper by Bonanos et al. (2019, AAp) for more details on the computation and meaning of these quantities.
This could be done in a SQL query through the CasJobs interface. Here we use the Astropy.table.join
function instead.
w475 = np.where(tab['Filter']=='ACS_F475W')
w814 = np.where(tab['Filter']=='ACS_F814W')
# the only key needed to do the join is MatchID, but we include other common columns
# so that join includes only one copy of them
jtab = join(tab[w475],tab[w814],
keys=['MatchID','GroupID','SubGroupID','RA','Dec','AutoClass','ExpertClass'],
table_names=['f475','f814'])
print(len(jtab),"matched F475W+F814W objects")
jtab[jtab['AutoClass']>0]
We mark the galaxy center as well. Note that this field is in the outskirts of IC 1613. The 0.5 degree search radius (which is the maximum allowed in the API) allows finding these objects.
pylab.rcParams.update({'font.size': 16})
pylab.figure(1,(10,10))
pylab.plot(jtab['RA'], jtab['Dec'], 'bo', markersize=1,
label='{:,} HCV measurements'.format(len(tab)))
pylab.plot(ra,dec,'rx',label=target,markersize=10)
pylab.gca().invert_xaxis()
pylab.gca().set_aspect('equal')
pylab.xlabel('RA [deg]')
pylab.ylabel('Dec [deg]')
pylab.legend(loc='best')
The median absolute deviation variability index is used by the HCV to identify variables. It measures the scatter among the multi-epoch measurements. Some scatter is expected from noise (which increases for fainter objects). Objects with MAD values that are high are likely to be variable.
This plots single-filter and multi-filter variable candidates (SFVC and MFVC) in different colors. Note that variable objects with low F475W MAD values are variable in a different filter (typically F814W in this field).
This plot is similar to the upper panel of Figure 4 in Bonanos et al. (2019, AAp).
wnovar = np.where(jtab['AutoClass']==0)[0]
wsfvc = np.where(jtab['AutoClass']==1)[0]
wmfvc = np.where(jtab['AutoClass']==2)[0]
x = jtab['MeanCorrMag_f475']
y = jtab['MAD_f475']
pylab.rcParams.update({'font.size': 16})
pylab.figure(1,(15,10))
pylab.plot(x[wnovar], y[wnovar], 'x', markersize=4, color='silver',
label='{:,} non-variable'.format(len(wnovar)))
pylab.plot(x[wsfvc], y[wsfvc], 'o', markersize=5, color='blue',
label='{:,} single-filter variable candidates'.format(len(wsfvc)))
pylab.plot(x[wmfvc], y[wmfvc], 'o', markersize=5, color='tab:cyan',
label='{:,} multi-filter variable candidates'.format(len(wmfvc)))
pylab.xlabel('ACS_F475W [mag]')
pylab.ylabel('ACS_F475W MAD [mag]')
pylab.legend(loc='best', title='{} HSC measurements near {}'.format(len(jtab),target))
Many of the candidate variables lie on the instability strip.
This plot is similar to the lower panel of Figure 4 in Bonanos et al. (2019, AAp).
wnovar = np.where(jtab['AutoClass']==0)[0]
wsfvc = np.where(jtab['AutoClass']==1)[0]
wmfvc = np.where(jtab['AutoClass']==2)[0]
x = jtab['MeanCorrMag_f475'] - jtab['MeanCorrMag_f814']
y = jtab['MeanCorrMag_f475']
pylab.rcParams.update({'font.size': 16})
pylab.figure(1,(15,10))
pylab.plot(x[wnovar], y[wnovar], 'x', markersize=4, color='silver',
label='{:,} non-variable'.format(len(wnovar)))
pylab.plot(x[wsfvc], y[wsfvc], 'o', markersize=5, color='blue',
label='{:,} single-filter variable candidates'.format(len(wsfvc)))
pylab.plot(x[wmfvc], y[wmfvc], 'o', markersize=5, color='tab:cyan',
label='{:,} multi-filter variable candidates'.format(len(wmfvc)))
pylab.gca().invert_yaxis()
pylab.xlim(-1.1, 4)
pylab.xlabel('ACS F475W - F814W [mag]')
pylab.ylabel('ACS F475W [mag]')
pylab.legend(loc='best', title='{} HSC measurements near {}'.format(len(jtab),target))
Note that the MatchID
could be determined by positional searches, filtering the catalog, etc. This object comes from the top left panel of Figure 9 in Bonanos et al. (2019, AAp).
matchid = 1905457
t0 = time.time()
# get light curves for F606W and F814W
nova_606 = ascii.read(hcvsearch(table='hcv',MatchID=matchid,Filter='ACS_F606W'))
print("{:.1f} sec: retrieved {} F606W measurements".format(time.time()-t0,len(nova_606)))
nova_814 = ascii.read(hcvsearch(table='hcv',MatchID=matchid,Filter='ACS_F814W'))
print("{:.1f} sec: retrieved {} F814W measurements".format(time.time()-t0,len(nova_814)))
# get the object RA and Dec as well
nova_tab = ascii.read(hcvsearch(table='hcvsummary',MatchID=matchid,Filter='ACS_F814W'))
print("{:.1f} sec: retrieved object info".format(time.time()-t0))
nova_606
pylab.rcParams.update({'font.size': 16})
pylab.figure(1,(15,10))
x = nova_606['MJD']
y = nova_606['CorrMag']
e = nova_606['MagErr']
pylab.errorbar(x,y,yerr=e,fmt='ob',ecolor='k',elinewidth=1,markersize=8,label='ACS F606W')
x = nova_814['MJD']
y = nova_814['CorrMag']
e = nova_814['MagErr']
pylab.errorbar(x,y,yerr=e,fmt='or',ecolor='k',elinewidth=1,markersize=8,label='ACS F814W')
pylab.gca().invert_yaxis()
pylab.xlabel('MJD [days]')
pylab.ylabel('magnitude')
pylab.legend(loc='best', title='Nova in M87 (MatchID: {})'.format(matchid))
The Hubble Legacy Archive (HLA) images were the source of the measurements in the HSC and HCV, and it can be useful to look at the images. Examination of the images can be useful to identified cosmic-ray contamination and other possible image artifacts. In this case, no issues are seen, so the light curve is reliable.
Note that the ACS F606W images of M87 have only a single exposure, so they do have cosmic ray contamination. The accompanying F814W images have multiple exposures, allowing CRs to be removed. In this case the F814W combined image is used to find objects, while the F606W exposure is used only for photometry. That reduces the effects of F606W CRs on the catalog but it is still a good idea to confirm the quality of the images.
The get_hla_cutout
function reads a single cutout image (as a JPEG grayscale image) and returns a PIL image object. See the documentation on the fitscut image cutout service for more information on the web service being used.
def get_hla_cutout(imagename,ra,dec,size=33,autoscale=99.5,asinh=True,zoom=1):
"""Get JPEG cutout for an image"""
url = "https://hla.stsci.edu/cgi-bin/fitscut.cgi"
r = requests.get(url, params=dict(ra=ra, dec=dec, size=size,
format="jpeg", red=imagename, autoscale=autoscale, asinh=asinh, zoom=zoom))
im = Image.open(BytesIO(r.content))
return im
# sort images by magnitude from brightest to faintest
phot = nova_606
isort = np.argsort(phot['CorrMag'])
# select the brightest, median and faintest magnitudes
ind = [isort[0], isort[len(isort)//2], isort[-1]]
# we plot zoomed-in and zoomed-out views side-by-side for each selected image
nim = len(ind)*2
ncols = 2 # images per row
nrows = (nim+ncols-1)//ncols
imsize1 = 19
imsize2 = 101
mra = nova_tab['RA'][0]
mdec = nova_tab['Dec'][0]
pylab.rcParams.update({"font.size":16})
pylab.figure(1,(12, (12/ncols)*nrows))
t0 = time.time()
ip = 0
for k in ind:
im1 = get_hla_cutout(phot['ImageName'][k],mra,mdec,size=imsize1)
ip += 1
pylab.subplot(nrows,ncols,ip)
pylab.imshow(im1,origin="upper",cmap="gray")
pylab.title('{} m={:.3f}'.format(phot['ImageName'][k],phot['CorrMag'][k]),fontsize=14)
im2 = get_hla_cutout(phot['ImageName'][k],mra,mdec,size=imsize2)
ip += 1
pylab.subplot(nrows,ncols,ip)
pylab.imshow(im2,origin="upper",cmap="gray")
xbox = np.array([-1,1])*imsize1/2 + (imsize2-1)//2
pylab.plot(xbox[[0,1,1,0,0]],xbox[[0,0,1,1,0]],'r-',linewidth=1)
pylab.title('{} m={:.3f}'.format(phot['ImageName'][k],phot['CorrMag'][k]),fontsize=14)
pylab.tight_layout()
print("{:.1f} s: got {} cutouts".format(time.time()-t0,ip))
The HCV includes an automatic classification AutoClass
for candidate variables as well as an expert validation for some fields that were selected for visual examination. For this example, we select all the objects in the HCV that have expert classification information.
t0 = time.time()
#get data for objects with an expert validation
constraints = {"ExpertClass.gte": 1}
tab = ascii.read(hcvsearch(table="hcvsummary",**constraints))
print("Retrieved {} rows in {:.1f} sec".format(len(tab),time.time()-t0))
# clean up the output format
tab['MeanMag'].format = "{:.3f}"
tab['MeanCorrMag'].format = "{:.3f}"
tab['MAD'].format = "{:.4f}"
tab['Chi2'].format = "{:.4f}"
tab['RA'].format = "{:.6f}"
tab['Dec'].format = "{:.6f}"
# tab includes 1 row for each filter (so multiple rows for objects with multiple filters)
# get an array that has only one row per object
mval, uindex = np.unique(tab['MatchID'],return_index=True)
utab = tab[uindex]
print("{} unique MatchIDs in table".format(len(utab)))
tab
An ExpertClass
value of 1 indicates that the object is confidently confirmed to be a variable; 2 means that the measurements do not have apparent problems and so the object is likely to be variable (usually the variability is too small to be obvious in the image); 4 means that the variability is likely to be the result of artifacts in the image (e.g., residual cosmic rays or diffraction spikes from nearby bright stars).
Compare the distributions for single-filter variable candidates (SFVC, AutoClass
=1) and multi-filter variable candidates (MFVC, AutoClass
=2). The fraction of artifacts is lower in the MFVC sample.
sfcount = np.bincount(utab['ExpertClass'][utab['AutoClass']==1])
mfcount = np.bincount(utab['ExpertClass'][utab['AutoClass']==2])
sfrat = sfcount/sfcount.sum()
mfrat = mfcount/mfcount.sum()
print("Type Variable Likely Artifact Total")
print("SFVC {:8d} {:6d} {:8d} {:5d} counts".format(sfcount[1],sfcount[2],sfcount[4],sfcount.sum()))
print("MFVC {:8d} {:6d} {:8d} {:5d} counts".format(mfcount[1],mfcount[2],mfcount[4],mfcount.sum()))
print("SFVC {:8.3f} {:6.3f} {:8.3f} {:5.3f} fraction".format(sfrat[1],sfrat[2],sfrat[4],sfrat.sum()))
print("MFVC {:8.3f} {:6.3f} {:8.3f} {:5.3f} fraction".format(mfrat[1],mfrat[2],mfrat[4],mfrat.sum()))
Note that only the filters identified as variable (FilterDetFlag
> 0) are included here.
This version of the plot shows the distributions for the various ExpertClass
values along with, for comparison, the distribution for all objects in gray (which is identical in each panel). Most objects are classified as confident or likely variables. Objects with lower MAD values (indicating a lower amplitude of variability) are less likely to be identified as confident variables because low-level variability is more difficult to confirm via visual examination.
w = np.where(tab['FilterDetFlag']>0)[0]
mad = tab['MAD'][w]
e = tab['ExpertClass'][w]
xrange = [7.e-3, 2.0]
bins = xrange[0]*(xrange[1]/xrange[0])**np.linspace(0.0,1.0,50)
pylab.rcParams.update({'font.size':16})
pylab.figure(1,(12,12))
pylab.subplot(311)
pylab.hist(mad,bins=bins,log=True,color='lightgray',label='All')
wp = np.where(e==1)[0]
pylab.hist(mad[wp],bins=bins,log=True,label='Confident',color='C2')
pylab.xscale('log')
pylab.ylabel('Count')
pylab.legend(loc='upper left')
pylab.title('HCV Expert Validation')
pylab.subplot(312)
pylab.hist(mad,bins=bins,log=True,color='lightgray',label='All')
wp = np.where(e==2)[0]
pylab.hist(mad[wp],bins=bins,log=True,label='Likely',color='C1')
pylab.xscale('log')
pylab.ylabel('Count')
pylab.legend(loc='upper left')
pylab.subplot(313)
pylab.hist(mad,bins=bins,log=True,color='lightgray',label='All')
wp = np.where(e==4)[0]
pylab.hist(mad[wp],bins=bins,log=True,label='Artifact',color='C0')
pylab.xscale('log')
pylab.ylabel('Count')
pylab.legend(loc='upper left')
pylab.xlabel('MAD Variability Index [mag]')
The plot below shows the same distributions, but plotted as stacked histograms. The top panel uses a linear scale on the y-axis and the bottom panel uses a log y scale.
w = np.where(tab['FilterDetFlag']>0)[0]
mad = tab['MAD'][w]
e = tab['ExpertClass'][w]
xrange = [7.e-3, 2.0]
bins = xrange[0]*(xrange[1]/xrange[0])**np.linspace(0.0,1.0,50)
pylab.rcParams.update({'font.size':16})
pylab.figure(1,(15,12))
pylab.subplot(211)
hlog = False
pylab.hist(mad,bins=bins,log=hlog,label='Artifact')
wp = np.where(e<4)[0]
pylab.hist(mad[wp],bins=bins,log=hlog,label='Likely Variable')
wp = np.where(e==1)[0]
pylab.hist(mad[wp],bins=bins,log=hlog,label='Confident Variable')
pylab.xscale('log')
pylab.xlabel('MAD Variability Index [mag]')
pylab.ylabel('Count')
pylab.legend(loc='upper right',title='HCV Expert Validation')
pylab.subplot(212)
hlog = True
pylab.hist(mad,bins=bins,log=hlog,label='Artifact')
wp = np.where(e<4)[0]
pylab.hist(mad[wp],bins=bins,log=hlog,label='Likely Variable')
wp = np.where(e==1)[0]
pylab.hist(mad[wp],bins=bins,log=hlog,label='Confident Variable')
pylab.xscale('log')
pylab.xlabel('MAD Variability Index [mag]')
pylab.ylabel('Count')
pylab.legend(loc='upper right',title='HCV Expert Validation')