Query PS1 catalog database for moving targets¶

This notebook is available for download. In the future it will be incorporated into the MAST notebooks repository.

This searches for Pan-STARRS1 catalog detections of the asteroid 174361 Rickwhite (2002 TV315). The CADC solar system image search API is used to get an initial list of observations. Then the JPL Horizons database is queried (using the astroquery.jplhorizons module) to get the position of the target as seen from the PS1 observatory on Mt Haleakala at the exact times of the PS1 observations. Those positions are passed to the PS1 CasJobs database where they are cross-matched with the PS1 detection catalog.

For asteroid rickwhite, there are 93 potential observations identified. Image cutouts are extracted from all those images using the ps1images image cutout service. 30 images are blank because the object position falls in a chip gap or off the edge of the camera. The 63 remaining images are displayed and show that the positions are indeed correct. Some of the non-blank images still have gaps or other artifacts at the object position.

Using CasJobs, the PS1 catalog positions for the good images are queried and all matching detections from those images are returned. 46 matches are found for the 63 images.

The time for all queries is about 1 minute (including the time to extract the image cutouts and to query the MAST database). The time and number of results returned should be roughly similar for most other solar system objects in the PS1 observing area.

Dependencies:

  • mastcasjobs module to query MAST CasJobs database from Python
  • astroquery.jplhorizons to access JPL Horizons ephemerides for objects
  • Pillow (PIL image module) to read JPEG images
  • Custom code is used to access the CADC SSOS service

R. White, 2025 February 3

Imports and configuration¶

In [1]:
import astropy
from astropy.table import Table, join, vstack
from astropy.time import Time
from astropy.io import fits

from astroquery.jplhorizons import Horizons
import mastcasjobs
from PIL import Image, UnidentifiedImageError

import numpy as np
import matplotlib.pyplot as plt
from io import BytesIO, StringIO
import time
import requests
import sys
import os
import getpass

saveplots = False # set to save some of the plots to png files
In [2]:
# code to customize display
from IPython.display import display, Markdown, Latex, HTML

# Set page width to fill browser for longer output lines
# and make tables left-aligned
display(HTML("""
<style>
.container { width:100% !important; }
table { align:left; display:block; }
</style>
"""))
# set width for pprint
astropy.conf.max_width = 150

Get and check CasJobs user login information¶

This uses the CASJOBS_USERID and CASJOBS_PW environment variables if they are defined. If prompts for them if the environment variables are not defined.

If you do not have a CasJobs account, create one using the MAST CasJobs website.

In [3]:
if not os.environ.get('CASJOBS_USERID'):
    os.environ['CASJOBS_USERID'] = input('Enter Casjobs UserID:')
if not os.environ.get('CASJOBS_PW'):
    os.environ['CASJOBS_PW'] = getpass.getpass('Enter Casjobs password:')

jobs = mastcasjobs.MastCasJobs(context="MyDB")
print("Confirmed MAST CasJobs account access")

Confirmed MAST CasJobs account access

Function to search moving object path using CADC¶

See the CADC moving target API documentation for more details about the moving target search.

In [4]:
def cadc_ssos_query(object_name, search="bynameCADC", epoch1=54985, epoch2=57079,
                    xyres="no", telinst="Pan-STARRS1", lang="en", format="tsv",
                    url="https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/cadcbin/ssos/ssosclf.pl"):
    """Use CADC moving object query to find PS1 observations by target name
    
    The epoch parameters give the range in MJD for the PS1 observations.
    The only parameter that might be usefully modified is the search.  
    E.g., use search="bynameHorizons" to search using the JPL Horizons emphemeris rather
    than the ephemeris cached at CADC.  That might be useful if the CADC cache has become
    out-of-date.
    
    This return an astropy table with the observations.
    """
    r = requests.get(url, params=dict(lang=lang, object=object_name, search=search,
                                      epoch1=epoch1, epoch2=epoch2, xyres=xyres,
                                      telinst=telinst, format=format))
    r.raise_for_status()
    # raise exception for object not found
    msg = r.text.split('\n',1)[0]
    if msg.lower().find("not found") >= 0:
        if msg.lower().endswith("<br/>"):
            msg = msg[:-5]
            raise ValueError(f"Error from CADC search: {msg}")
    tab = Table.read(r.text, format="ascii.csv", delimiter="\t",
                     converters=dict(Image_target=str))
    return tab

Functions to get images from PS1 cutout server¶

In [5]:
def getimages(tab, format="jpg", verbose=False, **kw):
    
    """Return list of images specified by the table with usual columns
    (projcell, subcell, image, ra and dec columns)
    
    Parameters:
    format = data format (options are "jpg", "png", "fits")
    verbose = if true, print progress info
    Other parameters passed on to tab2cutout:
    size = extracted image size in pixels (0.25 arcsec/pixel)
    output_size = output (display) image size in pixels (default = size).
                  output_size has no effect for fits format images.
    racol, deccol, imagecol = names of table columns
    Returns the image as a PIL Image or an astropy fits hdulist
    """
    
    if format not in ("jpg", "png", "fits"):
        raise ValueError("format must be jpg, png or fits")
    urls = tab2cutout(tab, format=format, verbose=verbose, **kw)
    images = []
    if verbose:
        t0 = time.time()
        print(f"{time.time()-t0:.1f} s: completed 0 of {len(urls)} images", end="")
    for i, url in enumerate(urls,start=1):
        r = requests.get(url)
        if format == "fits":
            images.append(fits.open(BytesIO(r.content)))
        else:
            try:
                images.append(Image.open(BytesIO(r.content)))
            except UnidentifiedImageError as e:
                print(f"\nImage not found {url}")
                images.append(None)
        if verbose:
            print(f"\r{time.time()-t0:.1f} s: completed {i} of {len(urls)} images", end="")
    if verbose:
        print()
    return images


def tab2path(tab, imagecol='Image', verbose=False):
    """Convert table with usual columns (projcell, subcell, filter, expStart) to warp filter path
    
    Convert expStart from TAI to UTC for the file label
    
    Returns a list of file paths
    """
    for c in ('projcell','subcell',imagecol):
        assert c in tab.colnames, f"'{c}' column not found in table"
    rv = []
    for i, row in enumerate(tab):
        projcell = f"{row['projcell']:04d}"
        subcell = f"{row['subcell']:03d}"
        image = row[imagecol]
        rv.append(f"/rings.v3.skycell/{projcell}/{subcell}/{image}.fits")
    return rv


def tab2cutout(tab, size=64, output_size=256, format="jpg", autoscale=99.5, asinh=True,
               racol="ra", deccol="dec", imagecol="Image", verbose=False,
               urlbase="https://ps1images.stsci.edu/cgi-bin/fitscut.cgi"):
    """Convert table with usual columns (projcell, subcell, image, ra and dec columns) 
    to URL for warp image cutout
    
    Returns a list of URLs
    """
    for c in ('projcell', 'subcell', imagecol, racol, deccol):
        assert c in tab.colnames, f"'{c}' column not found in table"
    paths = tab2path(tab, imagecol=imagecol, verbose=verbose)
    rv = []
    for row, path in zip(tab, paths):
        params = dict(red=path, ra=row[racol], dec=row[deccol], size=size, format=format)
        if filter != "fits":
            # these parameters are used only for non-fits formats
            params['autoscale'] = autoscale
            if output_size is not None:
                params['output_size'] = output_size
            if asinh:
                params['asinh'] = 1
        sparams = '&'.join([f"{key}={str(value)}" for key, value in params.items()])
        rv.append(f"{urlbase}?{sparams}")
    return rv

Functions to locate PanSTARRS skycell for a given RA/Dec¶

In [6]:
# pixel scale is 0.25 arcsec
pixscale = 0.25

