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_WSID and CASJOBS_PW environment variables with your Casjobs account information. You can get your WSID by going to https://mastweb.stsci.edu/ps1casjobs/changedetails.aspx after you login to Casjobs. Your password is what you enter when logging into Casjobs.

This script prompts for your Casjobs WSID and password if the environment variables are not defined.

You can also specify your wsid and password directly in the MastCasJobs initialization using the userid and password keyword parameters.

This notebook is available for download.

In [1]:
%matplotlib inline
import mastcasjobs
from astropy.io import ascii
from astropy.table import Table

import sys
import os
import re
import numpy as np
import pylab
import json

try: # Python 3.x
    from urllib.parse import quote as urlencode
    from urllib.request import urlretrieve
except ImportError:  # Python 2.x
    from urllib import pathname2url as urlencode
    from urllib import urlretrieve

try: # Python 3.x
    import http.client as httplib 
except ImportError:  # Python 2.x
    import httplib

# get the WSID and password if not already defined
import getpass
if not os.environ.get('CASJOBS_WSID'):
    os.environ['CASJOBS_WSID'] = input('Enter Casjobs WSID:')
if not os.environ.get('CASJOBS_PW'):
    os.environ['CASJOBS_PW'] = getpass.getpass('Enter Casjobs password:')

Useful functions

Resolve names using MAST query and fix up column names returned from Casjobs.

In [2]:
def mastQuery(request):
    """Perform a MAST query.

    Parameters
    ----------
    request (dictionary): The MAST request json object

    Returns head,content where head is the response HTTP headers, and content is the returned data"""
    
    server='mast.stsci.edu'

    # Grab Python Version 
    version = ".".join(map(str, sys.version_info[:3]))

    # Create Http Header Variables
    headers = {"Content-type": "application/x-www-form-urlencoded",
               "Accept": "text/plain",
               "User-agent":"python-requests/"+version}

    # Encoding the request as a json string
    requestString = json.dumps(request)
    requestString = urlencode(requestString)
    
    # opening the https connection
    conn = httplib.HTTPSConnection(server)

    # Making the query
    conn.request("POST", "/api/v0/invoke", "request="+requestString, headers)

    # Getting the response
    resp = conn.getresponse()
    head = resp.getheaders()
    content = resp.read().decode('utf-8')

    # Close the https connection
    conn.close()

    return head,content


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'
                                },
                      }
    headers,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 fixcolnames(tab):
    """Fix column names returned by the casjobs query
    
    Parameters
    ----------
    tab (astropy.table.Table): Input table

    Returns reference to original table with column names modified"""

    pat = re.compile(r'\[(?P<name>[^[]+)\]')
    for c in tab.colnames:
        m = pat.match(c)
        if not m:
            raise ValueError("Unable to parse column name '{}'".format(c))
        newname = m.group('name')
        tab.rename_column(c,newname)
    return tab

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.

