Query Pan-STARRS DR2 catalog using MAST API¶
The new MAST interface to the Pan-STARRS catalog supports queries to both the DR1 and DR2 PS1 catalogs. It also has an associated API, which is used in this script.
This script shows how to query the Pan-STARRS DR2 catalog using the PS1 search API. 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 notebook is available for download.
In [1]:
from astropy.io import ascii
from astropy.table import Table
import sys
import re
import numpy as np
import matplotlib.pyplot as plt
import json
import requests
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
Useful functions¶
Execute PS1 searches and resolve names using MAST query.
In [2]:
def ps1cone(ra,dec,radius,table="mean",release="dr1",format="csv",columns=None,
baseurl="https://catalogs.mast.stsci.edu/api/v0.1/panstarrs", verbose=False,
**kw):
"""Do a cone search of the PS1 catalog
Parameters
----------
ra (float): (degrees) J2000 Right Ascension
dec (float): (degrees) J2000 Declination
radius (float): (degrees) Search radius (<= 0.5 degrees)
table (string): mean, stack, or detection
release (string): dr1 or dr2
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., 'nDetections.min':2)
"""
data = kw.copy()
data['ra'] = ra
data['dec'] = dec
data['radius'] = radius
return ps1search(table=table,release=release,format=format,columns=columns,
baseurl=baseurl, verbose=verbose, **data)
def ps1search(table="mean",release="dr1",format="csv",columns=None,
baseurl="https://catalogs.mast.stsci.edu/api/v0.1/panstarrs", verbose=False,
**kw):
"""Do a general search of the PS1 catalog (possibly without ra/dec/radius)
Parameters
----------
table (string): mean, stack, or detection
release (string): dr1 or dr2
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., 'nDetections.min':2). Note this is required!
"""
data = kw.copy()
if not data:
raise ValueError("You must specify some parameters for search")
checklegal(table,release)
if format not in ("csv","votable","json"):
raise ValueError("Bad value for format")
url = f"{baseurl}/{release}/{table}.{format}"
if columns:
# check that column values are legal
# create a dictionary to speed this up
dcols = {}
for col in ps1metadata(table,release)['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 checklegal(table,release):
"""Checks if this combination of table and release is acceptable
Raises a VelueError exception if there is problem
"""
releaselist = ("dr1", "dr2")
if release not in ("dr1","dr2"):
raise ValueError("Bad value for release (must be one of {})".format(', '.join(releaselist)))
if release=="dr1":
tablelist = ("mean", "stack")
else:
tablelist = ("mean", "stack", "detection")
if table not in tablelist:
raise ValueError("Bad value for table (for {} must be one of {})".format(release, ", ".join(tablelist)))
def ps1metadata(table="mean",release="dr1",
baseurl="https://catalogs.mast.stsci.edu/api/v0.1/panstarrs"):
"""Return metadata for the specified catalog and table
Parameters
----------
table (string): mean, stack, or detection
release (string): dr1 or dr2
baseurl: base URL for the request
Returns an astropy table with columns name, type, description
"""
checklegal(table,release)
url = f"{baseurl}/{release}/{table}/metadata"
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 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)
Use metadata query to get information on available columns¶
This query works for any of the tables in the API (mean, stack, detection).
In [3]:
meta = ps1metadata("mean","dr2")
meta
Out[3]:
Simple positional query¶
Search mean object table with nDetections > 1.
This searches the mean object catalog for objects within 50 arcsec of M87 (RA=187.706, Dec=12.391 in degrees). 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.
In [4]:
ra = 187.706
dec = 12.391
radius = 50.0/3600.0
constraints = {'nDetections.gt':1}
# strip blanks and weed out blank and commented-out values
columns = """objID,raMean,decMean,nDetections,ng,nr,ni,nz,ny,
gMeanPSFMag,rMeanPSFMag,iMeanPSFMag,zMeanPSFMag,yMeanPSFMag""".split(',')
columns = [x.strip() for x in columns]
columns = [x for x in columns if x and not x.startswith('#')]
results = ps1cone(ra,dec,radius,release='dr2',columns=columns,verbose=True,**constraints)
# print first few lines
lines = results.split('\n')
print(len(lines),"rows in results -- first 5 rows:")
print('\n'.join(lines[:6]))
Convert the results to an astropy table¶
The CSV results string is easily converted to an astropy table. This table is easily manipulated to extract information on individual columns or rows.
In [5]:
tab = ascii.read(results)
# improve the format
for filter in 'grizy':
col = filter+'MeanPSFMag'
try:
tab[col].format = ".4f"
tab[col][tab[col] == -999.0] = np.nan
except KeyError:
print("{} not found".format(col))
tab
Out[5]:
Query for a single object by objID¶
It is possible to query the catalog directly on the object identifier without any RA/Dec restriction. This might not be particularly useful when search for objects, but it is very useful when searching the detection table for time-series data.
In [6]:
results1 = ps1search(release='dr2',columns=columns,verbose=True,objid=122851876947049923)
tab1 = ascii.read(results1)
# improve the format
for filter in 'grizy':
col = filter+'MeanPSFMag'
try:
tab1[col].format = ".4f"
tab1[col][tab1[col] == -999.0] = np.nan
except KeyError:
print("{} not found".format(col))
tab1
Out[6]:
Search stack objects at same position¶
There is no need for the nDetections limit for stack objects, which can in fact have nDetections = 0 for objects that are too faint to be detected on single-epoch images. But we require primaryDetection=1 in order to eliminate duplicate sources at the edges of the skycell regions used for processing. (There is another column bestDetection that would be better suited for this test but is currently not correct in the database.)
In [7]:
sconstraints = {'primaryDetection':1}
scolumns = """objID,raMean,decMean,nDetections,ng,nr,ni,nz,ny,
nStackDetections,primaryDetection,
gPSFMag,rPSFMag,iPSFMag,zPSFMag,yPSFMag""".split(',')
# strip blanks and weed out blank and commented-out values
scolumns = [x.strip() for x in scolumns]
scolumns = [x for x in scolumns if x and not x.startswith('#')]
sresults = ps1cone(ra,dec,radius,table="stack",release="dr2",columns=scolumns,verbose=True,**sconstraints)
stab = ascii.read(sresults)
for col in scolumns:
try:
stab[col]
except KeyError:
print(col,"not found")
# improve the format
for filter in 'grizy':
col = filter+'PSFMag'
try:
stab[col].format = ".4f"
stab[col][stab[col] == -999.0] = np.nan
except KeyError:
print("{} not found".format(col))
stab
Out[7]:
Explore the tables¶
Match the stack and mean tables and look at the subset of sources that are detected in both catalogs. The commented-out lines show how to restrict the joined table to only stack detections and to only stack non-detections.
In [8]:
from astropy.table import join
jtab = join(stab,tab,join_type='outer')
jtab.sort('objID')
jtab_both = jtab[(jtab['primaryDetection']==1) & (jtab['nDetections']>1)]
#jtab[jtab['nStackDetections'].mask].show_in_notebook()
#jtab[~jtab['nStackDetections'].mask].show_in_notebook()
#jtab.show_in_notebook()
jtab_both.show_in_notebook()
Out[8]:
Plot RA & Dec positions of mean and stack objects¶
Note that raMean and decMean are defined for all objects, including stack-only objects. For objects detected only on the stacked images, the raStack and decStack values are given in the raMean and decMean columns. That makes it simple to analyze the positions without testing to see what positions are available.
In [9]:
plt.rcParams.update({'font.size': 16})
plt.figure(1,(10,10))
plt.plot(tab['raMean'], tab['decMean'], 'ro', label='Mean (nDet>1)')
plt.plot(stab['raMean'], stab['decMean'], 'bo', label='Stack')
plt.plot(jtab_both['raMean'], jtab_both['decMean'], 'go', label='Both')
plt.xlabel('RA [deg]')
plt.ylabel('Dec [deg]')
plt.legend(loc='best')
Out[9]:
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. A small search radius is used so only a single object is returned.
In [10]:
objname = 'KQ UMa'
ra, dec = resolve(objname)
radius = 1.0/3600.0 # radius = 1 arcsec
results = ps1cone(ra,dec,radius,release='dr2',columns=columns,**constraints)
tab = ascii.read(results)
# improve the format
for filter in 'grizy':
col = filter+'MeanPSFMag'
tab[col].format = ".4f"
tab[col][tab[col] == -999.0] = np.nan
tab
Out[10]:
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.
In [11]:
def addfilter(dtab):
"""Add filter name as column in detection table by translating filterID
This modifies the table in place. If the 'filter' column already exists,
the table is returned unchanged.
"""
if 'filter' not in dtab.colnames:
# the filterID value goes from 1 to 5 for grizy
id2filter = np.array(list('grizy'))
dtab['filter'] = id2filter[(dtab['filterID']-1).data]
return dtab
objid = tab['objID'][0]
dconstraints = {'objID': objid}
dcolumns = ("""objID,detectID,filterID,obsTime,ra,dec,psfFlux,psfFluxErr,psfMajorFWHM,psfMinorFWHM,
psfQfPerfect,apFlux,apFluxErr,infoFlag,infoFlag2,infoFlag3""").split(',')
# strip blanks and weed out blank and commented-out values
dcolumns = [x.strip() for x in dcolumns]
dcolumns = [x for x in dcolumns if x and not x.startswith('#')]
dresults = ps1search(table='detection',release='dr2',columns=dcolumns,**dconstraints)
dtab = addfilter(ascii.read(dresults))
dtab.sort('obsTime')
dtab
Out[11]:
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 [12]:
# 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.
In [13]:
# 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.
In [14]:
dtab[ (dtab['filter']=='z') & (np.abs(mag-tab['zMeanPSFMag'][0]) > 2) ]
Out[14]:
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 [15]:
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[15]:
Repeat the plot with bad psfQfPerfect values excluded¶
Do the plot again but exclude low psfQfPerfect values.
In [16]:
# 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).
In [17]:
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¶
In [18]:
objname = 'KIC 2161623'
ra, dec = resolve(objname)
radius = 1.0/3600.0 # radius = 1 arcsec
results = ps1cone(ra,dec,radius,release='dr2',columns=columns,**constraints)
tab = ascii.read(results)
# improve the format
for filter in 'grizy':
col = filter+'MeanPSFMag'
tab[col].format = ".4f"
tab[col][tab[col] == -999.0] = np.nan
tab
Out[18]:
Get the detection information¶
This time include the psfQfPerfect limit directly in the database query.
In [19]:
objid = tab['objID'][0]
dconstraints = {'objID': objid, 'psfQfPerfect.min': 0.9}
dcolumns = ("""objID,detectID,filterID,obsTime,ra,dec,psfFlux,psfFluxErr,psfMajorFWHM,psfMinorFWHM,
psfQfPerfect,apFlux,apFluxErr,infoFlag,infoFlag2,infoFlag3""").split(',')
# strip blanks and weed out blank and commented-out values
dcolumns = [x.strip() for x in dcolumns]
dcolumns = [x for x in dcolumns if x and not x.startswith('#')]
dresults = ps1search(table='detection',release='dr2',columns=dcolumns,**dconstraints)
dtab = addfilter(ascii.read(dresults))
dtab.sort('obsTime')
# 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[19]:
In [20]:
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.
In [21]:
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¶
In [22]:
objname = 'KIC 8153568'
ra, dec = resolve(objname)
radius = 1.0/3600.0 # radius = 1 arcsec
results = ps1cone(ra,dec,radius,release='dr2',columns=columns,**constraints)
tab = ascii.read(results)
# improve the format
for filter in 'grizy':
col = filter+'MeanPSFMag'
tab[col].format = ".4f"
tab[col][tab[col] == -999.0] = np.nan
tab
Out[22]:
In [23]:
objid = tab['objID'][0]
dconstraints = {'objID': objid, 'psfQfPerfect.min': 0.9}
dcolumns = ("""objID,detectID,filterID,obsTime,ra,dec,psfFlux,psfFluxErr,psfMajorFWHM,psfMinorFWHM,
psfQfPerfect,apFlux,apFluxErr,infoFlag,infoFlag2,infoFlag3""").split(',')
# strip blanks and weed out blank and commented-out values
dcolumns = [x.strip() for x in dcolumns]
dcolumns = [x for x in dcolumns if x and not x.startswith('#')]
dresults = ps1search(table='detection',release='dr2',columns=dcolumns,**dconstraints)
dtab = addfilter(ascii.read(dresults))
dtab.sort('obsTime')
# 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[23]:
In [24]:
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.
In [25]:
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()
In [ ]: