Query Pan-STARRS DR2 catalog using CasJobs¶
This script shows how to query the Pan-STARRS DR2 catalog using a Python interface to Casjobs. The examples show how to do a simple cone search, how to manipulate the table of results, and how to get a light curve from the table of detections.
This relies on the mastcasjobs Python module. Follow the installation instructions given here: https://github.com/rlwastro/mastcasjobs.
You must have a MAST Casjobs account (see https://mastweb.stsci.edu/ps1casjobs to create one). Note that MAST Casjobs accounts are independent of SDSS Casjobs accounts.
For easy startup, set the CASJOBS_USERID and CASJOBS_PW environment variables with your Casjobs account information. These are the username and password you enter when logging into Casjobs.
This script prompts for your Casjobs username and password if the environment variables are not defined.
You can also specify your username and password directly in the MastCasJobs initialization using the username and password keyword parameters.
This notebook is available for download.
%matplotlib inline
import mastcasjobs
import requests
import os
import numpy as np
import matplotlib.pyplot as plt
import json
# get the WSID and password if not already defined
import getpass
if not os.environ.get('CASJOBS_USERID'):
os.environ['CASJOBS_USERID'] = input('Enter Casjobs username:')
if not os.environ.get('CASJOBS_PW'):
os.environ['CASJOBS_PW'] = getpass.getpass('Enter Casjobs password:')
Useful functions¶
Resolve names using MAST query. This uses the requests module to do the queries.
def mastQuery(request, json_return=False):
"""Perform a MAST query.
Parameters
----------
request (dictionary): The MAST request json object
Returns the text response or (if json_return=True) the json response
"""
url = "https://mast.stsci.edu/api/v0/invoke"
# Encoding the request as a json string
requestString = json.dumps(request)
# make the query
r = requests.post(url, data=dict(request=requestString))
# raise exception on error
r.raise_for_status()
if json_return:
return r.json()
else:
return r.text
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'
},
}
resolvedObject = mastQuery(resolverRequest, json_return=True)
# 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)
Simple positional query¶
This searches the mean object catalog for objects within 50 arcsec of M87 (RA=187.706, Dec=12.391 in degrees). The fGetNearbyObjEq function returns a list of object IDs, and the subsequent joins extract information from the ObjectThin table (which has information on object positions and the number of available measurements) and the MeanObject table (which has information on photometry averaged over the multiple epochs of observation).
Note that the results are restricted to objects with nDetections>1, where nDetections is the total number of times the object was detected on the single-epoch images in any filter at any time. Objects with nDetections=1 tend to be artifacts, so this is a quick way to eliminate most spurious objects from the catalog.
This query runs in the interactive "quick" Casjobs queue, where the function pauses until the query is complete and returns the results as a CSV string. If one of these queries times out (which can happen if the database server is heavily loaded), try simply running it again. Often the same query will succeed the second time.
The more robust way to handle long-running queries is to run them in batch mode, but that requires waiting in the Casjobs queue if the system is busy.
The quick query returns the results as an astropy table. This table is easily manipulated to extract information on individual columns or rows.
query = """select o.objID, o.raMean, o.decMean,
o.nDetections, o.ng, o.nr, o.ni, o.nz, o.ny,
m.gMeanPSFMag, m.rMeanPSFMag, m.iMeanPSFMag, m.zMeanPSFMag, m.yMeanPSFMag
from fGetNearbyObjEq(187.706,12.391,50.0/60.0) nb
inner join ObjectThin o on o.objid=nb.objid and o.nDetections>1
inner join MeanObject m on o.objid=m.objid and o.uniquePspsOBid=m.uniquePspsOBid
"""
jobs = mastcasjobs.MastCasJobs(context="PanSTARRS_DR2")
tab = jobs.quick(query, task_name="python cone search")
tab
Get DR2 light curve for RR Lyrae star KQ UMa¶
This time we start with the object name, use the MAST name resolver (which relies on Simbad and NED) to convert the name to RA and Dec, and then query the PS1 DR2 mean object catalog at that position. This time we use the fGetNearestObjEq function to select just the object closest to the search position.
objname = 'KQ UMa'
ra, dec = resolve(objname)
radius = 1.0/60.0 # radius = 1 arcsec
query = """select o.objID, o.raMean, o.decMean,
o.nDetections, o.ng, o.nr, o.ni, o.nz, o.ny,
m.gMeanPSFMag, m.rMeanPSFMag, m.iMeanPSFMag, m.zMeanPSFMag, m.yMeanPSFMag
from fGetNearestObjEq({},{},{}) nb
inner join ObjectThin o on o.objid=nb.objid and o.nDetections>1
inner join MeanObject m on o.objid=m.objid and o.uniquePspsOBid=m.uniquePspsOBid
""".format(ra,dec,radius)
print(query)
jobs = mastcasjobs.MastCasJobs(context="PanSTARRS_DR2")
tab = jobs.quick(query, task_name="python cone search")
tab
Get the detection information¶
Extract all the objects with the same object ID from the Detection table, which contains all the individual measurements for this source. The results are joined to the Filter table to convert the filter numbers to names. The somewhat odd structure for the SQL (with the inner query in parentheses) ensures that the select by objID occurs before the match to the Filter table. If the SQL optimizer gets confused and does the join between Detection and Filter before selecting the subset of objects, the query is very slow!
objid = tab['objID'][0]
query = """select
objID, detectID, filter=f.filterType, obsTime, ra, dec,
psfFlux, psfFluxErr, psfMajorFWHM, psfMinorFWHM, psfQfPerfect,
apFlux, apFluxErr, infoFlag, infoFlag2, infoFlag3
from (
select * from Detection where objID={}
) d
join Filter f on d.filterID=f.filterID
order by d.filterID, obsTime
""".format(objid)
print(query)
dtab = jobs.quick(query, task_name="python detection search")
dtab
Plot the light curves¶
The psfFlux values from the Detection table are converted from Janskys to AB magnitudes. Measurements in the 5 different filters are plotted separately.
# convert flux in Jy to magnitudes
t = dtab['obsTime']
mag = -2.5*np.log10(dtab['psfFlux']) + 8.90
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,10))
for i, filter in enumerate("grizy"):
plt.subplot(511+i)
w = np.where(dtab['filter']==filter)
plt.plot(t[w],mag[w],'-o')
plt.ylabel(filter+' [mag]')
plt.xlim(xlim)
plt.gca().invert_yaxis()
if i==0:
plt.title(objname)
plt.xlabel('Time [MJD]')
plt.tight_layout()
Plot differences from the mean magnitudes in the initial search.
# convert flux in Jy to magnitudes
t = dtab['obsTime']
mag = -2.5*np.log10(dtab['psfFlux']) + 8.90
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,10))
for i, filter in enumerate("grizy"):
plt.subplot(511+i)
w = np.where(dtab['filter']==filter)
magmean = tab[filter+'MeanPSFMag'][0]
plt.plot(t[w],mag[w] - magmean,'-o')
plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))
plt.xlim(xlim)
plt.gca().invert_yaxis()
if i==0:
plt.title(objname)
plt.xlabel('Time [MJD]')
plt.tight_layout()
Identify bad data¶
There is one clearly bad $z$ magnitude with a very large difference. Select the bad point and look at it in more detail.
Note that indexing a table (or numpy array) with a logical expression selects just the rows where that expression is true.
dtab[ (dtab['filter']=='z') & (np.abs(mag-tab['zMeanPSFMag'][0]) > 2) ]
From examining this table, it looks like psfQfPerfect is bad. This flag is the PSF-weighted fraction of unmasked pixels in the image (see the documentation for more details). Values near unity indicate good data that is not significantly affected by bad pixels.
Check all the psfQfPerfect values for the $z$ filter to see if this value really is unusual. The list of values below are sorted by magnitude. The bad point is the only value with psfQfPerfect < 0.95.
w = np.where(dtab['filter']=='z')
zdtab = dtab[w]
zdtab['mag'] = mag[w]
zdtab['dmag'] = zdtab['mag'] - tab['zMeanPSFMag'][0]
ii = np.argsort(-np.abs(zdtab['dmag']))
zdtab = zdtab[ii]
zdtab['objID','obsTime','mag','dmag','psfQfPerfect']
Repeat the plot with bad psfQfPerfect values excluded¶
Do the plot again but exclude low psfQfPerfect values.
# convert flux in Jy to magnitudes
t = dtab['obsTime']
mag = -2.5*np.log10(dtab['psfFlux']) + 8.90
magmean = 0.0*mag
for filter in "grizy":
magmean[dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]
dmag = mag - magmean
dmag1 = dmag[dtab['psfQfPerfect']>0.9]
# fix the x and y axis ranges
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
# flip axis direction for magnitude
ylim = np.array([dmag1.max(),dmag1.min()])
ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,10))
for i, filter in enumerate("grizy"):
plt.subplot(511+i)
w = np.where((dtab['filter']==filter) & (dtab['psfQfPerfect']>0.9))[0]
plt.plot(t[w],dmag[w],'-o')
plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))
plt.xlim(xlim)
plt.ylim(ylim)
if i==0:
plt.title(objname)
plt.xlabel('Time [MJD]')
plt.tight_layout()
Plot versus the periodic phase instead of epoch¶
Plot versus phase using known RR Lyr period from Simbad (table J/AJ/132/1202/table4).
period = 0.48636
# convert flux in Jy to magnitudes
t = (dtab['obsTime'] % period) / period
mag = -2.5*np.log10(dtab['psfFlux']) + 8.90
magmean = 0.0*mag
for filter in "grizy":
magmean[dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]
dmag = mag - magmean
dmag1 = dmag[dtab['psfQfPerfect']>0.9]
# fix the x and y axis ranges
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
# flip axis direction for magnitude
ylim = np.array([dmag1.max(),dmag1.min()])
ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,10))
for i, filter in enumerate("grizy"):
plt.subplot(511+i)
w = np.where((dtab['filter']==filter) & (dtab['psfQfPerfect']>0.9))[0]
w = w[np.argsort(t[w])]
plt.plot(t[w],dmag[w],'-o')
plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))
plt.xlim(xlim)
plt.ylim(ylim)
if i==0:
plt.title(objname)
plt.xlabel('Phase')
plt.tight_layout()
Repeat search using eclipsing binary KIC 2161623¶
objname = 'KIC 2161623'
ra, dec = resolve(objname)
radius = 1.0/60.0 # radius = 1 arcsec
query = """select o.objID, o.raMean, o.decMean,
o.nDetections, o.ng, o.nr, o.ni, o.nz, o.ny,
m.gMeanPSFMag, m.rMeanPSFMag, m.iMeanPSFMag, m.zMeanPSFMag, m.yMeanPSFMag
from fGetNearestObjEq({},{},{}) nb
inner join ObjectThin o on o.objid=nb.objid and o.nDetections>1
inner join MeanObject m on o.objid=m.objid and o.uniquePspsOBid=m.uniquePspsOBid
""".format(ra,dec,radius)
print(query)
jobs = mastcasjobs.MastCasJobs(context="PanSTARRS_DR2")
tab = jobs.quick(query, task_name="python cone search")
tab
Get the detection information¶
This time include the psfQfPerfect limit directly in the database query.
objid = tab['objID'][0]
query = """select
objID, detectID, filter=f.filterType, obsTime, ra, dec,
psfFlux, psfFluxErr, psfMajorFWHM, psfMinorFWHM, psfQfPerfect,
apFlux, apFluxErr, infoFlag, infoFlag2, infoFlag3
from (
select * from Detection where objID={}
) d
join Filter f on d.filterID=f.filterID
where psfQfPerfect>0.9
order by d.filterID, obsTime
""".format(objid)
print(query)
dtab = jobs.quick(query, task_name="python detection search")
# add magnitude and difference from mean
dtab['magmean'] = 0.0
for filter in "grizy":
dtab['magmean'][dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]
dtab['mag'] = -2.5*np.log10(dtab['psfFlux']) + 8.90
dtab['dmag'] = dtab['mag']-dtab['magmean']
dtab
t = dtab['obsTime']
dmag = dtab['dmag']
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
ylim = np.array([dmag.max(),dmag.min()])
ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,10))
for i, filter in enumerate("grizy"):
plt.subplot(511+i)
w = np.where(dtab['filter']==filter)[0]
plt.plot(t[w],dmag[w],'-o')
magmean = dtab['magmean'][w[0]]
plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))
plt.xlim(xlim)
plt.ylim(ylim)
if i==0:
plt.title(objname)
plt.xlabel('Time [MJD]')
plt.tight_layout()
Plot versus phase using known period¶
Eclipsing binaries basically vary by same amount in all filters since it is a geometrical effect, so combine the data into a single light curve. Wrap using known period and plot versus phase.
period = 2.2834698
bjd0 = 54999.599837
t = ((dtab['obsTime']-bjd0) % period) / period
dmag = dtab['dmag']
w = np.argsort(t)
t = t[w]
dmag = dmag[w]
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
ylim = np.array([dmag.max(),dmag.min()])
ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,6))
plt.plot(t,dmag,'-o')
plt.xlim(xlim)
plt.ylim(ylim)
plt.xlabel('Phase')
plt.ylabel('Delta magnitude from mean [mag]')
plt.title(objname)
plt.tight_layout()
Repeat search for another eclipsing binary KIC 8153568¶
objname = 'KIC 8153568'
ra, dec = resolve(objname)
radius = 1.0/60.0 # radius = 1 arcsec
query = """select o.objID, o.raMean, o.decMean,
o.nDetections, o.ng, o.nr, o.ni, o.nz, o.ny,
m.gMeanPSFMag, m.rMeanPSFMag, m.iMeanPSFMag, m.zMeanPSFMag, m.yMeanPSFMag
from fGetNearestObjEq({},{},{}) nb
inner join ObjectThin o on o.objid=nb.objid and o.nDetections>1
inner join MeanObject m on o.objid=m.objid and o.uniquePspsOBid=m.uniquePspsOBid
""".format(ra,dec,radius)
jobs = mastcasjobs.MastCasJobs(context="PanSTARRS_DR2")
tab = jobs.quick(query, task_name="python cone search")
tab
objid = tab['objID'][0]
query = """select
objID, detectID, filter=f.filterType, obsTime, ra, dec,
psfFlux, psfFluxErr, psfMajorFWHM, psfMinorFWHM, psfQfPerfect,
apFlux, apFluxErr, infoFlag, infoFlag2, infoFlag3
from (
select * from Detection where objID={}
) d
join Filter f on d.filterID=f.filterID
where psfQfPerfect>0.9
order by d.filterID, obsTime
""".format(objid)
dtab = jobs.quick(query, task_name="python detection search")
# add magnitude and difference from mean
dtab['magmean'] = 0.0
for filter in "grizy":
dtab['magmean'][dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]
dtab['mag'] = -2.5*np.log10(dtab['psfFlux']) + 8.90
dtab['dmag'] = dtab['mag']-dtab['magmean']
dtab
t = dtab['obsTime']
dmag = dtab['dmag']
xlim = np.array([t.min(),t.max()])
xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])
ylim = np.array([dmag.max(),dmag.min()])
ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(10,10))
for i, filter in enumerate("grizy"):
plt.subplot(511+i)
w = np.where(dtab['filter']==filter)[0]
plt.plot(t[w],dmag[w],'-o')
magmean = dtab['magmean'][w[0]]
plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))
plt.xlim(xlim)
plt.ylim(ylim)
if i==0:
plt.title(objname)
plt.xlabel('Time [MJD]')
plt.tight_layout()
Eclipsing binaries basically vary by same amount in all filters since it is a geometrical effect, so combine the data into a single light curve.
Wrap using known period and plot versus phase. Plot two periods of the light curve this time.
This nice light curve appears to show a secondary eclipse.
period = 3.6071431
bjd0 = 54999.289794
t = ((dtab['obsTime']-bjd0) % period) / period
dmag = dtab['dmag']
w = np.argsort(t)
# extend to two periods
nw = len(w)
w = np.append(w,w)
t = t[w]
# add one to second period
t[-nw:] += 1
dmag = dmag[w]
xlim = [0,2.0]
ylim = np.array([dmag.max(),dmag.min()])
ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])
plt.rcParams.update({'font.size': 14})
plt.figure(1,(12,6))
plt.plot(t,dmag,'-o')
plt.xlim(xlim)
plt.ylim(ylim)
plt.xlabel('Phase')
plt.ylabel('Delta magnitude from mean [mag]')
plt.title(objname)
plt.tight_layout()