In [3]:
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")
results = jobs.quick(query, task_name="python cone search")
print(results)
[objID]:Integer,[raMean]:Float,[decMean]:Float,[nDetections]:Integer,[ng]:Integer,[nr]:Integer,[ni]:Integer,[nz]:Integer,[ny]:Integer,[gMeanPSFMag]:Float,[rMeanPSFMag]:Float,[iMeanPSFMag]:Float,[zMeanPSFMag]:Float,[yMeanPSFMag]:Float
122861877070368181,187.70700504,12.38964881,2,0,2,0,0,0,-999,15.2261,-999,-999,-999
122861877090664549,187.70906689,12.38671491,2,0,0,0,0,2,-999,-999,-999,-999,-999
122861877111116495,187.71111183,12.38833651,2,1,1,0,0,0,-999,-999,-999,-999,-999
122861877112026930,187.71120247,12.38869951,2,0,1,1,0,0,-999,-999,-999,-999,-999
122861877184282370,187.7182975,12.38498455,2,0,0,2,0,0,-999,-999,18.1435,-999,-999
122871876925340070,187.69247065,12.39128192,6,2,1,2,1,0,19.3276,19.0171,19.5632,19.2438,-999
122871876941435711,187.6940799,12.3959326,4,0,1,2,1,0,-999,-999,19.0092,19.1839,-999
122871876945758919,187.69449749,12.39870483,2,0,1,1,0,0,-999,20.8721,19.6038,-999,-999
122871876946391889,187.6946297,12.39284071,2,0,2,0,0,0,-999,17.6493,-999,-999,-999
122871876951258032,187.69498862,12.39809102,2,1,0,0,0,1,-999,-999,-999,-999,17.2506
122871877012582071,187.70130031,12.39292937,3,2,0,0,0,1,17.3948,-999,-999,-999,-999
122871877026011277,187.70253784,12.39234062,62,13,12,16,10,11,16.9417,15.3886,15.3299,15.7681,16.9338
122871877062551436,187.70622365,12.39236505,2,0,0,0,0,2,-999,-999,-999,-999,15.1102
122871877063310741,187.70630472,12.39182315,3,0,2,0,1,0,-999,14.4229,-999,14.4249,-999
122871877068858435,187.70691585,12.39827062,7,3,1,2,1,0,18.6591,17.8651,17.7776,17.9392,-999
122871877075050430,187.70747912,12.3915839,2,0,0,0,2,0,-999,-999,-999,15.5752,-999
122881876951360506,187.69529229,12.40011853,2,0,0,1,0,1,-999,-999,19.2936,-999,17.7803
122881877085143466,187.7084996,12.40240408,37,9,11,7,7,3,18.8261,18.0068,17.7669,17.8751,17.6812
122881877095525205,187.70951703,12.40389466,2,0,0,2,0,0,-999,-999,18.8612,-999,-999
122881877105204765,187.71049731,12.4035168,2,0,0,2,0,0,-999,-999,19.1368,-999,-999
122881877112164157,187.71126262,12.40316813,2,0,0,2,0,0,-999,-999,19.0408,-999,-999
122851876947049923,187.69476028,12.382815,8,2,1,2,3,0,19.8674,18.5927,19.0692,18.1394,-999
122851877008084336,187.70080125,12.37828064,15,3,2,5,5,0,21.1629,-999,18.7543,18.4389,-999
122851877018527082,187.70183167,12.38053514,5,2,2,1,0,0,20.3146,-999,18.1687,-999,-999
122851877049294726,187.70495239,12.37851452,7,2,2,2,1,0,19.0306,18.4537,18.5761,-999,-999
122851877062673559,187.70624956,12.37754834,7,1,2,3,0,1,20.521,18.9063,18.9627,-999,18.1421
122851877090477761,187.70902824,12.38101615,29,6,7,11,5,0,19.7468,17.7461,17.4915,17.6221,-999
122851877096718451,187.70969564,12.38158628,2,1,0,1,0,0,20.5734,-999,-999,-999,-999
122851877109888891,187.71098809,12.38200012,2,0,0,2,0,0,-999,-999,-999,-999,-999
122851877115528846,187.71134341,12.38181909,2,1,1,0,0,0,23.3524,17.3302,-999,-999,-999
122851877124637737,187.71246077,12.38100609,50,9,10,15,9,7,19.3247,18.0296,17.7988,17.9304,17.5089
122861876929948166,187.69304884,12.38978617,7,2,0,3,2,0,19.7578,-999,18.5519,19.226,-999
122861876944428797,187.69448057,12.39022676,2,0,0,0,2,0,-999,-999,-999,20.1199,-999
122861876946372059,187.69472326,12.38464777,2,0,0,0,2,0,-999,-999,-999,18.8606,-999
122861876946987056,187.69442859,12.38863934,2,0,1,0,1,0,-999,17.1852,-999,18.3888,-999
122861876947177869,187.69469555,12.38944923,2,0,2,0,0,0,-999,18.3706,-999,-999,-999
122861876952960254,187.69534396,12.38300575,2,0,0,0,0,2,-999,-999,-999,-999,18.6671
122861876979515942,187.69789315,12.38779688,44,8,7,13,8,8,18.2713,17.4948,17.2884,17.4338,17.1059
122861877044629638,187.70451598,12.39086724,2,0,0,0,2,0,-999,-999,-999,15.5294,-999
122861877050678994,187.70512805,12.39043743,5,2,2,1,0,0,17.0287,15.0843,16.169,-999,-999
122861877056308967,187.70559948,12.39030211,2,0,0,0,0,2,-999,-999,-999,-999,14.9985
122861877056688054,187.70558996,12.38960002,2,0,0,0,2,0,-999,-999,-999,14.4064,-999
122861877058698594,187.70580919,12.39013311,2,0,0,0,0,2,-999,-999,-999,-999,14.7305
122861877059169881,187.70591961,12.39112604,54,10,10,14,8,12,10.295,10.732,11.076,11.342,11.557
122861877064767617,187.70646616,12.38914179,2,0,2,0,0,0,-999,15.1197,-999,-999,-999

Convert the results to an astropy table

The CSV results string is easily converted to an astropy table. The column names are messy (with embedded datatypes), so the fixcolnames function (defined above) is called to clean up the names. This table is easily manipulated to extract information on individual columns or rows.

In [4]:
tab = fixcolnames(ascii.read(results))
tab
Out[4]:
Table length=45
objIDraMeandecMeannDetectionsngnrninznygMeanPSFMagrMeanPSFMagiMeanPSFMagzMeanPSFMagyMeanPSFMag
int64float64float64int64int64int64int64int64int64float64float64float64float64float64
122861877070368181187.7070050412.38964881202000-999.015.2261-999.0-999.0-999.0
122861877090664549187.7090668912.38671491200002-999.0-999.0-999.0-999.0-999.0
122861877111116495187.7111118312.38833651211000-999.0-999.0-999.0-999.0-999.0
122861877112026930187.7112024712.38869951201100-999.0-999.0-999.0-999.0-999.0
122861877184282370187.718297512.38498455200200-999.0-999.018.1435-999.0-999.0
122871876925340070187.6924706512.3912819262121019.327619.017119.563219.2438-999.0
122871876941435711187.694079912.3959326401210-999.0-999.019.009219.1839-999.0
122871876945758919187.6944974912.39870483201100-999.020.872119.6038-999.0-999.0
122871876946391889187.694629712.39284071202000-999.017.6493-999.0-999.0-999.0
122871876951258032187.6949886212.39809102210001-999.0-999.0-999.0-999.017.2506
..........................................
122861876947177869187.6946955512.38944923202000-999.018.3706-999.0-999.0-999.0
122861876952960254187.6953439612.38300575200002-999.0-999.0-999.0-999.018.6671
122861876979515942187.6978931512.387796884487138818.271317.494817.288417.433817.1059
122861877044629638187.7045159812.39086724200020-999.0-999.0-999.015.5294-999.0
122861877050678994187.7051280512.3904374352210017.028715.084316.169-999.0-999.0
122861877056308967187.7055994812.39030211200002-999.0-999.0-999.0-999.014.9985
122861877056688054187.7055899612.38960002200020-999.0-999.0-999.014.4064-999.0
122861877058698594187.7058091912.39013311200002-999.0-999.0-999.0-999.014.7305
122861877059169881187.7059196112.391126045410101481210.29510.73211.07611.34211.557
122861877064767617187.7064661612.38914179202000-999.015.1197-999.0-999.0-999.0

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.

In [5]:
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")
results = jobs.quick(query, task_name="python cone search")
tab = fixcolnames(ascii.read(results))
tab
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(139.3344625554179,68.6350884268329,0.016666666666666666) 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

Out[5]:
Table length=1
objIDraMeandecMeannDetectionsngnrninznygMeanPSFMagrMeanPSFMagiMeanPSFMagzMeanPSFMagyMeanPSFMag
int64float64float64int64int64int64int64int64int64float64float64float64float64float64
190361393344112894139.3344497268.63506046681021131415.040214.55314.210914.281414.3041

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!

In [6]:
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)

dresults = jobs.quick(query, task_name="python detection search")
dtab = fixcolnames(ascii.read(dresults))
dtab
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=190361393344112894
    ) d
join Filter f on d.filterID=f.filterID
order by d.filterID, obsTime

Out[6]:
Table length=66
objIDdetectIDfilterobsTimeradecpsfFluxpsfFluxErrpsfMajorFWHMpsfMinorFWHMpsfQfPerfectapFluxapFluxErrinfoFlaginfoFlag2infoFlag3
int64int64str1float64float64float64float64float64float64float64float64float64float64int64int64int64
190361393344112894153347716310000010g55634.477414139.334520768.635035770.008261921.14074e-051.881851.767340.9929160.008617051.14233e-05102760517128124782656
190361393344112894153348968310000008g55634.4899457139.3344882168.635061460.00773731.10268e-051.810311.605180.9984610.007921721.09406e-05102760517128124782656
190361393344112894232228791560000017g56423.2881719139.3344944168.635039520.003351987.26745e-061.602041.450480.9985890.003386587.2363e-06102760517128108038208
190361393344112894255559866370000015g56656.5989211139.3344459268.635045890.003729097.37688e-061.828311.686920.9991910.003721637.52088e-06102760517128124815424
190361393344112894262040070370000016g56721.4009635139.3344514868.635048070.0035767.22894e-061.605851.503960.9986050.003596037.43463e-061027605171287374912
190361393344112894262040708370000014g56721.4073411139.3344530468.635044060.003505547.14809e-061.624151.488270.9994870.003596247.43434e-061027605171287374912
190361393344112894264231864260000022g56743.3189045139.334443168.635044970.003438677.81574e-061.539571.47610.9990490.003612397.85237e-061027605171287374912
190361393344112894264232516260000024g56743.3254198139.3344451268.635052850.003474037.87789e-061.527511.433250.9993670.003603497.8686e-061027605171287374912
190361393344112894153441340410000012r55635.4136426139.3344717368.63505340.009787581.10843e-051.678791.56530.9985090.009907631.11231e-05102760517128124815424
190361393344112894153442549410000013r55635.4257356139.3344773268.635052790.00936211.1005e-051.347241.209520.4074960.009506791.09949e-0510276051716034880
................................................
190361393344112894183254354520000014y55933.5437238139.334453468.635059040.00793942.4924e-050.9592310.9014580.9968530.008030782.37126e-051027605171287374912
190361393344112894183255163520000015y55933.5518128139.3344423268.635056690.009155732.64422e-051.113310.953680.9971150.009145152.48501e-0510276051712874483776
190361393344112894191726725200000011y56018.2674345139.3344607968.635064030.008296462.59461e-051.073321.015630.9979240.008337592.33303e-051027605171287374912
190361393344112894191727369200000012y56018.2738664139.3344518968.635061520.008199962.55524e-051.095041.059140.9978190.008207532.3153e-051027605171287374912
190361393344112894229122569150000014y56392.2258719139.3344274268.635062340.007804222.57689e-051.028780.8931950.9986180.007772792.16455e-051027605171287374912
190361393344112894229123273150000011y56392.2329084139.3344459568.635057440.008035032.55974e-050.9418010.8887880.998050.008070372.2092e-051027605171287374912
190361393344112894264322648150000024y56744.2269452139.33445468.635058640.005404461.32146e-051.388291.34770.9989220.005498671.13527e-05102760517128124782656
190361393344112894265623008150000027y56757.2305455139.3344431968.635063750.006366141.60784e-050.8950410.8033650.8525160.004128591.01303e-0510276058112834880
190361393344112894287165374630000023y56972.6539243139.3344422668.635060870.006878172.79326e-051.536861.288130.9916470.006999142.16585e-051027605171287374912
190361393344112894287265116630000028y56973.6513438139.3344499168.635057590.006653772.7704e-051.856191.703970.998640.006631062.13156e-051027605171287342144

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.

In [7]:
# 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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,10))
for i, filter in enumerate("grizy"):
    pylab.subplot(511+i)
    w = np.where(dtab['filter']==filter)
    pylab.plot(t[w],mag[w],'-o')
    pylab.ylabel(filter+' [mag]')
    pylab.xlim(xlim)
    pylab.gca().invert_yaxis()
    if i==0:
        pylab.title(objname)