# table of rings info
rings = Table.read("""zone projcell nband dec dec_min dec_max xcell ycell crpix1 crpix2
13 487 72 -38.0 -39.98569922229713 -35.98748871852416 6305 6279 239.5 238.0
14 559 76 -34.0 -35.98748871852416 -31.988947843740302 6265 6274 237.5 237.5
15 635 79 -30.0 -31.988947843740302 -27.990585819085737 6274 6269 239.5 241.5
16 714 82 -26.0 -27.990585819085737 -23.991856317278298 6255 6265 241.5 240.5
17 796 84 -22.0 -23.991856317278298 -19.99333106003115 6279 6261 241.0 242.0
18 880 86 -18.0 -19.99333106003115 -15.994464100667857 6274 6258 241.5 238.0
19 966 88 -14.0 -15.994464100667857 -11.995671545879146 6242 6254 240.5 239.5
20 1054 89 -10.0 -11.995671545879146 -7.997014822849273 6248 6250 240.0 242.0
21 1143 89 -6.0 -7.997014822849273 -3.998086931795508 6291 6247 237.5 240.0
22 1232 90 -2.0 -3.998086931795508 1.3877787807814457e-17 6240 6243 240.0 242.0
23 1322 90 2.0 1.3877787807814457e-17 3.998086931795508 6240 6243 240.0 242.0
24 1412 89 6.0 3.998086931795508 7.997014822849273 6291 6247 237.5 240.0
25 1501 89 10.0 7.997014822849273 11.995671545879146 6248 6250 240.0 242.0
26 1590 88 14.0 11.995671545879146 15.994464100667857 6242 6254 240.5 239.5
27 1678 86 18.0 15.994464100667857 19.99333106003115 6274 6258 241.5 238.0
28 1764 84 22.0 19.99333106003115 23.991856317278298 6279 6261 241.0 242.0
29 1848 82 26.0 23.991856317278298 27.990585819085737 6255 6265 241.5 240.5
30 1930 79 30.0 27.990585819085737 31.988947843740302 6274 6269 239.5 241.5
31 2009 76 34.0 31.988947843740302 35.98748871852416 6265 6274 237.5 237.5
32 2085 72 38.0 35.98748871852416 39.98569922229713 6305 6279 239.5 238.0
33 2157 68 42.0 39.98569922229713 43.98384988850235 6320 6284 239.5 240.0
34 2225 64 46.0 43.98384988850235 47.98179099058205 6307 6289 238.0 242.0
35 2289 60 50.0 47.98179099058205 51.97972846465295 6261 6295 241.0 239.5
36 2349 55 54.0 51.97972846465295 55.97676677254534 6283 6302 239.0 242.0
37 2404 50 58.0 55.97676677254534 59.97411834356914 6278 6311 238.5 237.5
38 2454 45 62.0 59.97411834356914 63.97005521318337 6240 6319 240.0 241.0
39 2499 39 66.0 63.97005521318337 67.96475284198297 6307 6333 239.5 240.0
40 2538 33 70.0 67.96475284198297 71.95817764373305 6365 6350 238.5 240.0
41 2571 27 74.0 71.95817764373305 75.94974140116577 6413 6371 240.5 242.0
42 2598 21 78.0 75.94974140116577 79.93968190571101 6452 6398 240.0 242.0
43 2619 15 82.0 79.93968190571101 83.93610604798963 6481 6429 241.0 242.0
44 2634 9 86.0 83.93610604798963 87.96940975734594 6501 6419 239.0 239.5
45 2643 1 90.0 87.96940975734594 90.0001 6240 6240 239.5 240.0
""", format="ascii.csv", delimiter=" ")

dec_limit = rings['dec_min'].min()

def findskycell(ra, dec):

    """Given input arrays RA, DEC (degrees), returns an astropy table with the skycell info

    RA and Dec can be scalars or arrays.
    This returns just the best sky cell for each ra/dec (the one where the position is closest to the center).
    This uses the rings table for info on the tessellation.
    The returned table has these columns:
    Column      Value
    ra          Input position
    dec         
    index       0-based index into original array
    projcell    Projection cell (0 if outside coverage)
    subcell     Subcell (0..99)
    crval1      Sky position of reference pixel for projection cell
    crval2      
    crpix1      Reference pixel in this skycell
    crpix2      
    x           Source pixel position in this skycell
    y           
    iring       Index into rings array
    """

    if np.isscalar(ra) and np.isscalar(dec):
        return _findskycell_array(np.array([ra]),np.array([dec]))
    if len(ra) == len(dec):
        return _findskycell_array(np.asarray(ra),np.asarray(dec))
    else:
        raise ValueError("ra and dec must both be scalars or be matching length arrays")


def _findskycell_array(ra, dec):

    """Internal function: ra and dec are known to be arrays"""

    if ra.ndim != 1 or dec.ndim != 1:
        raise ValueError("ra and dec must be 1-D arrays")
    index = np.arange(len(ra),dtype=int)
    # find dec zone where rings.dec_min <= dec < rings.dec_max
    idec = np.searchsorted(rings['dec_max'], dec)

    # special handling at pole where overlap is complicated
    # do extra checks for top 2 rings
    # always start with the ring just below the pole
    nearpole = np.where(idec >= len(rings)-2)
    idec[nearpole] = len(rings)-2

    nband = rings['nband'][idec]
    # get normalized RA in range 0..360
    nra = ra % 360.0
    ira = np.rint(nra*nband/360.0).astype(int) % nband

    projcell = rings['projcell'][idec] + ira
    dec_cen = rings['dec'][idec]
    ra_cen = ira*360.0/nband

    # locate subcell(s) within the projection cell

    # use tangent project to get pixel offsets
    x, y = sky2xy_tan(nra, dec, ra_cen, dec_cen)

    pad = 480

    if len(nearpole[0]) > 0:
        # handle the points near the pole (if any)
        # we know that this "ring" has only a single field
        projcell2 = rings['projcell'][-1]
        dec_cen2 = rings['dec'][-1]
        ra_cen2 = 0.0
        x2, y2 = sky2xy_tan(nra[nearpole], dec[nearpole], ra_cen2, dec_cen2)
        # compare x,y and x2,y2 to image sizes to select best image
        # returns a Boolean array with true for values where 2nd image is better
        use2 = poleselect(x[nearpole], y[nearpole], x2, y2, rings[-2], rings[-1], pad)
        if use2.any():
            # slightly obscure syntax here makes this work even if ra, dec are multi-dimensional
            wuse2 = np.where(use2)[0]
            w2 = tuple(xx[wuse2] for xx in nearpole)
            idec[w2] = len(rings)-1
            nband[w2] = 1
            ira[w2] = 0
            projcell[w2] = projcell2
            dec_cen[w2] = dec_cen2
            ra_cen[w2] = ra_cen2
            x[w2] = x2[wuse2]
            y[w2] = y2[wuse2]

    # compute the subcell from the pixel location
    px = rings['xcell'][idec]-pad
    py = rings['ycell'][idec]-pad
    k = np.rint(4.5+x/px).astype(int).clip(0,9)
    j = np.rint(4.5+y/py).astype(int).clip(0,9)
    subcell = 10*j + k

    # get pixel coordinates within the skycell image
    crpix1 = rings['crpix1'][idec] + px*(5-k)
    crpix2 = rings['crpix2'][idec] + py*(5-j)
    ximage = x + crpix1
    yimage = y + crpix2
    iring = idec

    # insert zeros where we are below lowest dec_min
    w = np.where(dec < dec_limit)
    projcell[w] = 0
    subcell[w] = 0
    crpix1[w] = 0
    crpix2[w] = 0
    ximage[w] = 0
    yimage[w] = 0

    # return an astropy Table
    return Table([ra,dec,index,projcell,subcell,ra_cen,dec_cen,crpix1,crpix2,ximage,yimage,iring],
                 names="ra,dec,index,projcell,subcell,crval1,crval2,crpix1,crpix2,x,y,iring".split(","))