pylab.xlabel('Time [MJD]')
pylab.tight_layout()

Plot differences from the mean magnitudes in the initial search.

In [8]:
# 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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,10))
for i, filter in enumerate("grizy"):
    pylab.subplot(511+i)
    w = np.where(dtab['filter']==filter)
    magmean = tab[filter+'MeanPSFMag'][0]
    pylab.plot(t[w],mag[w] - magmean,'-o')
    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))
    pylab.xlim(xlim)
    pylab.gca().invert_yaxis()
    if i==0:
        pylab.title(objname)
pylab.xlabel('Time [MJD]')
pylab.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.

In [9]:
dtab[ (dtab['filter']=='z') & (np.abs(mag-tab['zMeanPSFMag'][0]) > 2) ]
Out[9]:
Table length=1
objIDdetectIDfilterobsTimeradecpsfFluxpsfFluxErrpsfMajorFWHMpsfMinorFWHMpsfQfPerfectapFluxapFluxErrinfoFlaginfoFlag2infoFlag3
int64int64str1float64float64float64float64float64float64float64float64float64float64int64int64int64
190361393344112894183252627520000234z55933.5264577139.3348816868.635322730.0003179456.73008e-061.075371.01530.3229860.0002132172.36939e-0610276045312832768

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.

In [10]:
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']
Out[10]:
Table length=13
objIDobsTimemagdmagpsfQfPerfect
int64float64float64float64float64
19036139334411289455933.526457717.644120000843743.3627200008437390.322986
19036139334411289456289.615934613.890659532156441-0.390740467843558650.997811
19036139334411289456289.624111213.916806313581588-0.364593686418411170.988369
19036139334411289456351.416848313.998972905236004-0.282427094763995970.999257
19036139334411289455281.252828514.5379146386864730.25651463868647360.99754
19036139334411289456351.42407614.032501974717418-0.248898025282581870.999187
19036139334411289455527.650891914.512118206878710.230718206878711030.997265
19036139334411289456648.567601914.056429712068542-0.224970287931457240.997982
19036139334411289455527.638146914.4651409733243030.183740973324303170.99738
19036139334411289455281.262515114.4285539422482680.147153942248268170.955584
19036139334411289455933.53451514.3076397658521160.0262397658521162920.997489
19036139334411289456019.296878214.27863588374898-0.00276411625102035430.997654
19036139334411289456019.303801414.2828912657662080.0014912657662087270.997338

Repeat the plot with bad psfQfPerfect values excluded

Do the plot again but exclude low psfQfPerfect values.

In [11]:
# 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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,10))
for i, filter in enumerate("grizy"):
    pylab.subplot(511+i)
    w = np.where((dtab['filter']==filter) & (dtab['psfQfPerfect']>0.9))[0]
    pylab.plot(t[w],dmag[w],'-o')
    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))
    pylab.xlim(xlim)
    pylab.ylim(ylim)
    if i==0:
        pylab.title(objname)
pylab.xlabel('Time [MJD]')
pylab.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).

In [12]:
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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,10))
for i, filter in enumerate("grizy"):
    pylab.subplot(511+i)
    w = np.where((dtab['filter']==filter) & (dtab['psfQfPerfect']>0.9))[0]
    w = w[np.argsort(t[w])]
    pylab.plot(t[w],dmag[w],'-o')
    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))
    pylab.xlim(xlim)
    pylab.ylim(ylim)
    if i==0:
        pylab.title(objname)
pylab.xlabel('Phase')
pylab.tight_layout()

Repeat search using eclipsing binary KIC 2161623

From Villanova Kepler Eclipsing Binaries

In [13]:
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")
results = jobs.quick(query, task_name="python cone search")
tab = fixcolnames(ascii.read(results))
tab
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(291.7444617740266,37.5910012223162,0.016666666666666666) 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

Out[13]:
Table length=1
objIDraMeandecMeannDetectionsngnrninznygMeanPSFMagrMeanPSFMagiMeanPSFMagzMeanPSFMagyMeanPSFMag
int64float64float64int64int64int64int64int64int64float64float64float64float64float64
153102917444859851291.7444645137.591000167101612151414.599814.282114.158714.200414.0672

Get the detection information

This time include the psfQfPerfect limit directly in the database query.

In [14]:
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)

dresults = jobs.quick(query, task_name="python detection search")
dtab = fixcolnames(ascii.read(dresults))

# 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
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=153102917444859851
    ) d
join Filter f on d.filterID=f.filterID
where psfQfPerfect>0.9
order by d.filterID, obsTime

Out[14]:
Table length=45
objIDdetectIDfilterobsTimeradecpsfFluxpsfFluxErrpsfMajorFWHMpsfMinorFWHMpsfQfPerfectapFluxapFluxErrinfoFlaginfoFlag2infoFlag3magmeanmagdmag
int64int64str1float64float64float64float64float64float64float64float64float64float64int64int64int64float64float64float64
15310291744485985190150443710000088g55002.5047803291.7444639137.59100060.005010086.65687e-061.987441.812110.9994210.005174137.31225e-0610276051712812481542414.599814.6503883483652380.05058834836523829
15310291744485985190151959710000087g55002.5199475291.7444661137.591001560.004873346.48844e-061.984871.707930.9992340.005115427.27714e-0610276051712812481542414.599814.6804332199703220.08063321997032169
153102917444859851131534996560000129g55416.3502179291.7444608237.590996920.005256228.41646e-061.624031.578930.9985660.00529928.62578e-06102760517128737491214.599814.598316163992422-0.0014838360075781765
153102917444859851131536180650000138g55416.362065291.7444593737.590997330.005330288.62169e-061.320811.271860.9987570.005400828.91577e-0610276051712812481542414.599814.58312494213292-0.016675057867079346
153102917444859851204528930560000076g56146.2895546291.7444591437.590980730.005223998.43794e-061.333691.292110.9997580.005325068.88465e-06102760517128737491214.599814.6049941576962470.0051941576962466485
153102917444859851204530041560000085g56146.3006673291.7444648637.590976490.005235398.44281e-061.428831.226830.9646460.005324198.85594e-06102760517128737491214.599814.6026274021838770.0028274021838772256
153102917444859851241032265120000094g56511.322909291.7444565737.590997520.005286418.3427e-061.151231.107890.9987220.005374848.8135e-06102760517128737491214.599814.59209789291219-0.007702107087810717
153102917444859851241033321120000081g56511.3334693291.7444594537.590996430.005278898.31263e-061.078980.9873230.9952790.00533468.77771e-06102760517128737491214.599814.59364346902239-0.006156530977609265
15310291744485985191252543710000102r55013.5256677291.7444598537.590994540.007002438.7848e-061.901871.686890.9989620.007235139.82915e-06102760517128737491214.282114.2868780598010370.00477805980103696
153102917444859851131334556560000161r55414.3458067291.744462237.590998720.007148638.94809e-061.348021.239560.9978170.007153779.46631e-06102760517128737491214.282114.264442951598488-0.017657048401511943
.........................................................
153102917444859851159961616410000125y55700.6163536291.7444634637.590997680.008672862.66106e-050.7800350.6504660.9595020.008731062.5663e-05102760517128737491214.067214.054594160173803-0.012605839826196785
153102917444859851160460087410000136y55705.6010605291.7444637237.590996670.008389142.4518e-050.8922940.7857140.9974980.008507052.38093e-05102760517128737491214.067214.0907063948286750.02350639482867578
153102917444859851160461274410000134y55705.6129264291.7444647937.590994840.008463762.45964e-050.8435670.7193540.9821990.008459892.36938e-05102760517128737491214.067214.0810916501795480.01389165017954852
153102917444859851171734372010000127y55818.3439014291.7444641937.590999730.008995372.60801e-051.256050.9636240.9975090.008931252.37426e-051027605171285770656014.067214.014952421081372-0.05224757891862808
153102917444859851171735523010000129y55818.3554064291.7444605837.591002080.008938852.61524e-051.26840.9412590.9979290.008983012.38145e-051027605171284092934414.067214.02179587604703-0.04540412395297011
153102917444859851194462075640000160y56045.6209324291.744462737.590999140.009013122.46438e-051.411571.293440.9983750.009005782.28883e-051027605171285770656014.067214.01281211682762-0.05438788317237986
153102917444859851194462752640000158y56045.6277083291.7444598937.590997180.008816422.44413e-051.403771.198980.9995120.008887092.26963e-05102760517128737491214.067214.03676932228063-0.03043067771936947
153102917444859851213319324300000159y56234.1934155291.7444637337.59100090.008511512.42822e-051.045040.8496060.994970.008423522.41261e-05102760517128737491214.067214.0749834656833920.007783465683392166
153102917444859851213319880300000163y56234.1989869291.7444691537.591000790.008522162.36733e-051.038710.8439690.9949770.008332252.42534e-051027605171287448377614.067214.07362579090310.006425790903099582
153102917444859851245335559360000173y56554.3560604291.7444648837.590998740.008707051.611e-051.429061.389350.9596460.008868841.39697e-051027605171284092934414.067214.050322404022426-0.01687759597757399
In [15]:
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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,10))
for i, filter in enumerate("grizy"):
    pylab.subplot(511+i)
    w = np.where(dtab['filter']==filter)[0]
    pylab.plot(t[w],dmag[w],'-o')
    magmean = dtab['magmean'][w[0]]
    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))
    pylab.xlim(xlim)
    pylab.ylim(ylim)
    if i==0:
        pylab.title(objname)
pylab.xlabel('Time [MJD]')
pylab.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.

In [16]:
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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,6))
pylab.plot(t,dmag,'-o')
pylab.xlim(xlim)
pylab.ylim(ylim)
pylab.xlabel('Phase')
pylab.ylabel('Delta magnitude from mean [mag]')
pylab.title(objname)
pylab.tight_layout()

Repeat search for another eclipsing binary KIC 8153568

In [17]:
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")
results = jobs.quick(query, task_name="python cone search")
tab = fixcolnames(ascii.read(results))
tab
Out[17]:
Table length=1
objIDraMeandecMeannDetectionsngnrninznygMeanPSFMagrMeanPSFMagiMeanPSFMagzMeanPSFMagyMeanPSFMag
int64float64float64int64int64int64int64int64int64float64float64float64float64float64
160802869044447231286.9044440244.0054736988161531101615.182514.989914.890715.199914.8484
In [18]:
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)

dresults = jobs.quick(query, task_name="python detection search")
dtab = fixcolnames(ascii.read(dresults))

# 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
Out[18]:
Table length=60
objIDdetectIDfilterobsTimeradecpsfFluxpsfFluxErrpsfMajorFWHMpsfMinorFWHMpsfQfPerfectapFluxapFluxErrinfoFlaginfoFlag2infoFlag3magmeanmagdmag
int64int64str1float64float64float64float64float64float64float64float64float64float64int64int64int64float64float64float64
16080286904444723191336429430000113g55014.3646518286.9044337444.005475750.0031095.98065e-061.781.238110.9983350.003124475.69291e-06102760517128737491215.182515.16844819489202-0.014051805107978765
16080286904444723191337480430000116g55014.3751566286.9044389444.005475020.003130566.00144e-061.669391.327730.9981970.003156135.71358e-0610276051712812481542415.182515.160944920386187-0.021555079613811756
160802869044447231126057565120000092g55361.57592286.9044531744.005474510.003095736.72312e-061.74981.624220.9988260.003100456.69387e-06102760517128737491215.182515.173092310398793-0.00940768960120586
160802869044447231126058741130000111g55361.587678286.9044515744.005480650.003028156.57752e-061.628391.525790.9983870.00310626.76425e-06102760517128737491215.182515.1970565394434590.01455653944345947
160802869044447231128744126310000284g55388.4415151286.904450244.005478910.0008064623.51444e-061.101320.9925010.9976310.0008104163.45294e-0610276051712812481542415.182516.633540228713331.45104022871333
160802869044447231128745341310000249g55388.4536712286.9044504344.005479360.0009817623.82737e-061.10851.074410.9974910.0009855843.78859e-0610276051712812481542415.182516.41998445419421.2374844541942007
160802869044447231164544841720000097g55746.4486739286.9044531544.005476460.003066556.32653e-061.328511.166930.998310.003047336.63894e-061027605171284092934415.182515.1833748744026880.0008748744026885191
160802869044447231164744408530000111g55748.4443394286.9044520744.005478350.003040576.5461e-061.371531.231930.9985440.003078286.69725e-06102760517128737491215.182515.1926124845225450.010112484522545984
160802869044447231164745526530000113g55748.4555151286.9044527544.00547970.003047286.57104e-061.275311.264240.9985190.003085646.71944e-06102760517128737491215.182515.19021909674720.0077190967472002825
160802869044447231235756984470000099g56458.5701018286.9044537844.005474560.003086676.74978e-061.565411.496550.90410.003092766.7849e-061027605171287448377615.182515.176274497593653-0.00622550240634645
.........................................................
160802869044447231281422935550000135z56915.2295351286.904444144.005471620.004052511.10757e-051.311431.133850.9995880.004052121.03719e-0510276051712812481542415.199914.880689762025602-0.3192102379743975
160802869044447231136027009120000145y55461.2702758286.9044550344.005478270.004297691.96242e-051.486521.389680.9986860.004233951.58905e-05102760517128737491214.848414.816912285364621-0.031487714635378694
160802869044447231136028204120000148y55461.2822322286.9044599444.005482290.004174612.0278e-051.765061.7130.9987050.004197361.56164e-05102760517128737491214.848414.8484602272529576.0227252957290034e-05
160802869044447231159959993400000147y55700.6001108286.9044469244.005480860.004009831.95314e-051.040480.8986510.9980290.003961771.71112e-051027605171284092934414.848414.8921850981421540.0437850981421537
160802869044447231159961165400000142y55700.611833286.9044446344.005480940.004096671.92559e-050.8760260.7410520.9988060.004039491.75641e-05102760517128737491214.848414.8689225461527830.02052254615278315
160802869044447231160459626400000159y55705.5964454286.904448944.00548110.003313431.64053e-051.01180.848480.9974650.003279931.48499e-0510276051712812481542414.848415.0993054998595560.2509054998595559
160802869044447231160460795400000171y55705.6081345286.9044468444.005479640.003211941.64645e-051.081210.9909370.9978730.00318011.44718e-0510276051712812481542414.848415.1330814401710430.2846814401710436
160802869044447231171023586640000136y55811.2360427286.904452544.005483570.00421441.82169e-050.7231510.6624750.99780.004252461.73593e-05102760517128734214414.848414.838160616795776-0.01023938320422424
160802869044447231171024776640000151y55811.2479472286.9044527644.005482270.004250271.862e-050.7636580.6921140.9977030.004313361.74107e-05102760517128734214414.848414.828958700882717-0.01944129911728254
160802869044447231245335002470000153y56554.3504901286.9044508544.005477840.004270931.24367e-051.463481.403940.9726070.004264339.76066e-06102760517128737491214.848414.823693866365545-0.02470613363445473
In [19]:
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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(10,10))
for i, filter in enumerate("grizy"):
    pylab.subplot(511+i)
    w = np.where(dtab['filter']==filter)[0]
    pylab.plot(t[w],dmag[w],'-o')
    magmean = dtab['magmean'][w[0]]
    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))
    pylab.xlim(xlim)
    pylab.ylim(ylim)
    if i==0:
        pylab.title(objname)
pylab.xlabel('Time [MJD]')
pylab.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.

In [20]:
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])

pylab.rcParams.update({'font.size': 14})
pylab.figure(1,(12,6))
pylab.plot(t,dmag,'-o')
pylab.xlim(xlim)
pylab.ylim(ylim)
pylab.xlabel('Phase')
pylab.ylabel('Delta magnitude from mean [mag]')
pylab.title(objname)
pylab.tight_layout()