def poleselect(x1, y1, x2, y2, rings1, rings2, pad):

    """Compares x,y values from 2 images to determine which is best

    Returns boolean array with True where x2,y2 is best
    """

    nx1 = 10*(rings1['xcell']-pad)+pad
    ny1 = 10*(rings1['ycell']-pad)+pad
    nx2 = 10*(rings2['xcell']-pad)+pad
    ny2 = 10*(rings2['ycell']-pad)+pad
    # compute minimum distances to image edges
    # note negative values are off edge
    d1 = np.minimum(np.minimum(x1+nx1//2,nx1//2-1-x1), np.minimum(y1+ny1//2,ny1//2-1-y1))
    d2 = np.minimum(np.minimum(x2+nx2//2,nx2//2-1-x2), np.minimum(y2+ny2//2,ny2//2-1-y2))
    return (d1 < d2)


def sky2xy_tan(ra, dec, ra_cen, dec_cen, crpix=(0.0,0.0)):

    """Convert RA,Dec sky position (degrees) to X,Y pixel position

    ra[n], dec[n] are input arrays in degrees
    ra_cen[n], dec_cen[n] are image centers in degrees
    crpix is the reference pixel position (x,y)
    Returns tuple (x,y) where x and y are arrays with pixel position for each RA,Dec
    """

    dtor = np.pi/180
    cd00 = -pixscale*dtor/3600
    cd01 = 0.0
    cd10 = 0.0
    cd11 = -cd00
    determ = cd00*cd11-cd01*cd10
    cdinv00 =  cd11/determ
    cdinv01 = -cd01/determ
    cdinv10 = -cd10/determ
    cdinv11 =  cd00/determ

    cos_crval1 = np.cos(dtor*dec_cen)
    sin_crval1 = np.sin(dtor*dec_cen)

    radif = (ra - ra_cen)*dtor
    w = np.where(radif > np.pi)
    radif[w] -= 2*np.pi
    w = np.where(radif < -np.pi)
    radif[w] += 2*np.pi

    decrad = dec*dtor
    cos_dec = np.cos(decrad)
    sin_dec = np.sin(decrad)
    cos_radif = np.cos(radif)
    sin_radif = np.sin(radif)
    h = sin_dec*sin_crval1 + cos_dec*cos_crval1*cos_radif
    xsi = cos_dec*sin_radif/h
    eta = (sin_dec*cos_crval1 - cos_dec*sin_crval1*cos_radif)/h
    xdif = cdinv00*xsi + cdinv01*eta
    ydif = cdinv10*xsi + cdinv11*eta
    return (xdif+crpix[0], ydif+crpix[1])


def xy2sky_tan(x, y, ra_cen, dec_cen, crpix=(0.0,0.0)):

    """Convert X,Y pixel position to RA,Dec sky position (degrees)

    x[n], y[n] are input arrays in pixels
    ra_cen[n], dec_cen[n] are image centers in degrees
    crpix is the reference pixel position (x,y)
    Returns tuple (ra,dec) where ra and dec are arrays with pixel position for each x,y
    """

    dtor = np.pi/180
    cd00 = -pixscale*dtor/3600
    cd01 = 0.0
    cd10 = 0.0
    cd11 = -cd00

    cos_crval1 = np.cos(dtor*dec_cen)
    sin_crval1 = np.sin(dtor*dec_cen)

    xdif = x - crpix[0]
    ydif = y - crpix[1]
    xsi = cd00*xdif + cd01*ydif
    eta = cd10*xdif + cd11*ydif
    beta = cos_crval1 - eta*sin_crval1
    ra = np.arctan2(xsi, beta) + dtor*ra_cen
    gamma = np.sqrt(xsi**2 + beta**2)
    dec = np.arctan2(eta*cos_crval1+sin_crval1, gamma)
    return (ra/dtor, dec/dtor)

Get image list from CADC moving object query ¶

This query takes about 10 seconds.

In [7]:
%%time 

source = "rickwhite"

# some other sample objects
# source = "33798"       # fast-moving object, 10 arcsec/hour
# source = "134340"      # Pluto
# source = "2003 YT1"    # asteroid that passes close to North Celestial Pole (only 2 epochs in PS1)
# source = "Cruithne"    # fast-moving, unusual Earth-associated orbit
# source = "5261"        # Mars-trojan asteroid
# source = "243"         # 243 Ida, first asteroid known to have a moon
# source = "2024 YR4"

cadc_tab = cadc_ssos_query(source)
print(f"Returned table has {len(cadc_tab)} rows")
cadc_tab[:8]

Returned table has 127 rows
CPU times: user 16.1 ms, sys: 3 ms, total: 19.1 ms
Wall time: 2.81 s
Out[7]:
Table length=8
ImageMJDFilterExptimeObject_RAObject_DecImage_targetTelescope/InstrumentMetaDataDatalink
str43float64str1int64float64float64str8str11int64int64
rings.v3.skycell.0956.062.wrp.y.54992_6094854992.6100484111y30319.176751857141-17.48580135266580956.062Pan-STARRS1----
rings.v3.skycell.0956.062.wrp.y.54992_6215554992.6221213111y30319.176589325724-17.48656667391820956.062Pan-STARRS1----
rings.v3.skycell.0783.051.wrp.i.55089_3047355089.3054836223i60304.364116324646-25.84291714492830783.051Pan-STARRS1----
rings.v3.skycell.0783.051.wrp.i.55089_3122255089.3129660222i60304.363786381988-25.84308024046160783.051Pan-STARRS1----
rings.v3.skycell.1240.063.wrp.y.55416_6158355416.6164040111y3032.5693639363845-1.353319102650011240.063Pan-STARRS1----
rings.v3.skycell.1240.063.wrp.y.55416_6265555416.6271279111y3032.5706287757064-1.353780617601911240.063Pan-STARRS1----
rings.v3.skycell.1240.051.wrp.z.55427_6279955427.6285596111z3033.5734442378021-1.961510433104871240.051Pan-STARRS1----
rings.v3.skycell.1240.051.wrp.z.55427_6285055427.6290649111z3033.5734768149145-1.961544150651371240.051Pan-STARRS1----

Recompute positions using accurate epochs¶

The interpolated positions from CADC can be off by a few arcsec, so recompute them using JPL Horizons at the exact times of the PS1 observations.

The location code F51 corresponds to the Pan-STARRS 1 observatory on Mt. Haleakala. The object positions are converted to RA and Dec as seen from the observatory.

Note that Horizons uses times in UTC rather than TAI (which is what is in the PS1 database). Appending TT to the start time is supposed to cause Horizons to use TT (= TAI $\pm0.002$ s), but that triggers some kind of parsing bug in the astroquery.Horizons module. Do the queries in UTC instead.

Divide this query into groups of 50 JD values to avoid errors from the service.

This query takes about 1 second.

In [8]:
# convert input epochs from TAI to UTC
jdepochs = Time(cadc_tab["MJD"],scale="tai",format="mjd").utc.jd.tolist()
tlist = []
block = 50
for k in range(0,len(jdepochs),block):
    obj = Horizons(id=source, location='F51', epochs=jdepochs[k:k+block])
    tlist.append(obj.ephemerides(extra_precision=True))
tab2 = vstack(tlist)
fullname = tab2['targetname'][0]
print(f"Got positions from Horizons for {len(jdepochs)} epochs with {len(tab2)} rows for {fullname}")
assert np.allclose(jdepochs, tab2['datetime_jd'].data)
# convert output epochs from UTC back to TAI for PS1 database matching
t2 = Time(tab2['datetime_jd'].data,format='jd',scale='utc').tai.mjd
assert np.allclose(t2,cadc_tab["MJD"].data)
tab2[:8]

Got positions from Horizons for 127 epochs with 127 rows for 174361 Rickwhite (2002 TV315)
Out[8]:
Table length=8
targetnamedatetime_strdatetime_jdHGsolar_presencelunar_presenceRADECRA_appDEC_appRA_rateDEC_rateAZELAZ_rateEL_ratesat_Xsat_Ysat_PANGsiderealtimeairmassmagextinctVsurfbrightilluminationillum_defectsat_sepsat_visang_widthPDObsLonPDObsLatPDSunLonPDSunLatSubSol_angSubSol_distNPole_angNPole_distEclLonEclLatrr_ratedeltadelta_ratelighttimevel_sunvel_obselongelongFlagalphalunar_elonglunar_illumsat_alphasunTargetPAvelocityPAOrbPlaneAngconstellationTDB-UTObsEclLonObsEclLatNPole_RANPole_DECGlxLonGlxLatsolartimeearth_lighttimeRA_3sigmaDEC_3sigmaSMAA_3sigmaSMIA_3sigmaTheta_3sigmaArea_3sigmaRSS_3sigmar_3sigmar_rate_3sigmaSBand_3sigmaXBand_3sigmaDoppDelay_3sigmatrue_anomhour_anglealpha_truePABLonPABLat
------dmag---------degdegdegdegarcsec / harcsec / hdegdegarcsec / minarcsec / minarcsecarcsecdegh---magmagmag / arcsec2%arcsecarcsec---arcsecdegdegdegdegdegarcsecdegarcsecdegdegAUkm / sAUkm / sminkm / skm / sdeg---degdeg%degdegdegdeg---sdegdegdegdegdegdeghminarcsecarcsecarcsecarcsecdegarcsec2arcseckmkm / sHzHzsdeghdegdegdeg
str29str24float64float64float64str1str1float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64int64float64int64float64str1int64int64int64int64int64float64float64int64int64float64float64float64float64float64float64float64float64float64float64str2float64float64float64float64float64float64float64str3float64float64float64int64int64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64
174361 Rickwhite (2002 TV315)2009-Jun-10 14:37:54.1832454993.10965489616.10.15Am319.175183288-17.48702118319.314219565-17.445783407-2.12126-9.19912184.61032405351.734420933858.3-68.01-396330.66-145783.37263.72821.48701509941.2720.16621.087--97.86733--444507.1*----------73.490.0----299.428641-1.2138282.9311942573462.3478962.24609047951069-17.761664218.680164916.089515218.2682872123.4742/L16.797825.893.02139.5953253.527264.891-2.99843Cap66.184675316.3601942-1.5860015----32.199262-39.7748694.22317177410.0003550.0170.0040.0170.002839.160.00015550.0188.43659e-070.020.065.6e-05142.83070.19940046216.7946307.8228-1.4144
174361 Rickwhite (2002 TV315)2009-Jun-10 14:55:17.2812454993.12172779316.10.15Nm319.175003548-17.487760789319.31404156-17.446524501-2.12504-9.18129191.22332544551.171403567845.13-164.46-396370.18-145789.42263.73221.77755794621.2820.16721.087--97.86796--444550.3*----------73.490.0----299.430765-1.2142742.9312106290392.34777932.24596675493409-17.727551618.6791359216.089410818.2333167123.4862/L16.795325.792.986239.5858253.528264.891-2.99772Cap66.184675316.3598085-1.5866567----32.198282-39.7749764.51287551010.0003550.0170.0040.0170.002839.160.00015560.0188.43569e-070.020.065.6e-05142.83290.48995517616.7921307.8237-1.415
174361 Rickwhite (2002 TV315)2009-Sep-15 07:19:19.7852455089.80509010416.10.15304.363508871-25.843867625304.518350914-25.813120861-7.0487-3.30645184.27475200543.359211645810.09-62.99472138.654-103506.83120.88620.53131647641.4540.18921.168--98.26698--460876.1*----------252.180.0----315.844676-4.5596193.0352469573481.36444922.3106983796918620.920815719.2174924215.435845821.2322411128.0211/T15.1258171.516.089436.988972.243269.5584.40871Cap66.18243300.859432-5.9961301----17.096164-29.50064420.9849150580.0003550.0180.0020.0180.001942.0270.00010890.01810.73881e-060.020.077.2e-05159.41960.23009308215.1292308.2877-5.3229
174361 Rickwhite (2002 TV315)2009-Sep-15 07:30:06.2642455089.81257250316.10.15304.363118317-25.844031687304.517959995-25.813287168-7.03858-3.28514187.59317841843.097713638804.83-111.62472114.077-103497.04120.88520.71138573611.4610.1921.168--98.26662--460848.5*----------252.180.0----315.845911-4.5598593.0352528531611.36437032.3107888327482620.941209819.218244715.435809321.2509173128.0135/T15.1273171.415.998136.99572.244269.5584.40887Cap66.18243300.8590495-5.9962122----17.095858-29.5003621.16453203580.0003550.0180.0020.0180.001942.0270.00010890.01810.73961e-060.020.077.2e-05159.42090.41018840315.1308308.2881-5.323
174361 Rickwhite (2002 TV315)2010-Aug-08 14:47:03.3072455417.11601049316.10.15Am32.570101791-1.35337625832.710947476-1.30000022416.76048-6.50686154.38387980265.834330269824.91364.75-365711.08-62465.88273.1541.50092226581.0950.14321.444--97.19585--379616.3*----------73.70.0----10.247126-11.5955452.960357257372-2.12842712.52415715782073-24.297186820.9927749515.902570425.9604158105.449/L19.281585.43.128655.124253.669248.1431.78785Cet66.18310330.0509646-13.6335521----162.493783-58.0336064.27301607630.0003550.0170.0060.0180.0026616.7320.00014850.01812.2141e-060.010.038.1e-05213.1751-0.67980756619.279920.0304-12.7972
174361 Rickwhite (2002 TV315)2010-Aug-08 15:02:29.8522455417.12673439316.10.15Am32.571300154-1.35384088832.712147823-1.30046501616.74478-6.50908163.52625372267.137166441869.65239.21-365748.07-62456.67273.1561.75900054911.0850.14121.444--97.1961--379648.8*----------73.70.0----10.249052-11.5956352.960344073409-2.1285342.52400676624768-24.267049820.9915241815.902653825.9295398105.458/L19.280685.63.083655.1159253.672248.1431.78872Cet66.18310330.0519592-13.6344018----162.496147-58.0333094.53040860970.0003550.0170.0060.0180.0026616.7320.00014850.01812.21361e-060.010.038.1e-05213.177-0.42180930619.27920.0318-12.7976
174361 Rickwhite (2002 TV315)2010-Aug-19 15:04:33.5502455428.12816609316.10.15A33.574050097-1.96160450333.715730495-1.9083679648.650545-10.0534190.36977747367.042780558888.9-152.09-404100.52-52478.66273.8122.51626011531.0850.14121.265--97.50912--413448.7*----------76.910.0----12.235226-11.681072.94647175111-2.23746972.37437634172711-22.725250419.7470858915.990601223.5942426114.8469/L18.1638119.777.455846.8448256.884248.032.72162Cet66.18286130.8084463-14.5496064----164.702584-57.945994.59845985240.0003550.0180.0060.0190.0029516.7060.00017540.01911.67631e-060.010.057.8e-05215.10990.26854474918.161421.3878-13.2816
174361 Rickwhite (2002 TV315)2010-Aug-19 15:05:17.2082455428.12867139316.10.15A33.574079276-1.96163835533.715759703-1.9084018188.650356-10.0534190.82980970867.011366365887.74-158.75-404102.33-52478.19273.8122.5284205251.0860.14121.265--97.50913--413450.3*----------76.910.0----12.235317-11.6810732.946471098088-2.23747472.3743697099712-22.72378519.7470307315.990605323.5928166114.8473/L18.1637119.777.461246.8444256.884248.032.72166Cet66.18286130.8084627-14.5496482----164.702666-57.9464.610588830.0003550.0180.0060.0190.0029516.7060.00017540.01911.67631e-060.010.057.8e-05215.110.28070321118.161321.3879-13.2816

Add new positions and differences to cadc_tab¶

Also add the illumination (illumination percentage) and delta (range to object in AU), which are useful to understand apparent brightness.

In [9]:
for c in ('ra','dec'):
    uc = c.upper()
    if c in cadc_tab.colnames:
        cadc_tab[c] = tab2[uc]
    else:
        cadc_tab.add_column(tab2[uc],name=c,index=cadc_tab.colnames.index('Object_RA'))

cadc_tab['dra'] = (cadc_tab['ra']-cadc_tab['Object_RA'])*3600*np.cos(np.radians(cadc_tab['dec']))
cadc_tab['ddec'] = (cadc_tab['dec']-cadc_tab['Object_Dec'])*3600
for c in ('ra','dec','Object_RA','Object_Dec'):
    cadc_tab[c].format = ".6f"
for c in ('dra','ddec'):
    cadc_tab[c].format = ".2f"
    cadc_tab[c].unit = "arcsec"

cadc_tab['illumination'] = tab2['illumination']
cadc_tab['illumination'].format = ".2f"
cadc_tab['delta'] = tab2['delta']
cadc_tab['delta'].format = ".2f"
cadc_tab['delta'].unit = "AU"

print(f"dra range  {cadc_tab['dra'].min():.2f} {cadc_tab['dra'].max():.2f}"
      f" MAD {np.median(np.abs(cadc_tab['dra'].data.filled())):.2f} arcsec")
print(f"ddec range {cadc_tab['ddec'].min():.2f} {cadc_tab['ddec'].max():.2f}"
      f" MAD {np.median(np.abs(cadc_tab['ddec'].data.filled())):.2f} arcsec")

dra range  -9.42 10.06 MAD 2.48 arcsec
ddec range -13.03 4.91 MAD 3.05 arcsec

Compute PS1 skycell position for each ra, dec value¶

Add the projcell and subcell to the table.

I could get this from the ForcedWarpMeta table in the database, but it is easier just to compute it from the RA and Dec values.

In [10]:
rv = findskycell(cadc_tab['ra'],cadc_tab['dec'])
for c in rv.colnames:
    if c not in cadc_tab.colnames:
        cadc_tab[c] = rv[c]
cadc_tab["crval1"].format = ".6f"
for c in ("x","y"):
    cadc_tab[c].format = ".2f"
cadc_tab[:8]
Out[10]:
Table length=8
ImageMJDFilterExptimeradecObject_RAObject_DecImage_targetTelescope/InstrumentMetaDataDatalinkdraddecilluminationdeltaindexprojcellsubcellcrval1crval2crpix1crpix2xyiring
degdegarcsecarcsec%AU
str43float64str1int64float64float64float64float64str8str11int64int64float64float64float64float64int64int64int64float64float64float64float64float64float64int64
rings.v3.skycell.0956.062.wrp.y.54992_6094854992.6100484111y30319.175183-17.487021319.176752-17.4858010956.062Pan-STARRS1-----5.39-4.3997.872.25095662318.139535-18.017623.5-5540.03397.481808.465
rings.v3.skycell.0956.062.wrp.y.54992_6215554992.6221213111y30319.175004-17.487761319.176589-17.4865670956.062Pan-STARRS1-----5.44-4.3097.872.25195662318.139535-18.017623.5-5540.03400.011797.825
rings.v3.skycell.0783.051.wrp.i.55089_3047355089.3054836223i60304.363509-25.843868304.364116-25.8429170783.051Pan-STARRS1-----1.97-3.4298.272.31278351302.926829-26.023341.5240.54719.582387.033
rings.v3.skycell.0783.051.wrp.i.55089_3122255089.3129660222i60304.363118-25.844032304.363786-25.8430800783.051Pan-STARRS1-----2.16-3.4398.272.31378351302.926829-26.023341.5240.54724.672384.723
rings.v3.skycell.1240.063.wrp.y.55416_6158355416.6164040111y3032.570102-1.35337632.569364-1.3533191240.063Pan-STARRS1----2.66-0.2197.202.52412406332.000000-2.011760.0-5521.03552.033789.819
rings.v3.skycell.1240.063.wrp.y.55416_6265555416.6271279111y3032.571300-1.35384132.570629-1.3537811240.063Pan-STARRS1----2.42-0.2297.202.52512406332.000000-2.011760.0-5521.03534.783783.129
rings.v3.skycell.1240.051.wrp.z.55427_6279955427.6285596111z3033.574050-1.96160533.573444-1.9615101240.051Pan-STARRS1----2.18-0.3497.512.37612405132.000000-2.023280.0242.0621.27784.249
rings.v3.skycell.1240.051.wrp.z.55427_6285055427.6290649111z3033.574079-1.96163833.573477-1.9615441240.051Pan-STARRS1----2.17-0.3497.512.37712405132.000000-2.023280.0242.0620.85783.759

Remove duplicate warp image entries¶

Note that several skycells can include a position. Those skycells, which come from the same single exposure, actually have the same data; they are not independent images. The findskycell() function just picks the best skycell (farthest from the edge). That is fine since each exposure that overlaps that skycell will generate a warp image.

Find the duplicate entries in the table that are from overlapping skycells, and keep the one that matches the projcell and subcell values. This is a little tricky because the corrected positions can cause objects to move from one skycell to another. Handle that as well in getting a unique set of images.

In [11]:
umjd, umjd_counts = np.unique(cadc_tab['MJD'], return_counts=True)
nuniq = len(umjd)
print(f"cadc_tab has {len(cadc_tab)-nuniq} duplicate entries, leaving {nuniq} unique entries")
keep = np.array([image.startswith(f"rings.v3.skycell.{projcell:04d}.{subcell:03d}.") 
                  for image, projcell, subcell in zip(cadc_tab['Image'].data,cadc_tab['projcell'].data,cadc_tab['subcell'].data)])

if keep.sum() != nuniq:
    # find which MJD values will disappear
    jtest = join(cadc_tab['MJD','Image','projcell','subcell'], cadc_tab['MJD','Image_target'][keep], join_type='left')
    w = np.where(jtest['Image_target'].mask)[0]
    if len(w) > nuniq-keep.sum():
        # this happens when there were multiple skycells matched from CADC, but they are
        # replaced by a single best skycell that is not either of the CADC entries
        # e.g. I found a case where CADC returned skycells 1589.003 and 1589.004 (at the same MJD)
        # but the best skycell is 1589.013
        # just keep one of the entries for duplicates at the same MJD
        j = len("rings.v3.skycell.nnnn.nnn.")
        print(f"First iteration of keep has {len(keep)} entries, {len(w)=}")
        cmulti = {}
        nchanged = 0
        for i in w:
            tt = jtest['MJD'][i]
            cmulti[tt] = cmulti.get(tt,0) + 1
            if cmulti[tt] == 2:
                # change filename for this row
                projcell = jtest['projcell'][i]
                subcell = jtest['subcell'][i]
                image = jtest['Image'][i]
                kk = np.where((cadc_tab['MJD'] == tt) & 
                              (cadc_tab['projcell'] == projcell) &
                              (cadc_tab['subcell'] == subcell) &
                              (cadc_tab['Image'] == image))[0]
                assert len(kk)==1, f"Trouble for {jtest['Image']} MJD {tt} {len(kk)=}"
                kk = kk[0]
                oldname = cadc_tab['Image'][kk]
                if oldname.startswith(f"rings.v3.skycell.{projcell:04d}."):
                    cadc_tab['Image'][kk] = f"rings.v3.skycell.{projcell:04d}.{subcell:03d}.{oldname[j:]}"
                    print(f"Changed {oldname} to {cadc_tab['Image'][kk]}")
                    nchanged += 1
        if nchanged == 0:
            print("Warning: no changes in attempt to fix mismatch")
        else:
            # recreate the keep and jtest arrays
            keep = np.array([image.startswith(f"rings.v3.skycell.{projcell:04d}.{subcell:03d}.") 
                      for image, projcell, subcell in zip(cadc_tab['Image'].data,cadc_tab['projcell'].data,cadc_tab['subcell'].data)])
            jtest = join(cadc_tab['MJD','Image','projcell','subcell'], cadc_tab['MJD','Image_target'][keep], join_type='left')
            w = np.where(jtest['Image_target'].mask)[0]
            print(f"Second iteration of keep has {len(keep)} entries, {len(w)=} (changed {nchanged} names)")

    assert len(w) == nuniq-keep.sum(), f"Mismatch {len(w)=} {nuniq=} {keep.sum()=} {nuniq-keep.sum()=}"
    assert (cadc_tab['Image'][w]==jtest['Image'][w]).all()
    # keep the image name from CADC, and change the projcell/subcell columns to match
    new_projcell = [int(x.split('.')[0]) for x in cadc_tab['Image_target'][w].data]
    new_subcell = [int(x.split('.')[1]) for x in cadc_tab['Image_target'][w].data]
    print(f"Patching projcell and subcell for {len(w)} images")
    cadc_tab['projcell'][w] = new_projcell
    cadc_tab['subcell'][w] = new_subcell
    # recompute the keep array, which now should be complete
    keep = np.array([image.startswith(f"rings.v3.skycell.{projcell:04d}.{subcell:03d}.") 
                  for image, projcell, subcell in zip(cadc_tab['Image'].data,cadc_tab['projcell'].data,cadc_tab['subcell'].data)])
    assert keep.sum() == nuniq, f"ERROR: unable to fix discrepancy, keep.sum={keep.sum()} nuniq={nuniq}"
    print("Modified lines")
    print(cadc_tab['Image','MJD','Filter','Image_target','projcell','subcell'][w])

cadc_tab = cadc_tab[keep]
cadc_tab[:8]

cadc_tab has 34 duplicate entries, leaving 93 unique entries
Patching projcell and subcell for 3 images
Modified lines
                   Image                          MJD        Filter Image_target projcell subcell
------------------------------------------- ---------------- ------ ------------ -------- -------
rings.v3.skycell.1060.082.wrp.r.55498_38975 55498.3903793815      r     1060.082     1060      82
rings.v3.skycell.1060.082.wrp.r.55498_40076 55498.4013951815      r     1060.082     1060      82
rings.v3.skycell.0970.099.wrp.r.56971_28764 56971.2883144167      r     0970.099      970      99
Out[11]:
Table length=8
ImageMJDFilterExptimeradecObject_RAObject_DecImage_targetTelescope/InstrumentMetaDataDatalinkdraddecilluminationdeltaindexprojcellsubcellcrval1crval2crpix1crpix2xyiring
degdegarcsecarcsec%AU
str43float64str1int64float64float64float64float64str8str11int64int64float64float64float64float64int64int64int64float64float64float64float64float64float64int64
rings.v3.skycell.0956.062.wrp.y.54992_6094854992.6100484111y30319.175183-17.487021319.176752-17.4858010956.062Pan-STARRS1-----5.39-4.3997.872.25095662318.139535-18.017623.5-5540.03397.481808.465
rings.v3.skycell.0956.062.wrp.y.54992_6215554992.6221213111y30319.175004-17.487761319.176589-17.4865670956.062Pan-STARRS1-----5.44-4.3097.872.25195662318.139535-18.017623.5-5540.03400.011797.825
rings.v3.skycell.0783.051.wrp.i.55089_3047355089.3054836223i60304.363509-25.843868304.364116-25.8429170783.051Pan-STARRS1-----1.97-3.4298.272.31278351302.926829-26.023341.5240.54719.582387.033
rings.v3.skycell.0783.051.wrp.i.55089_3122255089.3129660222i60304.363118-25.844032304.363786-25.8430800783.051Pan-STARRS1-----2.16-3.4398.272.31378351302.926829-26.023341.5240.54724.672384.723
rings.v3.skycell.1240.063.wrp.y.55416_6158355416.6164040111y3032.570102-1.35337632.569364-1.3533191240.063Pan-STARRS1----2.66-0.2197.202.52412406332.000000-2.011760.0-5521.03552.033789.819
rings.v3.skycell.1240.063.wrp.y.55416_6265555416.6271279111y3032.571300-1.35384132.570629-1.3537811240.063Pan-STARRS1----2.42-0.2297.202.52512406332.000000-2.011760.0-5521.03534.783783.129
rings.v3.skycell.1240.051.wrp.z.55427_6279955427.6285596111z3033.574050-1.96160533.573444-1.9615101240.051Pan-STARRS1----2.18-0.3497.512.37612405132.000000-2.023280.0242.0621.27784.249
rings.v3.skycell.1240.051.wrp.z.55427_6285055427.6290649111z3033.574079-1.96163833.573477-1.9615441240.051Pan-STARRS1----2.17-0.3497.512.37712405132.000000-2.023280.0242.0620.85783.759

Extract and plot the full path of object from JPL Horizons¶

This is not needed and can be skipped to save time. This takes about 6 seconds. It does result in an interesting plot, however.

Query JPL Horizons to get the source position at one day intervals. Undo the RA-wrapping and plot the path for the 5 years of the PS1 survey. Mark the positions where PS1 observations were made.

In [12]:
%%time 

# PS1 epoch range with epoch edges set to round days
emin = Time(54985.0, format='mjd')
emax = Time(57079.0, format='mjd')

print(f"Searching for {source} from {emin.isot} to {emax.isot} ({emax.mjd-emin.mjd:.0f} days)")

# query using 1-day time steps
obj = Horizons(id=source, location='F51', 
               epochs=dict(start=emin.isot, stop=emax.isot, step='1d'))#, 

jpltab = obj.ephemerides(extra_precision=True)
fullname = jpltab['targetname'][0]
print(f"Got positions from Horizons with {len(jpltab)} rows for {fullname}")

# undo RA-wrapping from 0-360 degrees
# add or subtract multiples of 360 to create a smooth path
njumps = 0
ratab = jpltab['RA'].data.copy()
w = np.where(ratab[1:] > ratab[:-1]+180)[0]
for k in w:
    njumps += 1
    ratab[k+1:] -= 360
w = np.where(ratab[1:] < ratab[:-1]-180)[0]
for k in w:
    njumps += 1
    ratab[k+1:] += 360
while ratab.min() < 0:
    ratab += 360

plt.rcParams.update({"font.size":14})
plt.figure(1,(12,6))
plt.plot(ratab, jpltab['DEC'])

# plot every 30 days as a dot to see movement
mstep = int(30/(jpltab['datetime_jd'][1]-jpltab['datetime_jd'][0]) + 0.5)
plt.plot(ratab[::mstep], jpltab['DEC'][::mstep], 'o', color="tab:blue", label="Monthly positions")

# show positions at the times of PS1 observations
jpl_mjd = Time(jpltab['datetime_jd'],format='jd').mjd
ra_interp = np.interp(cadc_tab['MJD'], jpl_mjd, ratab)
dec_interp = np.interp(cadc_tab['MJD'], jpl_mjd, jpltab['DEC'])
plt.plot(ra_interp, dec_interp, 'o', label=f"{len(cadc_tab)} PS1 observations") 

plt.xlabel("Unwrapped RA [deg]")
plt.ylabel("Dec [deg]")
plt.legend(loc="best")
plt.title(f"{fullname} from {emin.isot[:10]} to {emax.isot[:10]}")
plt.tight_layout()
if saveplots:
    plt.savefig("sky_path.png", facecolor="white")

Searching for rickwhite from 2009-06-03T00:00:00.000 to 2015-02-26T00:00:00.000 (2094 days)
Got positions from Horizons with 2095 rows for 174361 Rickwhite (2002 TV315)
CPU times: user 3.01 s, sys: 14.2 ms, total: 3.03 s
Wall time: 3.03 s
No description has been provided for this image

Extract JPEG cutout images to see whether the object is there¶

Note some images are blank or have artifacts near the source positions. Approximately half of the positions have data at the center. That will be typical for these moving object queries.

Blank images have a minimum value of 255 (all white), so they are easily excluded.

There are ~100 cutout images of which ~60 are not completely blank. That will be true for any object, although the exact numbers will vary depending on the fraction of observations where the object winds up in a gap in the camera.

In [13]:
%%time

imsize = 65 # cutout size in 0.25 arcsec pixels
output_size = 130

images = getimages(cadc_tab, size=imsize, output_size=output_size, verbose=True)
print(f"Got {len(images)} images")

# omit blank images
goodimages = []
igood = []
for i, image in enumerate(images):
    if image is not None and np.array(image).min() != 255:
        goodimages.append(image)
        igood.append(i)
ngood = len(goodimages)
igood = np.array(igood)
goodtab = cadc_tab[igood]
# add SrcID
goodtab.add_column(np.arange(len(goodtab)), name="SrcID", index=0)
print(f"{len(goodimages)} of {len(images)} are not blank")

17.5 s: completed 93 of 93 images
Got 93 images
63 of 93 are not blank
CPU times: user 707 ms, sys: 200 ms, total: 907 ms
Wall time: 17.5 s

Show non-blank images in a grid¶

In [14]:
ncols = 7
nrows = (ngood+ncols-1)//ncols
hsize = ncols * 2.5
vsize = nrows * 2.5
extent = np.array([1,-1,-1,1])*0.25*imsize/2
plt.rcParams.update({"font.size":10})
plt.figure(1,(hsize,vsize))
for i in range(ngood):
    plt.subplot(nrows,ncols,i+1)
    plt.imshow(goodimages[i], origin="upper", cmap="gray", vmin=0, vmax=255, extent=extent)
    plt.title(f"{goodtab['Filter'][i]} {goodtab['MJD'][i]:.5f}")
plt.tight_layout()
plt.suptitle(f"{ngood} {imsize}x{imsize} pixel ({imsize/4:.0f}x{imsize/4:.0f} arcsec) PS1 images for {fullname}",
             y=1.01, fontsize=16, fontweight="bold")
if saveplots:
    plt.savefig("images.png",facecolor="white", bbox_inches="tight")
No description has been provided for this image

Test: Extract some images at the original CADC positions to compare with the JPL positions¶

Select the first 14 images with shifts small enough to see on the image. Then make a plot comparing the test images and the good images.

The first and third rows use the CADC position, and the second and fourth rows use the refined JPL Horizons position. The image center is marked. Clearly the updated JPL positions are much better when there is a shift.

This test is turned off by default since it is not needed for the example. Set testing = True to run it.

In [15]:
testing = False
if not testing:
    print("Not testing, skipped this step")
else:
    goodsep = np.sqrt(goodtab['dra'].data**2+goodtab['ddec'].data**2)
    itest = np.where(goodsep < 0.25*imsize/2)[0]
    itest = itest[:14]

    testimages = getimages(goodtab[itest], size=imsize, output_size=output_size, verbose=True,
                           racol='Object_RA', deccol='Object_Dec')
    print(f"Got {len(testimages)} test images")

    ncols = 7
    nrows = (len(itest)*2+ncols-1)//ncols
    hsize = ncols * 2.5
    vsize = nrows * 2.5
    extent = np.array([-1,1,-1,1])*0.25*imsize/2
    plt.rcParams.update({"font.size":10})
    plt.figure(1,(hsize,vsize))
    for i in range(len(itest)):
        j = itest[i]
        plt.subplot(nrows,ncols,i+1+(i//ncols)*ncols)
        plt.imshow(testimages[i], origin="upper", cmap="gray", vmin=0, vmax=255, extent=extent)
        plt.plot(0.0, 0.0, 'o', color="tab:orange", markerfacecolor="none", markersize=15)
        plt.title(f"{goodtab['SrcID'][j]}: CADC {goodtab['Filter'][j]} {goodtab['MJD'][j]:.5f}", color="tab:red")
        plt.subplot(nrows,ncols,i+1+(i//ncols+1)*ncols)
        plt.imshow(goodimages[j], origin="upper", cmap="gray", vmin=0, vmax=255, extent=extent)
        plt.plot(0.0, 0.0, 'o', color="tab:blue", markerfacecolor="none", markersize=15)
        plt.title(f"{goodtab['SrcID'][j]}: JPL {goodtab['Filter'][j]} {goodtab['MJD'][j]:.5f}", color="tab:blue")
    plt.tight_layout()
    _ = plt.suptitle(f"{ngood} {imsize}x{imsize} pixel ({imsize/4:.0f}x{imsize/4:.0f} arcsec) PS1 images for {fullname}",
                     y=1.01, fontsize=16, fontweight="bold")

Not testing, skipped this step

Upload the final match table and extract Detection table information¶

This creates a table MyDB.jplhorizons3 in CasJobs and then does a cross-match to the ObjectThin and Detection tables.

The constraints are:

  • Search radius 2 arcsec
  • MJD match within half of exposure time (typically 15-20 sec)

The matched objects that are found are never close to these limits.

Note the filterID is not included in the constraints. It could be included, but that does not seem necessary given the tight constraint on MJD. There is only one exposure in a single filter at any given time. The next code blocks checks for any filter mismatches, but should not ever find any problems unless there are errors in the database.

If you do not want to see the images, this crossmatch can work from the original cadc_tab rather than the good images list in goodtab. That will return the same number of matches.

In [16]:
%%time

# convert table to a csv string
cnames = ["SrcID", "ra", "dec", "Filter", "MJD"]
with StringIO() as fh:
    goodtab[cnames].write(fh, format="ascii.csv")
    hdata = fh.getvalue()

table = "jplhorizons3"

jobs = mastcasjobs.MastCasJobs(context="MyDB", request_type="POST")
jobs.drop_table_if_exists(table)
print(f"Dropped MyDB.{table}")

jobs.upload_table(table, hdata, exists=False)
print(f"Uploaded table to MyDB.{table}")
 
query = f"""
select t.SrcID, dbo.fDistanceArcminEq(d.ra,d.dec,t.ra,t.dec)*60 as darcsec,
    t.ra as hra, t.dec as hdec,
    t.filter as hfilter, t.mjd as hmjd,
    d.objID, d.detectID, d.filterID, d.imageID, d.obsTime, d.ra, d.dec,
    d.raErr, d.decErr, d.expTime, d.psfFlux, d.psfFluxErr, 
    d.psfMajorFWHM, d.psfMinorFWHM, d.psfTheta, d.psfCore, 
    d.psfQfPerfect, d.psfChiSq, d.apFlux, d.apFluxErr, 
    d.infoFlag, d.infoFlag2, d.infoFlag3
from MyDB.{table} t 
cross apply fGetNearbyObjEq(t.ra, t.dec, 2.0/60.0) nb
join Detection d on nb.objid=d.objid and d.obsTime between t.mjd-0.5*d.expTime/86400 and t.mjd+0.5*d.expTime/86400
order by d.obsTime, darcsec
"""

dtab = jobs.quick(query, context="PanSTARRS_DR2")
print(f"Retrieved {len(dtab)} row table")

dtab['filterID'] = dtab['filterID'].astype(np.uint8)
dtab['darcsec'].format = ".2f"
for c in ('hmjd','obsTime'):
    dtab[c].format = ".5f"
dtab

Dropped MyDB.jplhorizons3
Uploaded table to MyDB.jplhorizons3
Retrieved 46 row table
CPU times: user 45.6 ms, sys: 7.39 ms, total: 52.9 ms
Wall time: 3.01 s
Out[16]:
Table length=46
SrcIDdarcsechrahdechfilterhmjdobjIDdetectIDfilterIDimageIDobsTimeradecraErrdecErrexpTimepsfFluxpsfFluxErrpsfMajorFWHMpsfMinorFWHMpsfThetapsfCorepsfQfPerfectpsfChiSqapFluxapFluxErrinfoFlaginfoFlag2infoFlag3
int64float64float64float64str1float64int64int64uint8int64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64int64int64int64
20.06304.363509-25.843868i55089.3054876983043635077697988305121500036743969301555089.30548304.36352642-25.843862150.05195260.051952660.01.21332e-051.15707e-061.191191.0210189.8058-0.2920290.9981261.132681.17748e-053.2792e-073565158507374912
30.08304.363118-25.844032i55089.3129776983043631157517988312611500026953969381555089.31297304.36313469-25.844014680.04905940.049059460.01.29931e-051.17005e-061.07470.997377-87.965-0.01341260.9980071.13999.41763e-062.94723e-073565158507374912
40.0833.57405-1.961605z55427.6285610564033574016652613266283765000062142095286555427.6285633.57405744-1.961626760.06446870.064468730.01.89817e-052.24692e-060.7962770.73786912.24480.3103980.9989120.8455812.40206e-057.97542e-073565158507374912
60.2328.072054-7.870439r55485.556529855028072065591613845562857000023522414145755485.5565228.07208716-7.870493890.03548950.035489540.03.63915e-052.37473e-062.331192.237854.921650.08059580.998981.101224.16959e-057.27543e-073565158507342144
70.0424.781017-8.892393g55500.350939732024780989955313993506775000028412481277555500.3509324.78101246-8.892404150.03804770.038047743.02.14361e-051.50116e-061.548251.17630.7006-0.193640.9992430.8390431.67801e-054.9641e-073565158507374912
80.0824.778556-8.89289g55500.361959732024778535891313993617067000016412481446755500.3619524.77854484-8.892910420.03694150.036941543.02.10933e-051.43444e-061.427171.3378436.4017-0.08803460.9991140.7036922.32896e-055.87499e-073565158507374912
130.09175.24212912.105662y55934.6669712252175242122731818336667856000060254364705655934.66697175.2421059712.105669480.07394350.073943530.04.78354e-056.49878e-060.8880250.810951-12.71780.1769580.9935530.8585264.39961e-051.70125e-063565158507374912
140.13175.24389912.106305y55934.6752112252175243950809018336750356000039054364865655934.67521175.2439331312.106317940.08897050.088970530.04.46878e-057.3053e-061.020890.899402-3.55538-0.04036980.9981440.8249093.27531e-051.45738e-063565158507374912
160.02174.36623.198033r56000.3160913583174366090823018993158524000009124644392456000.31609174.3659955423.198034620.004228010.0042410140.00.000209861.84846e-060.9757250.87791548.6828-0.1289160.9988791.093630.0002146121.70012e-061027605171287374912
.......................................................................................
520.2314.730025-11.989597i56969.381639361014730084284328683813636000070838140483656969.3816314.73008502-11.989621050.07161650.071616545.03.2669e-054.28825e-062.609892.4563551.31020.07952890.9992751.060012.1796e-055.19359e-073565158507342144
530.2414.728101-11.98955i56969.393529361014728129285628683932536000083338140673656969.3935214.72813339-11.989608540.07160590.071605945.02.65778e-053.48817e-061.948981.7654736.2210.05250980.9989211.031532.33346e-055.44176e-073565158507342144
550.0414.447885-11.976463r56971.263319362014447883867328702630414000380828153021456971.2633114.44787964-11.976452720.0295110.02951145.02.3693e-051.28286e-061.39221.2086663.89510.032360.9981991.206642.54061e-055.51953e-073565158507374912
560.1214.445983-11.97634r56971.275819362014445965877728702755414000312028153211456971.2758114.44596473-11.976368270.0432190.04321945.01.58947e-051.26072e-061.263241.1708871.33110.1836830.9965051.209841.96222e-054.88875e-073565158507374912
570.1714.444077-11.976215r56971.288319362014444119894528702880414000354028153391456971.2883114.44412545-11.97622290.04766430.047664345.01.62038e-051.41711e-061.291831.1846246.70370.174620.99891.116412.30559e-055.23975e-073566796907374912
580.0914.442139-11.976086r56971.301039362014442106912128703007614000385128153601456971.3010314.44211462-11.976080840.03655340.036553445.02.30187e-051.54391e-061.320591.213747.90010.06874880.9986971.136812.8126e-055.77588e-073565158507374912
590.0914.154488-11.946825r56973.347639366014154298414228723473647000026828169954756973.3476314.15446308-11.946818420.03068490.030684945.03.00238e-051.69391e-061.301611.26482-74.44510.2056270.1601020.8526783.09146e-056.03717e-0735651585327374912
600.1914.152793-11.946608r56973.359559366014152861438828723592847000039328170144756973.3595514.15284202-11.946628070.04021830.040218345.02.75852e-052.03998e-061.685481.53239-62.29210.1412240.9984920.8667362.68547e-055.57764e-073565158507374912
610.1714.151106-11.946389r56973.371449366014151071472028723711627000045328170312756973.3714414.15106183-11.946371870.05057050.050570545.02.00454e-051.86387e-061.372541.31645-68.14070.2997330.9991820.8277961.75901e-054.52477e-073565158507374912
620.1014.149334-11.946158r56973.383949366014149331494528723836627000045228170512756973.3839414.14932258-11.946183840.05621570.056215745.01.66604e-051.72211e-061.299441.2293-49.40460.2066590.9984130.8608521.7375e-054.51461e-073565158507374912

Check for mismatches in filter or time¶

Also checks to see whether any SrcID entries have duplicate matches. That can happen if the object happens to be moving through a very crowded region such as the Galactic plane. If there are multiple matches, the closest one is retained. (You might want to look more carefully at such cases, however.)

In [17]:
print((np.abs(dtab['hmjd']-dtab['obsTime']) > 0.00002).sum(), "sources with MJD discrepancy")
print((dtab['filterID'] != [("grizy".find(x)+1) for x in dtab['hfilter']]).sum(), "sources with filter discrepancy")
ndup = len(dtab)-len(np.unique(dtab['SrcID']))
print(f"{ndup} duplicate matches to same SrcID")
if ndup > 0:
    # pick closer of each duplicate entry
    dtab.sort(['SrcID', 'darcsec'])
    h = np.bincount(dtab['SrcID'])
    w = np.where(h > 1)[0]
    keep = np.ones(len(dtab), dtype=bool)
    for s in w:
        ww = np.where(dtab['SrcID']==s)[0]
        keep[ww[1:]] = False
    dtab = dtab[keep]
    print(f"Retaining {len(dtab)} entries with only closest for dups")
        
print(f"Maximum epoch difference: {(dtab['hmjd']-dtab['obsTime']).max()*86400:.4f} s")

0 sources with MJD discrepancy
0 sources with filter discrepancy
0 duplicate matches to same SrcID
Maximum epoch difference: 0.0037 s

Scatterplot showing position differences between PS1 Detection and JPL Horizons¶

The position differences will be a combination of PS1 position errors and uncertainties in the object's ephemeris. The psfQfPerfect flag is used to highlight likely poor quality detections that have bad pixels in the core of the PSF. They account for the points that have large position errors.

In [18]:
plt.rcParams.update({"font.size":14})
plt.figure(1,(12,12))
wgood = np.where(dtab['psfQfPerfect'] > 0.95)[0]
wbad = np.where(dtab['psfQfPerfect'] <= 0.95)[0]
dra = (dtab['ra']-dtab['hra'])*3600*np.cos(np.radians(dtab['hdec']))
ddec = (dtab['dec']-dtab['hdec'])*3600
plt.plot(dra[wgood], ddec[wgood], 'o', label=f"{len(wgood)} good points")
plt.plot(dra[wbad], ddec[wbad], 'o', label=f"{len(wbad)} bad points")
plt.xlabel("Delta RA [arcsec]")
plt.ylabel("Delta Dec [arcsec]")
plt.gca().set_aspect(1.0)
plt.legend(loc="best", title="psfQfPerfect flags")
plt.title(f"PS1 catalog detections for {fullname}")
if saveplots:
    plt.savefig("position-diffs.png", facecolor="white", bbox_inches="tight")
No description has been provided for this image

Show images with sources in catalog marked¶

The shifts are usually too small to see, but this does show which of the images have PS1 detections and which do not.

Magnitudes are shown in the plot titles for detected objects. A magnitude followed by * is dubious because it has psfQfPerfect <= 0.95, indicating that there are bad pixels in the core of the object. The marker uses a different color to highlight those cases.

Note: The plotted positions are not strictly accurate. The cutouts are centered on the nearest pixel from the PS1 image, which means they can be shifted by up to half a pixel (0.125 arcsec) from the requested coordinate system. It would be more accurate to use the WCS for the cutout (which is available for FITS cutouts or in the comment section of JPEG cutouts) to compute the marker position.
In [19]:
goodmatch = join(goodtab, dtab['SrcID','ra','dec','psfFlux','psfQfPerfect'], keys='SrcID', join_type='left')
# offset of source from center
dx = (goodmatch['ra_2']-goodmatch['ra_1'])*3600*np.cos(np.radians(goodmatch['dec_1']))
dy = (goodmatch['dec_2']-goodmatch['dec_1'])*3600

# convert the fluxes to AB magnitudes
mag = -2.5*np.log10(goodmatch['psfFlux']) + 8.90
goodmatch['mag'] = mag

ncols = 7
nrows = (ngood+ncols-1)//ncols
hsize = ncols * 2.5
vsize = nrows * 2.5
extent = np.array([1,-1,-1,1])*0.25*imsize/2

plt.rcParams.update({"font.size":10})
plt.figure(1,(hsize,vsize))
for i in range(ngood):
    plt.subplot(nrows,ncols,i+1)
    plt.imshow(goodimages[i], origin="upper", cmap="gray", vmin=0, vmax=255, extent=extent)
    if hasattr(goodmatch['ra_2'],'mask') and goodmatch['ra_2'].mask[i]:
        # no detection
        plt.title(f"{goodmatch['Filter'][i]} {goodmatch['MJD'][i]:.5f}")
    else:
        # show detection
        if goodmatch['psfQfPerfect'][i]<=0.95:
            # dubious magnitude
            warning = '*'
            color = "tab:orange"
            tcolor = color
        else:
            warning = ''
            color = "tab:cyan"
            tcolor = "tab:blue"
        plt.plot(dx[i], dy[i], 'o', color=color, markerfacecolor="none", markersize=15, markeredgewidth=2)
        plt.title(f"{goodmatch['Filter'][i]}={mag[i]:.2f}{warning} {goodmatch['MJD'][i]:.5f}",
                  color=tcolor)
plt.tight_layout()
# print(f"{ngood} {imsize}x{imsize} pixel ({imsize/4:.0f}x{imsize/4:.0f} arcsec) PS1 images for {fullname}")
plt.suptitle(f"{ngood} {imsize}x{imsize} pixel ({imsize/4:.0f}x{imsize/4:.0f} arcsec) PS1 images for {fullname}",
             y=1.01, fontsize=16, fontweight="bold")
if saveplots:
    plt.savefig("images-catalogs.png", facecolor="white", bbox_inches="tight")
No description has been provided for this image

Plot light curve¶

This removes poor quality measurements using the psfQfPerfect > 0.95 test. It also shows a plot with the magnitudes corrected using a simple model determined by the illumination fraction and distance to the object.

Note that the model does account for most of the variation in apparent brightness.

In [20]:
jdtab = join(dtab, goodtab['SrcID','illumination','delta'])
assert len(jdtab) == len(dtab)
jdtab['mag'] = -2.5*np.log10(jdtab['psfFlux']) + 8.90

# illumination magnitude factor
jdtab['illumfac'] = -2.5*np.log10(jdtab['illumination']/jdtab['delta']**2)
jdtab['illumfac'] -= jdtab['illumfac'].mean()

fdtab = jdtab.group_by('filterID')

plt.rcParams.update({"font.size":14})
plt.figure(1,(10,10))
sub1 = plt.subplot(2,1,1)
for g in fdtab.groups:
    hfilter = g['hfilter'][0]
    ww = np.where(g['psfQfPerfect']>0.95)[0]
    lines = plt.plot(g['hmjd'][ww], g['mag'][ww], 'o-',label=f"{hfilter} ({len(ww)} pts)")
# plt.xlabel("MJD [days]")
plt.ylabel("Brightness [mag]")
ylim = plt.ylim()[::-1]
plt.legend(loc="upper left", fontsize=12)
plt.title(f"PS1 light curve for {fullname}")

sub2 = plt.subplot(2,1,2)
for g in fdtab.groups:
    hfilter = g['hfilter'][0]
    ww = np.where(g['psfQfPerfect']>0.95)[0]
    plt.plot(g['hmjd'][ww], g['mag'][ww]-g['illumfac'][ww], 'o-', label=f"{hfilter} ({len(ww)} pts)")
plt.xlabel("MJD [days]")
ylim2 = plt.ylim()[::-1]
ylim = (max(ylim[0],ylim2[0]), min(ylim[1],ylim2[1]))
sub2.set_ylim(ylim)
sub1.set_ylim(ylim)
plt.ylabel("Brightness - Prediction [mag]")
plt.legend(loc="upper left", fontsize=12)
plt.title("Light curve with distance and illumination effects removed")
plt.tight_layout()
if saveplots:
    plt.savefig("light-curves.png", facecolor="white")
No description has been provided for this image
In [ ]: